LU decomposition of a matrix with complete pivoting, and related features. More...
Public Types | |
enum | { MaxSmallDimAtCompileTime } |
typedef Matrix< Scalar, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1 > | ColVectorType |
typedef Matrix< typename MatrixType::Scalar, MatrixType::RowsAtCompileTime, Dynamic, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime > | ImageResultType |
typedef Matrix< int, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1 > | IntColVectorType |
typedef Matrix< int, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime > | IntRowVectorType |
typedef Matrix< typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, Dynamic, MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxColsAtCompileTime > | KernelResultType |
typedef NumTraits< typename MatrixType::Scalar >::Real | RealScalar |
typedef Matrix< Scalar, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime > | RowVectorType |
typedef MatrixType::Scalar | Scalar |
Public Member Functions | |
template<typename ImageMatrixType > | |
void | computeImage (ImageMatrixType *result) const |
template<typename ResultType > | |
void | computeInverse (ResultType *result) const |
template<typename KernelMatrixType > | |
void | computeKernel (KernelMatrixType *result) const |
ei_traits< MatrixType >::Scalar | determinant () const |
int | dimensionOfKernel () const |
const ImageResultType | image () const |
MatrixType | inverse () const |
bool | isInjective () const |
bool | isInvertible () const |
bool | isSurjective () const |
const KernelResultType | kernel () const |
LU (const MatrixType &matrix) | |
const MatrixType & | matrixLU () const |
const IntColVectorType & | permutationP () const |
const IntRowVectorType & | permutationQ () const |
int | rank () const |
template<typename OtherDerived , typename ResultType > | |
bool | solve (const MatrixBase< OtherDerived > &b, ResultType *result) const |
Protected Attributes | |
int | m_det_pq |
MatrixType | m_lu |
const MatrixType & | m_originalMatrix |
IntColVectorType | m_p |
RealScalar | m_precision |
IntRowVectorType | m_q |
int | m_rank |
LU decomposition of a matrix with complete pivoting, and related features.
MatrixType | the type of the matrix of which we are computing the LU decomposition |
This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q are permutation matrices. This is a rank-revealing LU decomposition. The eigenvalues (diagonal coefficients) of U are sorted in such a way that any zeros are at the end, so that the rank of A is the index of the first zero on the diagonal of U (with indices starting at 0) if any.
This decomposition provides the generic approach to solving systems of linear equations, computing the rank, invertibility, inverse, kernel, and determinant.
This LU decomposition is very stable and well tested with large matrices. Even exact rank computation works at sizes larger than 1000x1000. However there are use cases where the SVD decomposition is inherently more stable when dealing with numerically damaged input. For example, computing the kernel is more stable with SVD because the SVD can determine which singular values are negligible while LU has to work at the level of matrix coefficients that are less meaningful in this respect.
The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP(), permutationQ().
As an exemple, here is how the original matrix can be retrieved:
typedef Matrix<double, 5, 3> Matrix5x3; typedef Matrix<double, 5, 5> Matrix5x5; Matrix5x3 m = Matrix5x3::Random(); cout << "Here is the matrix m:" << endl << m << endl; Eigen::LU<Matrix5x3> lu(m); cout << "Here is, up to permutations, its LU decomposition matrix:" << endl << lu.matrixLU() << endl; cout << "Here is the L part:" << endl; Matrix5x5 l = Matrix5x5::Identity(); l.block<5,3>(0,0).part<StrictlyLowerTriangular>() = lu.matrixLU(); cout << l << endl; cout << "Here is the U part:" << endl; Matrix5x3 u = lu.matrixLU().part<UpperTriangular>(); cout << u << endl; cout << "Let us now reconstruct the original matrix m:" << endl; Matrix5x3 x = l * u; Matrix5x3 y; for(int i = 0; i < 5; i++) for(int j = 0; j < 3; j++) y(i, lu.permutationQ()[j]) = x(lu.permutationP()[i], j); cout << y << endl; // should be equal to the original matrix m
Output:
Here is the matrix m: 0.68 -0.605 -0.0452 -0.211 -0.33 0.258 0.566 0.536 -0.27 0.597 -0.444 0.0268 0.823 0.108 0.904 Here is, up to permutations, its LU decomposition matrix: 0.904 0.823 0.108 -0.299 0.812 0.569 -0.05 0.888 -1.1 0.0296 0.705 0.768 0.285 -0.549 0.0436 Here is the L part: 1 0 0 0 0 -0.299 1 0 0 0 -0.05 0.888 1 0 0 0.0296 0.705 0.768 1 0 0.285 -0.549 0.0436 0 1 Here is the U part: 0.904 0.823 0.108 0 0.812 0.569 0 0 -1.1 0 0 0 0 0 0 Let us now reconstruct the original matrix m: 0.68 -0.605 -0.0452 -0.211 -0.33 0.258 0.566 0.536 -0.27 0.597 -0.444 0.0268 0.823 0.108 0.904
LU | ( | const MatrixType & | matrix | ) |
Constructor.
matrix | the matrix of which to compute the LU decomposition. |
void computeImage | ( | ImageMatrixType * | result | ) | const |
Computes a basis of the image of the matrix, also called the column-space or range of he matrix.
result | a pointer to the matrix in which to store the image. The columns of this matrix will be set to form a basis of the image (it will be resized if necessary). |
Example:
MatrixXd m(3,3); m << 1,1,0, 1,3,2, 0,1,1; cout << "Here is the matrix m:" << endl << m << endl; LU<Matrix3d> lu(m); // allocate the matrix img with the correct size to avoid reallocation MatrixXd img(m.rows(), lu.rank()); lu.computeImage(&img); cout << "Notice that the middle column is the sum of the two others, so the " << "columns are linearly dependent." << endl; cout << "Here is a matrix whose columns have the same span but are linearly independent:" << endl << img << endl;
Output:
Here is the matrix m: 1 1 0 1 3 2 0 1 1 Notice that the middle column is the sum of the two others, so the columns are linearly dependent. Here is a matrix whose columns have the same span but are linearly independent: 1 1 3 1 1 0
void computeInverse | ( | ResultType * | result | ) | const [inline] |
Computes the inverse of the matrix of which *this is the LU decomposition.
result | a pointer to the matrix into which to store the inverse. Resized if needed. |
void computeKernel | ( | KernelMatrixType * | result | ) | const |
Computes a basis of the kernel of the matrix, also called the null-space of the matrix.
result | a pointer to the matrix in which to store the kernel. The columns of this matrix will be set to form a basis of the kernel (it will be resized if necessary). |
Example:
MatrixXf m = MatrixXf::Random(3,5); cout << "Here is the matrix m:" << endl << m << endl; LU<MatrixXf> lu(m); // allocate the matrix ker with the correct size to avoid reallocation MatrixXf ker(m.rows(), lu.dimensionOfKernel()); lu.computeKernel(&ker); cout << "Here is a matrix whose columns form a basis of the kernel of m:" << endl << ker << endl; cout << "By definition of the kernel, m*ker is zero:" << endl << m*ker << endl;
Output:
Here is the matrix m: 0.68 0.597 -0.33 0.108 -0.27 -0.211 0.823 0.536 -0.0452 0.0268 0.566 -0.605 -0.444 0.258 0.904 Here is a matrix whose columns form a basis of the kernel of m: -0.219 0.763 0.00335 -0.447 0 1 1 0 -0.145 -0.285 By definition of the kernel, m*ker is zero: -9.46e-09 1.96e-08 2.53e-10 3.82e-08 6.06e-10 1.45e-08
ei_traits< MatrixType >::Scalar determinant | ( | ) | const |
int dimensionOfKernel | ( | ) | const [inline] |
const LU< MatrixType >::ImageResultType image | ( | ) | const |
Example:
Matrix3d m; m << 1,1,0, 1,3,2, 0,1,1; cout << "Here is the matrix m:" << endl << m << endl; cout << "Notice that the middle column is the sum of the two others, so the " << "columns are linearly dependent." << endl; cout << "Here is a matrix whose columns have the same span but are linearly independent:" << endl << m.lu().image() << endl;
Output:
Here is the matrix m: 1 1 0 1 3 2 0 1 1 Notice that the middle column is the sum of the two others, so the columns are linearly dependent. Here is a matrix whose columns have the same span but are linearly independent: 1 1 3 1 1 0
MatrixType inverse | ( | void | ) | const [inline] |
bool isInjective | ( | ) | const [inline] |
bool isInvertible | ( | ) | const [inline] |
bool isSurjective | ( | ) | const [inline] |
const LU< MatrixType >::KernelResultType kernel | ( | ) | const |
Example:
MatrixXf m = MatrixXf::Random(3,5); cout << "Here is the matrix m:" << endl << m << endl; MatrixXf ker = m.lu().kernel(); cout << "Here is a matrix whose columns form a basis of the kernel of m:" << endl << ker << endl; cout << "By definition of the kernel, m*ker is zero:" << endl << m*ker << endl;
Output:
Here is the matrix m: 0.68 0.597 -0.33 0.108 -0.27 -0.211 0.823 0.536 -0.0452 0.0268 0.566 -0.605 -0.444 0.258 0.904 Here is a matrix whose columns form a basis of the kernel of m: -0.219 0.763 0.00335 -0.447 0 1 1 0 -0.145 -0.285 By definition of the kernel, m*ker is zero: -9.46e-09 1.96e-08 2.53e-10 3.82e-08 6.06e-10 1.45e-08
const MatrixType& matrixLU | ( | ) | const [inline] |
const IntColVectorType& permutationP | ( | ) | const [inline] |
const IntRowVectorType& permutationQ | ( | ) | const [inline] |
int rank | ( | ) | const [inline] |
bool solve | ( | const MatrixBase< OtherDerived > & | b, |
ResultType * | result | ||
) | const |
This method finds a solution x to the equation Ax=b, where A is the matrix of which *this is the LU decomposition, if any exists.
b | the right-hand-side of the equation to solve. Can be a vector or a matrix, the only requirement in order for the equation to make sense is that b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition. |
result | a pointer to the vector or matrix in which to store the solution, if any exists. Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). If no solution exists, *result is left with undefined coefficients. |
Example:
typedef Matrix<float,2,3> Matrix2x3; typedef Matrix<float,3,2> Matrix3x2; Matrix2x3 m = Matrix2x3::Random(); Matrix2f y = Matrix2f::Random(); cout << "Here is the matrix m:" << endl << m << endl; cout << "Here is the matrix y:" << endl << y << endl; Matrix3x2 x; if(m.lu().solve(y, &x)) { assert(y.isApprox(m*x)); cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; } else cout << "The equation mx=y does not have any solution." << endl;
Output:
Here is the matrix m: 0.68 0.566 0.823 -0.211 0.597 -0.605 Here is the matrix y: -0.33 -0.444 0.536 0.108 Here is a solution x to the equation mx=y: 0 0 0.291 -0.216 -0.6 -0.391