SVD.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
SVD.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$svd.cpp                           singular value decomposition
00003 
00004 // Copyright (C) 1991,2,3,4,5: R B Davies
00005 // Updated 17 July, 1995
00006 
00007 #define WANT_MATH
00008 
00009 #include "include.h"
00010 #include "newmat.h"
00011 #include "newmatrm.h"
00012 #include "precisio.h"
00013 
00014 #ifdef use_namespace
00015 namespace NEWMAT {
00016 #endif
00017 
00018     static Real pythag(Real f, Real g, Real& c, Real& s) {
00019         // return z=sqrt(f*f+g*g), c=f/z, s=g/z
00020         // set c=1,s=0 if z==0
00021         // avoid floating point overflow or divide by zero
00022         if (f==0 && g==0) { c=1.0; s=0.0; return 0.0; }
00023         Real af = f>=0 ? f : -f;
00024         Real ag = g>=0 ? g : -g;
00025         if (ag<af) {
00026             Real h = g/f; Real sq = sqrt(1.0+h*h);
00027             if (f<0) sq = -sq;                          // make return value non-negative
00028             c = 1.0/sq; s = h/sq; return sq*f;
00029         }
00030         else {
00031             Real h = f/g; Real sq = sqrt(1.0+h*h);
00032             if (g<0) sq = -sq;
00033             s = 1.0/sq; c = h/sq; return sq*g;
00034         }
00035     }
00036 
00037     void SVD(const Matrix& A, DiagonalMatrix& Q, Matrix& U, Matrix& V,
00038              bool withU, bool withV) {
00039         // from Wilkinson and Reinsch: "Handbook of Automatic Computation"
00040         Tracer trace("SVD");
00041         Real eps = FloatingPointPrecision::Epsilon();
00042         Real tol = FloatingPointPrecision::Minimum()/eps;
00043 
00044         int m = A.Nrows(); int n = A.Ncols();
00045         if (m<n)
00046             Throw(ProgramException("Want no. Rows >= no. Cols", A));
00047 
00048         U = A; Real g = 0.0; Real f,h; Real x = 0.0; int i;
00049         RowVector E(n); RectMatrixRow EI(E,0); Q.ReSize(n);
00050         RectMatrixCol UCI(U,0); RectMatrixRow URI(U,0,1,n-1);
00051 
00052         for (i=0; i<n; i++) {
00053             EI.First() = g; Real ei = g; EI.Right(); Real s = UCI.SumSquare();
00054             if (s<tol) Q.element(i) = 0.0;
00055             else {
00056                 f = UCI.First(); g = -sign(sqrt(s), f); h = f*g-s; UCI.First() = f-g;
00057                 Q.element(i) = g; RectMatrixCol UCJ = UCI; int j=n-i;
00058                 while (--j) { UCJ.Right(); UCJ.AddScaled(UCI, (UCI*UCJ)/h); }
00059             }
00060 
00061             s = URI.SumSquare();
00062             if (s<tol) g = 0.0;
00063             else {
00064                 f = URI.First(); g = -sign(sqrt(s), f); URI.First() = f-g;
00065                 EI.Divide(URI,f*g-s); RectMatrixRow URJ = URI; int j=m-i;
00066                 while (--j) { URJ.Down(); URJ.AddScaled(EI, URI*URJ); }
00067             }
00068 
00069             Real y = fabs(Q.element(i)) + fabs(ei); if (x<y) x = y;
00070             UCI.DownDiag(); URI.DownDiag();
00071         }
00072 
00073         if (withV) {
00074             V.ReSize(n,n); V = 0.0; RectMatrixCol VCI(V,n,n,0);
00075             for (i=n-1; i>=0; i--) {
00076                 URI.UpDiag(); VCI.Left();
00077                 if (g!=0.0) {
00078                     VCI.Divide(URI, URI.First()*g); int j = n-i;
00079                     RectMatrixCol VCJ = VCI;
00080                     while (--j) { VCJ.Right(); VCJ.AddScaled( VCI, (URI*VCJ) ); }
00081                 }
00082                 VCI.Zero(); VCI.Up(); VCI.First() = 1.0; g=E.element(i);
00083             }
00084         }
00085 
00086         if (withU) {
00087             for (i=n-1; i>=0; i--) {
00088                 UCI.UpDiag(); g = Q.element(i); URI.Reset(U,i,i+1,n-i-1); URI.Zero();
00089                 if (g!=0.0) {
00090                     h=UCI.First()*g; int j=n-i; RectMatrixCol UCJ = UCI;
00091                     while (--j) {
00092                         UCJ.Right(); UCI.Down(); UCJ.Down(); Real s = UCI*UCJ;
00093                         UCI.Up(); UCJ.Up(); UCJ.AddScaled(UCI,s/h);
00094                     }
00095                     UCI.Divide(g);
00096                 }
00097                 else UCI.Zero();
00098                 UCI.First() += 1.0;
00099             }
00100         }
00101 
00102         eps *= x;
00103         for (int k=n-1; k>=0; k--) {
00104             Real z; Real y; int limit = 50; int l;
00105             while (limit--) {
00106                 Real c, s; int i; int l1=k; bool tfc=false;
00107                 for (l=k; l>=0; l--) {
00108                     //          if (fabs(E.element(l))<=eps) goto test_f_convergence;
00109                     if (fabs(E.element(l))<=eps) { tfc=true; break; }
00110                     if (fabs(Q.element(l-1))<=eps) { l1=l; break; }
00111                 }
00112                 if (!tfc) {
00113                     l=l1; l1=l-1; s = -1.0; c = 0.0;
00114                     for (i=l; i<=k; i++) {
00115                         f = - s * E.element(i); E.element(i) *= c;
00116                         //             if (fabs(f)<=eps) goto test_f_convergence;
00117                         if (fabs(f)<=eps) break;
00118                         g = Q.element(i); h = pythag(g,f,c,s); Q.element(i) = h;
00119                         if (withU) {
00120                             RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,l1);
00121                             ComplexScale(UCJ, UCI, c, s);
00122                         }
00123                     }
00124                 }
00125                 //    test_f_convergence: z = Q.element(k); if (l==k) goto convergence;
00126                 z = Q.element(k);  if (l==k) break;
00127 
00128                 x = Q.element(l); y = Q.element(k-1);
00129                 g = E.element(k-1); h = E.element(k);
00130                 f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y);
00131                 if (f>1)         g = f * sqrt(1 + square(1/f));
00132                 else if (f<-1)   g = -f * sqrt(1 + square(1/f));
00133                 else             g = sqrt(f*f + 1);
00134                 f = ((x-z)*(x+z) + h*(y / ((f<0.0) ? f-g : f+g)-h)) / x;
00135 
00136                 c = 1.0; s = 1.0;
00137                 for (i=l+1; i<=k; i++) {
00138                     g = E.element(i); y = Q.element(i); h = s*g; g *= c;
00139                     z = pythag(f,h,c,s); E.element(i-1) = z;
00140                     f = x*c + g*s; g = -x*s + g*c; h = y*s; y *= c;
00141                     if (withV) {
00142                         RectMatrixCol VCI(V,i); RectMatrixCol VCJ(V,i-1);
00143                         ComplexScale(VCI, VCJ, c, s);
00144                     }
00145                     z = pythag(f,h,c,s); Q.element(i-1) = z;
00146                     f = c*g + s*y; x = -s*g + c*y;
00147                     if (withU) {
00148                         RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,i-1);
00149                         ComplexScale(UCI, UCJ, c, s);
00150                     }
00151                 }
00152                 E.element(l) = 0.0; E.element(k) = f; Q.element(k) = x;
00153             }
00154             if (l!=k) { Throw(ConvergenceException(A)); }
00155             // convergence:
00156             if (z < 0.0) {
00157                 Q.element(k) = -z;
00158                 if (withV) { RectMatrixCol VCI(V,k); VCI.Negate(); }
00159             }
00160         }
00161     }
00162 
00163     void SVD(const Matrix& A, DiagonalMatrix& D)
00164     { Matrix U; SVD(A, D, U, U, false, false); }
00165 
00166 #ifdef use_namespace
00167 }
00168 #endif
Generated on Fri Aug 23 01:38:38 2002 for VirtualU by
1.2.17