calculate sparse subset of inverse of sparse matrix More...

#include <SparseInverse.h>

Public Types
typedef SparseMatrix< Scalar, ColMajor >	MatrixType

typedef SparseMatrix< Scalar, RowMajor >	RowMatrixType

Public Member Functions
	SparseInverse ()

	SparseInverse (const SparseLU< MatrixType > &slu)
	This Constructor is for if you already have a factored SparseLU and would like to use it to calculate a sparse inverse. More...

SparseInverse &	compute (const SparseMatrix< Scalar > &A)
	Calculate the sparse inverse from a given sparse input. More...

const MatrixType &	inverse () const
	return the already-calculated sparse inverse, or a 0x0 matrix if it could not be computed More...

Static Public Member Functions
static MatrixType	computeInverse (const SparseLU< MatrixType > &slu)
	Internal function to calculate the sparse inverse in a functional way. More...

static MatrixType	computeInverse (const RowMatrixType &Upper, const Matrix< Scalar, Dynamic, 1 > &inverseDiagonal, const MatrixType &Lower)
	Internal function to calculate the inverse from strictly upper, diagonal and strictly lower components. More...

Private Attributes
MatrixType	_result

Detailed Description

template<typename Scalar>
class Eigen::SparseInverse< Scalar >

calculate sparse subset of inverse of sparse matrix

This class returns a sparse subset of the inverse of the input matrix. The nonzeros correspond to the nonzeros of the input, plus any additional elements required due to fill-in of the internal LU factorization. This is is minimized via a applying a fill-reducing permutation as part of the LU factorization.

If there are specific entries of the input matrix which you need inverse values for, which are zero for the input, you need to insert entries into the input sparse matrix for them to be calculated.

Due to the sensitive nature of matrix inversion, particularly on large matrices which are made possible via sparsity, high accuracy dot products based on Kahan summation are used to reduce numerical error. If you still encounter numerical errors you may with to equilibrate your matrix before calculating the inverse, as well as making sure it is actually full rank.

Member Typedef Documentation

◆ MatrixType

template<typename Scalar >

typedef SparseMatrix<Scalar, ColMajor> Eigen::SparseInverse< Scalar >::MatrixType

◆ RowMatrixType

template<typename Scalar >

typedef SparseMatrix<Scalar, RowMajor> Eigen::SparseInverse< Scalar >::RowMatrixType

Constructor & Destructor Documentation

◆ SparseInverse() [1/2]

template<typename Scalar >

Eigen::SparseInverse< Scalar >::SparseInverse ( )

inline

133 {}

◆ SparseInverse() [2/2]

template<typename Scalar >

Eigen::SparseInverse< Scalar >::SparseInverse ( const SparseLU< MatrixType > & slu )

inline

This Constructor is for if you already have a factored SparseLU and would like to use it to calculate a sparse inverse.

Just call this constructor with your already factored SparseLU class and you can directly call the .inverse() method to get the result.

142 { _result = computeInverse(slu); }

Eigen::SparseInverse::computeInverse

static MatrixType computeInverse(const SparseLU< MatrixType > &slu)

Internal function to calculate the sparse inverse in a functional way.

Definition: SparseInverse.h:163

Eigen::SparseInverse::_result

MatrixType _result

Definition: SparseInverse.h:228

References Eigen::SparseInverse< Scalar >::_result, and Eigen::SparseInverse< Scalar >::computeInverse().

Member Function Documentation

◆ compute()

template<typename Scalar >

SparseInverse& Eigen::SparseInverse< Scalar >::compute ( const SparseMatrix< Scalar > & A )

inline

Calculate the sparse inverse from a given sparse input.

                                                         {
     SparseLU<MatrixType> slu;
     slu.compute(A);
     _result = computeInverse(slu);
     return *this;
   }

References Eigen::SparseInverse< Scalar >::_result, Eigen::SparseLU< MatrixType_, OrderingType_ >::compute(), and Eigen::SparseInverse< Scalar >::computeInverse().

Referenced by check_sparse_inverse().

◆ computeInverse() [1/2]

template<typename Scalar >

static MatrixType Eigen::SparseInverse< Scalar >::computeInverse	(	const RowMatrixType &	Upper,
		const Matrix< Scalar, Dynamic, 1 > &	inverseDiagonal,
		const MatrixType &	Lower
	)

inlinestatic

Internal function to calculate the inverse from strictly upper, diagonal and strictly lower components.

                                                             {
     // Calculate the 'minimal set', which is the nonzeros of (L+U).transpose()
     // It could be zeroed, but we will overwrite all non-zeros anyways.
     MatrixType colInv = Lower.transpose().template triangularView<UnitUpper>();
     colInv += Upper.transpose();
  
     // We also need rowmajor representation in order to do efficient row-wise dot products
     RowMatrixType rowInv = Upper.transpose().template triangularView<UnitLower>();
     rowInv += Lower.transpose();
  
     // Use the Takahashi algorithm to build the supporting elements of the inverse
     // upwards and to the left, from the bottom right element, 1 col/row at a time
     for (Index recurseLevel = Upper.cols() - 1; recurseLevel >= 0; recurseLevel--) {
       const auto& col = Lower.col(recurseLevel);
       const auto& row = Upper.row(recurseLevel);
  
       // Calculate the inverse values for the nonzeros in this column
       typename MatrixType::ReverseInnerIterator colIter(colInv, recurseLevel);
       for (; recurseLevel < colIter.index(); --colIter) {
         const Scalar element = -accurateDot(col, rowInv.row(colIter.index()));
         colIter.valueRef() = element;
         rowInv.coeffRef(colIter.index(), recurseLevel) = element;
       }
  
       // Calculate the inverse values for the nonzeros in this row
       typename RowMatrixType::ReverseInnerIterator rowIter(rowInv, recurseLevel);
       for (; recurseLevel < rowIter.index(); --rowIter) {
         const Scalar element = -accurateDot(row, colInv.col(rowIter.index()));
         rowIter.valueRef() = element;
         colInv.coeffRef(recurseLevel, rowIter.index()) = element;
       }
  
       // And finally the diagonal, which corresponds to both row and col iterator now
       const Scalar diag = inverseDiagonal(recurseLevel) - accurateDot(row, colInv.col(recurseLevel));
       rowIter.valueRef() = diag;
       colIter.valueRef() = diag;
     }
  
     return colInv;
   }

References Eigen::accurateDot(), Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::coeffRef(), col(), diag, Eigen::Lower, row(), and Eigen::Upper.

◆ computeInverse() [2/2]

template<typename Scalar >

static MatrixType Eigen::SparseInverse< Scalar >::computeInverse ( const SparseLU< MatrixType > & slu )

inlinestatic

Internal function to calculate the sparse inverse in a functional way.

Returns: A sparse inverse representation, or, if the decomposition didn't complete, a 0x0 matrix.

                                                                     {
     if (slu.info() != Success) {
       return MatrixType(0, 0);
     }
  
     // Extract from SparseLU and decompose into L, inverse D and U terms
     Matrix<Scalar, Dynamic, 1> invD;
     RowMatrixType Upper;
     {
       RowMatrixType DU = slu.matrixU().toSparse();
       invD = DU.diagonal().cwiseInverse();
       Upper = (invD.asDiagonal() * DU).template triangularView<StrictlyUpper>();
     }
     MatrixType Lower = slu.matrixL().toSparse().template triangularView<StrictlyLower>();
  
     // Compute the inverse and reapply the permutation matrix from the LU decomposition
     return slu.colsPermutation().transpose() * computeInverse(Upper, invD, Lower) * slu.rowsPermutation();
   }