LinAlg Library

Classes
class	Imagine::Matrix< T >
	Matrix of variable size. More...
class	Imagine::SymMatrix< T >
	Symmetric Matrix of variable size. More...
class	Imagine::Vector< T >
	Vector of variable size. More...

Functions
template<typename T , int N>
FMatrix< T, N, N >	Imagine::cholesky (const FMatrix< T, N, N > &A, bool low=true)
	Cholesky decomposition. More...
template<typename T , int N>
T	Imagine::detFMatrix (const FMatrix< T, N, N > &A)
	Determinant. More...
template<typename T >
Matrix< T >	Imagine::Diagonal (const Vector< T > &d)
	Diagonal. More...
template<typename T >
void	Imagine::eigenvalues (const Matrix< T > &A, Vector< T > &wr, Vector< T > &wi)
	Eigenvalues calculation Given square matrix, only eigenvalues are obtained. More...
template<typename T , int N>
FMatrix< T, N, N >	Imagine::inverseFMatrix (const FMatrix< T, N, N > &A)
	Inverse. More...
template<typename T , int M, int N>
FVector< T, N >	Imagine::linSolve (const FMatrix< T, M, N > &A, const FVector< T, M > &b)
	Linear system. More...
template<typename T >
Matrix< T >	Imagine::outerProduct (const Vector< T > &v1, const Vector< T > &v2)
	Outer product matrix. More...
template<typename T , int M, int N>
FMatrix< T, N, M >	Imagine::pseudoInverse (const FMatrix< T, M, N > &A, T tolrel=0)
	Pseudo inverse. More...
template<typename T , int M, int N>
bool	Imagine::QR (const FMatrix< T, M, N > &A, FMatrix< T, M, N > &Q, FMatrix< T, N, N > &R)
	QR decomposition. More...
template<typename T , int M, int N>
bool	Imagine::QRAll (const FMatrix< T, M, N > &A, FMatrix< T, M, M > &Q, FMatrix< T, M, N > &R)
	QR decomposition. More...
template<typename T , int M, int N>
void	Imagine::svd (const FMatrix< T, M, N > &A, FMatrix< T, M, M > &U, FVector< T, M > &S, FMatrix< T, N, N > &Vt, bool all=false)
	SVD. More...

Detailed Description

Function Documentation

cholesky()

template<typename T , int N>

FMatrix<T,N,N> Imagine::cholesky	(	const FMatrix< T, N, N > &	A,
		bool	low = true
	)

Cholesky decomposition of symmetric positive-definite matrix.

Parameters

A	matrix
low	Choose lower (default) or upper true: A=L*L^T with L lower false: A=U^T*U with U upper

Returns

L or U

    FMatrix<T,5,5> L=cholesky(E);                   // Cholesky decomposition:
    cout << "Lower Cholesky " 
         << norm(E-L*transpose(L)) << endl;         // -lower
    FMatrix<T,5,5> V=cholesky(E,false);
    cout << "Upper Cholesky " 
         << norm(E-transpose(V)*V) << endl;         // -upper...

Examples:

LinAlg/test/test.cpp.

detFMatrix()

template<typename T , int N>

T Imagine::detFMatrix

(

const FMatrix< T, N, N > &

A

)

Determinant. Try Cholesky decomposition just in case. If not sym def pos, uses SVD.

Parameters

A

matrix

Returns

determinant

    FMatrix<T,5,5> E=B*transpose(B);                // sym pos-def
    cout << "Determinants: "     // Determinant
         << detFMatrix(B)*detFMatrix(inverseFMatrix(B))-1 << endl;
    cout << "Determinants: "                        
         << detFMatrix(E)*detFMatrix(inverseFMatrix(E))-1 << endl;
                                 // Check sym pos-def acceleration...

Examples:

LinAlg/test/test.cpp.

Diagonal()

template<typename T >

Matrix<T> Imagine::Diagonal ( const Vector< T > &  d )

inline

Diagonal matrix

Parameters

d	diagonal vector

Returns

matrix

    T tt[]={1,2,3,4,5,6};                // pre-allocated
    Vector<T> d(tt,6);                   // diagonal
    B=Diagonal(d);                       // ...

