LinearAlgebra.C File Reference

#include "Image/LinearAlgebra.H"
#include "Image/Image.H"
#include "Image/MatrixOps.H"
#include "Image/gsl.H"
#include "Image/lapack.H"
#include "rutz/trace.h"
#include "inst/Image/LinearAlgebra.I"

Include dependency graph for LinearAlgebra.C:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
void	svd (const Image< double > &A, Image< double > &U, Image< double > &S, Image< double > &V, SvdAlgo algo, SvdFlag flags)
	Compute the singular decomposition of M: M = U * S * V^T.
void	svdf (const Image< float > &A, Image< float > &U, Image< float > &S, Image< float > &V, SvdAlgo algo, SvdFlag flags)
	Single-precision counterpart of svd().
Image< double >	svdPseudoInv (const Image< double > &x, const SvdAlgo algo, int *const rank, const double tolfactor)
	Compute the pseudo-inverse using the method described in svd().
Image< float >	svdPseudoInvf (const Image< float > &x, const SvdAlgo algo, int *const rank, const float tolfactor)
	Single-precision counterpart of svdPseudoInv().
template<class T >
Image< typename promote_trait < T, float >::TP >	naiveUnstablePseudoInv (const Image< T > &M)
	NOTE that this procedure is called "naive, unstable" for a reason!
template Image< promote_trait < double, float >::TP >	naiveUnstablePseudoInv (const Image< double > &)

---

Detailed Description

Definition in file LinearAlgebra.C.

---

Function Documentation

template<class T >

Image<typename promote_trait<T, float>::TP> naiveUnstablePseudoInv

(

const Image< T > &

M

)

[inline]

NOTE that this procedure is called "naive, unstable" for a reason!

This procedure computes the pseudo-inverse using the direct formula pinv(M) = (M^T * M)^-1 * M^T. The formula is analytically correct, but in practice it can be numerically unstable. Indeed it might work just fine for your matrix, but to be safe you'll need to check the identity pinv(M)*M ~= I -- for example, do RMSerr(matrixMult(M,pinv(M)), eye<float>(M.getHeight())) and make sure that it's small (say, less than 0.01 or so). If you do get a large RMS error, then you can consider using svdPseudoInv() instead, which offers more numerical stability.

Definition at line 217 of file LinearAlgebra.C.

References matrixInv(), matrixMult(), and transpose().

void svd	(	const Image< double > &	M,
		Image< double > &	U,
		Image< double > &	S,
		Image< double > &	V,
		SvdAlgo	algo,
		const SvdFlag	flags = 0
	)

Compute the singular decomposition of M: M = U * S * V^T.

The U, S, V matrices have special properties:

• U is column-orthogonal, such that (U^T * U) = I
• S is square, diagonal, so that S is easily invertible
• V is square and orthogonal, so (V^T * V) = (V * V^T) = I

Note that these properties mean that the svd can be used to compute the pseudo-inverse of a general matrix M:

• M = U*S*V^T
• U^T*M = S*V^T (since U is column-orthogonal)
• S^-1*U^T*M = V^T (since S is diagonal and invertible)
• V*S^-1*U^T*M = I (since V is square and orthogonal)

by definition, pinv(M)*M = I, therefore

• V*S^-1*U^T = pinv(M)

(Note ^T is used to designate the matrix transpose, and ^-1 is used to designate the matrix inverse).

You can use the routine svdPseudoInv to compute the pseudo-inverse using just this method.

Definition at line 108 of file LinearAlgebra.C.

References SVD_GSL, SVD_LAPACK, and SVD_TALL.

void svdf	(	const Image< float > &	A,
		Image< float > &	U,
		Image< float > &	S,
		Image< float > &	V,
		SvdAlgo	algo,
		SvdFlag	flags
	)

Single-precision counterpart of svd().

Definition at line 141 of file LinearAlgebra.C.

References SVD_GSL, SVD_LAPACK, and SVD_TALL.

Image<double> svdPseudoInv	(	const Image< double > &	x,
		const SvdAlgo	algo,
		int *const	rank = 0,
		const double	tolfactor = 1.0e-8
	)

Compute the pseudo-inverse using the method described in svd().

Optionally throw out small values from the diagonal S matrix -- this helps numerical stability.

The number of eigenvalues retained is controlled by tolfactor. If tolfactor is less than 1.0, then it serves as a multiplier of the largest eigenvalue, below which any smaller eigenvalues are thrown out. If tolfactor is greater than 1.0, then the integer portion of tolfactor specifies the number of eigenvalues to be retained.

If rank is a non-null pointer, then the rank of the diagonal matrix (after setting small values to zero) will be returned through that pointer.

Definition at line 174 of file LinearAlgebra.C.

References Image< T >::getHeight(), Image< T >::getWidth(), matrixMult(), and transpose().

Image<float> svdPseudoInvf	(	const Image< float > &	x,
		const SvdAlgo	algo,
		int *const	rank,
		const float	tolfactor
	)

Single-precision counterpart of svdPseudoInv().

Definition at line 196 of file LinearAlgebra.C.

References Image< T >::getHeight(), Image< T >::getWidth(), matrixMult(), and transpose().