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

bug #1493: Make representation of HouseholderSequence consistent and working for complex numbers. Made corresponding unit test actually test that. Also simplify implementation of QR decompositions

Browse Source
This commit is contained in:
Christoph Hertzberg 2018-04-15 10:15:28 +02:00
parent c9ecfff2e6
commit 42715533f1
9 changed files with 60 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Eigenvalues/HessenbergDecomposition.h
View File

@ -315,7 +315,7 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector
// A = A H'
matA.rightCols(remainingSize)
.applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0));
.applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1), numext::conj(h), &temp.coeffRef(0));
}
}

8
Eigen/src/Householder/Householder.h
View File

@ -103,7 +103,7 @@ void MatrixBase<Derived>::makeHouseholder(
* \param essential the essential part of the vector \c v
* \param tau the scaling factor of the Householder transformation
* \param workspace a pointer to working space with at least
* this->cols() * essential.size() entries
* this->cols() entries
*
* \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
* MatrixBase::applyHouseholderOnTheRight()
@ -140,7 +140,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft(
* \param essential the essential part of the vector \c v
* \param tau the scaling factor of the Householder transformation
* \param workspace a pointer to working space with at least
* this->cols() * essential.size() entries
* this->rows() entries
*
* \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
* MatrixBase::applyHouseholderOnTheLeft()
@ -160,10 +160,10 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight(
{
Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
tmp.noalias() = right * essential.conjugate();
tmp.noalias() = right * essential;
tmp += this->col(0);
this->col(0) -= tau * tmp;
right.noalias() -= tau * tmp * essential.transpose();
right.noalias() -= tau * tmp * essential.adjoint();
}
}

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

@ -140,6 +140,22 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
Side
> ConjugateReturnType;
typedef HouseholderSequence<
VectorsType,
typename internal::conditional<NumTraits<Scalar>::IsComplex,
typename internal::remove_all<typename CoeffsType::ConjugateReturnType>::type,
CoeffsType>::type,
Side
> AdjointReturnType;
typedef HouseholderSequence<
typename internal::conditional<NumTraits<Scalar>::IsComplex,
typename internal::remove_all<typename VectorsType::ConjugateReturnType>::type,
VectorsType>::type,
CoeffsType,
Side
> TransposeReturnType;
/** \brief Constructor.
* \param[in] v %Matrix containing the essential parts of the Householder vectors
* \param[in] h Vector containing the Householder coefficients
@ -206,9 +222,12 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
}
/** \brief %Transpose of the Householder sequence. */
HouseholderSequence transpose() const
TransposeReturnType transpose() const
{
return HouseholderSequence(*this).setTrans(!m_trans);
return TransposeReturnType(m_vectors.conjugate(), m_coeffs)
.setTrans(!m_trans)
.setLength(m_length)
.setShift(m_shift);
}
/** \brief Complex conjugate of the Householder sequence. */
@ -221,13 +240,16 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
}
/** \brief Adjoint (conjugate transpose) of the Householder sequence. */
ConjugateReturnType adjoint() const
AdjointReturnType adjoint() const
{
return conjugate().setTrans(!m_trans);
return AdjointReturnType(m_vectors, m_coeffs.conjugate())
.setTrans(!m_trans)
.setLength(m_length)
.setShift(m_shift);
}
/** \brief Inverse of the Householder sequence (equals the adjoint). */
ConjugateReturnType inverse() const { return adjoint(); }
AdjointReturnType inverse() const { return adjoint(); }
/** \internal */
template<typename DestType> inline void evalTo(DestType& dst) const
@ -273,10 +295,10 @@ template<typename VectorsType, typename CoeffsType, int Side> class HouseholderS
Index cornerSize = rows() - k - m_shift;
if(m_trans)
dst.bottomRightCorner(cornerSize, cornerSize)
.applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
.applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
else
dst.bottomRightCorner(cornerSize, cornerSize)
.applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), &workspace.coeffRef(0));
.applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
}
}
}

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

@ -595,11 +595,7 @@ void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &
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()
);
c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() );
m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
.template triangularView<Upper>()

7
Eigen/src/QR/CompleteOrthogonalDecomposition.h
View File

@ -452,7 +452,7 @@ void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
// Apply Z(k) to the first k rows of X_k
m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
.applyHouseholderOnTheRight(
m_cpqr.m_qr.row(k).tail(cols - rank).transpose(), m_zCoeffs(k),
m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k),
&m_temp(0));
}
if (k != rank - 1) {
@ -500,11 +500,8 @@ void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
}
// Compute c = Q^* * rhs
// Note that the matrix Q = H_0^* H_1^*... so its inverse is
// Q^* = (H_0 H_1 ...)^T
typename RhsType::PlainObject c(rhs);
c.applyOnTheLeft(
householderSequence(matrixQTZ(), hCoeffs()).setLength(rank).transpose());
c.applyOnTheLeft(matrixQ().setLength(rank).adjoint());
// Solve T z = c(1:rank, :)
dst.topRows(rank) = matrixT()

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

