source: src/main/java/agents/Jama/SingularValueDecomposition.java@ 126

Last change on this file since 126 was 1, checked in by Wouter Pasman, 6 years ago

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