mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-02-17 18:09:55 +08:00
* rename QR to HouseholderQR because really that impacts the API, not just the impl.
* rename qr() to householderQr(), for same reason. * clarify that it's non-pivoting, non-rank-revealing, so remove all the rank API, make solve() be void instead of bool, update the docs/test, etc. * fix warning in SVD
This commit is contained in:
parent
0c2232e5d9
commit
e093b43b2c
@ -653,7 +653,7 @@ template<typename Derived> class MatrixBase
|
||||
|
||||
/////////// QR module ///////////
|
||||
|
||||
const QR<PlainMatrixType> qr() const;
|
||||
const HouseholderQR<PlainMatrixType> householderQr() const;
|
||||
|
||||
EigenvaluesReturnType eigenvalues() const;
|
||||
RealScalar operatorNorm() const;
|
||||
|
@ -112,7 +112,7 @@ template<typename MatrixType, int Direction = BothDirections> class Reverse;
|
||||
|
||||
template<typename MatrixType> class LU;
|
||||
template<typename MatrixType> class PartialLU;
|
||||
template<typename MatrixType> class QR;
|
||||
template<typename MatrixType> class HouseholderQR;
|
||||
template<typename MatrixType> class SVD;
|
||||
template<typename MatrixType> class LLT;
|
||||
template<typename MatrixType> class LDLT;
|
||||
|
@ -22,24 +22,26 @@
|
||||
// License and a copy of the GNU General Public License along with
|
||||
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
#ifndef EIGEN_QR_H
|
||||
#define EIGEN_QR_H
|
||||
#ifndef EIGEN_HouseholderQR_H
|
||||
#define EIGEN_HouseholderQR_H
|
||||
|
||||
/** \ingroup QR_Module
|
||||
/** \ingroup HouseholderQR_Module
|
||||
* \nonstableyet
|
||||
*
|
||||
* \class QR
|
||||
* \class HouseholderQR
|
||||
*
|
||||
* \brief QR decomposition of a matrix
|
||||
* \brief Householder QR decomposition of a matrix
|
||||
*
|
||||
* \param MatrixType the type of the matrix of which we are computing the QR decomposition
|
||||
*
|
||||
* This class performs a QR decomposition using Householder transformations. The result is
|
||||
* stored in a compact way compatible with LAPACK.
|
||||
*
|
||||
* Note that no pivoting is performed. This is \b not a rank-revealing decomposition.
|
||||
*
|
||||
* \sa MatrixBase::qr()
|
||||
*/
|
||||
template<typename MatrixType> class QR
|
||||
template<typename MatrixType> class HouseholderQR
|
||||
{
|
||||
public:
|
||||
|
||||
@ -53,88 +55,23 @@ template<typename MatrixType> class QR
|
||||
* \brief Default Constructor.
|
||||
*
|
||||
* The default constructor is useful in cases in which the user intends to
|
||||
* perform decompositions via QR::compute(const MatrixType&).
|
||||
* perform decompositions via HouseholderQR::compute(const MatrixType&).
|
||||
*/
|
||||
QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
|
||||
HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
|
||||
|
||||
QR(const MatrixType& matrix)
|
||||
HouseholderQR(const MatrixType& matrix)
|
||||
: m_qr(matrix.rows(), matrix.cols()),
|
||||
m_hCoeffs(matrix.cols()),
|
||||
m_isInitialized(false)
|
||||
{
|
||||
compute(matrix);
|
||||
}
|
||||
|
||||
/** \deprecated use isInjective()
|
||||
* \returns whether or not the matrix is of full rank
|
||||
*
|
||||
* \note Since the rank is computed only once, i.e. the first time it is needed, this
|
||||
* method almost does not perform any further computation.
|
||||
*/
|
||||
EIGEN_DEPRECATED bool isFullRank() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
return rank() == m_qr.cols();
|
||||
}
|
||||
|
||||
/** \returns the rank of the matrix of which *this is the QR decomposition.
|
||||
*
|
||||
* \note Since the rank is computed only once, i.e. the first time it is needed, this
|
||||
* method almost does not perform any further computation.
|
||||
*/
|
||||
int rank() const;
|
||||
|
||||
/** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
|
||||
*
|
||||
* \note Since the rank is computed only once, i.e. the first time it is needed, this
|
||||
* method almost does not perform any further computation.
|
||||
*/
|
||||
inline int dimensionOfKernel() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
return m_qr.cols() - rank();
|
||||
}
|
||||
|
||||
/** \returns true if the matrix of which *this is the QR decomposition represents an injective
|
||||
* linear map, i.e. has trivial kernel; false otherwise.
|
||||
*
|
||||
* \note Since the rank is computed only once, i.e. the first time it is needed, this
|
||||
* method almost does not perform any further computation.
|
||||
*/
|
||||
inline bool isInjective() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
return rank() == m_qr.cols();
|
||||
}
|
||||
|
||||
/** \returns true if the matrix of which *this is the QR decomposition represents a surjective
|
||||
* linear map; false otherwise.
|
||||
*
|
||||
* \note Since the rank is computed only once, i.e. the first time it is needed, this
|
||||
* method almost does not perform any further computation.
|
||||
*/
|
||||
inline bool isSurjective() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
return rank() == m_qr.rows();
|
||||
}
|
||||
|
||||
/** \returns true if the matrix of which *this is the QR decomposition is invertible.
|
||||
*
|
||||
* \note Since the rank is computed only once, i.e. the first time it is needed, this
|
||||
* method almost does not perform any further computation.
|
||||
*/
|
||||
inline bool isInvertible() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
return isInjective() && isSurjective();
|
||||
}
|
||||
|
||||
|
||||
/** \returns a read-only expression of the matrix R of the actual the QR decomposition */
|
||||
const Part<NestByValue<MatrixRBlockType>, UpperTriangular>
|
||||
matrixR(void) const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
|
||||
int cols = m_qr.cols();
|
||||
return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>();
|
||||
}
|
||||
@ -148,22 +85,14 @@ template<typename MatrixType> class QR
|
||||
* 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().
|
||||
*
|
||||
* \note The case where b is a matrix is not yet implemented. Also, this
|
||||
* code is space inefficient.
|
||||
*
|
||||
* Example: \include QR_solve.cpp
|
||||
* Output: \verbinclude QR_solve.out
|
||||
*
|
||||
* \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()
|
||||
* Example: \include HouseholderQR_solve.cpp
|
||||
* Output: \verbinclude HouseholderQR_solve.out
|
||||
*/
|
||||
template<typename OtherDerived, typename ResultType>
|
||||
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
|
||||
void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
|
||||
|
||||
MatrixType matrixQ(void) const;
|
||||
|
||||
@ -172,34 +101,14 @@ template<typename MatrixType> class QR
|
||||
protected:
|
||||
MatrixType m_qr;
|
||||
VectorType m_hCoeffs;
|
||||
mutable int m_rank;
|
||||
mutable bool m_rankIsUptodate;
|
||||
bool m_isInitialized;
|
||||
};
|
||||
|
||||
/** \returns the rank of the matrix of which *this is the QR decomposition. */
|
||||
template<typename MatrixType>
|
||||
int QR<MatrixType>::rank() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
if (!m_rankIsUptodate)
|
||||
{
|
||||
RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff();
|
||||
int n = m_qr.cols();
|
||||
m_rank = 0;
|
||||
while(m_rank<n && !ei_isMuchSmallerThan(m_qr.diagonal().coeff(m_rank), maxCoeff))
|
||||
++m_rank;
|
||||
m_rankIsUptodate = true;
|
||||
}
|
||||
return m_rank;
|
||||
}
|
||||
|
||||
#ifndef EIGEN_HIDE_HEAVY_CODE
|
||||
|
||||
template<typename MatrixType>
|
||||
void QR<MatrixType>::compute(const MatrixType& matrix)
|
||||
void HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
|
||||
{
|
||||
m_rankIsUptodate = false;
|
||||
m_qr = matrix;
|
||||
m_hCoeffs.resize(matrix.cols());
|
||||
|
||||
@ -262,12 +171,12 @@ void QR<MatrixType>::compute(const MatrixType& matrix)
|
||||
|
||||
template<typename MatrixType>
|
||||
template<typename OtherDerived, typename ResultType>
|
||||
bool QR<MatrixType>::solve(
|
||||
void HouseholderQR<MatrixType>::solve(
|
||||
const MatrixBase<OtherDerived>& b,
|
||||
ResultType *result
|
||||
) const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
|
||||
const int rows = m_qr.rows();
|
||||
ei_assert(b.rows() == rows);
|
||||
result->resize(rows, b.cols());
|
||||
@ -276,27 +185,17 @@ bool QR<MatrixType>::solve(
|
||||
// Q^T without explicitly forming matrixQ(). Investigate.
|
||||
*result = matrixQ().transpose()*b;
|
||||
|
||||
if(!isSurjective())
|
||||
{
|
||||
// is result is in the image of R ?
|
||||
RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff();
|
||||
for(int col = 0; col < result->cols(); ++col)
|
||||
for(int row = m_rank; row < result->rows(); ++row)
|
||||
if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res))
|
||||
return false;
|
||||
}
|
||||
m_qr.corner(TopLeft, m_rank, m_rank)
|
||||
const int rank = std::min(result->rows(), result->cols());
|
||||
m_qr.corner(TopLeft, rank, rank)
|
||||
.template marked<UpperTriangular>()
|
||||
.solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols()));
|
||||
|
||||
return true;
|
||||
.solveTriangularInPlace(result->corner(TopLeft, rank, result->cols()));
|
||||
}
|
||||
|
||||
/** \returns the matrix Q */
|
||||
template<typename MatrixType>
|
||||
MatrixType QR<MatrixType>::matrixQ() const
|
||||
MatrixType HouseholderQR<MatrixType>::matrixQ() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "QR is not initialized.");
|
||||
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
|
||||
// compute the product Q_0 Q_1 ... Q_n-1,
|
||||
// where Q_k is the k-th Householder transformation I - h_k v_k v_k'
|
||||
// and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
|
||||
@ -319,16 +218,16 @@ MatrixType QR<MatrixType>::matrixQ() const
|
||||
|
||||
#endif // EIGEN_HIDE_HEAVY_CODE
|
||||
|
||||
/** \return the QR decomposition of \c *this.
|
||||
/** \return the Householder QR decomposition of \c *this.
|
||||
*
|
||||
* \sa class QR
|
||||
* \sa class HouseholderQR
|
||||
*/
|
||||
template<typename Derived>
|
||||
const QR<typename MatrixBase<Derived>::PlainMatrixType>
|
||||
MatrixBase<Derived>::qr() const
|
||||
const HouseholderQR<typename MatrixBase<Derived>::PlainMatrixType>
|
||||
MatrixBase<Derived>::householderQr() const
|
||||
{
|
||||
return QR<PlainMatrixType>(eval());
|
||||
return HouseholderQR<PlainMatrixType>(eval());
|
||||
}
|
||||
|
||||
|
||||
#endif // EIGEN_QR_H
|
||||
#endif // EIGEN_HouseholderQR_H
|
||||
|
@ -145,7 +145,6 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
|
||||
{
|
||||
const int m = matrix.rows();
|
||||
const int n = matrix.cols();
|
||||
const int nu = std::min(m,n);
|
||||
|
||||
m_matU.resize(m, m);
|
||||
m_matU.setZero();
|
||||
|
@ -6,5 +6,5 @@
|
||||
</head><body>
|
||||
<a name="top"></a>
|
||||
<a class="logo" href="http://eigen.tuxfamily.org/">
|
||||
<img src="Eigen_Silly_Professor_64x64.png" width=64 height=64 alt="Eiegn's silly professor"
|
||||
<img src="Eigen_Silly_Professor_64x64.png" width=64 height=64 alt="Eigen's silly professor"
|
||||
style="position:absolute; border:none" /></a>
|
@ -4,11 +4,6 @@ Matrix3f y = Matrix3f::Random();
|
||||
cout << "Here is the matrix m:" << endl << m << endl;
|
||||
cout << "Here is the matrix y:" << endl << y << endl;
|
||||
Matrix3f x;
|
||||
if(m.qr().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;
|
||||
|
||||
m.householderQr().solve(y, &x));
|
||||
assert(y.isApprox(m*x));
|
||||
cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
|
@ -64,7 +64,7 @@ template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane)
|
||||
// transform
|
||||
if (!NumTraits<Scalar>::IsComplex)
|
||||
{
|
||||
MatrixType rot = MatrixType::Random(dim,dim).qr().matrixQ();
|
||||
MatrixType rot = MatrixType::Random(dim,dim).householderQr().matrixQ();
|
||||
DiagonalMatrix<Scalar,HyperplaneType::AmbientDimAtCompileTime> scaling(VectorType::Random());
|
||||
Translation<Scalar,HyperplaneType::AmbientDimAtCompileTime> translation(VectorType::Random());
|
||||
|
||||
|
@ -241,8 +241,8 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri
|
||||
const int diag_size = std::min(d.rows(),d.cols());
|
||||
d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank);
|
||||
|
||||
QR<MatrixType> qra(a);
|
||||
QR<MatrixType> qrb(b);
|
||||
HouseholderQR<MatrixType> qra(a);
|
||||
HouseholderQR<MatrixType> qrb(b);
|
||||
m = (qra.matrixQ() * d * qrb.matrixQ()).lazy();
|
||||
}
|
||||
|
||||
|
62
test/qr.cpp
62
test/qr.cpp
@ -36,7 +36,7 @@ template<typename MatrixType> void qr(const MatrixType& m)
|
||||
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
|
||||
|
||||
MatrixType a = MatrixType::Random(rows,cols);
|
||||
QR<MatrixType> qrOfA(a);
|
||||
HouseholderQR<MatrixType> qrOfA(a);
|
||||
VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR());
|
||||
VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR());
|
||||
|
||||
@ -55,40 +55,6 @@ template<typename MatrixType> void qr(const MatrixType& m)
|
||||
VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint());
|
||||
}
|
||||
|
||||
template<typename MatrixType> void qr_non_invertible()
|
||||
{
|
||||
/* this test covers the following files: QR.h */
|
||||
int rows = ei_random<int>(20,200), cols = ei_random<int>(20,rows), cols2 = ei_random<int>(20,rows);
|
||||
int rank = ei_random<int>(1, std::min(rows, cols)-1);
|
||||
|
||||
MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1);
|
||||
createRandomMatrixOfRank(rank, rows, cols, m1);
|
||||
|
||||
QR<MatrixType> lu(m1);
|
||||
// typename LU<MatrixType>::KernelResultType m1kernel = lu.kernel();
|
||||
// typename LU<MatrixType>::ImageResultType m1image = lu.image();
|
||||
std::cerr << rows << "x" << cols << " " << rank << " " << lu.rank() << "\n";
|
||||
if (rank != lu.rank())
|
||||
std::cerr << lu.matrixR().diagonal().transpose() << "\n";
|
||||
VERIFY(rank == lu.rank());
|
||||
VERIFY(cols - lu.rank() == lu.dimensionOfKernel());
|
||||
VERIFY(!lu.isInjective());
|
||||
VERIFY(!lu.isInvertible());
|
||||
VERIFY(lu.isSurjective() == (lu.rank() == rows));
|
||||
// VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
|
||||
// VERIFY(m1image.lu().rank() == rank);
|
||||
// MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
|
||||
// sidebyside << m1, m1image;
|
||||
// VERIFY(sidebyside.lu().rank() == rank);
|
||||
m2 = MatrixType::Random(cols,cols2);
|
||||
m3 = m1*m2;
|
||||
m2 = MatrixType::Random(cols,cols2);
|
||||
lu.solve(m3, &m2);
|
||||
VERIFY_IS_APPROX(m3, m1*m2);
|
||||
m3 = MatrixType::Random(rows,cols2);
|
||||
VERIFY(!lu.solve(m3, &m2));
|
||||
}
|
||||
|
||||
template<typename MatrixType> void qr_invertible()
|
||||
{
|
||||
/* this test covers the following files: QR.h */
|
||||
@ -105,33 +71,18 @@ template<typename MatrixType> void qr_invertible()
|
||||
m1 += a * a.adjoint();
|
||||
}
|
||||
|
||||
QR<MatrixType> lu(m1);
|
||||
VERIFY(0 == lu.dimensionOfKernel());
|
||||
VERIFY(size == lu.rank());
|
||||
VERIFY(lu.isInjective());
|
||||
VERIFY(lu.isSurjective());
|
||||
VERIFY(lu.isInvertible());
|
||||
// VERIFY(lu.image().lu().isInvertible());
|
||||
HouseholderQR<MatrixType> qr(m1);
|
||||
m3 = MatrixType::Random(size,size);
|
||||
lu.solve(m3, &m2);
|
||||
qr.solve(m3, &m2);
|
||||
//std::cerr << m3 - m1*m2 << "\n\n";
|
||||
VERIFY_IS_APPROX(m3, m1*m2);
|
||||
// VERIFY_IS_APPROX(m2, lu.inverse()*m3);
|
||||
m3 = MatrixType::Random(size,size);
|
||||
VERIFY(lu.solve(m3, &m2));
|
||||
}
|
||||
|
||||
template<typename MatrixType> void qr_verify_assert()
|
||||
{
|
||||
MatrixType tmp;
|
||||
|
||||
QR<MatrixType> qr;
|
||||
VERIFY_RAISES_ASSERT(qr.isFullRank())
|
||||
VERIFY_RAISES_ASSERT(qr.rank())
|
||||
VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
|
||||
VERIFY_RAISES_ASSERT(qr.isInjective())
|
||||
VERIFY_RAISES_ASSERT(qr.isSurjective())
|
||||
VERIFY_RAISES_ASSERT(qr.isInvertible())
|
||||
HouseholderQR<MatrixType> qr;
|
||||
VERIFY_RAISES_ASSERT(qr.matrixR())
|
||||
VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
|
||||
VERIFY_RAISES_ASSERT(qr.matrixQ())
|
||||
@ -149,11 +100,6 @@ void test_qr()
|
||||
}
|
||||
|
||||
for(int i = 0; i < g_repeat; i++) {
|
||||
CALL_SUBTEST( qr_non_invertible<MatrixXf>() );
|
||||
CALL_SUBTEST( qr_non_invertible<MatrixXd>() );
|
||||
// TODO fix issue with complex
|
||||
// CALL_SUBTEST( qr_non_invertible<MatrixXcf>() );
|
||||
// CALL_SUBTEST( qr_non_invertible<MatrixXcd>() );
|
||||
CALL_SUBTEST( qr_invertible<MatrixXf>() );
|
||||
CALL_SUBTEST( qr_invertible<MatrixXd>() );
|
||||
// TODO fix issue with complex
|
||||
|
Loading…
Reference in New Issue
Block a user