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

include the fixes of the third edition

Browse Source
This commit is contained in:
Gael Guennebaud 2009-07-06 15:01:30 +02:00
parent 0c2232e5d9
commit f84bd3e7b1
1 changed files with 66 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

91
Eigen/src/SVD/SVD.h
View File

@ -111,7 +111,7 @@ template<typename MatrixType> class SVD
protected: protected:
// Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow. // Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow.
inline static Scalar pythagora(Scalar a, Scalar b) inline static Scalar pythag(Scalar a, Scalar b)
{ {
Scalar abs_a = ei_abs(a); Scalar abs_a = ei_abs(a);
Scalar abs_b = ei_abs(b); Scalar abs_b = ei_abs(b);
@ -138,14 +138,13 @@ template<typename MatrixType> class SVD
/** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
* *
* \note this code has been adapted from Numerical Recipes, second edition. * \note this code has been adapted from Numerical Recipes, third edition.
*/ */
template<typename MatrixType> template<typename MatrixType>
void SVD<MatrixType>::compute(const MatrixType& matrix) void SVD<MatrixType>::compute(const MatrixType& matrix)
{ {
const int m = matrix.rows(); const int m = matrix.rows();
const int n = matrix.cols(); const int n = matrix.cols();
const int nu = std::min(m,n);
m_matU.resize(m, m); m_matU.resize(m, m);
m_matU.setZero(); m_matU.setZero();
@ -158,22 +157,24 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
MatrixType A = matrix; MatrixType A = matrix;
SingularValuesType& W = m_sigma; SingularValuesType& W = m_sigma;
int flag,i,its,j,jj,k,l,nm; bool flag;
int i,its,j,jj,k,l,nm;
Scalar anorm, c, f, g, h, s, scale, x, y, z; Scalar anorm, c, f, g, h, s, scale, x, y, z;
bool convergence = true; bool convergence = true;
Scalar eps = precision<Scalar>();
Matrix<Scalar,Dynamic,1> rv1(n); Matrix<Scalar,Dynamic,1> rv1(n);
g = scale = anorm = 0; g = scale = anorm = 0;
// Householder reduction to bidiagonal form. // Householder reduction to bidiagonal form.
for (i=0; i<n; i++) for (i=0; i<n; i++)
{ {
l = i+1; l = i+2;
rv1[i] = scale*g; rv1[i] = scale*g;
g = s = scale = 0.0; g = s = scale = 0.0;
if (i < m) if (i < m)
{ {
scale = A.col(i).end(m-i).cwise().abs().sum(); scale = A.col(i).end(m-i).cwise().abs().sum();
if (scale) if (scale != Scalar(0))
{ {
for (k=i; k<m; k++) for (k=i; k<m; k++)
{ {
@ -184,7 +185,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
g = -sign( ei_sqrt(s), f ); g = -sign( ei_sqrt(s), f );
h = f*g - s; h = f*g - s;
A(i, i)=f-g; A(i, i)=f-g;
for (j=l; j<n; j++) for (j=l-1; j<n; j++)
{ {
s = A.col(i).end(m-i).dot(A.col(j).end(m-i)); s = A.col(i).end(m-i).dot(A.col(j).end(m-i));
f = s/h; f = s/h;
@ -195,39 +196,38 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
} }
W[i] = scale * g; W[i] = scale * g;
g = s = scale = 0.0; g = s = scale = 0.0;
if (i < m && i != (n-1)) if (i+1 <= m && i+1 != n)
{ {
scale = A.row(i).end(n-l).cwise().abs().sum(); scale = A.row(i).end(n-l+1).cwise().abs().sum();
if (scale) if (scale != Scalar(0))
{ {
for (k=l; k<n; k++) for (k=l-1; k<n; k++)
{ {
A(i, k) /= scale; A(i, k) /= scale;
s += A(i, k)*A(i, k); s += A(i, k)*A(i, k);
} }
f = A(i, l); f = A(i,l-1);
g = -sign(ei_sqrt(s),f); g = -sign(ei_sqrt(s),f);
h = f*g - s; h = f*g - s;
A(i, l) = f-g; A(i,l-1) = f-g;
for (k=l; k<n; k++) rv1.end(n-l+1) = A.row(i).end(n-l+1)/h;
rv1[k] = A(i, k)/h; for (j=l-1; j<m; j++)
for (j=l; j<m; j++)
{ {
s = A.row(j).end(n-l).dot(A.row(i).end(n-l)); s = A.row(j).end(n-l+1).dot(A.row(i).end(n-l+1));
A.row(j).end(n-l) += s*rv1.end(n-l).transpose(); A.row(j).end(n-l+1) += s*rv1.end(n-l+1).transpose();
} }
A.row(i).end(n-l) *= scale; A.row(i).end(n-l+1) *= scale;
} }
} }
anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) ); anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) );
} }
// Accumulation of right-hand transformations. // Accumulation of right-hand transformations.
for (i=(n-1); i>=0; i--) for (i=n-1; i>=0; i--)
{ {
//Accumulation of right-hand transformations. //Accumulation of right-hand transformations.
if (i < (n-1)) if (i < n-1)
{ {
if (g) if (g != Scalar(0.0))
{ {
for (j=l; j<n; j++) //Double division to avoid possible underflow. for (j=l; j<n; j++) //Double division to avoid possible underflow.
V(j, i) = (A(i, j)/A(i, l))/g; V(j, i) = (A(i, j)/A(i, l))/g;
@ -249,14 +249,14 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
{ {
l = i+1; l = i+1;
g = W[i]; g = W[i];
for (j=l; j<n; j++) if (n-l>0)
A(i, j)=0.0; A.row(i).end(n-l).setZero();
if (g) if (g != Scalar(0.0))
{ {
g = (Scalar)1.0/g; g = Scalar(1.0)/g;
for (j=l; j<n; j++) for (j=l; j<n; j++)
{ {
s = A.col(i).end(m-i).dot(A.col(j).end(m-i)); s = A.col(i).end(m-l).dot(A.col(j).end(m-l));
f = (s/A(i,i))*g; f = (s/A(i,i))*g;
A.col(j).end(m-i) += f * A.col(i).end(m-i); A.col(j).end(m-i) += f * A.col(i).end(m-i);
} }
@ -268,38 +268,41 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
} }
// Diagonalization of the bidiagonal form: Loop over // Diagonalization of the bidiagonal form: Loop over
// singular values, and over allowed iterations. // singular values, and over allowed iterations.
for (k=(n-1); k>=0; k--) for (k=n-1; k>=0; k--)
{ {
for (its=1; its<=max_iters; its++) for (its=0; its<max_iters; its++)
{ {
flag=1; flag = true;
for (l=k; l>=0; l--) for (l=k; l>=0; l--)
{ {
// Test for splitting. // Test for splitting.
nm = l-1; nm = l-1;
// Note that rv1[1] is always zero. // Note that rv1[1] is always zero.
if ((double)(ei_abs(rv1[l])+anorm) == anorm) //if ((double)(ei_abs(rv1[l])+anorm) == anorm)
if (l==0 || ei_abs(rv1[l]) <= eps*anorm)
{ {
flag=0; flag = false;
break; break;
} }
if ((double)(ei_abs(W[nm])+anorm) == anorm) //if ((double)(ei_abs(W[nm])+anorm) == anorm)
if (ei_abs(W[nm]) <= eps*anorm)
break; break;
} }
if (flag) if (flag)
{ {
c=0.0; //Cancellation of rv1[l], if l > 1. c = 0.0; //Cancellation of rv1[l], if l > 0.
s = 1.0; s = 1.0;
for (i=l ;i<=k; i++) for (i=l ;i<k+1; i++)
{ {
f = s*rv1[i]; f = s*rv1[i];
rv1[i] = c*rv1[i]; rv1[i] = c*rv1[i];
if ((double)(ei_abs(f)+anorm) == anorm) //if ((double)(ei_abs(f)+anorm) == anorm)
if (ei_abs(f) <= eps*anorm)
break; break;
g = W[i]; g = W[i];
h = pythagora(f,g); h = pythag(f,g);
W[i] = h; W[i] = h;
h = (Scalar)1.0/h; h = Scalar(1.0)/h;
c = g*h; c = g*h;
s = -f*h; s = -f*h;
for (j=0; j<m; j++) for (j=0; j<m; j++)
@ -320,7 +323,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
} }
break; break;
} }
if (its == max_iters) if (its+1 == max_iters)
{ {
convergence = false; convergence = false;
} }
@ -329,8 +332,8 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
y = W[nm]; y = W[nm];
g = rv1[nm]; g = rv1[nm];
h = rv1[k]; h = rv1[k];
f = ((y-z)*(y+z) + (g-h)*(g+h))/((Scalar)2.0*h*y); f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y);
g = pythagora(f,1.0); g = pythag(f,1.0);
f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x; f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x;
c = s = 1.0; c = s = 1.0;
//Next QR transformation: //Next QR transformation:
@ -341,7 +344,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
y = W[i]; y = W[i];
h = s*g; h = s*g;
g = c*g; g = c*g;
z = pythagora(f,h); z = pythag(f,h);
rv1[j] = z; rv1[j] = z;
c = f/z; c = f/z;
s = h/z; s = h/z;
@ -356,10 +359,10 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
V(jj,j) = x*c + z*s; V(jj,j) = x*c + z*s;
V(jj,i) = z*c - x*s; V(jj,i) = z*c - x*s;
} }
z = pythagora(f,h); z = pythag(f,h);
W[j] = z; W[j] = z;
// Rotation can be arbitrary if z = 0. // Rotation can be arbitrary if z = 0.
if (z) if (z!=Scalar(0))
{ {
z = Scalar(1.0)/z; z = Scalar(1.0)/z;
c = f*z; c = f*z;