linear algebra algorithms - see the c++ source of an example

This file provides a short description of algorithms implemented in file "capd/matrixAlgorithms/floatMatrixAlgorithms.hpp" from CAPD package - for more details see the source file.

See also:

Content of this file:

Solving of linear equations by Gauss elimination

The following function solves system of linear equations by using Gauss elimination algorithm

   template<typename MatrixType, typename VectorType>
   VectorType gauss(MatrixType a, VectorType b);

Since both arguments are modified in this function, they are copied to the actual arguments a,b

Example of usage:

   // we create a matrix
   double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15};
   Matrix<double,0,0> P(4,4,data);
   Vector<double,0> v = P.row(0);
   // here we call the function gauss
   Vector<double,0> solution = gauss(P,v);
   std::cout << "gauss(P,v)=" << solution << std::endl;
   // this should be close to v
   std::cout << "verification of the result:\nP*gauss(P,v)=" << P*solution << std::endl;

Orthonormalization of a matrix

The following function orthonormalizes columns of a nonsingular matrix. The function throws an exception when the matrix is singular.

   template<typename MatrixType>
   void orthonormalize(MatrixType& Q);

The above function performs Gramm-Schmidt orthonormalization of columns of matrix Q sorted by decreasing norm. The argument Q is modified by the function and after calling orthonormalize(Q), the matrix Q contains result of this function.

Example of usage:

   double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15};
   Matrix<double,0,0> Q(4,4,data);
   std::cout << "orthonormalization of Q=" << Q << std::endl;
   orthonormalize(Q);
   std::cout << "is equal to " << Q << std::endl;

QR decomposition

The following function computes QR-decomposition of a square nonsingular matrix.

   template<typename MatrixType>
   void QR_decompose(const MatrixType& A, MatrixType& Q, MatrixType& R);

As a result we obtain orthogonal and upper triangular matrix. First argument of the function is a matrix to be decomposed. Matrices Q and R are modified by the function and contain the result.

Example of usage:

   double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15};
   Matrix<double,0,0> P(4,4,data), Q(4,4), R(4,4);
   QR_decompose(P,Q,R); 
   // P should be equal to Q*R
   std::cout << "QR decomposition of P:" << std::endl;
   std::cout << "Q=" << Q << std::endl;
   std::cout << "R=" << R << std::endl;
   std::cout << "verification if P-Q*R is close to zero\n";
   std::cout << "P-Q*R=" << P-Q*R << std::endl;

The function throws an exception when the matrix is singular.

Diagonalization of a symmetric matrix

The following function performs Jacobi rotation algorithm for computing diagonalization of a square symmetric matrix.

   template<typename MatrixType>
   int symMatrixDiagonalize(const MatrixType& A, MatrixType& D, typename MatrixType::ScalarType diagonalizingRelTolerance = typename MatrixType::ScalarType(1e-5));

The argument D contains the result of the function. The last argument is a diagonalizing tolerance we want to obtain, namely if by 'd' and 's' we denote

the algorithm stops when s < d*diagonalizingRelTolerance or the next iteration would worsen the estimation. It should be noted that the procedure does not guarantee that such tolerance is achieved. This procedure is the first step in functions which compute spectral radius and maximal eigenvalue of a symmetric matrix.

Example of usage:

   double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15};
   Matrix<double,0,0> P(4,4,data), D(4,4);
   Matrix<double,0,0> S = 0.5*(P+Transpose(P));
   // diagonalization of S=0.5*(P+P^T)
   symMatrixDiagonalize(S,D);
   std::cout << "Diagonalization of 0.5*(P+P^T)=" << D << std::endl;

Spectral radius of a symmetric matrix

The following function computes an upper bound for spectral radius of a symmetric matrix.

   template<typename MatrixType>
   typename MatrixType::ScalarType spectralRadiusOfSymMatrix(const MatrixType &A, typename MatrixType::ScalarType diagonalizingRelTolerance = typename MatrixType::ScalarType(1e-5));

First, the function computes matrix which has the same eigenvalues and which is close to diagonal, next upper bound is computed from Gerschgorin theorem

Example of usage:

   double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15};
   Matrix<double,0,0> P(4,4,data);
   Matrix<double,0,0> S = 0.5*(P+Transpose(P));
   std::cout << "spectralRadiusOfSymMatrix(S)=" << spectralRadiusOfSymMatrix(S) << std::endl;

Maximal eigenvalue of a symmetric matrix

The following function computes an upper bound for maximal eigenvalue of a symmetric matrix.

   template<typename MatrixType>
   typename MatrixType::ScalarType maxEigenValueOfSymMatrix(const MatrixType &A, typename MatrixType::ScalarType diagonalizingRelTolerance = typename MatrixType::ScalarType(1e-5));

First, the function computes matrix which has the same eigenvalues and which is close to diagonal, next upper bound is computed from Gerschgorin theorem

Example of usage:

   double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15};
   Matrix<double,0,0> P(4,4,data);
   Matrix<double,0,0> S = 0.5*(P+Transpose(P));
   std::cout << "maxEigenValueOfSymMatrix(S)=" << maxEigenValueOfSymMatrix(S) << std::endl;

Computation of inverse matrix

The matrixAlgorithm module implements two functions for computing the inverse of a matrix. First function decomposes a matrix A=Q*R, where Q-orthogonal, R-upper diagonal and compute A-1 = R-1*QT

   template<typename MatrixType>
   MatrixType inverseMatrix(const MatrixType &A);

The second function (recommended) computes the inverse matrix using Gauss elimination algorithm

   template<typename MatrixType>
   MatrixType gaussInverseMatrix(const MatrixType &A);

Example of usage:

   double data[] = {1,2,1,2,2,2,2,0,1,97,3,4,1,19,1,15};
   Matrix<double,0,0> P(4,4,data);
   Matrix<double,0,0>  invP = inverseMatrix(P);
   std::cout << "Inverse of P by using QR-decomposition:" << invP << std::endl;
   std::cout << "verification: P*inverseMatrix(P)=" << P*invP << std::endl;

   invP = gaussInverseMatrix(P);
   std::cout << "Inverse of P by using Gauss elimination:" << invP << std::endl;
   std::cout << "verification: P*gaussInverseMatrix(P)=" << P*invP << std::endl;

Author: Daniel Wilczak, last modified on July 7, 2008.