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

* HouseholderSequence:

Browse Source 
* be aware of number of actual householder vectors
    (optimization in non-full-rank case, no behavior change)
  * fix applyThisOnTheRight, it was using k instead of actual_k
* QR: rename matrixQ() to householderQ() where applicable
This commit is contained in:
Benoit Jacob 2009-12-02 11:11:09 -05:00
parent 3279e39340
commit 3e73f6036c
7 changed files with 43 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
Eigen/src/Householder/HouseholderSequence.h
View File

@ -69,31 +69,35 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
typedef HouseholderSequence<VectorsType,
typename ei_meta_if<NumTraits<Scalar>::IsComplex,
typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type,
NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >,
CoeffsType>::ret> ConjugateReturnType;
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
: m_vectors(v), m_coeffs(h), m_trans(trans)
: m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize())
{}
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
: m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors)
{}
int rows() const { return m_vectors.rows(); }
int cols() const { return m_vectors.rows(); }
HouseholderSequence transpose() const
{ return HouseholderSequence(m_vectors, m_coeffs, !m_trans); }
{ return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors); }
ConjugateReturnType conjugate() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors); }
ConjugateReturnType adjoint() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); }
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors); }
ConjugateReturnType inverse() const { return adjoint(); }
/** \internal */
template<typename DestType> void evalTo(DestType& dst) const
{
int vecs = std::min(m_vectors.cols(),m_vectors.rows());
int vecs = m_actualVectors;
int length = m_vectors.rows();
dst.setIdentity();
Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows());
@ -111,22 +115,22 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
/** \internal */
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
{
int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
int length = m_vectors.rows(); // size of the largest householder vector
Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows());
int vecs = m_actualVectors; // number of householder vectors
int length = m_vectors.rows(); // size of the largest householder vector
Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows());
for(int k = 0; k < vecs; ++k)
{
int actual_k = m_trans ? vecs-k-1 : k;
dst.corner(BottomRight, dst.rows(), length-k)
.applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
dst.corner(BottomRight, dst.rows(), length-actual_k)
.applyHouseholderOnTheRight(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
}
}
/** \internal */
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
{
int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
int length = m_vectors.rows(); // size of the largest householder vector
int vecs = m_actualVectors; // number of householder vectors
int length = m_vectors.rows(); // size of the largest householder vector
Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
for(int k = 0; k < vecs; ++k)
{
@ -156,6 +160,7 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
typename VectorsType::Nested m_vectors;
typename CoeffsType::Nested m_coeffs;
bool m_trans;
int m_actualVectors;
private:
HouseholderSequence& operator=(const HouseholderSequence&);
@ -167,4 +172,10 @@ HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsTyp
return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans);
}
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
{
return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors);
}
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H

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

@ -105,7 +105,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
return ei_solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
}
HouseholderSequenceType matrixQ(void) const;
HouseholderSequenceType householderQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored
*/
@ -447,7 +447,8 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
c.applyOnTheLeft(householderSequence(
dec().matrixQR(),
dec().hCoeffs(),
true
true,
dec().nonzeroPivots()
));
dec().matrixQR()
@ -470,10 +471,11 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
/** \returns the matrix Q as a sequence of householder transformations */
template<typename MatrixType>
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>::matrixQ() const
typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
::householderQ() const
{
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false);
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots);
}
#endif // EIGEN_HIDE_HEAVY_CODE

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

