1 | /// \ingroup newmat |
---|
2 | ///@{ |
---|
3 | |
---|
4 | /// \file svd.cpp |
---|
5 | /// Singular value decomposition. |
---|
6 | |
---|
7 | // Copyright (C) 1991,2,3,4,5: R B Davies |
---|
8 | // Updated 17 July, 1995 |
---|
9 | |
---|
10 | #define WANT_MATH |
---|
11 | |
---|
12 | #include "include.h" |
---|
13 | #include "newmatap.h" |
---|
14 | #include "newmatrm.h" |
---|
15 | #include "precisio.h" |
---|
16 | |
---|
17 | #ifdef use_namespace |
---|
18 | namespace NEWMAT { |
---|
19 | #endif |
---|
20 | |
---|
21 | #ifdef DO_REPORT |
---|
22 | #define REPORT { static ExeCounter ExeCount(__LINE__,15); ++ExeCount; } |
---|
23 | #else |
---|
24 | #define REPORT {} |
---|
25 | #endif |
---|
26 | |
---|
27 | |
---|
28 | |
---|
29 | |
---|
30 | void SVD(const Matrix& A, DiagonalMatrix& Q, Matrix& U, Matrix& V, |
---|
31 | bool withU, bool withV) |
---|
32 | // from Wilkinson and Reinsch: "Handbook of Automatic Computation" |
---|
33 | { |
---|
34 | REPORT |
---|
35 | Tracer trace("SVD"); |
---|
36 | Real eps = FloatingPointPrecision::Epsilon(); |
---|
37 | Real tol = FloatingPointPrecision::Minimum()/eps; |
---|
38 | |
---|
39 | int m = A.Nrows(); int n = A.Ncols(); |
---|
40 | if (m<n) |
---|
41 | Throw(ProgramException("Want no. Rows >= no. Cols", A)); |
---|
42 | if (withV && &U == &V) |
---|
43 | Throw(ProgramException("Need different matrices for U and V", U, V)); |
---|
44 | U = A; Real g = 0.0; Real f,h; Real x = 0.0; int i; |
---|
45 | RowVector E(n); RectMatrixRow EI(E,0); Q.ReSize(n); |
---|
46 | RectMatrixCol UCI(U,0); RectMatrixRow URI(U,0,1,n-1); |
---|
47 | |
---|
48 | if (n) for (i=0;;) |
---|
49 | { |
---|
50 | EI.First() = g; Real ei = g; EI.Right(); Real s = UCI.SumSquare(); |
---|
51 | if (s<tol) { REPORT Q.element(i) = 0.0; } |
---|
52 | else |
---|
53 | { |
---|
54 | REPORT |
---|
55 | f = UCI.First(); g = -sign(sqrt(s), f); h = f*g-s; UCI.First() = f-g; |
---|
56 | Q.element(i) = g; RectMatrixCol UCJ = UCI; int j=n-i; |
---|
57 | while (--j) { UCJ.Right(); UCJ.AddScaled(UCI, (UCI*UCJ)/h); } |
---|
58 | } |
---|
59 | |
---|
60 | s = URI.SumSquare(); |
---|
61 | if (s<tol) { REPORT g = 0.0; } |
---|
62 | else |
---|
63 | { |
---|
64 | REPORT |
---|
65 | f = URI.First(); g = -sign(sqrt(s), f); URI.First() = f-g; |
---|
66 | EI.Divide(URI,f*g-s); RectMatrixRow URJ = URI; int j=m-i; |
---|
67 | while (--j) { URJ.Down(); URJ.AddScaled(EI, URI*URJ); } |
---|
68 | } |
---|
69 | |
---|
70 | Real y = fabs(Q.element(i)) + fabs(ei); if (x<y) { REPORT x = y; } |
---|
71 | if (++i == n) { REPORT break; } |
---|
72 | UCI.DownDiag(); URI.DownDiag(); |
---|
73 | } |
---|
74 | |
---|
75 | if (withV) |
---|
76 | { |
---|
77 | REPORT |
---|
78 | V.ReSize(n,n); V = 0.0; RectMatrixCol VCI(V,n-1,n-1,1); |
---|
79 | if (n) { VCI.First() = 1.0; g=E.element(n-1); if (n!=1) URI.UpDiag(); } |
---|
80 | for (i=n-2; i>=0; i--) |
---|
81 | { |
---|
82 | VCI.Left(); |
---|
83 | if (g!=0.0) |
---|
84 | { |
---|
85 | VCI.Divide(URI, URI.First()*g); int j = n-i; |
---|
86 | RectMatrixCol VCJ = VCI; |
---|
87 | while (--j) { VCJ.Right(); VCJ.AddScaled( VCI, (URI*VCJ) ); } |
---|
88 | } |
---|
89 | VCI.Zero(); VCI.Up(); VCI.First() = 1.0; g=E.element(i); |
---|
90 | if (i==0) break; |
---|
91 | URI.UpDiag(); |
---|
92 | } |
---|
93 | } |
---|
94 | |
---|
95 | if (withU) |
---|
96 | { |
---|
97 | REPORT |
---|
98 | for (i=n-1; i>=0; i--) |
---|
99 | { |
---|
100 | g = Q.element(i); URI.Reset(U,i,i+1,n-i-1); URI.Zero(); |
---|
101 | if (g!=0.0) |
---|
102 | { |
---|
103 | h=UCI.First()*g; int j=n-i; RectMatrixCol UCJ = UCI; |
---|
104 | while (--j) |
---|
105 | { |
---|
106 | UCJ.Right(); UCI.Down(); UCJ.Down(); Real s = UCI*UCJ; |
---|
107 | UCI.Up(); UCJ.Up(); UCJ.AddScaled(UCI,s/h); |
---|
108 | } |
---|
109 | UCI.Divide(g); |
---|
110 | } |
---|
111 | else UCI.Zero(); |
---|
112 | UCI.First() += 1.0; |
---|
113 | if (i==0) break; |
---|
114 | UCI.UpDiag(); |
---|
115 | } |
---|
116 | } |
---|
117 | |
---|
118 | eps *= x; |
---|
119 | for (int k=n-1; k>=0; k--) |
---|
120 | { |
---|
121 | Real z = -FloatingPointPrecision::Maximum(); // to keep Gnu happy |
---|
122 | Real y; int limit = 50; int l = 0; |
---|
123 | while (limit--) |
---|
124 | { |
---|
125 | Real c, s; int i; int l1=k; bool tfc=false; |
---|
126 | for (l=k; l>=0; l--) |
---|
127 | { |
---|
128 | // if (fabs(E.element(l))<=eps) goto test_f_convergence; |
---|
129 | if (fabs(E.element(l))<=eps) { REPORT tfc=true; break; } |
---|
130 | if (fabs(Q.element(l-1))<=eps) { REPORT l1=l; break; } |
---|
131 | REPORT |
---|
132 | } |
---|
133 | if (!tfc) |
---|
134 | { |
---|
135 | REPORT |
---|
136 | l=l1; l1=l-1; s = -1.0; c = 0.0; |
---|
137 | for (i=l; i<=k; i++) |
---|
138 | { |
---|
139 | f = - s * E.element(i); E.element(i) *= c; |
---|
140 | // if (fabs(f)<=eps) goto test_f_convergence; |
---|
141 | if (fabs(f)<=eps) { REPORT break; } |
---|
142 | g = Q.element(i); h = pythag(g,f,c,s); Q.element(i) = h; |
---|
143 | if (withU) |
---|
144 | { |
---|
145 | REPORT |
---|
146 | RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,l1); |
---|
147 | ComplexScale(UCJ, UCI, c, s); |
---|
148 | } |
---|
149 | } |
---|
150 | } |
---|
151 | // test_f_convergence: z = Q.element(k); if (l==k) goto convergence; |
---|
152 | z = Q.element(k); if (l==k) { REPORT break; } |
---|
153 | |
---|
154 | x = Q.element(l); y = Q.element(k-1); |
---|
155 | g = E.element(k-1); h = E.element(k); |
---|
156 | f = ((y-z)*(y+z) + (g-h)*(g+h)) / (2*h*y); |
---|
157 | if (f>1) { REPORT g = f * sqrt(1 + square(1/f)); } |
---|
158 | else if (f<-1) { REPORT g = -f * sqrt(1 + square(1/f)); } |
---|
159 | else { REPORT g = sqrt(f*f + 1); } |
---|
160 | { REPORT f = ((x-z)*(x+z) + h*(y / ((f<0.0) ? f-g : f+g)-h)) / x; } |
---|
161 | |
---|
162 | c = 1.0; s = 1.0; |
---|
163 | for (i=l+1; i<=k; i++) |
---|
164 | { |
---|
165 | g = E.element(i); y = Q.element(i); h = s*g; g *= c; |
---|
166 | z = pythag(f,h,c,s); E.element(i-1) = z; |
---|
167 | f = x*c + g*s; g = -x*s + g*c; h = y*s; y *= c; |
---|
168 | if (withV) |
---|
169 | { |
---|
170 | REPORT |
---|
171 | RectMatrixCol VCI(V,i); RectMatrixCol VCJ(V,i-1); |
---|
172 | ComplexScale(VCI, VCJ, c, s); |
---|
173 | } |
---|
174 | z = pythag(f,h,c,s); Q.element(i-1) = z; |
---|
175 | f = c*g + s*y; x = -s*g + c*y; |
---|
176 | if (withU) |
---|
177 | { |
---|
178 | REPORT |
---|
179 | RectMatrixCol UCI(U,i); RectMatrixCol UCJ(U,i-1); |
---|
180 | ComplexScale(UCI, UCJ, c, s); |
---|
181 | } |
---|
182 | } |
---|
183 | E.element(l) = 0.0; E.element(k) = f; Q.element(k) = x; |
---|
184 | } |
---|
185 | if (l!=k) { Throw(ConvergenceException(A)); } |
---|
186 | // convergence: |
---|
187 | if (z < 0.0) |
---|
188 | { |
---|
189 | REPORT |
---|
190 | Q.element(k) = -z; |
---|
191 | if (withV) { RectMatrixCol VCI(V,k); VCI.Negate(); } |
---|
192 | } |
---|
193 | } |
---|
194 | if (withU & withV) SortSV(Q, U, V); |
---|
195 | else if (withU) SortSV(Q, U); |
---|
196 | else if (withV) SortSV(Q, V); |
---|
197 | else sort_descending(Q); |
---|
198 | } |
---|
199 | |
---|
200 | void SVD(const Matrix& A, DiagonalMatrix& D) |
---|
201 | { REPORT Matrix U; SVD(A, D, U, U, false, false); } |
---|
202 | |
---|
203 | |
---|
204 | |
---|
205 | #ifdef use_namespace |
---|
206 | } |
---|
207 | #endif |
---|
208 | |
---|
209 | ///@} |
---|