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