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

make ColPivotingQR use HouseholderSequence

Browse Source
This commit is contained in:
Gael Guennebaud 2009-09-16 15:56:20 +02:00
parent 49dd5d7847
commit 24950bdfcb
6 changed files with 68 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Eigen/src/Core/VectorBlock.h
View File

@ -77,11 +77,12 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock
typedef Block<VectorType,
ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size,
ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size,
PacketAccess> Base;
PacketAccess> _Base;
enum {
IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1
};
public:
_EIGEN_GENERIC_PUBLIC_INTERFACE(VectorBlock, _Base)
using Base::operator=;
using Base::operator+=;

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

@ -80,10 +80,10 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
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); }
ConjugateReturnType conjugate() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
ConjugateReturnType adjoint() const
{ return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); }
@ -136,6 +136,22 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
}
}
template<typename OtherDerived>
typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other) const
{
typename OtherDerived::PlainMatrixType res(other);
applyThisOnTheLeft(res);
return res;
}
template<typename OtherDerived> friend
typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence& h)
{
typename OtherDerived::PlainMatrixType res(other);
h.applyThisOnTheRight(res);
return res;
}
protected:
typename VectorsType::Nested m_vectors;
@ -143,4 +159,10 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
bool m_trans;
};
template<typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType,CoeffsType> makeHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
{
return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans);
}
#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H

59
Eigen/src/QR/ColPivotingHouseholderQR.h
View File

@ -45,14 +45,14 @@
template<typename MatrixType> class ColPivotingHouseholderQR
{
public:
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@ -62,6 +62,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType;
typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
/**
* \brief Default Constructor.
@ -99,7 +100,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
template<typename OtherDerived, typename ResultType>
bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
MatrixQType matrixQ(void) const;
HouseholderSequenceType matrixQ(void) const;
/** \returns a reference to the matrix where the Householder QR decomposition is stored
*/
@ -110,13 +111,13 @@ template<typename MatrixType> class ColPivotingHouseholderQR
}
ColPivotingHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
return m_cols_permutation;
}
/** \returns the absolute value of the determinant of the matrix of which
* *this is the QR decomposition. It has only linear complexity
* (that is, O(n) where n is the dimension of the square matrix)
@ -145,7 +146,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR
* \sa absDeterminant(), MatrixBase::determinant()
*/
typename MatrixType::RealScalar logAbsDeterminant() const;
/** \returns the rank of the matrix of which *this is the QR decomposition.
*
* \note This is computed at the time of the construction of the QR decomposition. This
@ -268,7 +269,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
int cols = matrix.cols();
int size = std::min(rows,cols);
m_rank = size;
m_qr = matrix;
m_hCoeffs.resize(size);
@ -279,18 +280,18 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
IntRowVectorType cols_transpositions(matrix.cols());
m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0;
RealRowVectorType colSqNorms(cols);
for(int k = 0; k < cols; ++k)
colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
RealScalar biggestColSqNorm = colSqNorms.maxCoeff();
for (int k = 0; k < size; ++k)
{
int biggest_col_in_corner;
RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner);
biggest_col_in_corner += k;
// if the corner is negligible, then we have less than full rank, and we can finish early
if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision))
{
@ -302,7 +303,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
}
break;
}
cols_transpositions.coeffRef(k) = biggest_col_in_corner;
if(k != biggest_col_in_corner) {
m_qr.col(k).swap(m_qr.col(biggest_col_in_corner));
@ -316,7 +317,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
m_qr.corner(BottomRight, rows-k, cols-k-1)
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2();
}
@ -326,7 +327,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
return *this;
}
@ -352,16 +353,11 @@ bool ColPivotingHouseholderQR<MatrixType>::solve(
const int rows = m_qr.rows();
const int cols = b.cols();
ei_assert(b.rows() == rows);
typename OtherDerived::PlainMatrixType c(b);
Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
for (int k = 0; k < m_rank; ++k)
{
int remainingSize = rows-k;
c.corner(BottomRight, remainingSize, cols)
.applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
}
// Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose());
if(!isSurjective())
{
@ -381,25 +377,12 @@ bool ColPivotingHouseholderQR<MatrixType>::solve(
return true;
}
/** \returns the matrix Q */
/** \returns the matrix Q as a sequence of householder transformations */
template<typename MatrixType>
typename ColPivotingHouseholderQR<MatrixType>::MatrixQType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
typename ColPivotingHouseholderQR<MatrixType>::HouseholderSequenceType ColPivotingHouseholderQR<MatrixType>::matrixQ() const
{
ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized.");
// compute the product H'_0 H'_1 ... H'_n-1,
// where H_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), ...]
int rows = m_qr.rows();
int cols = m_qr.cols();
int size = std::min(rows,cols);
MatrixQType res = MatrixQType::Identity(rows, rows);
Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
for (int k = size-1; k >= 0; k--)
{
res.block(k, k, rows-k, rows-k)
.applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
}
return res;
return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
}
#endif // EIGEN_HIDE_HEAVY_CODE

