View Javadoc

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