Statistics
| Revision:

root / tmp / org.txm.analec.rcp / src / JamaPlus / SingularValueDecomposition.java0 @ 1726

History | View | Annotate | Download (15.8 kB)

1
package JamaPlus;
2
import JamaPlus.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

    
18
public 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
   @param A    Rectangular matrix
47
   @return     Structure to access U, S and V.
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] = Math.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] = Math.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 = Math.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 = Math.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 = Math.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 = Math.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
}