mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

finally, the good approach was two-sided Jacobi. Indeed, it allows

Browse Source 
to guarantee the precision of the output, which is very valuable.
Here, we guarantee that the diagonal matrix returned by the SVD is
actually diagonal, to machine precision.

Performance isn't bad at all at 50% of the current householder SVD
performance for a 200x200 matrix (no vectorization) and we have
lots of room for improvement.
This commit is contained in:
Benoit Jacob 2009-08-12 02:35:07 -04:00
parent ce033ebdfe
commit 22d65d47d0
4 changed files with 63 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -782,8 +782,9 @@ template<typename Derived> class MatrixBase
void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s); void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s);
void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s); void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s);
bool makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s); bool makeJacobi(int p, int q, Scalar *c, Scalar *s) const;
bool makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s); bool makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const;
bool makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const;
#ifdef EIGEN_MATRIXBASE_PLUGIN #ifdef EIGEN_MATRIXBASE_PLUGIN
#include EIGEN_MATRIXBASE_PLUGIN #include EIGEN_MATRIXBASE_PLUGIN

6
Eigen/src/Jacobi/CMakeLists.txt Normal file
View File

@ -0,0 +1,6 @@
FILE(GLOB Eigen_Jacobi_SRCS "*.h")
INSTALL(FILES
${Eigen_Jacobi_SRCS}
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Jacobi COMPONENT Devel
)

30
Eigen/src/Jacobi/Jacobi.h
View File

@ -48,13 +48,13 @@ void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s
} }
template<typename Scalar> template<typename Scalar>
bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar max_coeff, Scalar *c, Scalar *s) bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s)
{ {
if(ei_abs(y) < max_coeff * 0.5 * machine_epsilon<Scalar>()) if(ei_abs(y) < ei_abs(z-x) * 0.5 * machine_epsilon<Scalar>())
{ {
*c = Scalar(1); *c = Scalar(1);
*s = Scalar(0); *s = Scalar(0);
return true; return false;
} }
else else
{ {
@ -67,23 +67,31 @@ bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar max_coeff, Scalar *c, Sc
t = Scalar(1) / (tau - w); t = Scalar(1) / (tau - w);
*c = Scalar(1) / ei_sqrt(1 + ei_abs2(t)); *c = Scalar(1) / ei_sqrt(1 + ei_abs2(t));
*s = *c * t; *s = *c * t;
return false; return true;
} }
} }
template<typename Derived> template<typename Derived>
inline bool MatrixBase<Derived>::makeJacobi(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s) inline bool MatrixBase<Derived>::makeJacobi(int p, int q, Scalar *c, Scalar *s) const
{ {
return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), max_coeff, c, s); return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), c, s);
} }
template<typename Derived> template<typename Derived>
inline bool MatrixBase<Derived>::makeJacobiForAtA(int p, int q, Scalar max_coeff, Scalar *c, Scalar *s) inline bool MatrixBase<Derived>::makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const
{ {
return ei_makeJacobi(col(p).squaredNorm(), return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(q,p)),
col(p).dot(col(q)), ei_conj(coeff(p,p))*coeff(p,q) + ei_conj(coeff(q,p))*coeff(q,q),
col(q).squaredNorm(), ei_abs2(coeff(p,q)) + ei_abs2(coeff(q,q)),
max_coeff, c,s);
}
template<typename Derived>
inline bool MatrixBase<Derived>::makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const
{
return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(p,q)),
ei_conj(coeff(q,p))*coeff(p,p) + ei_conj(coeff(q,q))*coeff(p,q),
ei_abs2(coeff(q,p)) + ei_abs2(coeff(q,q)),
c,s); c,s);
} }

74
Eigen/src/SVD/JacobiSquareSVD.h
View File

