[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

| Advanced Matrix Algebra [Linear Algebra] |

Solution of linear systems, eigen systems, linear least squares etc. More...

## Functions

template<class T, class C1, class C2>
bool | choleskyDecomposition (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &L)
template<class T, class C1>
T | determinant (MultiArrayView< 2, T, C1 > const &a, std::string method="LU")
template<class T, class C>
TemporaryMatrix< T > | inverse (const MultiArrayView< 2, T, C > &v)
template<class T, class C1, class C2>
bool | inverse (const MultiArrayView< 2, T, C1 > &v, MultiArrayView< 2, T, C2 > &res)
template<class T, class C1, class C2, class C3>
bool | leastSquares (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, std::string method="QR")
template<class T, class C1, class C2, class C3>
bool | linearSolve (const MultiArrayView< 2, T, C1 > &A, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &res, std::string method="QR")
template<class T, class C1, class C2, class C3>
bool | linearSolveLowerTriangular (const MultiArrayView< 2, T, C1 > &l, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
template<class T, class C1, class C2, class C3>
bool | linearSolveUpperTriangular (const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
template<class T, class C1>
T | logDeterminant (MultiArrayView< 2, T, C1 > const &a)
template<class T, class C1, class C2, class C3>
bool | nonsymmetricEigensystem (MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, std::complex< T >, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
template<class POLYNOMIAL, class VECTOR>
bool | polynomialRealRootsEigenvalueMethod (POLYNOMIAL const &p, VECTOR &roots, bool polishRoots)
template<class POLYNOMIAL, class VECTOR>
bool | polynomialRootsEigenvalueMethod (POLYNOMIAL const &poly, VECTOR &roots, bool polishRoots)
template<class T, class C1, class C2, class C3>
bool | qrDecomposition (MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, T, C2 > &q, MultiArrayView< 2, T, C3 > &r, double epsilon=0.0)
template<class T, class C1, class C2, class C3>
bool | reverseElimination (const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
template<class T, class C1, class C2, class C3>
bool | ridgeRegression (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, double lambda)
template<class T, class C1, class C2, class C3, class Array>
bool | ridgeRegressionSeries (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, Array const &lambda)
template<class T, class C1, class C2, class C3, class C4>
unsigned int | singularValueDecomposition (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &U, MultiArrayView< 2, T, C3 > &S, MultiArrayView< 2, T, C4 > &V)
template<class T, class C1, class C2, class C3>
bool | symmetricEigensystem (MultiArrayView< 2, T, C1 > const &a, MultiArrayView< 2, T, C2 > &ew, MultiArrayView< 2, T, C3 > &ev)
template<class T, class C1, class C2, class C3, class C4>
bool | weightedLeastSquares (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > const &weights, MultiArrayView< 2, T, C4 > &x, std::string method="QR")
template<class T, class C1, class C2, class C3, class C4>
bool | weightedRidgeRegression (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > const &weights, MultiArrayView< 2, T, C4 > &x, double lambda)

## Detailed Description

Solution of linear systems, eigen systems, linear least squares etc.

## Function Documentation

template<class T, class C1, class C2, class C3>
bool vigra::linalg::symmetricEigensystem | ( | MultiArrayView< 2, T, C1 > const & | a,
|  | MultiArrayView< 2, T, C2 > & | ew,
|  | MultiArrayView< 2, T, C3 > & | ev |
| ) |  |  |

Compute the eigensystem of a symmetric matrix.

a is a real symmetric matrix, ew is a single-column matrix holding the eigenvalues, and ev is a matrix of the same size as a whose columns are the corresponding eigenvectors. Eigenvalues will be sorted from largest to smallest magnitude. The algorithm returns `false` when it doesn't converge. It can be applied in-place, i.e. `&a == &ev` is allowed. The code of this function was adapted from JAMA.

#include <vigra/eigensystem.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3>
bool vigra::linalg::nonsymmetricEigensystem | ( | MultiArrayView< 2, T, C1 > const & | a,
|  | MultiArrayView< 2, std::complex< T >, C2 > & | ew,
|  | MultiArrayView< 2, T, C3 > & | ev |
| ) |  |  |

Compute the eigensystem of a square, but not necessarily symmetric matrix.

a is a real square matrix, ew is a single-column matrix holding the possibly complex eigenvalues, and ev is a matrix of the same size as a whose columns are the corresponding eigenvectors. Eigenvalues will be sorted from largest to smallest magnitude. The algorithm returns `false` when it doesn't converge. It can be applied in-place, i.e. `&a == &ev` is allowed. The code of this function was adapted from JAMA.

#include <vigra/eigensystem.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class POLYNOMIAL, class VECTOR>
bool vigra::linalg::polynomialRootsEigenvalueMethod | ( | POLYNOMIAL const & | poly,
|  | VECTOR & | roots,
|  | bool | polishRoots |
| ) |  |  |

Compute the roots of a polynomial using the eigenvalue method.

poly is a real polynomial (compatible to vigra::PolynomialView), and roots a complex valued vector (compatible to `std::vector` with a `value_type` compatible to the type `POLYNOMIAL::Complex`) to which the roots are appended. The function calls nonsymmetricEigensystem() with the standard companion matrix yielding the roots as eigenvalues. It returns `false` if it fails to converge.

#include <vigra/eigensystem.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

See also:
  polynomialRoots(), vigra::Polynomial

template<class POLYNOMIAL, class VECTOR>
bool vigra::linalg::polynomialRealRootsEigenvalueMethod | ( | POLYNOMIAL const & | p,
|  | VECTOR & | roots,
|  | bool | polishRoots |
| ) |  |  |

Compute the real roots of a real polynomial using the eigenvalue method.

poly is a real polynomial (compatible to vigra::PolynomialView), and roots a real valued vector (compatible to `std::vector` with a `value_type` compatible to the type `POLYNOMIAL::Real`) to which the roots are appended. The function calls polynomialRootsEigenvalueMethod() and throws away all complex roots. It returns `false` if it fails to converge.

#include <vigra/eigensystem.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

See also:
  polynomialRealRoots(), vigra::Polynomial

template<class T, class C1, class C2>
bool vigra::linalg::inverse | ( | const MultiArrayView< 2, T, C1 > & | v,
|  | MultiArrayView< 2, T, C2 > & | res |
| ) |  |  |

Create the inverse or pseudo-inverse of matrix v.

If the matrix v is square, res must have the same shape and will contain the inverse of v. If v is rectangular, it must have more rows than columns, and res must have the transposed shape of v. The inverse is then computed in the least-squares sense, i.e. res will be the pseudo-inverse (Moore-Penrose inverse). The function returns `true` upon success, and `false` if v is not invertible (has not full rank). The inverse is computed by means of QR decomposition. This function can be applied in-place.

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C>
TemporaryMatrix<T> vigra::linalg::inverse | ( | const MultiArrayView< 2, T, C > & | v | ) |

Create the inverse or pseudo-inverse of matrix v.

The result is returned as a temporary matrix. If the matrix v is square, the result will have the same shape and contains the inverse of v. If v is rectangular, it must have more rows than columns, and the result will have the transposed shape of v. The inverse is then computed in the least-squares sense, i.e. res will be the pseudo-inverse (Moore-Penrose inverse). The inverse is computed by means of QR decomposition. If v is not invertible, `vigra::PreconditionViolation` exception is thrown. Usage:

            vigra::Matrix<double> v(n, n);
            v = ...;

            vigra::Matrix<double> m = inverse(v);

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1>
T vigra::linalg::determinant | ( | MultiArrayView< 2, T, C1 > const & | a,
|  | std::string | method = `"LU"` |
| ) |  |  |

Compute the determinant of a square matrix.

method must be one of the following:

"Cholesky"
  Compute the solution by means of Cholesky decomposition. This method is faster than "LU", but requires the matrix a to be symmetric positive definite. If this is not the case, a `ContractViolation` exception is thrown.

"LU"
  (default) Compute the solution by means of LU decomposition.

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1>
T vigra::linalg::logDeterminant | ( | MultiArrayView< 2, T, C1 > const & | a | ) |

Compute the logarithm of the determinant of a symmetric positive definite matrix.

This is useful to avoid multiplication of very large numbers in big matrices. It is implemented by means of Cholesky decomposition.

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2>
bool vigra::linalg::choleskyDecomposition | ( | MultiArrayView< 2, T, C1 > const & | A,
|  | MultiArrayView< 2, T, C2 > & | L |
| ) |  |  |

Cholesky decomposition.

A must be a symmetric positive definite matrix, and L will be a lower triangular matrix, such that (up to round-off errors):

            A == L * transpose(L);

This implementation cannot be applied in-place, i.e. `&L == &A` is an error. If A is not symmetric, a `ContractViolation` exception is thrown. If it is not positive definite, the function returns `false`.

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3>
bool vigra::linalg::qrDecomposition | ( | MultiArrayView< 2, T, C1 > const & | a,
|  | MultiArrayView< 2, T, C2 > & | q,
|  | MultiArrayView< 2, T, C3 > & | r,
|  | double | epsilon = `0.0` |
| ) |  |  |

QR decomposition.

a contains the original matrix, results are returned in q and r, where q is a orthogonal matrix, and r is an upper triangular matrix, such that (up to round-off errors):

            a == q * r;

If a dosn't have full rank, the function returns `false`. The decomposition is computed by householder transformations. It can be applied in-place, i.e. `&a == &q` or `&a == &r` are allowed.

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3>
bool vigra::linalg::reverseElimination | ( | const MultiArrayView< 2, T, C1 > & | r,
|  | const MultiArrayView< 2, T, C2 > & | b,
|  | MultiArrayView< 2, T, C3 > | x |
| ) |  |  |

Deprecated, use linearSolveUpperTriangular().

template<class T, class C1, class C2, class C3>
bool vigra::linalg::linearSolveUpperTriangular | ( | const MultiArrayView< 2, T, C1 > & | r,
|  | const MultiArrayView< 2, T, C2 > & | b,
|  | MultiArrayView< 2, T, C3 > | x |
| ) |  |  |

Solve a linear system with upper-triangular coefficient matrix.

The square matrix r must be an upper-triangular coefficient matrix as can, for example, be obtained by means of QR decomposition. If r doesn't have full rank the function fails and returns `false`, otherwise it returns `true`. The lower triangular part of matrix r will not be touched, so it doesn't need to contain zeros.

The column vectors of matrix b are the right-hand sides of the equation (several equations with the same coefficients can thus be solved in one go). The result is returned int x, whose columns contain the solutions for the corresponding columns of b. This implementation can be applied in-place, i.e. `&b == &x` is allowed. The following size requirements apply:

            rowCount(r) == columnCount(r);
            rowCount(r) == rowCount(b);
            columnCount(r) == rowCount(x);
            columnCount(b) == columnCount(x);

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3>
bool vigra::linalg::linearSolveLowerTriangular | ( | const MultiArrayView< 2, T, C1 > & | l,
|  | const MultiArrayView< 2, T, C2 > & | b,
|  | MultiArrayView< 2, T, C3 > | x |
| ) |  |  |

Solve a linear system with lower-triangular coefficient matrix.

The square matrix l must be a lower-triangular coefficient matrix. If l doesn't have full rank the function fails and returns `false`, otherwise it returns `true`. The upper triangular part of matrix l will not be touched, so it doesn't need to contain zeros.

The column vectors of matrix b are the right-hand sides of the equation (several equations with the same coefficients can thus be solved in one go). The result is returned in x, whose columns contain the solutions for the correspoinding columns of b. This implementation can be applied in-place, i.e. `&b == &x` is allowed. The following size requirements apply:

            rowCount(l) == columnCount(l);
            rowCount(l) == rowCount(b);
            columnCount(l) == rowCount(x);
            columnCount(b) == columnCount(x);

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3>
bool vigra::linalg::linearSolve | ( | const MultiArrayView< 2, T, C1 > & | A,
|  | const MultiArrayView< 2, T, C2 > & | b,
|  | MultiArrayView< 2, T, C3 > & | res,
|  | std::string | method = `"QR"` |
| ) |  |  |

Solve a linear system.

A is the coefficient matrix, and the column vectors in b are the right-hand sides of the equation (so, several equations with the same coefficients can be solved in one go). The result is returned in res, whose columns contain the solutions for the corresponding columns of b. The number of columns of A must equal the number of rows of both b and res, and the number of columns of b and res must match.

method must be one of the following:

"Cholesky"
  Compute the solution by means of Cholesky decomposition. The coefficient matrix A must by symmetric positive definite. If this is not the case, the function returns `false`.

"QR"
  (default) Compute the solution by means of QR decomposition. The coefficient matrix A can be square or rectangular. In the latter case, it must have more rows than columns, and the solution will be computed in the least squares sense. If A doesn't have full rank, the function returns `false`.

"SVD"
  Compute the solution by means of singular value decomposition. The coefficient matrix A can be square or rectangular. In the latter case, it must have more rows than columns, and the solution will be computed in the least squares sense. If A doesn't have full rank, the function returns `false`.

"NE"
  Compute the solution by means of the normal equations, i.e. by applying Cholesky decomposition to the equivalent problem `A'*A*x = A'*b`. This only makes sense when the equation is to be solved in the least squares sense, i.e. when A is a rectangular matrix with more rows than columns. If A doesn't have full column rank, the function returns `false`.

This function can be applied in-place, i.e. `&b == &res` or `&A == &res` are allowed (provided they have the required shapes).

The following size requirements apply:

            rowCount(r) == rowCount(b);
            columnCount(r) == rowCount(x);
            columnCount(b) == columnCount(x);

#include <vigra/linear_solve.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3>
bool vigra::linalg::leastSquares | ( | MultiArrayView< 2, T, C1 > const & | A,
|  | MultiArrayView< 2, T, C2 > const & | b,
|  | MultiArrayView< 2, T, C3 > & | x,
|  | std::string | method = `"QR"` |
| ) |  |  |

Ordinary Least Squares Regression.

Given a matrix A with `m` rows and `n` columns (with `m >= n`), and a column vector b of length `m` rows, this function computes a column vector x of length `n` rows such that the residual

is minimized. When b is a matrix with `k` columns, x must also have `k` columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

This function is just another name for linearSolve(), perhaps leading to more readable code when A is a rectangular matrix. It returns `false` when the rank of A is less than `n`. See linearSolve() for more documentation.

#include <vigra/regression_8hxx.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3, class C4>
bool vigra::linalg::weightedLeastSquares | ( | MultiArrayView< 2, T, C1 > const & | A,
|  | MultiArrayView< 2, T, C2 > const & | b,
|  | MultiArrayView< 2, T, C3 > const & | weights,
|  | MultiArrayView< 2, T, C4 > & | x,
|  | std::string | method = `"QR"` |
| ) |  |  |

Weighted Least Squares Regression.

Given a matrix A with `m` rows and `n` columns (with `m >= n`), a vector b of length `m`, and a weight vector weights of length `m` with non-negative entries, this function computes a vector x of length `n` such that the weighted residual

is minimized, where `diag(weights)` creates a diagonal matrix from weights. The algorithm calls leastSquares() on the equivalent problem

where the square root of weights is just taken element-wise.

When b is a matrix with `k` columns, x must also have `k` columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

The function returns `false` when the rank of the weighted matrix A is less than `n`.

#include <vigra/regression.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3>
bool vigra::linalg::ridgeRegression | ( | MultiArrayView< 2, T, C1 > const & | A,
|  | MultiArrayView< 2, T, C2 > const & | b,
|  | MultiArrayView< 2, T, C3 > & | x,
|  | double | lambda |
| ) |  |  |

Ridge Regression.

Given a matrix A with `m` rows and `n` columns (with `m >= n`), a vector b of length `m`, and a regularization parameter `lambda >= 0.0`, this function computes a vector x of length `n` such that the residual

is minimized. This is implemented by means of singularValueDecomposition().

When b is a matrix with `k` columns, x must also have `k` columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

The function returns `false` if the rank of A is less than `n` and `lambda == 0.0`.

#include <vigra/regression.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3, class C4>
bool vigra::linalg::weightedRidgeRegression | ( | MultiArrayView< 2, T, C1 > const & | A,
|  | MultiArrayView< 2, T, C2 > const & | b,
|  | MultiArrayView< 2, T, C3 > const & | weights,
|  | MultiArrayView< 2, T, C4 > & | x,
|  | double | lambda |
| ) |  |  |

Weighted ridge Regression.

Given a matrix A with `m` rows and `n` columns (with `m >= n`), a vector b of length `m`, a weight vector weights of length `m` with non-negative entries, and a regularization parameter `lambda >= 0.0` this function computes a vector x of length `n` such that the weighted residual

is minimized, where `diag(weights)` creates a diagonal matrix from weights. The algorithm calls ridgeRegression() on the equivalent problem

where the square root of weights is just taken element-wise. This solution is computed by means of singularValueDecomposition().

When b is a matrix with `k` columns, x must also have `k` columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

The function returns `false` if the rank of A is less than `n` and `lambda == 0.0`.

#include <vigra/regression.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3, class Array>
bool vigra::linalg::ridgeRegressionSeries | ( | MultiArrayView< 2, T, C1 > const & | A,
|  | MultiArrayView< 2, T, C2 > const & | b,
|  | MultiArrayView< 2, T, C3 > & | x,
|  | Array const & | lambda |
| ) |  |  |

Ridge Regression with many lambdas.

This executes ridgeRegression() for a sequence of regularization parameters. This is implemented so that the singularValueDecomposition() has to be executed only once. lambda must be an array conforming to the `std::vector` interface, i.e. must support `lambda.size()` and `lambda[k]`. The columns of the matrix x will contain the solutions for the corresponding lambda, so the number of columns of the matrix x must be equal to `lambda.size()`, and b must be a columns vector, i.e. cannot contain several right hand sides at once.

The function returns `false` when the matrix A is rank deficient. If this happens, and one of the lambdas is zero, the corresponding column of x will be skipped.

#include <vigra/regression.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

template<class T, class C1, class C2, class C3, class C4>
unsigned int vigra::linalg::singularValueDecomposition | ( | MultiArrayView< 2, T, C1 > const & | A,
|  | MultiArrayView< 2, T, C2 > & | U,
|  | MultiArrayView< 2, T, C3 > & | S,
|  | MultiArrayView< 2, T, C4 > & | V |
| ) |  |  |

Singular Value Decomposition.

For an m-by-n matrix A with m >= n, the singular value decomposition is an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and an n-by-n orthogonal matrix V so that A = U*S*V'.

To save memory, this functions stores the matrix S in a column vector of appropriate length (a diagonal matrix can be obtained by `diagonalMatrix(S)`). The singular values, sigma[k] = S(k, 0), are ordered so that sigma[0] >= sigma[1] >= ... >= sigma[n-1].

The singular value decomposition always exists, so this function will never fail (except if the shapes of the argument matrices don't match). The effective numerical rank of A is returned.

(Adapted from JAMA, a Java Matrix Library, developed jointly by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).

#include <vigra/singular_value_decomposition.hxx> or
#include <vigra/linear_algebra.hxx>
Namespaces: vigra and vigra::linalg

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python VIGRA 1.6.0 (13 Aug 2008)