Imagine++
Classes | Functions
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>
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

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
Amatrix
lowChoose 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.
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
Amatrix
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.
template<typename T >
Matrix<T> Imagine::Diagonal ( const Vector< T > &  d)
inline

Diagonal matrix

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

template<typename T >
void Imagine::eigenvalues ( const Matrix< T > &  A,
Vector< T > &  wr,
Vector< T > &  wi 
)
Parameters
Asquare matrix to decompose
wrreal parts of eigenvalues
wiimaginary parts of eigenvalues
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
Amatrix 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.
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
Amatrix
bright 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.
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,v2arguments
Returns
Result
Matrix<T> O=outerProduct(v1,v2); // outer product

Examples:
LinAlg/test/test.cpp.
template<typename T , int M, int N>
FMatrix<T,N,M> Imagine::pseudoInverse ( const FMatrix< T, M, N > &  A,
tolrel = 0 
)

Pseudo-inverse using SVD

Parameters
Amatrix
tolrelrelative 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.
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
    Amatrix
    Q,Routput
    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.
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
    Amatrix
    Q,Routput
    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.
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 min(M,N) with decreasing singular values. 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, expect all is true.

Parameters
Amatrix
U,S,VtSVD output
allcomplete U and V is true
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: " // Check Vt...

Examples:
LinAlg/test/test.cpp.