eigenvalues()

template<typename T >

void Imagine::eigenvalues	(	const Matrix< T > &	A,
		Vector< T > &	wr,
		Vector< T > &	wi
	)

Parameters

A	square matrix to decompose
wr	real parts of eigenvalues
wi	imaginary parts of eigenvalues

    T coeff[] = {1,2,2,1};               // Eigenvalues
    Matrix<T> E(coeff,2,2);
    Vector<T> lambdaR, lambdaI;
    eigenvalues(E,lambdaR, lambdaI);
    cout << "Eigen values of matrix E=" << E;
    cout << "Real parts: "      << lambdaR << " (should be 3, -1)" << endl;
    cout << "Imaginary parts: " << lambdaI << " (should be 0)" << endl; // ...

Examples:

LinAlg/test/test.cpp.

inverseFMatrix()

template<typename T , int N>

FMatrix<T,N,N> Imagine::inverseFMatrix

(

const FMatrix< T, N, N > &

A

)

Inverse of an FMatrix. If non invertible, ouptuts a message to cerr and returns a matrix with zeroed elements. This works for any FMatrix size, on the contrary to inverse(), the latter being limited to N<=3.

Parameters

A	matrix to inverse

Returns

the inverse matrix

    FMatrix<T,3,3> A2=A*transpose(A);
    A2=inverseFMatrix(A2);               // test fast call to direct inversion
    FMatrix<T,4,4> A3=transpose(A)*A;
    FMatrix<T,4,4> A4 = inverseFMatrix(A3);               // inverse

Examples:

LinAlg/test/test.cpp.

linSolve()

template<typename T , int M, int N>

FVector<T,N> Imagine::linSolve	(	const FMatrix< T, M, N > &	A,
		const FVector< T, M > &	b
	)

Solve linear system. Returns x such that:

• M=N, square system: Ax=b
• M>N, over determined system: x=argmin_y|Ay-b|
• M<N, under determined system: x=argmin_{y,Ay=b}|y|

Parameters

A	matrix
b	right term

Returns

solution (zeroed vector if no solution)

    initRandom(1);                       // Linear solve:
    FMatrix<T,5,5> B;
    for (int i=0;i<5;i++)
        for (int j=0;j<5;j++)
            B(i,j)=T(doubleRandom());
    FVector<T,5> b,x;
    for (int j=0;j<5;j++)
        b[j]=T(doubleRandom());
    x=linSolve(B,b);
    cout << "Square system: " << norm(B*x-b) << endl; // -square
    FMatrix<T,5,4> C( B.data() );
    FVector<T,4> x2=linSolve(C,b);
    cout << "Over determined system: " 
         << norm(C*x2-b) << endl;                   // -over-determined (min possible but not 0)
    FMatrix<T,4,5> D( B.data() );
    FVector<T,4> d(b.data());
    FVector<T,5> x3=linSolve(D,d);
    cout << "Minimum norm under determined system: " 
        << norm(D*x3-d) << " " << norm(x3) << endl; // -under-determined (min |x| possible but not 0)...

Examples:

LinAlg/test/test.cpp.

outerProduct()

template<typename T >

Matrix<T> Imagine::outerProduct	(	const Vector< T > &	v1,
		const Vector< T > &	v2
	)

Vector outer product, i.e. the matrix given by v1 * v2^T

Parameters

v1,v2

arguments

Returns

Result

    Matrix<T> O=outerProduct(v1,v2); // outer product

Examples:

LinAlg/test/test.cpp.

pseudoInverse()

template<typename T , int M, int N>

FMatrix<T,N,M> Imagine::pseudoInverse	(	const FMatrix< T, M, N > &	A,
		T	tolrel = 0
	)

Pseudo-inverse using SVD

Parameters

A	matrix
tolrel	relative tolerance (under which a sing. value is considered nul)

Returns