@ -353,11 +353,7 @@ void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) c
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()
);
c.applyOnTheLeft(householderQ().setLength(rank).adjoint() );
m_qr.topLeftCorner(rank, rank)
.template triangularView<Upper>()

2
Eigen/src/SVD/UpperBidiagonalization.h
View File

@ -127,7 +127,7 @@ void upperbidiagonalization_inplace_unblocked(MatrixType& mat,
.makeHouseholderInPlace(mat.coeffRef(k,k+1), upper_diagonal[k]);
// apply householder transform to remaining part of mat on the left
mat.bottomRightCorner(remainingRows-1, remainingCols)
.applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).transpose(), mat.coeff(k,k+1), tempData);
.applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).adjoint(), mat.coeff(k,k+1), tempData);
}
}

29
test/householder.cpp
View File

@ -49,6 +49,17 @@ template<typename MatrixType> void householder(const MatrixType& m)
v1.applyHouseholderOnTheLeft(essential,beta,tmp);
VERIFY_IS_APPROX(v1.norm(), v2.norm());
// reconstruct householder matrix:
SquareMatrixType id, H1, H2;
id.setIdentity(rows, rows);
H1 = H2 = id;
VectorType vv(rows);
vv << Scalar(1), essential;
H1.applyHouseholderOnTheLeft(essential, beta, tmp);
H2.applyHouseholderOnTheRight(essential, beta, tmp);
VERIFY_IS_APPROX(H1, H2);
VERIFY_IS_APPROX(H1, id - beta * vv*vv.adjoint());
MatrixType m1(rows, cols),
m2(rows, cols);
@ -69,7 +80,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
m3.rowwise() = v1.transpose();
m4 = m3;
m3.row(0).makeHouseholder(essential, beta, alpha);
m3.applyHouseholderOnTheRight(essential,beta,tmp);
m3.applyHouseholderOnTheRight(essential.conjugate(),beta,tmp);
VERIFY_IS_APPROX(m3.norm(), m4.norm());
if(rows>=2) VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
VERIFY_IS_MUCH_SMALLER_THAN(numext::imag(m3(0,0)), numext::real(m3(0,0)));
@ -104,14 +115,14 @@ template<typename MatrixType> void householder(const MatrixType& m)
VERIFY_IS_APPROX(hseq_mat.adjoint(), hseq_mat_adj);
VERIFY_IS_APPROX(hseq_mat.conjugate(), hseq_mat_conj);
VERIFY_IS_APPROX(hseq_mat.transpose(), hseq_mat_trans);
VERIFY_IS_APPROX(hseq_mat * m6, hseq_mat * m6);
VERIFY_IS_APPROX(hseq_mat.adjoint() * m6, hseq_mat_adj * m6);
VERIFY_IS_APPROX(hseq_mat.conjugate() * m6, hseq_mat_conj * m6);
VERIFY_IS_APPROX(hseq_mat.transpose() * m6, hseq_mat_trans * m6);
VERIFY_IS_APPROX(m6 * hseq_mat, m6 * hseq_mat);
VERIFY_IS_APPROX(m6 * hseq_mat.adjoint(), m6 * hseq_mat_adj);
VERIFY_IS_APPROX(m6 * hseq_mat.conjugate(), m6 * hseq_mat_conj);
VERIFY_IS_APPROX(m6 * hseq_mat.transpose(), m6 * hseq_mat_trans);
VERIFY_IS_APPROX(hseq * m6, hseq_mat * m6);
VERIFY_IS_APPROX(hseq.adjoint() * m6, hseq_mat_adj * m6);
VERIFY_IS_APPROX(hseq.conjugate() * m6, hseq_mat_conj * m6);
VERIFY_IS_APPROX(hseq.transpose() * m6, hseq_mat_trans * m6);
VERIFY_IS_APPROX(m6 * hseq, m6 * hseq_mat);
VERIFY_IS_APPROX(m6 * hseq.adjoint(), m6 * hseq_mat_adj);
VERIFY_IS_APPROX(m6 * hseq.conjugate(), m6 * hseq_mat_conj);
VERIFY_IS_APPROX(m6 * hseq.transpose(), m6 * hseq_mat_trans);
// test householder sequence on the right with a shift

2
test/qr.cpp
View File

@ -85,7 +85,7 @@ template<typename MatrixType> void qr_invertible()
qr.compute(m1);
VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
// This test is tricky if the determinant becomes too small.
// Since we generate random numbers with magnitude rrange [0,1], the average determinant is 0.5^size
// Since we generate random numbers with magnitude range [0,1], the average determinant is 0.5^size
VERIFY_IS_MUCH_SMALLER_THAN( abs(absdet-qr.absDeterminant()), numext::maxi(RealScalar(pow(0.5,size)),numext::maxi<RealScalar>(abs(absdet),abs(qr.absDeterminant()))) );
}