Public Types | Public Member Functions | Protected Attributes
LU< MatrixType > Class Template Reference

LU decomposition of a matrix with complete pivoting, and related features. More...

List of all members.

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 IntColVectorTypepermutationP () const
const IntRowVectorTypepermutationQ () 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

Detailed Description

template<typename MatrixType>
class Eigen::LU< MatrixType >

LU decomposition of a matrix with complete pivoting, and related features.

Parameters:
MatrixTypethe 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
See also:
MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse()

Constructor & Destructor Documentation

LU ( const MatrixType &  matrix)

Constructor.

Parameters:
matrixthe matrix of which to compute the LU decomposition.

Member Function Documentation

void computeImage ( ImageMatrixType *  result) const

Computes a basis of the image of the matrix, also called the column-space or range of he matrix.

Note:
Calling this method on the zero matrix will make an assertion fail.
Parameters:
resulta 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
See also:
image(), computeKernel(), kernel()
void computeInverse ( ResultType *  result) const [inline]

Computes the inverse of the matrix of which *this is the LU decomposition.

Parameters:
resulta pointer to the matrix into which to store the inverse. Resized if needed.
Note:
If this matrix is not invertible, *result is left with undefined coefficients. Use isInvertible() to first determine whether this matrix is invertible.
See also:
MatrixBase::computeInverse(), inverse()
void computeKernel ( KernelMatrixType *  result) const

Computes a basis of the kernel of the matrix, also called the null-space of the matrix.

Note:
This method is only allowed on non-invertible matrices, as determined by isInvertible(). Calling it on an invertible matrix will make an assertion fail.
Parameters:
resulta 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:
-1.12e-08  2.98e-08
 -1.4e-09 -3.86e-08
 1.49e-08         0
See also:
kernel(), computeImage(), image()
ei_traits< MatrixType >::Scalar determinant ( ) const
Returns:
the determinant of the matrix of which *this is the LU decomposition. It has only linear complexity (that is, O(n) where n is the dimension of the square matrix) as the LU decomposition has already been computed.
Note:
This is only for square matrices.
For fixed-size matrices of size up to 4, MatrixBase::determinant() offers optimized paths.
Warning:
a determinant can be very big or small, so for matrices of large enough dimension, there is a risk of overflow/underflow.
See also:
MatrixBase::determinant()
int dimensionOfKernel ( ) const [inline]
Returns:
the dimension of the kernel of the matrix of which *this is the LU decomposition.
Note:
Since the rank is computed at the time of the construction of the LU decomposition, this method almost does not perform any further computation.
const LU< MatrixType >::ImageResultType image ( ) const
Returns:
the image of the matrix, also called its column-space. The columns of the returned matrix will form a basis of the kernel.
Note:
: Calling this method on the zero matrix will make an assertion fail.
: this method returns a matrix by value, which induces some inefficiency. If you prefer to avoid this overhead, use computeImage() instead.

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
See also:
computeImage(), kernel()
MatrixType inverse ( void  ) const [inline]
Returns:
the inverse of the matrix of which *this is the LU decomposition.
Note:
If this matrix is not invertible, the returned matrix has undefined coefficients. Use isInvertible() to first determine whether this matrix is invertible.
See also:
computeInverse(), MatrixBase::inverse()
bool isInjective ( ) const [inline]
Returns:
true if the matrix of which *this is the LU decomposition represents an injective linear map, i.e. has trivial kernel; false otherwise.
Note:
Since the rank is computed at the time of the construction of the LU decomposition, this method almost does not perform any further computation.
bool isInvertible ( ) const [inline]
Returns:
true if the matrix of which *this is the LU decomposition is invertible.
Note:
Since the rank is computed at the time of the construction of the LU decomposition, this method almost does not perform any further computation.
bool isSurjective ( ) const [inline]
Returns:
true if the matrix of which *this is the LU decomposition represents a surjective linear map; false otherwise.
Note:
Since the rank is computed at the time of the construction of the LU decomposition, this method almost does not perform any further computation.
const LU< MatrixType >::KernelResultType kernel ( ) const
Returns:
the kernel of the matrix, also called its null-space. The columns of the returned matrix will form a basis of the kernel.
Note:
: this method is only allowed on non-invertible matrices, as determined by isInvertible(). Calling it on an invertible matrix will make an assertion fail.
: this method returns a matrix by value, which induces some inefficiency. If you prefer to avoid this overhead, use computeKernel() instead.

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:
-1.12e-08  2.98e-08
 -1.4e-09 -3.86e-08
 1.49e-08         0
See also:
computeKernel(), image()
const MatrixType& matrixLU ( ) const [inline]
Returns:
the LU decomposition matrix: the upper-triangular part is U, the unit-lower-triangular part is L (at least for square matrices; in the non-square case, special care is needed, see the documentation of class LU).
See also:
matrixL(), matrixU()
const IntColVectorType& permutationP ( ) const [inline]
Returns:
a vector of integers, whose size is the number of rows of the matrix being decomposed, representing the P permutation i.e. the permutation of the rows. For its precise meaning, see the examples given in the documentation of class LU.
See also:
permutationQ()
const IntRowVectorType& permutationQ ( ) const [inline]
Returns:
a vector of integers, whose size is the number of columns of the matrix being decomposed, representing the Q permutation i.e. the permutation of the columns. For its precise meaning, see the examples given in the documentation of class LU.
See also:
permutationP()
int rank ( ) const [inline]
Returns:
the rank of the matrix of which *this is the LU decomposition.
Note:
This is computed at the time of the construction of the LU decomposition. This method does not perform any further computation.
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.

Parameters:
bthe 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.
resulta 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.
Returns:
true if any solution exists, false if no solution exists.
Note:
If there exist more than one solution, this method will arbitrarily choose one. If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel().

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
See also:
MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()

The documentation for this class was generated from the following file: