GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-05 03:30:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Port QR module to Solve/Inverse

Browse Source
This commit is contained in:
Gael Guennebaud 2014-03-11 11:47:32 +01:00
parent ae40583965
commit 9be72cda2a
3 changed files with 229 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -13,6 +13,15 @@
namespace Eigen {
namespace internal {
template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
: traits<_MatrixType>
{
enum { Flags = 0 };
};
} // end namespace internal
/** \ingroup QR_Module
*
* \class ColPivHouseholderQR
@ -56,6 +65,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
typedef typename MatrixType::PlainObject PlainObject;
private:
@ -137,6 +147,15 @@ template<typename _MatrixType> class ColPivHouseholderQR
* Example: \include ColPivHouseholderQR_solve.cpp
* Output: \verbinclude ColPivHouseholderQR_solve.out
*/
#ifdef EIGEN_TEST_EVALUATORS
template<typename Rhs>
inline const Solve<ColPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return Solve<ColPivHouseholderQR, Rhs>(*this, b.derived());
}
#else
template<typename Rhs>
inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
@ -144,9 +163,10 @@ template<typename _MatrixType> class ColPivHouseholderQR
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
}
#endif
HouseholderSequenceType householderQ(void) const;
HouseholderSequenceType matrixQ(void) const
HouseholderSequenceType householderQ() const;
HouseholderSequenceType matrixQ() const
{
return householderQ();
}
@ -284,6 +304,13 @@ template<typename _MatrixType> class ColPivHouseholderQR
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
*/
#ifdef EIGEN_TEST_EVALUATORS
inline const Inverse<ColPivHouseholderQR> inverse() const
{
eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return Inverse<ColPivHouseholderQR>(*this);
}
#else
inline const
internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
inverse() const
@ -292,6 +319,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
(*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
}
#endif
inline Index rows() const { return m_qr.rows(); }
inline Index cols() const { return m_qr.cols(); }
@ -382,6 +410,12 @@ template<typename _MatrixType> class ColPivHouseholderQR
eigen_assert(m_isInitialized && "Decomposition is not initialized.");
return Success;
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename RhsType, typename DstType>
EIGEN_DEVICE_FUNC
void _solve_impl(const RhsType &rhs, DstType &dst) const;
#endif
protected:
MatrixType m_qr;
@ -514,8 +548,41 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
return *this;
}
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename _MatrixType>
template<typename RhsType, typename DstType>
void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
eigen_assert(rhs.rows() == rows());
const Index nonzero_pivots = nonzeroPivots();
if(nonzero_pivots == 0)
{
dst.setZero();
return;
}
typename RhsType::PlainObject c(rhs);
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(householderSequence(m_qr, m_hCoeffs)
.setLength(nonzero_pivots)
.transpose()
);
m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
.template triangularView<Upper>()
.solveInPlace(c.topRows(nonzero_pivots));
for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
}
#endif
namespace internal {
#ifndef EIGEN_TEST_EVALUATORS
template<typename _MatrixType, typename Rhs>
struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
: solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
@ -524,34 +591,23 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
eigen_assert(rhs().rows() == dec().rows());
const Index cols = dec().cols(),
nonzero_pivots = dec().nonzeroPivots();
if(nonzero_pivots == 0)
{
dst.setZero();
return;
}
typename Rhs::PlainObject c(rhs());
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
.setLength(dec().nonzeroPivots())
.transpose()
);
dec().matrixR()
.topLeftCorner(nonzero_pivots, nonzero_pivots)
.template triangularView<Upper>()
.solveInPlace(c.topRows(nonzero_pivots));
for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
dec()._solve_impl(rhs(), dst);
}
};
#endif
#ifdef EIGEN_TEST_EVALUATORS
template<typename DstXprType, typename MatrixType, typename Scalar>
struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
{
typedef ColPivHouseholderQR<MatrixType> QrType;
typedef Inverse<QrType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
}
};
#endif
} // end namespace internal

135
Eigen/src/QR/FullPivHouseholderQR.h
View File