26
Eigen/src/QR/FullPivotingHouseholderQR.h
View File

@ -45,14 +45,14 @@
template<typename MatrixType> class FullPivotingHouseholderQR
{
public:
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
Options = MatrixType::Options,
DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@ -106,13 +106,13 @@ template<typename MatrixType> class FullPivotingHouseholderQR
}
FullPivotingHouseholderQR& compute(const MatrixType& matrix);
const IntRowVectorType& colsPermutation() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
return m_cols_permutation;
}
const IntColVectorType& rowsTranspositions() const
{
ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized.");
@ -147,7 +147,7 @@ template<typename MatrixType> class FullPivotingHouseholderQR
* \sa absDeterminant(), MatrixBase::determinant()
*/
typename MatrixType::RealScalar logAbsDeterminant() const;
/** \returns the rank of the matrix of which *this is the QR decomposition.
*
* \note This is computed at the time of the construction of the QR decomposition. This
@ -271,7 +271,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co
int cols = matrix.cols();
int size = std::min(rows,cols);
m_rank = size;
m_qr = matrix;
m_hCoeffs.resize(size);
@ -283,9 +283,9 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co
IntRowVectorType cols_transpositions(matrix.cols());
m_cols_permutation.resize(matrix.cols());
int number_of_transpositions = 0;
RealScalar biggest(0);
for (int k = 0; k < size; ++k)
{
int row_of_biggest_in_corner, col_of_biggest_in_corner;
@ -297,7 +297,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co
row_of_biggest_in_corner += k;
col_of_biggest_in_corner += k;
if(k==0) biggest = biggest_in_corner;
// if the corner is negligible, then we have less than full rank, and we can finish early
if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
{
@ -336,7 +336,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co
m_det_pq = (number_of_transpositions%2) ? -1 : 1;
m_isInitialized = true;
return *this;
}
@ -358,13 +358,13 @@ bool FullPivotingHouseholderQR<MatrixType>::solve(
}
else return false;
}
const int rows = m_qr.rows();
const int cols = b.cols();
ei_assert(b.rows() == rows);
typename OtherDerived::PlainMatrixType c(b);
Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
for (int k = 0; k < m_rank; ++k)
{

8
test/qr_colpivoting.cpp
View File

@ -42,18 +42,18 @@ template<typename MatrixType> void qr()
VERIFY(!qr.isInjective());
VERIFY(!qr.isInvertible());
VERIFY(!qr.isSurjective());
MatrixType r = qr.matrixQR();
// 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 b = qr.matrixQ() * r;
MatrixType c = MatrixType::Zero(rows,cols);
for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i);
VERIFY_IS_APPROX(m1, c);
MatrixType m2 = MatrixType::Random(cols,cols2);
MatrixType m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);

8
test/qr_fullpivoting.cpp
View File

@ -46,14 +46,14 @@ template<typename MatrixType> void qr()
MatrixType r = qr.matrixQR();
// 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 b = qr.matrixQ() * r;
MatrixType c = MatrixType::Zero(rows,cols);
for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i);
VERIFY_IS_APPROX(m1, c);
MatrixType m2 = MatrixType::Random(cols,cols2);
MatrixType m3 = m1*m2;
m2 = MatrixType::Random(cols,cols2);
@ -88,7 +88,7 @@ template<typename MatrixType> void qr_invertible()
m3 = MatrixType::Random(size,size);
VERIFY(qr.solve(m3, &m2));
VERIFY_IS_APPROX(m3, m1*m2);
// now construct a matrix with prescribed determinant
m1.setZero();
for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();