#ifndef _SVD_H #define _SVD_H #include "dense_matrix.h" #include "dense_vector.h" // #include "tnt_array1d.h" // #include "tnt_array1d_utils.h" // #include "tnt_array2d.h" // #include "tnt_array2d_utils.h" /// #include "tnt_math_utils.h" /** Singular Value Decomposition.

For an m-by-n matrix A with m >= n, the singular value decomposition is an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and an n-by-n orthogonal matrix V so that A = U*S*V'.

The singular values, sigma[k] = S[k][k], are ordered so that sigma[0] >= sigma[1] >= ... >= sigma[n-1].

The singular value decompostion always exists, so the constructor will never fail. The matrix condition number and the effective numerical rank can be computed from this decomposition.

(Adapted from JAMA, a Java Matrix Library, developed by jointly by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). */ /** @returns hypotenuse of real (non-complex) Reals a and b by avoiding underflow/overflow using (a * sqrt( 1 + (b/a) * (b/a))), rather than sqrt(a*a + b*b). */ template Real hypot(const Real &a, const Real &b) { if (a== 0) return abs(b); else { Real c = b/a; return abs(a) * sqrt(1 + c*c); } } /** @returns the minimum of Reals a and b. */ template Real min(const Real &a, const Real &b) { return a < b ? a : b; } /** @returns the maximum of Reals a and b. */ template Real max(const Real &a, const Real &b) { return a > b ? a : b; } /** @returns the absolute value of a real (no-complex) Real. */ template Real abs(const Real &a) { return (a > 0 ? a : -a); } template class SVD { DenseMatrix U; DenseMatrix V; DenseVector s; int m, n; public: SVD (const DenseMatrix &Arg); void solve ( const DenseVector &b, DenseVector &x); void getU (DenseMatrix &A); /* Return the right singular vectors */ void getV (DenseMatrix &A); /* Return the one-dimensional array of singular values */ void getSingularValues (DenseVector &x); /* Return the diagonal matrix of singular values @return S */ void getS (DenseMatrix &A); /* Return true if convergent */ bool is_convergent (); /* Two norm (max(S)) */ double norm2 (); /* Two norm of condition number (max(S)/min(S)) */ double cond (); /* Effective numerical matrix rank @return Number of nonnegligible singular values. */ int rank (); }; //END OF Class definitions /* Return the left singular vectors */ template void SVD::getU (DenseMatrix &A) { int minm = min(m+1,n); A.resize(m, minm); for (int i=0; i void SVD::getV (DenseMatrix &A) { A = V; } template bool SVD::is_convergent () { return conv; } /** Return the one-dimensional array of singular values */ template void SVD::getSingularValues (DenseVector &x) { x = s; } /** Return the diagonal matrix of singular values @return S */ template void SVD::getS (DenseMatrix &A) { A.resize(n,n); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { A(i,j) = 0.0; } A(i,i) = s(i); } } /** Two norm (max(S)) */ template double SVD::norm2 () { return s(0); } /** Two norm of condition number (max(S)/min(S)) */ template double SVD::cond () { return s(0)/s(min(m,n)-1); } /** Effective numerical matrix rank @return Number of nonnegligible singular values. */ template int SVD::rank () { double eps = pow(2.0,-52.0); double tol = max(m,n)*s(0)*eps; int r = 0; for (int i = 0; i < s.size(); i++) { if (s(i) > tol) { r++; } } return r; } //////////////////////////////////////////////////////// template SVD::SVD (const DenseMatrix &Arg) { m = Arg.m(); n = Arg.n(); int nu = min(m,n); s = DenseVector(min(m+1,n)); s.zero(); U = DenseMatrix(m, nu); U.zero(); V = DenseMatrix(n,n); V.zero(); DenseVector e(n); DenseVector work(m); DenseMatrix A(Arg); int wantu = 1; /* boolean */ int wantv = 1; /* boolean */ int i=0, j=0, k=0; // Reduce A to bidiagonal form, storing the diagonal elements // in s and the super-diagonal elements in e. int nct = min(m-1,n); int nrt = max(0,min(n-2,m)); for (k = 0; k < max(nct,nrt); k++) { if (k < nct) { // Compute the transformation for the k-th column and // place the k-th diagonal in s(k). // Compute 2-norm of k-th column without under/overflow. s(k) = 0; for (i = k; i < m; i++) { s(k) = hypot(s(k),A(i,k)); } if (s(k) != 0.0) { if (A(k,k) < 0.0) { s(k) = -s(k); } for (i = k; i < m; i++) { A(i,k) /= s(k); } A(k,k) += 1.0; } s(k) = -s(k); } for (j = k+1; j < n; j++) { if ((k < nct) && (s(k) != 0.0)) { // Apply the transformation. double t = 0; for (i = k; i < m; i++) { t += A(i,k)*A(i,j); } t = -t/A(k,k); for (i = k; i < m; i++) { A(i,j) += t*A(i,k); } } // Place the k-th row of A into e for the // subsequent calculation of the row transformation. e(j) = A(k,j); } if (wantu & (k < nct)) { // Place the transformation in U for subsequent back // multiplication. for (i = k; i < m; i++) { U(i,k) = A(i,k); } } if (k < nrt) { // Compute the k-th row transformation and place the // k-th super-diagonal in e(k). // Compute 2-norm without under/overflow. e(k) = 0; for (i = k+1; i < n; i++) { e(k) = hypot(e(k),e(i)); } if (e(k) != 0.0) { if (e(k+1) < 0.0) { e(k) = -e(k); } for (i = k+1; i < n; i++) { e(i) /= e(k); } e(k+1) += 1.0; } e(k) = -e(k); if ((k+1 < m) & (e(k) != 0.0)) { // Apply the transformation. for (i = k+1; i < m; i++) { work(i) = 0.0; } for (j = k+1; j < n; j++) { for (i = k+1; i < m; i++) { work(i) += e(j)*A(i,j); } } for (j = k+1; j < n; j++) { double t = -e(j)/e(k+1); for (i = k+1; i < m; i++) { A(i,j) += t*work(i); } } } if (wantv) { // Place the transformation in V for subsequent // back multiplication. for (i = k+1; i < n; i++) { V(i,k) = e(i); } } } } // Set up the final bidiagonal matrix or order p. int p = min(n,m+1); if (nct < n) { s(nct) = A(nct,nct); } if (m < p) { s(p-1) = 0.0; } if (nrt+1 < p) { e(nrt) = A(nrt,p-1); } e(p-1) = 0.0; // If required, generate U. if (wantu) { for (j = nct; j < nu; j++) { for (i = 0; i < m; i++) { U(i,j) = 0.0; } U(j,j) = 1.0; } for (k = nct-1; k >= 0; k--) { if (s(k) != 0.0) { for (j = k+1; j < nu; j++) { double t = 0; for (i = k; i < m; i++) { t += U(i,k)*U(i,j); } t = -t/U(k,k); for (i = k; i < m; i++) { U(i,j) += t*U(i,k); } } for (i = k; i < m; i++ ) { U(i,k) = -U(i,k); } U(k,k) = 1.0 + U(k,k); for (i = 0; i < k-1; i++) { U(i,k) = 0.0; } } else { for (i = 0; i < m; i++) { U(i,k) = 0.0; } U(k,k) = 1.0; } } } // If required, generate V. if (wantv) { for (k = n-1; k >= 0; k--) { if ((k < nrt) & (e(k) != 0.0)) { for (j = k+1; j < nu; j++) { double t = 0; for (i = k+1; i < n; i++) { t += V(i,k)*V(i,j); } t = -t/V(k+1,k); for (i = k+1; i < n; i++) { V(i,j) += t*V(i,k); } } } for (i = 0; i < n; i++) { V(i,k) = 0.0; } V(k,k) = 1.0; } } // Main iteration loop for the singular values. int pp = p-1; int iter = 0; double eps = pow(2.0,-52.0); while (p > 0) { int k=0; int kase=0; // This section of the program inspects for // negligible elements in the s and e arrays. On // completion the variables kase and k are set as follows. // kase = 1 if s(p) and e(k-1) are negligible and k