@ -15,6 +15,12 @@ namespace Eigen {
namespace internal {
template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> >
: traits<_MatrixType>
{
enum { Flags = 0 };
};
template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
template<typename MatrixType>
@ -23,7 +29,7 @@ struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
typedef typename MatrixType::PlainObject ReturnType;
};
}
} // end namespace internal
/** \ingroup QR_Module
*
@ -69,6 +75,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
typedef typename MatrixType::PlainObject PlainObject;
/** \brief Default Constructor.
*
@ -144,6 +151,15 @@ template<typename _MatrixType> class FullPivHouseholderQR
* Example: \include FullPivHouseholderQR_solve.cpp
* Output: \verbinclude FullPivHouseholderQR_solve.out
*/
#ifdef EIGEN_TEST_EVALUATORS
template<typename Rhs>
inline const Solve<FullPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
}
#else
template<typename Rhs>
inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
@ -151,6 +167,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
}
#endif
/** \returns Expression object representing the matrix Q
*/
@ -280,7 +297,15 @@ template<typename _MatrixType> class FullPivHouseholderQR
*
* \note If this matrix is not invertible, the returned matrix has undefined coefficients.
* Use isInvertible() to first determine whether this matrix is invertible.
*/ inline const
*/
#ifdef EIGEN_TEST_EVALUATORS
inline const Inverse<FullPivHouseholderQR> inverse() const
{
eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
return Inverse<FullPivHouseholderQR>(*this);
}
#else
inline const
internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
inverse() const
{
@ -288,6 +313,7 @@ template<typename _MatrixType> class FullPivHouseholderQR
return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
(*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
}
#endif
inline Index rows() const { return m_qr.rows(); }
inline Index cols() const { return m_qr.cols(); }
@ -366,6 +392,12 @@ template<typename _MatrixType> class FullPivHouseholderQR
* diagonal coefficient of U.
*/
RealScalar maxPivot() const { return m_maxpivot; }
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename RhsType, typename DstType>
EIGEN_DEVICE_FUNC
void _solve_impl(const RhsType &rhs, DstType &dst) const;
#endif
protected:
MatrixType m_qr;
@ -485,8 +517,46 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
return *this;
}
namespace internal {
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename _MatrixType>
template<typename RhsType, typename DstType>
void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
eigen_assert(rhs.rows() == rows());
const Index l_rank = rank();
// FIXME introduce nonzeroPivots() and use it here. and more generally,
// make the same improvements in this dec as in FullPivLU.
if(l_rank==0)
{
dst.setZero();
return;
}
typename RhsType::PlainObject c(rhs);
Matrix<Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
for (Index k = 0; k < l_rank; ++k)
{
Index remainingSize = rows()-k;
c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
c.bottomRightCorner(remainingSize, rhs.cols())
.applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
m_hCoeffs.coeff(k), &temp.coeffRef(0));
}
m_qr.topLeftCorner(l_rank, l_rank)
.template triangularView<Upper>()
.solveInPlace(c.topRows(l_rank));
for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
}
#endif
namespace internal {
#ifndef EIGEN_TEST_EVALUATORS
template<typename _MatrixType, typename Rhs>
struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
: solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
@ -495,38 +565,23 @@ struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
const Index rows = dec().rows(), cols = dec().cols();
eigen_assert(rhs().rows() == rows);
// FIXME introduce nonzeroPivots() and use it here. and more generally,
// make the same improvements in this dec as in FullPivLU.
if(dec().rank()==0)
{
dst.setZero();
return;
}
typename Rhs::PlainObject c(rhs());
Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
for (Index k = 0; k < dec().rank(); ++k)
{
Index remainingSize = rows-k;
c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
c.bottomRightCorner(remainingSize, rhs().cols())
.applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
dec().hCoeffs().coeff(k), &temp.coeffRef(0));
}
dec().matrixQR()
.topLeftCorner(dec().rank(), dec().rank())
.template triangularView<Upper>()
.solveInPlace(c.topRows(dec().rank()));
for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
dec()._solve_impl(rhs(), dst);
}
};
#endif // EIGEN_TEST_EVALUATORS
#ifdef EIGEN_TEST_EVALUATORS
template<typename DstXprType, typename MatrixType, typename Scalar>
struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
{
typedef FullPivHouseholderQR<MatrixType> QrType;
typedef Inverse<QrType> SrcXprType;
static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
{
dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
}
};
#endif
/** \ingroup QR_Module
*
@ -534,6 +589,7 @@ struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
*
* \tparam MatrixType type of underlying dense matrix
*/
// #ifndef EIGEN_TEST_EVALUATORS
template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
: public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
{
@ -550,7 +606,7 @@ public:
: m_qr(qr),
m_hCoeffs(hCoeffs),
m_rowsTranspositions(rowsTranspositions)
{}
{}
template <typename ResultType>
void evalTo(ResultType& result) const
@ -580,8 +636,8 @@ public:
}
}
Index rows() const { return m_qr.rows(); }
Index cols() const { return m_qr.rows(); }
Index rows() const { return m_qr.rows(); }
Index cols() const { return m_qr.rows(); }
protected:
typename MatrixType::Nested m_qr;
@ -589,6 +645,13 @@ protected:
typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
};
// #ifdef EIGEN_TEST_EVALUATORS
// template<typename MatrixType>
// struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
// : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
// {};
// #endif
} // end namespace internal
template<typename MatrixType>

