1 |
#ifndef JAMA_SVD_H |
2 |
#define JAMA_SVD_H |
3 |
|
4 |
#include "math/DynamicRectMatrix.hpp" |
5 |
|
6 |
#include <algorithm> |
7 |
// for min(), max() below |
8 |
#include <cmath> |
9 |
// for abs() below |
10 |
|
11 |
using namespace oopse; |
12 |
using namespace std; |
13 |
|
14 |
namespace JAMA |
15 |
{ |
16 |
/** Singular Value Decomposition. |
17 |
<P> |
18 |
For an m-by-n matrix A with m >= n, the singular value decomposition is |
19 |
an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and |
20 |
an n-by-n orthogonal matrix V so that A = U*S*V'. |
21 |
<P> |
22 |
The singular values, sigma(k) = S(k,k), are ordered so that |
23 |
sigma(0) >= sigma(1) >= ... >= sigma(n-1). |
24 |
<P> |
25 |
The singular value decompostion always exists, so the constructor will |
26 |
never fail. The matrix condition number and the effective numerical |
27 |
rank can be computed from this decomposition. |
28 |
|
29 |
<p> |
30 |
(Adapted from JAMA, a Java Matrix Library, developed by jointly |
31 |
by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). |
32 |
*/ |
33 |
template <class Real> |
34 |
class SVD |
35 |
{ |
36 |
|
37 |
DynamicRectMatrix<Real> U, V; |
38 |
DynamicVector<Real> s; |
39 |
int m, n; |
40 |
|
41 |
public: |
42 |
|
43 |
|
44 |
SVD (const DynamicRectMatrix<Real> &Arg) { |
45 |
|
46 |
m = Arg.getNRow(); |
47 |
n = Arg.getNCol(); |
48 |
int nu = min(m,n); |
49 |
s = DynamicVector<Real>(min(m+1,n)); |
50 |
U = DynamicRectMatrix<Real>(m, nu, Real(0)); |
51 |
V = DynamicRectMatrix<Real>(n,n); |
52 |
DynamicVector<Real> e(n); |
53 |
DynamicVector<Real> work(m); |
54 |
DynamicRectMatrix<Real> A(Arg); |
55 |
int wantu = 1; /* boolean */ |
56 |
int wantv = 1; /* boolean */ |
57 |
int i=0, j=0, k=0; |
58 |
|
59 |
// Reduce A to bidiagonal form, storing the diagonal elements |
60 |
// in s and the super-diagonal elements in e. |
61 |
|
62 |
int nct = min(m-1,n); |
63 |
int nrt = max(0,min(n-2,m)); |
64 |
for (k = 0; k < max(nct,nrt); k++) { |
65 |
if (k < nct) { |
66 |
|
67 |
// Compute the transformation for the k-th column and |
68 |
// place the k-th diagonal in s(k). |
69 |
// Compute 2-norm of k-th column without under/overflow. |
70 |
s(k) = 0; |
71 |
for (i = k; i < m; i++) { |
72 |
s(k) = hypot(s(k),A(i,k)); |
73 |
} |
74 |
if (s(k) != 0.0) { |
75 |
if (A(k,k) < 0.0) { |
76 |
s(k) = -s(k); |
77 |
} |
78 |
for (i = k; i < m; i++) { |
79 |
A(i,k) /= s(k); |
80 |
} |
81 |
A(k,k) += 1.0; |
82 |
} |
83 |
s(k) = -s(k); |
84 |
} |
85 |
for (j = k+1; j < n; j++) { |
86 |
if ((k < nct) && (s(k) != 0.0)) { |
87 |
|
88 |
// Apply the transformation. |
89 |
|
90 |
Real t(0.0); |
91 |
for (i = k; i < m; i++) { |
92 |
t += A(i,k)*A(i,j); |
93 |
} |
94 |
t = -t/A(k,k); |
95 |
for (i = k; i < m; i++) { |
96 |
A(i,j) += t*A(i,k); |
97 |
} |
98 |
} |
99 |
|
100 |
// Place the k-th row of A into e for the |
101 |
// subsequent calculation of the row transformation. |
102 |
|
103 |
e(j) = A(k,j); |
104 |
} |
105 |
if (wantu & (k < nct)) { |
106 |
|
107 |
// Place the transformation in U for subsequent back |
108 |
// multiplication. |
109 |
|
110 |
for (i = k; i < m; i++) { |
111 |
U(i,k) = A(i,k); |
112 |
} |
113 |
} |
114 |
if (k < nrt) { |
115 |
|
116 |
// Compute the k-th row transformation and place the |
117 |
// k-th super-diagonal in e(k). |
118 |
// Compute 2-norm without under/overflow. |
119 |
e(k) = 0; |
120 |
for (i = k+1; i < n; i++) { |
121 |
e(k) = hypot(e(k),e(i)); |
122 |
} |
123 |
if (e(k) != 0.0) { |
124 |
if (e(k+1) < 0.0) { |
125 |
e(k) = -e(k); |
126 |
} |
127 |
for (i = k+1; i < n; i++) { |
128 |
e(i) /= e(k); |
129 |
} |
130 |
e(k+1) += 1.0; |
131 |
} |
132 |
e(k) = -e(k); |
133 |
if ((k+1 < m) & (e(k) != 0.0)) { |
134 |
|
135 |
// Apply the transformation. |
136 |
|
137 |
for (i = k+1; i < m; i++) { |
138 |
work(i) = 0.0; |
139 |
} |
140 |
for (j = k+1; j < n; j++) { |
141 |
for (i = k+1; i < m; i++) { |
142 |
work(i) += e(j)*A(i,j); |
143 |
} |
144 |
} |
145 |
for (j = k+1; j < n; j++) { |
146 |
Real t(-e(j)/e(k+1)); |
147 |
for (i = k+1; i < m; i++) { |
148 |
A(i,j) += t*work(i); |
149 |
} |
150 |
} |
151 |
} |
152 |
if (wantv) { |
153 |
|
154 |
// Place the transformation in V for subsequent |
155 |
// back multiplication. |
156 |
|
157 |
for (i = k+1; i < n; i++) { |
158 |
V(i,k) = e(i); |
159 |
} |
160 |
} |
161 |
} |
162 |
} |
163 |
|
164 |
// Set up the final bidiagonal matrix or order p. |
165 |
|
166 |
int p = min(n,m+1); |
167 |
if (nct < n) { |
168 |
s(nct) = A(nct,nct); |
169 |
} |
170 |
if (m < p) { |
171 |
s(p-1) = 0.0; |
172 |
} |
173 |
if (nrt+1 < p) { |
174 |
e(nrt) = A(nrt,p-1); |
175 |
} |
176 |
e(p-1) = 0.0; |
177 |
|
178 |
// If required, generate U. |
179 |
|
180 |
if (wantu) { |
181 |
for (j = nct; j < nu; j++) { |
182 |
for (i = 0; i < m; i++) { |
183 |
U(i,j) = 0.0; |
184 |
} |
185 |
U(j,j) = 1.0; |
186 |
} |
187 |
for (k = nct-1; k >= 0; k--) { |
188 |
if (s(k) != 0.0) { |
189 |
for (j = k+1; j < nu; j++) { |
190 |
Real t(0.0); |
191 |
for (i = k; i < m; i++) { |
192 |
t += U(i,k)*U(i,j); |
193 |
} |
194 |
t = -t/U(k,k); |
195 |
for (i = k; i < m; i++) { |
196 |
U(i,j) += t*U(i,k); |
197 |
} |
198 |
} |
199 |
for (i = k; i < m; i++ ) { |
200 |
U(i,k) = -U(i,k); |
201 |
} |
202 |
U(k,k) = 1.0 + U(k,k); |
203 |
for (i = 0; i < k-1; i++) { |
204 |
U(i,k) = 0.0; |
205 |
} |
206 |
} else { |
207 |
for (i = 0; i < m; i++) { |
208 |
U(i,k) = 0.0; |
209 |
} |
210 |
U(k,k) = 1.0; |
211 |
} |
212 |
} |
213 |
} |
214 |
|
215 |
// If required, generate V. |
216 |
|
217 |
if (wantv) { |
218 |
for (k = n-1; k >= 0; k--) { |
219 |
if ((k < nrt) & (e(k) != 0.0)) { |
220 |
for (j = k+1; j < nu; j++) { |
221 |
Real t(0.0); |
222 |
for (i = k+1; i < n; i++) { |
223 |
t += V(i,k)*V(i,j); |
224 |
} |
225 |
t = -t/V(k+1,k); |
226 |
for (i = k+1; i < n; i++) { |
227 |
V(i,j) += t*V(i,k); |
228 |
} |
229 |
} |
230 |
} |
231 |
for (i = 0; i < n; i++) { |
232 |
V(i,k) = 0.0; |
233 |
} |
234 |
V(k,k) = 1.0; |
235 |
} |
236 |
} |
237 |
|
238 |
// Main iteration loop for the singular values. |
239 |
|
240 |
int pp = p-1; |
241 |
int iter = 0; |
242 |
Real eps(pow(2.0,-52.0)); |
243 |
while (p > 0) { |
244 |
int k=0; |
245 |
int kase=0; |
246 |
|
247 |
// Here is where a test for too many iterations would go. |
248 |
|
249 |
// This section of the program inspects for |
250 |
// negligible elements in the s and e arrays. On |
251 |
// completion the variables kase and k are set as follows. |
252 |
|
253 |
// kase = 1 if s(p) and e(k-1) are negligible and k<p |
254 |
// kase = 2 if s(k) is negligible and k<p |
255 |
// kase = 3 if e(k-1) is negligible, k<p, and |
256 |
// s(k), ..., s(p) are not negligible (qr step). |
257 |
// kase = 4 if e(p-1) is negligible (convergence). |
258 |
|
259 |
for (k = p-2; k >= -1; k--) { |
260 |
if (k == -1) { |
261 |
break; |
262 |
} |
263 |
if (abs(e(k)) <= eps*(abs(s(k)) + abs(s(k+1)))) { |
264 |
e(k) = 0.0; |
265 |
break; |
266 |
} |
267 |
} |
268 |
if (k == p-2) { |
269 |
kase = 4; |
270 |
} else { |
271 |
int ks; |
272 |
for (ks = p-1; ks >= k; ks--) { |
273 |
if (ks == k) { |
274 |
break; |
275 |
} |
276 |
Real t( (ks != p ? abs(e(ks)) : 0.) + |
277 |
(ks != k+1 ? abs(e(ks-1)) : 0.)); |
278 |
if (abs(s(ks)) <= eps*t) { |
279 |
s(ks) = 0.0; |
280 |
break; |
281 |
} |
282 |
} |
283 |
if (ks == k) { |
284 |
kase = 3; |
285 |
} else if (ks == p-1) { |
286 |
kase = 1; |
287 |
} else { |
288 |
kase = 2; |
289 |
k = ks; |
290 |
} |
291 |
} |
292 |
k++; |
293 |
|
294 |
// Perform the task indicated by kase. |
295 |
|
296 |
switch (kase) { |
297 |
|
298 |
// Deflate negligible s(p). |
299 |
|
300 |
case 1: { |
301 |
Real f(e(p-2)); |
302 |
e(p-2) = 0.0; |
303 |
for (j = p-2; j >= k; j--) { |
304 |
Real t( hypot(s(j),f)); |
305 |
Real cs(s(j)/t); |
306 |
Real sn(f/t); |
307 |
s(j) = t; |
308 |
if (j != k) { |
309 |
f = -sn*e(j-1); |
310 |
e(j-1) = cs*e(j-1); |
311 |
} |
312 |
if (wantv) { |
313 |
for (i = 0; i < n; i++) { |
314 |
t = cs*V(i,j) + sn*V(i,p-1); |
315 |
V(i,p-1) = -sn*V(i,j) + cs*V(i,p-1); |
316 |
V(i,j) = t; |
317 |
} |
318 |
} |
319 |
} |
320 |
} |
321 |
break; |
322 |
|
323 |
// Split at negligible s(k). |
324 |
|
325 |
case 2: { |
326 |
Real f(e(k-1)); |
327 |
e(k-1) = 0.0; |
328 |
for (j = k; j < p; j++) { |
329 |
Real t(hypot(s(j),f)); |
330 |
Real cs( s(j)/t); |
331 |
Real sn(f/t); |
332 |
s(j) = t; |
333 |
f = -sn*e(j); |
334 |
e(j) = cs*e(j); |
335 |
if (wantu) { |
336 |
for (i = 0; i < m; i++) { |
337 |
t = cs*U(i,j) + sn*U(i,k-1); |
338 |
U(i,k-1) = -sn*U(i,j) + cs*U(i,k-1); |
339 |
U(i,j) = t; |
340 |
} |
341 |
} |
342 |
} |
343 |
} |
344 |
break; |
345 |
|
346 |
// Perform one qr step. |
347 |
|
348 |
case 3: { |
349 |
|
350 |
// Calculate the shift. |
351 |
|
352 |
Real scale = max(max(max(max( |
353 |
abs(s(p-1)),abs(s(p-2))),abs(e(p-2))), |
354 |
abs(s(k))),abs(e(k))); |
355 |
Real sp = s(p-1)/scale; |
356 |
Real spm1 = s(p-2)/scale; |
357 |
Real epm1 = e(p-2)/scale; |
358 |
Real sk = s(k)/scale; |
359 |
Real ek = e(k)/scale; |
360 |
Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; |
361 |
Real c = (sp*epm1)*(sp*epm1); |
362 |
Real shift = 0.0; |
363 |
if ((b != 0.0) || (c != 0.0)) { |
364 |
shift = sqrt(b*b + c); |
365 |
if (b < 0.0) { |
366 |
shift = -shift; |
367 |
} |
368 |
shift = c/(b + shift); |
369 |
} |
370 |
Real f = (sk + sp)*(sk - sp) + shift; |
371 |
Real g = sk*ek; |
372 |
|
373 |
// Chase zeros. |
374 |
|
375 |
for (j = k; j < p-1; j++) { |
376 |
Real t = hypot(f,g); |
377 |
Real cs = f/t; |
378 |
Real sn = g/t; |
379 |
if (j != k) { |
380 |
e(j-1) = t; |
381 |
} |
382 |
f = cs*s(j) + sn*e(j); |
383 |
e(j) = cs*e(j) - sn*s(j); |
384 |
g = sn*s(j+1); |
385 |
s(j+1) = cs*s(j+1); |
386 |
if (wantv) { |
387 |
for (i = 0; i < n; i++) { |
388 |
t = cs*V(i,j) + sn*V(i,j+1); |
389 |
V(i,j+1) = -sn*V(i,j) + cs*V(i,j+1); |
390 |
V(i,j) = t; |
391 |
} |
392 |
} |
393 |
t = hypot(f,g); |
394 |
cs = f/t; |
395 |
sn = g/t; |
396 |
s(j) = t; |
397 |
f = cs*e(j) + sn*s(j+1); |
398 |
s(j+1) = -sn*e(j) + cs*s(j+1); |
399 |
g = sn*e(j+1); |
400 |
e(j+1) = cs*e(j+1); |
401 |
if (wantu && (j < m-1)) { |
402 |
for (i = 0; i < m; i++) { |
403 |
t = cs*U(i,j) + sn*U(i,j+1); |
404 |
U(i,j+1) = -sn*U(i,j) + cs*U(i,j+1); |
405 |
U(i,j) = t; |
406 |
} |
407 |
} |
408 |
} |
409 |
e(p-2) = f; |
410 |
iter = iter + 1; |
411 |
} |
412 |
break; |
413 |
|
414 |
// Convergence. |
415 |
|
416 |
case 4: { |
417 |
|
418 |
// Make the singular values positive. |
419 |
|
420 |
if (s(k) <= 0.0) { |
421 |
s(k) = (s(k) < 0.0 ? -s(k) : 0.0); |
422 |
if (wantv) { |
423 |
for (i = 0; i <= pp; i++) { |
424 |
V(i,k) = -V(i,k); |
425 |
} |
426 |
} |
427 |
} |
428 |
|
429 |
// Order the singular values. |
430 |
|
431 |
while (k < pp) { |
432 |
if (s(k) >= s(k+1)) { |
433 |
break; |
434 |
} |
435 |
Real t = s(k); |
436 |
s(k) = s(k+1); |
437 |
s(k+1) = t; |
438 |
if (wantv && (k < n-1)) { |
439 |
for (i = 0; i < n; i++) { |
440 |
t = V(i,k+1); V(i,k+1) = V(i,k); V(i,k) = t; |
441 |
} |
442 |
} |
443 |
if (wantu && (k < m-1)) { |
444 |
for (i = 0; i < m; i++) { |
445 |
t = U(i,k+1); U(i,k+1) = U(i,k); U(i,k) = t; |
446 |
} |
447 |
} |
448 |
k++; |
449 |
} |
450 |
iter = 0; |
451 |
p--; |
452 |
} |
453 |
break; |
454 |
} |
455 |
} |
456 |
} |
457 |
|
458 |
|
459 |
void getU (DynamicRectMatrix<Real> &A) { |
460 |
|
461 |
int minm = min(m+1,n); |
462 |
|
463 |
A = DynamicRectMatrix<Real>(m, minm); |
464 |
|
465 |
for (int i=0; i<m; i++) |
466 |
for (int j=0; j<minm; j++) |
467 |
A(i,j) = U(i,j); |
468 |
} |
469 |
|
470 |
/* Return the right singular vectors */ |
471 |
void getV (DynamicRectMatrix<Real> &A) { |
472 |
A = V; |
473 |
} |
474 |
|
475 |
/** Return the one-dimensional array of singular values */ |
476 |
void getSingularValues (DynamicVector<Real> &x) { |
477 |
x = s; |
478 |
} |
479 |
|
480 |
/** Return the diagonal matrix of singular values |
481 |
@return S |
482 |
*/ |
483 |
void getS (DynamicRectMatrix<Real> &A) { |
484 |
A = DynamicRectMatrix<Real>(n,n); |
485 |
for (int i = 0; i < n; i++) { |
486 |
for (int j = 0; j < n; j++) { |
487 |
A(i,j) = 0.0; |
488 |
} |
489 |
A(i,i) = s(i); |
490 |
} |
491 |
} |
492 |
|
493 |
/** Two norm (max(S)) */ |
494 |
Real norm2 () { |
495 |
return s(0); |
496 |
} |
497 |
|
498 |
/** Two norm of condition number (max(S)/min(S)) */ |
499 |
Real cond () { |
500 |
return s(0)/s(min(m,n)-1); |
501 |
} |
502 |
|
503 |
/** Effective numerical matrix rank |
504 |
@return Number of nonnegligible singular values. |
505 |
*/ |
506 |
int rank () { |
507 |
Real eps = pow(2.0,-52.0); |
508 |
Real tol = max(m,n)*s(0)*eps; |
509 |
int r = 0; |
510 |
for (int i = 0; i < s.dim(); i++) { |
511 |
if (s(i) > tol) { |
512 |
r++; |
513 |
} |
514 |
} |
515 |
return r; |
516 |
} |
517 |
}; |
518 |
} |
519 |
#endif |
520 |
// JAMA_SVD_H |