@ -104,11 +104,10 @@ template<typename _MatrixType> class HouseholderQR
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
return ei_solve_retval<HouseholderQR, Rhs>(*this, b.derived());
}
MatrixQType matrixQ() const;
HouseholderSequenceType matrixQAsHouseholderSequence() const
HouseholderSequenceType householderQ() const
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
}
@ -240,14 +239,6 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs>
}
};
/** \returns the matrix Q */
template<typename MatrixType>
typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
return matrixQAsHouseholderSequence();
}
#endif // EIGEN_HIDE_HEAVY_CODE
/** \return the Householder QR decomposition of \c *this.

2
test/geo_hyperplane.cpp
View File

@ -64,7 +64,7 @@ template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane)
// transform
if (!NumTraits<Scalar>::IsComplex)
{
MatrixType rot = MatrixType::Random(dim,dim).householderQr().matrixQ();
MatrixType rot = MatrixType::Random(dim,dim).householderQr().householderQ();
DiagonalMatrix<Scalar,HyperplaneType::AmbientDimAtCompileTime> scaling(VectorType::Random());
Translation<Scalar,HyperplaneType::AmbientDimAtCompileTime> translation(VectorType::Random());

2
test/main.h
View File

@ -360,7 +360,7 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType&
HouseholderQR<MatrixAType> qra(a);
HouseholderQR<MatrixBType> qrb(b);
m = qra.matrixQ() * d * qrb.matrixQ();
m = qra.householderQ() * d * qrb.householderQ();
}
} // end namespace Eigen

10
test/qr.cpp
View File

@ -38,13 +38,13 @@ template<typename MatrixType> void qr(const MatrixType& m)
HouseholderQR<MatrixType> qrOfA(a);
MatrixType r = qrOfA.matrixQR();
MatrixQType q = qrOfA.matrixQ();
MatrixQType q = qrOfA.householderQ();
VERIFY_IS_UNITARY(q);
// FIXME need better way to construct trapezoid
for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
VERIFY_IS_APPROX(a, qrOfA.matrixQ() * r);
VERIFY_IS_APPROX(a, qrOfA.householderQ() * r);
}
template<typename MatrixType, int Cols2> void qr_fixedsize()
@ -58,7 +58,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
// FIXME need better way to construct trapezoid
for(int i = 0; i < Rows; i++) for(int j = 0; j < Cols; j++) if(i>j) r(i,j) = Scalar(0);
VERIFY_IS_APPROX(m1, qr.matrixQ() * r);
VERIFY_IS_APPROX(m1, qr.householderQ() * r);
Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
@ -93,7 +93,7 @@ template<typename MatrixType> void qr_invertible()
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
RealScalar absdet = ei_abs(m1.diagonal().prod());
m3 = qr.matrixQ(); // get a unitary
m3 = qr.householderQ(); // get a unitary
m1 = m3 * m1 * m3;
qr.compute(m1);
VERIFY_IS_APPROX(absdet, qr.absDeterminant());
@ -107,7 +107,7 @@ template<typename MatrixType> void qr_verify_assert()
HouseholderQR<MatrixType> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.matrixQ())
VERIFY_RAISES_ASSERT(qr.householderQ())
VERIFY_RAISES_ASSERT(qr.absDeterminant())
VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
}

8
test/qr_colpivoting.cpp
View File

@ -44,7 +44,7 @@ template<typename MatrixType> void qr()
VERIFY(!qr.isInvertible());
VERIFY(!qr.isSurjective());
MatrixQType q = qr.matrixQ();
MatrixQType q = qr.householderQ();
VERIFY_IS_UNITARY(q);
MatrixType r = qr.matrixQR().template triangularView<UpperTriangular>();
@ -73,7 +73,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
VERIFY(!qr.isSurjective());
Matrix<Scalar,Rows,Cols> r = qr.matrixQR().template triangularView<UpperTriangular>();
Matrix<Scalar,Rows,Cols> c = qr.matrixQ() * r * qr.colsPermutation().inverse();
Matrix<Scalar,Rows,Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse();
VERIFY_IS_APPROX(m1, c);
Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
@ -109,7 +109,7 @@ template<typename MatrixType> void qr_invertible()
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
RealScalar absdet = ei_abs(m1.diagonal().prod());
m3 = qr.matrixQ(); // get a unitary
m3 = qr.householderQ(); // get a unitary
m1 = m3 * m1 * m3;
qr.compute(m1);
VERIFY_IS_APPROX(absdet, qr.absDeterminant());
@ -123,7 +123,7 @@ template<typename MatrixType> void qr_verify_assert()
ColPivHouseholderQR<MatrixType> qr;
VERIFY_RAISES_ASSERT(qr.matrixQR())
VERIFY_RAISES_ASSERT(qr.solve(tmp))
VERIFY_RAISES_ASSERT(qr.matrixQ())
VERIFY_RAISES_ASSERT(qr.householderQ())
VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
VERIFY_RAISES_ASSERT(qr.isInjective())
VERIFY_RAISES_ASSERT(qr.isSurjective())