GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Enable inverse() method for computing pseudo-inverse.

Browse Source
This commit is contained in:
Rasmus Munk Larsen 2016-02-09 20:35:20 -08:00
parent 53f60e0afc
commit bb8811c655
2 changed files with 26 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
Eigen/src/QR/CompleteOrthogonalDecomposition.h
View File

@ -70,6 +70,7 @@ class CompleteOrthogonalDecomposition {
MatrixType, typename internal::remove_all< MatrixType, typename internal::remove_all<
typename HCoeffsType::ConjugateReturnType>::type> typename HCoeffsType::ConjugateReturnType>::type>
HouseholderSequenceType; HouseholderSequenceType;
typedef typename MatrixType::PlainObject PlainObject;
private: private:
typedef typename PermutationType::Index PermIndexType; typedef typename PermutationType::Index PermIndexType;
@ -246,26 +247,15 @@ class CompleteOrthogonalDecomposition {
*/ */
inline bool isInvertible() const { return m_cpqr.isInvertible(); } inline bool isInvertible() const { return m_cpqr.isInvertible(); }
/** \returns the inverse of the matrix of which *this is the complete /** \returns the pseudo-inverse of the matrix of which *this is the complete
* orthogonal decomposition. * orthogonal decomposition.
* * \warning: Do not compute \c this->inverse()*rhs to solve a linear systems.
* \note If this matrix is not invertible, the returned matrix has undefined * It is more efficient and numerically stable to call \c this->solve(rhs).
* coefficients. Use isInvertible() to first determine whether this matrix is
* invertible.
*/ */
inline const Inverse<CompleteOrthogonalDecomposition> inverse() const
// TODO(rmlarsen): Add method for pseudo-inverse. {
// inline const return Inverse<CompleteOrthogonalDecomposition>(*this);
// internal::solve_retval<CompleteOrthogonalDecomposition, typename }
// MatrixType::IdentityReturnType>
// inverse() const
// {
// eigen_assert(m_isInitialized && "CompleteOrthogonalDecomposition is not
// initialized.");
// return internal::solve_retval<CompleteOrthogonalDecomposition,typename
// MatrixType::IdentityReturnType>
// (*this, MatrixType::Identity(m_cpqr.rows(), m_cpqr.cols()));
// }
inline Index rows() const { return m_cpqr.rows(); } inline Index rows() const { return m_cpqr.rows(); }
inline Index cols() const { return m_cpqr.cols(); } inline Index cols() const { return m_cpqr.cols(); }
@ -513,6 +503,21 @@ void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
} }
#endif #endif
namespace internal {
template<typename DstXprType, typename MatrixType, typename Scalar>
struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
{
typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
typedef Inverse<CodType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows()));
}
};
} // end namespace internal
/** \returns the matrix Q as a sequence of householder transformations */ /** \returns the matrix Q as a sequence of householder transformations */
template <typename MatrixType> template <typename MatrixType>
typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType

3
test/qr_colpivoting.cpp
View File

@ -57,6 +57,9 @@ void cod() {
JacobiSVD<MatrixType> svd(matrix, ComputeThinU | ComputeThinV); JacobiSVD<MatrixType> svd(matrix, ComputeThinU | ComputeThinV);
MatrixType svd_solution = svd.solve(rhs); MatrixType svd_solution = svd.solve(rhs);
VERIFY_IS_APPROX(cod_solution, svd_solution); VERIFY_IS_APPROX(cod_solution, svd_solution);
MatrixType pinv = cod.inverse();
VERIFY_IS_APPROX(cod_solution, pinv * rhs);
} }
template <typename MatrixType, int Cols2> template <typename MatrixType, int Cols2>