MatrixType | the type of the matrix of which we are computing the LU decomposition |
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
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
Public Types | |
enum | { MaxSmallDimAtCompileTime } |
typedef Matrix< Scalar, MatrixType::RowsAtCompileTime, 1 > | ColVectorType |
typedef Matrix< typename MatrixType::Scalar, MatrixType::RowsAtCompileTime, Dynamic, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime > | ImageResultType |
typedef Matrix< int, MatrixType::RowsAtCompileTime, 1 > | IntColVectorType |
typedef Matrix< int, 1, MatrixType::ColsAtCompileTime > | 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 > | RowVectorType |
typedef MatrixType::Scalar | Scalar |
Public Member Functions | |
template<typename ImageMatrixType> | |
void | computeImage (ImageMatrixType *result) const |
void | computeInverse (MatrixType *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 |
IntRowVectorType | m_q |
int | m_rank |
LU | ( | const MatrixType & | matrix | ) | [inline] |
void computeImage | ( | ImageMatrixType * | result | ) | const [inline] |
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). |
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;
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 | ( | MatrixType * | 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 [inline] |
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). |
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;
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: -1.12e-08 2.98e-08 -1.4e-09 -3.86e-08 1.49e-08 0
ei_traits< MatrixType >::Scalar determinant | ( | ) | const [inline] |
For fixed-size matrices of size up to 4, MatrixBase::determinant() offers optimized paths.
int dimensionOfKernel | ( | ) | const [inline] |
const LU< MatrixType >::ImageResultType image | ( | ) | const [inline] |
: this method returns a matrix by value, which induces some inefficiency. If you prefer to avoid this overhead, use computeImage() instead.
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;
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 [inline] |
: this method returns a matrix by value, which induces some inefficiency. If you prefer to avoid this overhead, use computeKernel() instead.
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;
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: -1.12e-08 2.98e-08 -1.4e-09 -3.86e-08 1.49e-08 0
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 [inline] |
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. |
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;
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