= -1; k--) { if (k == -1) { break; } if (abs(e(k)) <= eps*(abs(s(k)) + abs(s(k+1)))) { e(k) = 0.0; break; } } if (k == p-2) { kase = 4; } else { int ks; for (ks = p-1; ks >= k; ks--) { if (ks == k) { break; } double t = (ks != p ? abs(e(ks)) : 0.) + (ks != k+1 ? abs(e(ks-1)) : 0.); if (abs(s(ks)) <= eps*t) { s(ks) = 0.0; break; } } if (ks == k) { kase = 3; } else if (ks == p-1) { kase = 1; } else { kase = 2; k = ks; } } k++; // Perform the task indicated by kase. switch (kase) { // Deflate negligible s(p). case 1: { double f = e(p-2); e(p-2) = 0.0; for (j = p-2; j >= k; j--) { double t = hypot(s(j),f); double cs = s(j)/t; double sn = f/t; s(j) = t; if (j != k) { f = -sn*e(j-1); e(j-1) = cs*e(j-1); } if (wantv) { for (i = 0; i < n; i++) { t = cs*V(i,j) + sn*V(i,p-1); V(i,p-1) = -sn*V(i,j) + cs*V(i,p-1); V(i,j) = t; } } } } break; // Split at negligible s(k). case 2: { double f = e(k-1); e(k-1) = 0.0; for (j = k; j < p; j++) { double t = hypot(s(j),f); double cs = s(j)/t; double sn = f/t; s(j) = t; f = -sn*e(j); e(j) = cs*e(j); if (wantu) { for (i = 0; i < m; i++) { t = cs*U(i,j) + sn*U(i,k-1); U(i,k-1) = -sn*U(i,j) + cs*U(i,k-1); U(i,j) = t; } } } } break; // Perform one qr step. case 3: { // Calculate the shift. double scale = max(max(max(max( abs(s(p-1)),abs(s(p-2))),abs(e(p-2))), abs(s(k))),abs(e(k))); double sp = s(p-1)/scale; double spm1 = s(p-2)/scale; double epm1 = e(p-2)/scale; double sk = s(k)/scale; double ek = e(k)/scale; double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; double c = (sp*epm1)*(sp*epm1); double shift = 0.0; if ((b != 0.0) | (c != 0.0)) { shift = sqrt(b*b + c); if (b < 0.0) { shift = -shift; } shift = c/(b + shift); } double f = (sk + sp)*(sk - sp) + shift; double g = sk*ek; // Chase zeros. for (j = k; j < p-1; j++) { double t = hypot(f,g); double cs = f/t; double sn = g/t; if (j != k) { e(j-1) = t; } f = cs*s(j) + sn*e(j); e(j) = cs*e(j) - sn*s(j); g = sn*s(j+1); s(j+1) = cs*s(j+1); if (wantv) { for (i = 0; i < n; i++) { t = cs*V(i,j) + sn*V(i,j+1); V(i,j+1) = -sn*V(i,j) + cs*V(i,j+1); V(i,j) = t; } } t = hypot(f,g); cs = f/t; sn = g/t; s(j) = t; f = cs*e(j) + sn*s(j+1); s(j+1) = -sn*e(j) + cs*s(j+1); g = sn*e(j+1); e(j+1) = cs*e(j+1); if (wantu && (j < m-1)) { for (i = 0; i < m; i++) { t = cs*U(i,j) + sn*U(i,j+1); U(i,j+1) = -sn*U(i,j) + cs*U(i,j+1); U(i,j) = t; } } } e(p-2) = f; iter = iter + 1; } break; // Convergence. case 4: { // Make the singular values positive. if (s(k) <= 0.0) { s(k) = (s(k) < 0.0 ? -s(k) : 0.0); if (wantv) { for (i = 0; i <= pp; i++) { V(i,k) = -V(i,k); } } } // Order the singular values. while (k < pp) { if (s(k) >= s(k+1)) { break; } double t = s(k); s(k) = s(k+1); s(k+1) = t; if (wantv && (k < n-1)) { for (i = 0; i < n; i++) { t = V(i,k+1); V(i,k+1) = V(i,k); V(i,k) = t; } } if (wantu && (k < m-1)) { for (i = 0; i < m; i++) { t = U(i,k+1); U(i,k+1) = U(i,k); U(i,k) = t; } } k++; } iter = 0; p--; } break; } } } ///////// End of function SVD (const DenseMatrix &Arg) //////////////////////////// /** * Solves a system of equations using SVD components. * @param u Orthogonal SVD matrix U. * @param w Singular value vector W. * @param v Orthogonal SVD matrix V. * @param b Solution vector. * @param x Variable vector. */ template void SVD::solve ( const DenseVector &b, DenseVector &x) { // Solution of A.x = b // A^-1 * b = V * (1/s) * U' *b int m = U.m(); // rows int n = U.m(); // columns int i,j; Real w; DenseVector tmp(n); for (i=0;i