65
Eigen/src/QR/HouseholderQR.h
View File

@ -117,6 +117,15 @@ template<typename _MatrixType> class HouseholderQR
* Example: \include HouseholderQR_solve.cpp
* Output: \verbinclude HouseholderQR_solve.out
*/
#ifdef EIGEN_TEST_EVALUATORS
template<typename Rhs>
inline const Solve<HouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
return Solve<HouseholderQR, Rhs>(*this, b.derived());
}
#else
template<typename Rhs>
inline const internal::solve_retval<HouseholderQR, Rhs>
solve(const MatrixBase<Rhs>& b) const
@ -124,6 +133,7 @@ template<typename _MatrixType> class HouseholderQR
eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
}
#endif
/** This method returns an expression of the unitary matrix Q as a sequence of Householder transformations.
*
@ -187,6 +197,12 @@ template<typename _MatrixType> class HouseholderQR
* For advanced uses only.
*/
const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename RhsType, typename DstType>
EIGEN_DEVICE_FUNC
void _solve_impl(const RhsType &rhs, DstType &dst) const;
#endif
protected:
MatrixType m_qr;
@ -308,6 +324,35 @@ struct householder_qr_inplace_blocked
}
};
} // end namespace internal
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename _MatrixType>
template<typename RhsType, typename DstType>
void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
{
const Index rank = (std::min)(rows(), cols());
eigen_assert(rhs.rows() == rows());
typename RhsType::PlainObject c(rhs);
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(householderSequence(
m_qr.leftCols(rank),
m_hCoeffs.head(rank)).transpose()
);
m_qr.topLeftCorner(rank, rank)
.template triangularView<Upper>()
.solveInPlace(c.topRows(rank));
dst.topRows(rank) = c.topRows(rank);
dst.bottomRows(cols()-rank).setZero();
}
#endif
namespace internal {
template<typename _MatrixType, typename Rhs>
struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
: solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
@ -316,25 +361,7 @@ struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
template<typename Dest> void evalTo(Dest& dst) const
{
const Index rows = dec().rows(), cols = dec().cols();
const Index rank = (std::min)(rows, cols);
eigen_assert(rhs().rows() == rows);
typename Rhs::PlainObject c(rhs());
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(householderSequence(
dec().matrixQR().leftCols(rank),
dec().hCoeffs().head(rank)).transpose()
);
dec().matrixQR()
.topLeftCorner(rank, rank)
.template triangularView<Upper>()
.solveInPlace(c.topRows(rank));
dst.topRows(rank) = c.topRows(rank);
dst.bottomRows(cols-rank).setZero();
dec()._solve_impl(rhs(), dst);
}
};