mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-11-27 06:30:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

change the make householder algorithm so that the remaining coefficient

Browse Source 
is real, and make Tridiagonalization use it
This commit is contained in:
Gael Guennebaud 2009-08-17 17:04:32 +02:00
parent e125c199bb
commit ff0f005d4c
7 changed files with 75 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Eigen/src/Core/MatrixBase.h
View File

@ -786,17 +786,17 @@ template<typename Derived> class MatrixBase
////////// Householder module ///////////
void makeHouseholderInPlace(RealScalar *tau, Scalar *beta);
void makeHouseholderInPlace(Scalar *tau, RealScalar *beta);
template<typename EssentialPart>
void makeHouseholder(EssentialPart *essential,
RealScalar *tau, Scalar *beta) const;
Scalar *tau, RealScalar *beta) const;
template<typename EssentialPart>
void applyHouseholderOnTheLeft(const EssentialPart& essential,
const RealScalar& tau,
const Scalar& tau,
Scalar* workspace);
template<typename EssentialPart>
void applyHouseholderOnTheRight(const EssentialPart& essential,
const RealScalar& tau,
const Scalar& tau,
Scalar* workspace);
///////// Jacobi module /////////

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

@ -42,7 +42,7 @@ void makeTrivialHouseholder(
}
template<typename Derived>
void MatrixBase<Derived>::makeHouseholderInPlace(RealScalar *tau, Scalar *beta)
void MatrixBase<Derived>::makeHouseholderInPlace(Scalar *tau, RealScalar *beta)
{
VectorBlock<Derived, ei_decrement_size<SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1);
makeHouseholder(&essentialPart, tau, beta);
@ -67,33 +67,35 @@ template<typename Derived>
template<typename EssentialPart>
void MatrixBase<Derived>::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<Scalar>()) * _squaredNorm)
VectorBlock<Derived, EssentialPart::SizeAtCompileTime> 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<Derived, EssentialPart::SizeAtCompileTime> 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<typename Derived>
template<typename EssentialPart>
void MatrixBase<Derived>::applyHouseholderOnTheLeft(
const EssentialPart& essential,
const RealScalar& tau,
const Scalar& tau,
Scalar* workspace)
{
Map<Matrix<Scalar, 1, ColsAtCompileTime, PlainMatrixType::Options, 1, MaxColsAtCompileTime> > tmp(workspace,cols());
@ -108,7 +110,7 @@ template<typename Derived>
template<typename EssentialPart>
void MatrixBase<Derived>::applyHouseholderOnTheRight(
const EssentialPart& essential,
const RealScalar& tau,
const Scalar& tau,
Scalar* workspace)
{
Map<Matrix<Scalar, RowsAtCompileTime, 1, PlainMatrixType::Options, MaxRowsAtCompileTime, 1> > tmp(workspace,rows());

23
Eigen/src/QR/QR.h
View File

@ -53,7 +53,7 @@ template<typename MatrixType> class HouseholderQR
typedef typename MatrixType::RealScalar RealScalar;
typedef Block<MatrixType, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixRBlockType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
typedef Matrix<RealScalar, MinSizeAtCompileTime, 1> HCoeffsType;
typedef Matrix<Scalar, MinSizeAtCompileTime, 1> HCoeffsType;
/**
* \brief Default Constructor.
@ -132,11 +132,13 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::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<MatrixType>::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<typename MatrixType>
MatrixType HouseholderQR<MatrixType>::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<MatrixType>::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;
}

80
Eigen/src/QR/Tridiagonalization.h
View File

@ -198,65 +198,29 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
{
assert(matA.rows()==matA.cols());
int n = matA.rows();
for (int i = 0; i<n-2; ++i)
Matrix<Scalar,1,Dynamic> aux(n);
for (int i = 0; i<n-1; ++i)
{
// let's consider the vector v = i-th column starting at position i+1
int remainingSize = n-i-1;
RealScalar beta;
Scalar h;
matA.col(i).end(remainingSize).makeHouseholderInPlace(&h, &beta);
// start of the householder transformation
// squared norm of the vector v skipping the first element
RealScalar v1norm2 = matA.col(i).end(n-(i+2)).squaredNorm();
// 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;
// FIXME comparing against 1
if (ei_isMuchSmallerThan(v1norm2,static_cast<Scalar>(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<LowerTriangular>()
* (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<LowerTriangular>()
* (h * matA.col(i).end(n-i-1)));
matA.corner(BottomRight, remainingSize, remainingSize).template selfadjointView<LowerTriangular>()
.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<LowerTriangular>()
.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<Scalar>::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<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con
Matrix<Scalar,1,Dynamic> 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<n;++j)
matQ.col(j).end(n-i-1) -= aux.coeff(j) * m_matrix.col(i).end(n-i-1);
// matQ.corner(BottomRight,n-i-1,n-i-1) -= (m_matrix.col(i).end(n-i-1) * aux.end(n-i-1)).lazy();
m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp;
matQ.corner(BottomRight,n-i-1,n-i-1)
.applyHouseholderOnTheLeft(m_matrix.col(i).end(n-i-2), ei_conj(m_hCoeffs.coeff(i)), &aux.coeffRef(0,0));
}
}

2
test/eigensolver_selfadjoint.cpp
View File

@ -120,8 +120,8 @@ void test_eigensolver_selfadjoint()
CALL_SUBTEST( selfadjointeigensolver(Matrix3f()) );
CALL_SUBTEST( selfadjointeigensolver(Matrix4d()) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXf(10,10)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(17,17)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) );
CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(17,17)) );
// some trivial but implementation-wise tricky cases
CALL_SUBTEST( selfadjointeigensolver(MatrixXd(1,1)) );

31
test/householder.cpp
View File

@ -27,6 +27,8 @@
template<typename MatrixType> void householder(const MatrixType& m)
{
static bool even = true;
even = !even;
/* this test covers the following files:
Householder.h
*/
@ -38,46 +40,55 @@ 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;
Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _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<double,2,2>()) );
CALL_SUBTEST( householder(Matrix<float,2,3>()) );
CALL_SUBTEST( householder(Matrix<double,3,5>()) );

2
test/qr.cpp
View File

@ -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++) {