pseudo-inverse

    FMatrix<T,5,5> B;
    FVector<T,5> b,x;
    FMatrix<T,5,4> C( B.data() );
    FMatrix<T,4,5> D( B.data() );
    FVector<T,4> d(b.data());
    x=pseudoInverse(B)*b;                         // Pseudo-inverse:
    cout << "Pseudo inverse, square: " 
         << norm(B*x-b) << endl;                  // -square case
    x2=pseudoInverse(C)*b;     
    cout << "Pseudo inverse, over determined: " 
         << norm(C*x2-b) << endl;                 // -over determined case
    x3=pseudoInverse(D)*d; // Erreur quand passage en double
    cout << "Pseudo inverse, under determined: " 
        << norm(D*x3-d) << " " << norm(x3) << endl; // -under determined case...

Examples:

LinAlg/test/test.cpp.

QR()

template<typename T , int M, int N>

bool Imagine::QR	(	const FMatrix< T, M, N > &	A,
		FMatrix< T, M, N > &	Q,
		FMatrix< T, N, N > &	R
	)

QR decomposition. A is MxN with M>=N. A=QR where:

• Q is MxN with orthogonal columns
• R is NxN upper triangular

  Parameters

  A	matrix
  Q,R	output

      FMatrix<T,3,2> Q; FMatrix<T,2,2> R;             // QR decomposition:
      QR(D2,Q,R);
      cout << "QR " << norm(D2-Q*R) << endl;          // -thin QR
      FMatrix<T,3,3> Q2; FMatrix<T,3,2> R2;
      QRAll(D2,Q2,R2);
      cout << "QR " << norm(D2-Q2*R2) << endl;        // -full QR...

Examples:

LinAlg/test/test.cpp.

QRAll()

template<typename T , int M, int N>

bool Imagine::QRAll	(	const FMatrix< T, M, N > &	A,
		FMatrix< T, M, M > &	Q,
		FMatrix< T, M, N > &	R
	)

QR decomposition. A is MxN with M>=N. A=QR where:

• Q is MxM orthogonal
• R is MxN with N first rows upper diagonal and M-N last rows nul.

  Parameters

  A	matrix
  Q,R	output

      FMatrix<T,3,2> Q; FMatrix<T,2,2> R;             // QR decomposition:
      QR(D2,Q,R);
      cout << "QR " << norm(D2-Q*R) << endl;          // -thin QR
      FMatrix<T,3,3> Q2; FMatrix<T,3,2> R2;
      QRAll(D2,Q2,R2);
      cout << "QR " << norm(D2-Q2*R2) << endl;        // -full QR...

Examples:

LinAlg/test/test.cpp.

svd()

template<typename T , int M, int N>

void Imagine::svd	(	const FMatrix< T, M, N > &	A,
		FMatrix< T, M, M > &	U,
		FVector< T, M > &	S,
		FMatrix< T, N, N > &	Vt,
		bool	all = false
	)

Singular value decomposition. A is MxN. S is of size M with decreasing singular values (only min(M,N) first coefficients are filled). A=U*diag(S)*Vt where U (resp. Vt) is MxM (resp NxN) orthonormal, diag(S) is MxN with S as diagonal. If M<N (resp. M>N) then Vt (resp. U) is incompletely filled, except if all is true.

Parameters

A	matrix
U,S,Vt	SVD output
all	completely fill U and V, even columns multiplied by zero

    FVector<T,3> S;                      // Singular value decomposition:
    FMatrix<T,3,3> U, Vt;
    svd(A2,U,S,Vt);                      // -square
    cout << "SVD check 1: " 
         << norm(A2-U*Diagonal(S)*Vt) << endl;
    FMatrix<T,4,4> Vt2;
    svd(A,U,S,Vt2);                      // -non square
    FMatrix<T,3,4> S2(T(0));
    for(int i=0; i < 3; i++)
        S2(i,i) = S[i];
    cout << "SVD check 2: " 
         << norm(A-U*S2*Vt2) << endl;
    svd(A,U,S,Vt2,true);                 // non square matrix, complete U and Vt
    cout << "SVD check 3: "              
         << norm(Vt2*transpose(Vt2)-FMatrix<T,4,4>::Identity()) << endl;
                                         // Check Vt...

See equivalent function for Matrix.

Examples:

LinAlg/test/test.cpp.