From ff0f005d4c0bd46e88d050b9f147eab810f4814d Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 17 Aug 2009 17:04:32 +0200 Subject: [PATCH] change the make householder algorithm so that the remaining coefficient is real, and make Tridiagonalization use it --- Eigen/src/Core/MatrixBase.h | 8 +-- Eigen/src/Householder/Householder.h | 34 ++++++------ Eigen/src/QR/QR.h | 23 +++++---- Eigen/src/QR/Tridiagonalization.h | 80 +++++++---------------------- test/eigensolver_selfadjoint.cpp | 2 +- test/householder.cpp | 31 +++++++---- test/qr.cpp | 2 +- 7 files changed, 75 insertions(+), 105 deletions(-) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 9e92c043f..1f4c6bf7a 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -786,17 +786,17 @@ template class MatrixBase ////////// Householder module /////////// - void makeHouseholderInPlace(RealScalar *tau, Scalar *beta); + void makeHouseholderInPlace(Scalar *tau, RealScalar *beta); template void makeHouseholder(EssentialPart *essential, - RealScalar *tau, Scalar *beta) const; + Scalar *tau, RealScalar *beta) const; template void applyHouseholderOnTheLeft(const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace); template void applyHouseholderOnTheRight(const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace); ///////// Jacobi module ///////// diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h index a7bbe96ce..769ba3d43 100644 --- a/Eigen/src/Householder/Householder.h +++ b/Eigen/src/Householder/Householder.h @@ -42,7 +42,7 @@ void makeTrivialHouseholder( } template -void MatrixBase::makeHouseholderInPlace(RealScalar *tau, Scalar *beta) +void MatrixBase::makeHouseholderInPlace(Scalar *tau, RealScalar *beta) { VectorBlock::ret> essentialPart(derived(), 1, size()-1); makeHouseholder(&essentialPart, tau, beta); @@ -67,33 +67,35 @@ template template void MatrixBase::makeHouseholder( EssentialPart *essential, - RealScalar *tau, - Scalar *beta) const + Scalar *tau, + RealScalar *beta) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) - RealScalar _squaredNorm = squaredNorm(); - Scalar c0; - if(ei_abs2(coeff(0)) <= ei_abs2(precision()) * _squaredNorm) + VectorBlock tail(derived(), 1, size()-1); + + RealScalar tailSqNorm; + Scalar c0 = coeff(0); + + if( (size()==1 || (tailSqNorm=tail.squaredNorm()) == RealScalar(0)) && ei_imag(c0)==RealScalar(0)) { - c0 = ei_sqrt(_squaredNorm); + *tau = 0; + *beta = ei_real(c0); } else { - Scalar sign = coeff(0) / ei_abs(coeff(0)); - c0 = coeff(0) + sign * ei_sqrt(_squaredNorm); + *beta = ei_sqrt(ei_abs2(c0) + tailSqNorm); + if (ei_real(c0)>=0.) + *beta = -*beta; + *essential = tail / (c0 - *beta); + *tau = ei_conj((*beta - c0) / *beta); } - VectorBlock tail(derived(), 1, size()-1); - *essential = tail / c0; - const RealScalar c0abs2 = ei_abs2(c0); - *tau = RealScalar(2) * c0abs2 / (c0abs2 + _squaredNorm - ei_abs2(coeff(0))); - *beta = coeff(0) - c0; } template template void MatrixBase::applyHouseholderOnTheLeft( const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace) { Map > tmp(workspace,cols()); @@ -108,7 +110,7 @@ template template void MatrixBase::applyHouseholderOnTheRight( const EssentialPart& essential, - const RealScalar& tau, + const Scalar& tau, Scalar* workspace) { Map > tmp(workspace,rows()); diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 0c6142129..90e6f8132 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -53,7 +53,7 @@ template class HouseholderQR typedef typename MatrixType::RealScalar RealScalar; typedef Block MatrixRBlockType; typedef Matrix MatrixTypeR; - typedef Matrix HCoeffsType; + typedef Matrix HCoeffsType; /** * \brief Default Constructor. @@ -132,11 +132,13 @@ HouseholderQR& HouseholderQR::compute(const MatrixType& int remainingRows = rows - k; int remainingCols = cols - k -1; - m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &m_qr.coeffRef(k,k)); - - if (remainingRows>1 && remainingCols>0) - m_qr.corner(BottomRight, remainingRows, remainingCols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + RealScalar beta; + m_qr.col(k).end(remainingRows).makeHouseholderInPlace(&m_hCoeffs.coeffRef(k), &beta); + m_qr.coeffRef(k,k) = beta; + + // apply H to remaining part of m_qr from the left + m_qr.corner(BottomRight, remainingRows, remainingCols) + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); } m_isInitialized = true; return *this; @@ -163,7 +165,7 @@ void HouseholderQR::solve( int remainingSize = rows-k; result->corner(BottomRight, remainingSize, cols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(0)); + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); } const int rank = std::min(result->rows(), result->cols()); @@ -177,8 +179,8 @@ template MatrixType HouseholderQR::matrixQ() const { 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' + // 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(); @@ -187,9 +189,8 @@ MatrixType HouseholderQR::matrixQ() const for (int k = cols-1; k >= 0; k--) { int remainingSize = rows-k; - res.corner(BottomRight, remainingSize, cols-k) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_real(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); + .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); } return res; } diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h index d98322fac..4cc7433c1 100644 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/QR/Tridiagonalization.h @@ -198,65 +198,29 @@ void Tridiagonalization::_compute(MatrixType& matA, CoeffVectorType& { assert(matA.rows()==matA.cols()); int n = matA.rows(); - for (int i = 0; i aux(n); + for (int i = 0; i(1))) - { - hCoeffs.coeffRef(i) = 0.; - } - else - { - Scalar v0 = matA.col(i).coeff(i+1); - RealScalar beta = ei_sqrt(ei_abs2(v0)+v1norm2); - if (ei_real(v0)>=0.) - beta = -beta; - matA.col(i).end(n-(i+2)) *= (Scalar(1)/(v0-beta)); - matA.col(i).coeffRef(i+1) = beta; - Scalar h = (beta - v0) / beta; - // end of the householder transformation + hCoeffs.end(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView() + * (ei_conj(h) * matA.col(i).end(remainingSize))); - // Apply similarity transformation to remaining columns, - // i.e., A = H' A H where H = I - h v v' and v = matA.col(i).end(n-i-1) - matA.col(i).coeffRef(i+1) = 1; + hCoeffs.end(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.end(remainingSize).dot(matA.col(i).end(remainingSize)))) * matA.col(i).end(n-i-1); - hCoeffs.end(n-i-1) = (matA.corner(BottomRight,n-i-1,n-i-1).template selfadjointView() - * (h * matA.col(i).end(n-i-1))); + matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView() + .rankUpdate(matA.col(i).end(remainingSize), hCoeffs.end(remainingSize), -1); - hCoeffs.end(n-i-1) += (h*Scalar(-0.5)*(hCoeffs.end(n-i-1).dot(matA.col(i).end(n-i-1)))) * matA.col(i).end(n-i-1); - - matA.corner(BottomRight, n-i-1, n-i-1).template selfadjointView() - .rankUpdate(matA.col(i).end(n-i-1), hCoeffs.end(n-i-1), -1); - - // note: at that point matA(i+1,i+1) is the (i+1)-th element of the final diagonal - // note: the sequence of the beta values leads to the subdiagonal entries - matA.col(i).coeffRef(i+1) = beta; - - hCoeffs.coeffRef(i) = h; - } - } - if (NumTraits::IsComplex) - { - // Householder transformation on the remaining single scalar - int i = n-2; - Scalar v0 = matA.col(i).coeff(i+1); - RealScalar beta = ei_abs(v0); - if (ei_real(v0)>=RealScalar(0)) - beta = -beta; matA.col(i).coeffRef(i+1) = beta; - // FIXME comparing against 1 - if(ei_isMuchSmallerThan(beta, Scalar(1))) hCoeffs.coeffRef(i) = Scalar(0); - else hCoeffs.coeffRef(i) = (beta - v0) / beta; - } - else - { - hCoeffs.coeffRef(n-2) = 0; + hCoeffs.coeffRef(i) = h; } } @@ -280,16 +244,8 @@ void Tridiagonalization::matrixQInPlace(MatrixBase* q) con Matrix aux(n); for (int i = n-2; i>=0; i--) { - Scalar tmp = m_matrix.coeff(i+1,i); - m_matrix.const_cast_derived().coeffRef(i+1,i) = 1; - - aux.end(n-i-1) = (m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy(); - // rank one update, TODO ! make it works efficiently as expected - for (int j=i+1;j void householder(const MatrixType& m) { + static bool even = true; + even = !even; /* this test covers the following files: Householder.h */ @@ -38,46 +40,55 @@ template void householder(const MatrixType& m) typedef Matrix VectorType; typedef Matrix::ret, 1> EssentialVectorType; typedef Matrix SquareMatrixType; + + Matrix _tmp(std::max(rows,cols)); + Scalar* tmp = &_tmp.coeffRef(0,0); - RealScalar beta; + Scalar beta; + RealScalar alpha; EssentialVectorType essential; VectorType v1 = VectorType::Random(rows), v2; v2 = v1; - v1.makeHouseholder(&essential, &beta); - v1.applyHouseholderOnTheLeft(essential,beta); - + v1.makeHouseholder(&essential, &beta, &alpha); + v1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(v1.norm(), v2.norm()); VERIFY_IS_MUCH_SMALLER_THAN(v1.end(rows-1).norm(), v1.norm()); v1 = VectorType::Random(rows); v2 = v1; - v1.applyHouseholderOnTheLeft(essential,beta); + v1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(v1.norm(), v2.norm()); MatrixType m1(rows, cols), m2(rows, cols); v1 = VectorType::Random(rows); + if(even) v1.end(rows-1).setZero(); m1.colwise() = v1; m2 = m1; - m1.col(0).makeHouseholder(&essential, &beta); - m1.applyHouseholderOnTheLeft(essential,beta); + m1.col(0).makeHouseholder(&essential, &beta, &alpha); + m1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(m1.norm(), m2.norm()); VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm()); + VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m1(0,0)), ei_real(m1(0,0))); + VERIFY_IS_APPROX(ei_real(m1(0,0)), alpha); v1 = VectorType::Random(rows); + if(even) v1.end(rows-1).setZero(); SquareMatrixType m3(rows,rows), m4(rows,rows); m3.rowwise() = v1.transpose(); m4 = m3; - m3.row(0).makeHouseholder(&essential, &beta); - m3.applyHouseholderOnTheRight(essential,beta); + m3.row(0).makeHouseholder(&essential, &beta, &alpha); + m3.applyHouseholderOnTheRight(essential,beta,tmp); VERIFY_IS_APPROX(m3.norm(), m4.norm()); VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm()); + VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0))); + VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha); } void test_householder() { - for(int i = 0; i < g_repeat; i++) { + for(int i = 0; i < 2*g_repeat; i++) { CALL_SUBTEST( householder(Matrix()) ); CALL_SUBTEST( householder(Matrix()) ); CALL_SUBTEST( householder(Matrix()) ); diff --git a/test/qr.cpp b/test/qr.cpp index 99df1d89b..f004a36ca 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -93,8 +93,8 @@ void test_qr() // FIXME : very weird bug here // CALL_SUBTEST( qr(Matrix2f()) ); CALL_SUBTEST( qr(Matrix4d()) ); - CALL_SUBTEST( qr(MatrixXcd(17,7)) ); CALL_SUBTEST( qr(MatrixXf(47,40)) ); + CALL_SUBTEST( qr(MatrixXcd(17,7)) ); } for(int i = 0; i < g_repeat; i++) {