GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 03:39:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

make jacobi SVD more robust after experimenting with very nasty matrices...

Browse Source 
it turns out to be better to repeat the jacobi steps on a given (p,q) pair until it
is diagonal to machine precision, before going to the next (p,q) pair. it's also
an optimization as experiments show that in a majority of cases this allows to find out
that the (p,q) pair is already diagonal to machine precision.
This commit is contained in:
Benoit Jacob 2009-08-12 18:23:39 -04:00
parent 309d540d4a
commit 2b618a2c16
2 changed files with 26 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -50,7 +50,7 @@ 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 *c, Scalar *s) bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s)
{ {
if(ei_abs(y) < ei_abs(z-x) * 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);

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

@ -102,14 +102,16 @@ 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);
while(true)
sweep_again:
for(int p = 1; p < size; ++p)
{ {
bool finished = true; for(int q = 0; q < p; ++q)
for(int p = 1; p < size; ++p)
{ {
for(int q = 0; q < p; ++q) Scalar c, s;
{ while(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
Scalar c, s; > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*machine_epsilon<Scalar>())
{
if(work_matrix.makeJacobiForAtA(p,q,&c,&s)) if(work_matrix.makeJacobiForAtA(p,q,&c,&s))
{ {
work_matrix.applyJacobiOnTheRight(p,q,c,s); work_matrix.applyJacobiOnTheRight(p,q,c,s);
@ -117,9 +119,13 @@ void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType&
} }
if(work_matrix.makeJacobiForAAt(p,q,&c,&s)) if(work_matrix.makeJacobiForAAt(p,q,&c,&s))
{ {
Scalar x = ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(p,q));
Scalar y = ei_conj(work_matrix.coeff(q,p))*work_matrix.coeff(p,p) + ei_conj(work_matrix.coeff(q,q))*work_matrix.coeff(p,q);
Scalar z = ei_abs2(work_matrix.coeff(q,p)) + ei_abs2(work_matrix.coeff(q,q));
work_matrix.applyJacobiOnTheLeft(p,q,c,s); work_matrix.applyJacobiOnTheLeft(p,q,c,s);
if(ComputeU) m_matrixU.applyJacobiOnTheRight(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)) )) if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
> std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q))) )
{ {
work_matrix.row(p).swap(work_matrix.row(q)); work_matrix.row(p).swap(work_matrix.row(q));
if(ComputeU) m_matrixU.col(p).swap(m_matrixU.col(q)); if(ComputeU) m_matrixU.col(p).swap(m_matrixU.col(q));
@ -127,15 +133,18 @@ void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType&
} }
} }
} }
RealScalar biggest = work_matrix.diagonal().cwise().abs().maxCoeff(); }
for(int p = 0; p < size; ++p)
{ RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff();
for(int q = 0; q < size; ++q) RealScalar maxAllowedOffDiag = biggestOnDiag * machine_epsilon<Scalar>();
{ for(int p = 0; p < size; ++p)
if(p!=q && ei_abs(work_matrix.coeff(p,q)) > biggest * machine_epsilon<Scalar>()) finished = false; {
} for(int q = 0; q < p; ++q)
} if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
if(finished) break; goto sweep_again;
for(int q = p+1; q < size; ++q)
if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
goto sweep_again;
} }
m_singularValues = work_matrix.diagonal().cwise().abs(); m_singularValues = work_matrix.diagonal().cwise().abs();