@ -102,69 +102,65 @@ void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType&
if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); if(ComputeU) m_matrixU = MatrixUType::Identity(size,size);
if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); if(ComputeV) m_matrixV = MatrixUType::Identity(size,size);
m_singularValues.resize(size); m_singularValues.resize(size);
RealScalar max_coeff = work_matrix.cwise().abs().maxCoeff(); while(true)
for(int k = 1; k < 40; ++k) { {
bool finished = true; bool finished = true;
for(int p = 1; p < size; ++p) for(int p = 1; p < size; ++p)
{ {
for(int q = 0; q < p; ++q) for(int q = 0; q < p; ++q)
{ {
Scalar c, s; Scalar c, s;
finished &= work_matrix.makeJacobiForAtA(p,q,max_coeff,&c,&s); if(work_matrix.makeJacobiForAtA(p,q,&c,&s))
work_matrix.applyJacobiOnTheRight(p,q,c,s); {
if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s); work_matrix.applyJacobiOnTheRight(p,q,c,s);
if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s);
}
if(work_matrix.makeJacobiForAAt(p,q,&c,&s))
{
work_matrix.applyJacobiOnTheLeft(p,q,c,s);
if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s);
if(std::max(ei_abs(work_matrix.coeff(p,q)), ei_abs(work_matrix.coeff(q,p))) > std::max(ei_abs(work_matrix.coeff(q,q)), ei_abs(work_matrix.coeff(p,p)) ))
{
work_matrix.row(p).swap(work_matrix.row(q));
if(ComputeU) m_matrixU.col(p).swap(m_matrixU.col(q));
}
}
}
}
RealScalar biggest = work_matrix.diagonal().cwise().abs().maxCoeff();
for(int p = 0; p < size; ++p)
{
for(int q = 0; q < size; ++q)
{
if(p!=q && ei_abs(work_matrix.coeff(p,q)) > biggest * machine_epsilon<Scalar>()) finished = false;
} }
} }
if(finished) break; if(finished) break;
} }
m_singularValues = work_matrix.diagonal().cwise().abs();
RealScalar biggestSingularValue = m_singularValues.maxCoeff();
for(int i = 0; i < size; ++i) for(int i = 0; i < size; ++i)
{ {
m_singularValues.coeffRef(i) = work_matrix.col(i).norm(); RealScalar a = ei_abs(work_matrix.coeff(i,i));
m_singularValues.coeffRef(i) = a;
if(ComputeU && !ei_isMuchSmallerThan(a, biggestSingularValue)) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a;
} }
int first_zero = size;
RealScalar biggest = m_singularValues.maxCoeff();
for(int i = 0; i < size; i++) for(int i = 0; i < size; i++)
{ {
int pos; int pos;
RealScalar biggest_remaining = m_singularValues.end(size-i).maxCoeff(&pos); m_singularValues.end(size-i).maxCoeff(&pos);
if(first_zero == size && ei_isMuchSmallerThan(biggest_remaining, biggest)) first_zero = pos + i;
if(pos) if(pos)
{ {
pos += i; pos += i;
std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos));
if(ComputeU) work_matrix.col(pos).swap(work_matrix.col(i)); if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i));
if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i));
} }
} }
if(ComputeU)
{
for(int i = 0; i < first_zero; ++i)
{
m_matrixU.col(i) = work_matrix.col(i) / m_singularValues.coeff(i);
}
if(first_zero < size)
{
for(int i = first_zero; i < size; ++i)
{
for(int j = 0; j < size; ++j)
{
m_matrixU.col(i).setZero();
m_matrixU.coeffRef(j,i) = Scalar(1);
for(int k = 0; k < first_zero; ++k)
m_matrixU.col(i) -= m_matrixU.col(i).dot(m_matrixU.col(k)) * m_matrixU.col(k);
RealScalar n = m_matrixU.col(i).norm();
if(!ei_isMuchSmallerThan(n, biggest))
{
m_matrixU.col(i) /= n;
break;
}
}
}
}
}
m_isInitialized = true; m_isInitialized = true;
} }
#endif // EIGEN_JACOBISQUARESVD_H #endif // EIGEN_JACOBISQUARESVD_H