mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-08-11 11:19:02 +08:00
* Fix a bug in HouseholderQR with mixed fixed/dynamic size: must use EIGEN_SIZE_MIN instead of EIGEN_ENUM_MIN, and there are many other occurences throughout Eigen!
* HouseholderSequence: - add shift parameter - add essentialVector() method to start abstracting the direction - add unit test in householder.cpp
This commit is contained in:
Benoit Jacob
parent
325da2ea3c
commit
24a09ceae8
@ -65,32 +65,44 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
|
||||
: public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> >
|
||||
{
|
||||
typedef typename VectorsType::Scalar Scalar;
|
||||
typedef Block<VectorsType, Dynamic, 1> EssentialVectorType;
|
||||
|
||||
public:
|
||||
|
||||
typedef HouseholderSequence<VectorsType,
|
||||
typename ei_meta_if<NumTraits<Scalar>::IsComplex,
|
||||
NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >,
|
||||
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_actualVectors(v.diagonalSize())
|
||||
{}
|
||||
: m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()),
|
||||
m_shift(0)
|
||||
{
|
||||
}
|
||||
|
||||
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
|
||||
: m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors)
|
||||
{}
|
||||
HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift)
|
||||
: m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors), m_shift(shift)
|
||||
{
|
||||
}
|
||||
|
||||
int rows() const { return m_vectors.rows(); }
|
||||
int cols() const { return m_vectors.rows(); }
|
||||
|
||||
const EssentialVectorType essentialVector(int k) const
|
||||
{
|
||||
ei_assert(k>= 0 && k < m_actualVectors);
|
||||
const int start = k+1+m_shift;
|
||||
return Block<VectorsType,Dynamic,1>(m_vectors, start, k, rows()-start, 1);
|
||||
}
|
||||
|
||||
HouseholderSequence transpose() const
|
||||
{ return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors); }
|
||||
{ return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors, m_shift); }
|
||||
|
||||
ConjugateReturnType conjugate() const
|
||||
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors); }
|
||||
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors, m_shift); }
|
||||
|
||||
ConjugateReturnType adjoint() const
|
||||
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors); }
|
||||
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors, m_shift); }
|
||||
|
||||
ConjugateReturnType inverse() const { return adjoint(); }
|
||||
|
||||
@ -98,45 +110,41 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
|
||||
template<typename DestType> void evalTo(DestType& dst) const
|
||||
{
|
||||
int vecs = m_actualVectors;
|
||||
int length = m_vectors.rows();
|
||||
dst.setIdentity();
|
||||
Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows());
|
||||
dst.setIdentity(rows(), rows());
|
||||
Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(rows());
|
||||
for(int k = vecs-1; k >= 0; --k)
|
||||
{
|
||||
int cornerSize = rows() - k - m_shift;
|
||||
if(m_trans)
|
||||
dst.corner(BottomRight, length-k, length-k)
|
||||
.applyHouseholderOnTheRight(m_vectors.col(k).tail(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
|
||||
dst.corner(BottomRight, cornerSize, cornerSize)
|
||||
.applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
|
||||
else
|
||||
dst.corner(BottomRight, length-k, length-k)
|
||||
.applyHouseholderOnTheLeft(m_vectors.col(k).tail(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(k));
|
||||
dst.corner(BottomRight, cornerSize, cornerSize)
|
||||
.applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &temp.coeffRef(0));
|
||||
}
|
||||
}
|
||||
|
||||
/** \internal */
|
||||
template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
|
||||
{
|
||||
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)
|
||||
for(int k = 0; k < m_actualVectors; ++k)
|
||||
{
|
||||
int actual_k = m_trans ? vecs-k-1 : k;
|
||||
dst.corner(BottomRight, dst.rows(), length-actual_k)
|
||||
.applyHouseholderOnTheRight(m_vectors.col(actual_k).tail(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
|
||||
int actual_k = m_trans ? m_actualVectors-k-1 : k;
|
||||
dst.corner(BottomRight, dst.rows(), rows()-m_shift-actual_k)
|
||||
.applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
|
||||
}
|
||||
}
|
||||
|
||||
/** \internal */
|
||||
template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
|
||||
{
|
||||
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)
|
||||
for(int k = 0; k < m_actualVectors; ++k)
|
||||
{
|
||||
int actual_k = m_trans ? k : vecs-k-1;
|
||||
dst.corner(BottomRight, length-actual_k, dst.cols())
|
||||
.applyHouseholderOnTheLeft(m_vectors.col(actual_k).tail(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
|
||||
int actual_k = m_trans ? k : m_actualVectors-k-1;
|
||||
dst.corner(BottomRight, rows()-m_shift-actual_k, dst.cols())
|
||||
.applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
|
||||
}
|
||||
}
|
||||
|
||||
@ -161,6 +169,7 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
|
||||
typename CoeffsType::Nested m_coeffs;
|
||||
bool m_trans;
|
||||
int m_actualVectors;
|
||||
int m_shift;
|
||||
};
|
||||
|
||||
template<typename VectorsType, typename CoeffsType>
|
||||
@ -170,9 +179,9 @@ HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsTyp
|
||||
}
|
||||
|
||||
template<typename VectorsType, typename CoeffsType>
|
||||
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
|
||||
HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift)
|
||||
{
|
||||
return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors);
|
||||
return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors, shift);
|
||||
}
|
||||
|
||||
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
|
||||
|
||||
@ -448,7 +448,8 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
|
||||
dec().matrixQR(),
|
||||
dec().hCoeffs(),
|
||||
true,
|
||||
dec().nonzeroPivots()
|
||||
dec().nonzeroPivots(),
|
||||
0
|
||||
));
|
||||
|
||||
dec().matrixQR()
|
||||
@ -475,7 +476,7 @@ typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHousehol
|
||||
::householderQ() const
|
||||
{
|
||||
ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
|
||||
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots);
|
||||
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots, 0);
|
||||
}
|
||||
|
||||
#endif // EIGEN_HIDE_HEAVY_CODE
|
||||
|
||||
@ -55,11 +55,11 @@ template<typename _MatrixType> class HouseholderQR
|
||||
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
|
||||
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
|
||||
Options = MatrixType::Options,
|
||||
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
|
||||
DiagSizeAtCompileTime = EIGEN_SIZE_MIN(ColsAtCompileTime,RowsAtCompileTime)
|
||||
};
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType;
|
||||
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor> MatrixQType;
|
||||
typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
|
||||
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
|
||||
typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
|
||||
@ -191,7 +191,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
|
||||
|
||||
RowVectorType temp(cols);
|
||||
|
||||
for (int k = 0; k < size; ++k)
|
||||
for(int k = 0; k < size; ++k)
|
||||
{
|
||||
int remainingRows = rows - k;
|
||||
int remainingCols = cols - k - 1;
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
|
||||
// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
|
||||
//
|
||||
// Eigen is free software; you can redistribute it and/or
|
||||
// modify it under the terms of the GNU Lesser General Public
|
||||
@ -23,7 +23,7 @@
|
||||
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
#include "main.h"
|
||||
#include <Eigen/Householder>
|
||||
#include <Eigen/QR>
|
||||
|
||||
template<typename MatrixType> void householder(const MatrixType& m)
|
||||
{
|
||||
@ -40,7 +40,13 @@ template<typename MatrixType> void householder(const MatrixType& m)
|
||||
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
|
||||
typedef Matrix<Scalar, ei_decrement_size<MatrixType::RowsAtCompileTime>::ret, 1> EssentialVectorType;
|
||||
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
|
||||
typedef Matrix<Scalar, Dynamic, MatrixType::ColsAtCompileTime> HBlockMatrixType;
|
||||
typedef Matrix<Scalar, Dynamic, 1> HCoeffsVectorType;
|
||||
|
||||
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> RightSquareMatrixType;
|
||||
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic> VBlockMatrixType;
|
||||
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime> TMatrixType;
|
||||
|
||||
Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols));
|
||||
Scalar* tmp = &_tmp.coeffRef(0,0);
|
||||
|
||||
@ -85,8 +91,42 @@ template<typename MatrixType> void householder(const MatrixType& m)
|
||||
VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0)));
|
||||
VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha);
|
||||
|
||||
// test householder sequence
|
||||
// TODO test HouseholderSequence
|
||||
// test householder sequence on the left with a shift
|
||||
|
||||
int shift = ei_random(0, std::max(rows-2,0));
|
||||
int brows = rows - shift;
|
||||
m1.setRandom(rows, cols);
|
||||
HBlockMatrixType hbm = m1.block(shift,0,brows,cols);
|
||||
HouseholderQR<HBlockMatrixType> qr(hbm);
|
||||
m2 = m1;
|
||||
m2.block(shift,0,brows,cols) = qr.matrixQR();
|
||||
HCoeffsVectorType hc = qr.hCoeffs().conjugate();
|
||||
HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc, false, hc.size(), shift);
|
||||
MatrixType m5 = m2;
|
||||
m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
|
||||
VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
|
||||
m3 = hseq;
|
||||
VERIFY_IS_APPROX(m3*m5, m1); // test evaluating hseq to a dense matrix, then applying
|
||||
|
||||
#if 0
|
||||
// test householder sequence on the right with a shift
|
||||
|
||||
TMatrixType tm1 = m1.transpose();
|
||||
TMatrixType tm2 = m2.transpose();
|
||||
|
||||
int bcols = cols - shift;
|
||||
VBlockMatrixType vbm =
|
||||
HouseholderQR<HBlockMatrixType> qr(hbm);
|
||||
m2 = m1;
|
||||
m2.block(shift,0,brows,cols) = qr.matrixQR();
|
||||
HCoeffsVectorType hc = qr.hCoeffs().conjugate();
|
||||
HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc, false, hc.size(), shift);
|
||||
MatrixType m5 = m2;
|
||||
m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
|
||||
VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
|
||||
m3 = hseq;
|
||||
VERIFY_IS_APPROX(m3*m5, m1); // test evaluating hseq to a dense matrix, then applying
|
||||
#endif
|
||||
}
|
||||
|
||||
void test_householder()
|
||||
@ -98,6 +138,6 @@ void test_householder()
|
||||
CALL_SUBTEST_4( householder(Matrix<float,4,4>()) );
|
||||
CALL_SUBTEST_5( householder(MatrixXd(10,12)) );
|
||||
CALL_SUBTEST_6( householder(MatrixXcf(16,17)) );
|
||||
CALL_SUBTEST_7( householder(MatrixXf(25,7)) );
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -36,14 +36,11 @@ template<typename MatrixType> void qr(const MatrixType& m)
|
||||
|
||||
MatrixType a = MatrixType::Random(rows,cols);
|
||||
HouseholderQR<MatrixType> qrOfA(a);
|
||||
MatrixType r = qrOfA.matrixQR();
|
||||
|
||||
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);
|
||||
|
||||
MatrixType r = qrOfA.matrixQR().template triangularView<Upper>();
|
||||
VERIFY_IS_APPROX(a, qrOfA.householderQ() * r);
|
||||
}
|
||||
|
||||
@ -114,7 +111,7 @@ template<typename MatrixType> void qr_verify_assert()
|
||||
|
||||
void test_qr()
|
||||
{
|
||||
for(int i = 0; i < 1; i++) {
|
||||
for(int i = 0; i < g_repeat; i++) {
|
||||
CALL_SUBTEST_1( qr(MatrixXf(47,40)) );
|
||||
CALL_SUBTEST_2( qr(MatrixXcd(17,7)) );
|
||||
CALL_SUBTEST_3(( qr_fixedsize<Matrix<float,3,4>, 2 >() ));
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user