001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.structure.jama;
022
023        /** Singular Value Decomposition.
024        <P>
025        For an m-by-n matrix A with m >= n, the singular value decomposition is
026        an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
027        an n-by-n orthogonal matrix V so that A = U*S*V'.
028        <P>
029        The singular values, sigma[k] = S[k][k], are ordered so that
030        sigma[0] >= sigma[1] >= ... >= sigma[n-1].
031        <P>
032        The singular value decompostion always exists, so the constructor will
033        never fail.  The matrix condition number and the effective numerical
034        rank can be computed from this decomposition.
035        */
036
037public class SingularValueDecomposition implements java.io.Serializable {
038
039         static final long serialVersionUID = 640239472093534756l;
040
041/* ------------------------
042        Class variables
043 * ------------------------ */
044
045        /** Arrays for internal storage of U and V.
046        @serial internal storage of U.
047        @serial internal storage of V.
048        */
049        private double[][] U, V;
050
051        /** Array for internal storage of singular values.
052        @serial internal storage of singular values.
053        */
054        private double[] s;
055
056        /** Row and column dimensions.
057        @serial row dimension.
058        @serial column dimension.
059        */
060        private int m, n;
061
062/* ------------------------
063        Constructor
064 * ------------------------ */
065
066        /** Construct the singular value decomposition. Provides a data structure to access U, S and V.
067        @param Arg    Rectangular matrix
068        */
069
070        public SingularValueDecomposition (Matrix Arg) {
071
072                // Derived from LINPACK code.
073                // Initialize.
074                double[][] A = Arg.getArrayCopy();
075                m = Arg.getRowDimension();
076                n = Arg.getColumnDimension();
077
078                /* Apparently the failing cases are only a proper subset of (m<n),
079         so let's not throw error.  Correct fix to come later?
080                if (m<n) {
081          throw new IllegalArgumentException("Jama SVD only works for m >= n"); }
082                */
083                int nu = Math.min(m,n);
084                s = new double [Math.min(m+1,n)];
085                U = new double [m][nu];
086                V = new double [n][n];
087                double[] e = new double [n];
088                double[] work = new double [m];
089                boolean wantu = true;
090                boolean wantv = true;
091
092                // Reduce A to bidiagonal form, storing the diagonal elements
093                // in s and the super-diagonal elements in e.
094
095                int nct = Math.min(m-1,n);
096                int nrt = Math.max(0,Math.min(n-2,m));
097                for (int k = 0; k < Math.max(nct,nrt); k++) {
098                        if (k < nct) {
099
100                                // Compute the transformation for the k-th column and
101                                // place the k-th diagonal in s[k].
102                                // Compute 2-norm of k-th column without under/overflow.
103                                s[k] = 0;
104                                for (int i = k; i < m; i++) {
105                                        s[k] = Maths.hypot(s[k],A[i][k]);
106                                }
107                                if (s[k] != 0.0) {
108                                        if (A[k][k] < 0.0) {
109                                                s[k] = -s[k];
110                                        }
111                                        for (int i = k; i < m; i++) {
112                                                A[i][k] /= s[k];
113                                        }
114                                        A[k][k] += 1.0;
115                                }
116                                s[k] = -s[k];
117                        }
118                        for (int j = k+1; j < n; j++) {
119                                if ((k < nct) && (s[k] != 0.0))  {
120
121                                // Apply the transformation.
122
123                                        double t = 0;
124                                        for (int i = k; i < m; i++) {
125                                                t += A[i][k]*A[i][j];
126                                        }
127                                        t = -t/A[k][k];
128                                        for (int i = k; i < m; i++) {
129                                                A[i][j] += t*A[i][k];
130                                        }
131                                }
132
133                                // Place the k-th row of A into e for the
134                                // subsequent calculation of the row transformation.
135
136                                e[j] = A[k][j];
137                        }
138                        if (wantu && (k < nct)) {
139
140                                // Place the transformation in U for subsequent back
141                                // multiplication.
142
143                                for (int i = k; i < m; i++) {
144                                        U[i][k] = A[i][k];
145                                }
146                        }
147                        if (k < nrt) {
148
149                                // Compute the k-th row transformation and place the
150                                // k-th super-diagonal in e[k].
151                                // Compute 2-norm without under/overflow.
152                                e[k] = 0;
153                                for (int i = k+1; i < n; i++) {
154                                        e[k] = Maths.hypot(e[k],e[i]);
155                                }
156                                if (e[k] != 0.0) {
157                                        if (e[k+1] < 0.0) {
158                                                e[k] = -e[k];
159                                        }
160                                        for (int i = k+1; i < n; i++) {
161                                                e[i] /= e[k];
162                                        }
163                                        e[k+1] += 1.0;
164                                }
165                                e[k] = -e[k];
166                                if ((k+1 < m) && (e[k] != 0.0)) {
167
168                                // Apply the transformation.
169
170                                        for (int i = k+1; i < m; i++) {
171                                                work[i] = 0.0;
172                                        }
173                                        for (int j = k+1; j < n; j++) {
174                                                for (int i = k+1; i < m; i++) {
175                                                        work[i] += e[j]*A[i][j];
176                                                }
177                                        }
178                                        for (int j = k+1; j < n; j++) {
179                                                double t = -e[j]/e[k+1];
180                                                for (int i = k+1; i < m; i++) {
181                                                        A[i][j] += t*work[i];
182                                                }
183                                        }
184                                }
185                                if (wantv) {
186
187                                // Place the transformation in V for subsequent
188                                // back multiplication.
189
190                                        for (int i = k+1; i < n; i++) {
191                                                V[i][k] = e[i];
192                                        }
193                                }
194                        }
195                }
196
197                // Set up the final bidiagonal matrix or order p.
198
199                int p = Math.min(n,m+1);
200                if (nct < n) {
201                        s[nct] = A[nct][nct];
202                }
203                if (m < p) {
204                        s[p-1] = 0.0;
205                }
206                if (nrt+1 < p) {
207                        e[nrt] = A[nrt][p-1];
208                }
209                e[p-1] = 0.0;
210
211                // If required, generate U.
212
213                if (wantu) {
214                        for (int j = nct; j < nu; j++) {
215                                for (int i = 0; i < m; i++) {
216                                        U[i][j] = 0.0;
217                                }
218                                U[j][j] = 1.0;
219                        }
220                        for (int k = nct-1; k >= 0; k--) {
221                                if (s[k] != 0.0) {
222                                        for (int j = k+1; j < nu; j++) {
223                                                double t = 0;
224                                                for (int i = k; i < m; i++) {
225                                                        t += U[i][k]*U[i][j];
226                                                }
227                                                t = -t/U[k][k];
228                                                for (int i = k; i < m; i++) {
229                                                        U[i][j] += t*U[i][k];
230                                                }
231                                        }
232                                        for (int i = k; i < m; i++ ) {
233                                                U[i][k] = -U[i][k];
234                                        }
235                                        U[k][k] = 1.0 + U[k][k];
236                                        for (int i = 0; i < k-1; i++) {
237                                                U[i][k] = 0.0;
238                                        }
239                                } else {
240                                        for (int i = 0; i < m; i++) {
241                                                U[i][k] = 0.0;
242                                        }
243                                        U[k][k] = 1.0;
244                                }
245                        }
246                }
247
248                // If required, generate V.
249
250                if (wantv) {
251                        for (int k = n-1; k >= 0; k--) {
252                                if ((k < nrt) && (e[k] != 0.0)) {
253                                        for (int j = k+1; j < nu; j++) {
254                                                double t = 0;
255                                                for (int i = k+1; i < n; i++) {
256                                                        t += V[i][k]*V[i][j];
257                                                }
258                                                t = -t/V[k+1][k];
259                                                for (int i = k+1; i < n; i++) {
260                                                        V[i][j] += t*V[i][k];
261                                                }
262                                        }
263                                }
264                                for (int i = 0; i < n; i++) {
265                                        V[i][k] = 0.0;
266                                }
267                                V[k][k] = 1.0;
268                        }
269                }
270
271                // Main iteration loop for the singular values.
272
273                int pp = p-1;
274                int iter = 0;
275                double eps = Math.pow(2.0,-52.0);
276                double tiny = Math.pow(2.0,-966.0);
277                while (p > 0) {
278                        int k,kase;
279
280                        // Here is where a test for too many iterations would go.
281
282                        // This section of the program inspects for
283                        // negligible elements in the s and e arrays.  On
284                        // completion the variables kase and k are set as follows.
285
286                        // kase = 1     if s(p) and e[k-1] are negligible and k<p
287                        // kase = 2     if s(k) is negligible and k<p
288                        // kase = 3     if e[k-1] is negligible, k<p, and
289                        //              s(k), ..., s(p) are not negligible (qr step).
290                        // kase = 4     if e(p-1) is negligible (convergence).
291
292                        for (k = p-2; k >= -1; k--) {
293                                if (k == -1) {
294                                        break;
295                                }
296                                if (Math.abs(e[k]) <=
297                                                tiny + eps*(Math.abs(s[k]) + Math.abs(s[k+1]))) {
298                                        e[k] = 0.0;
299                                        break;
300                                }
301                        }
302                        if (k == p-2) {
303                                kase = 4;
304                        } else {
305                                int ks;
306                                for (ks = p-1; ks >= k; ks--) {
307                                        if (ks == k) {
308                                                break;
309                                        }
310                                        double t = (ks != p ? Math.abs(e[ks]) : 0.) +
311                                                                  (ks != k+1 ? Math.abs(e[ks-1]) : 0.);
312                                        if (Math.abs(s[ks]) <= tiny + eps*t)  {
313                                                s[ks] = 0.0;
314                                                break;
315                                        }
316                                }
317                                if (ks == k) {
318                                        kase = 3;
319                                } else if (ks == p-1) {
320                                        kase = 1;
321                                } else {
322                                        kase = 2;
323                                        k = ks;
324                                }
325                        }
326                        k++;
327
328                        // Perform the task indicated by kase.
329
330                        switch (kase) {
331
332                                // Deflate negligible s(p).
333
334                                case 1: {
335                                        double f = e[p-2];
336                                        e[p-2] = 0.0;
337                                        for (int j = p-2; j >= k; 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                                                if (j != k) {
343                                                        f = -sn*e[j-1];
344                                                        e[j-1] = cs*e[j-1];
345                                                }
346                                                if (wantv) {
347                                                        for (int i = 0; i < n; i++) {
348                                                                t = cs*V[i][j] + sn*V[i][p-1];
349                                                                V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1];
350                                                                V[i][j] = t;
351                                                        }
352                                                }
353                                        }
354                                }
355                                break;
356
357                                // Split at negligible s(k).
358
359                                case 2: {
360                                        double f = e[k-1];
361                                        e[k-1] = 0.0;
362                                        for (int j = k; j < p; j++) {
363                                                double t = Maths.hypot(s[j],f);
364                                                double cs = s[j]/t;
365                                                double sn = f/t;
366                                                s[j] = t;
367                                                f = -sn*e[j];
368                                                e[j] = cs*e[j];
369                                                if (wantu) {
370                                                        for (int i = 0; i < m; i++) {
371                                                                t = cs*U[i][j] + sn*U[i][k-1];
372                                                                U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1];
373                                                                U[i][j] = t;
374                                                        }
375                                                }
376                                        }
377                                }
378                                break;
379
380                                // Perform one qr step.
381
382                                case 3: {
383
384                                        // Calculate the shift.
385
386                                        double scale = Math.max(Math.max(Math.max(Math.max(
387                                                          Math.abs(s[p-1]),Math.abs(s[p-2])),Math.abs(e[p-2])),
388                                                          Math.abs(s[k])),Math.abs(e[k]));
389                                        double sp = s[p-1]/scale;
390                                        double spm1 = s[p-2]/scale;
391                                        double epm1 = e[p-2]/scale;
392                                        double sk = s[k]/scale;
393                                        double ek = e[k]/scale;
394                                        double b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0;
395                                        double c = (sp*epm1)*(sp*epm1);
396                                        double shift = 0.0;
397                                        if ((b != 0.0) || (c != 0.0)) {
398                                                shift = Math.sqrt(b*b + c);
399                                                if (b < 0.0) {
400                                                        shift = -shift;
401                                                }
402                                                shift = c/(b + shift);
403                                        }
404                                        double f = (sk + sp)*(sk - sp) + shift;
405                                        double g = sk*ek;
406
407                                        // Chase zeros.
408
409                                        for (int j = k; j < p-1; j++) {
410                                                double t = Maths.hypot(f,g);
411                                                double cs = f/t;
412                                                double sn = g/t;
413                                                if (j != k) {
414                                                        e[j-1] = t;
415                                                }
416                                                f = cs*s[j] + sn*e[j];
417                                                e[j] = cs*e[j] - sn*s[j];
418                                                g = sn*s[j+1];
419                                                s[j+1] = cs*s[j+1];
420                                                if (wantv) {
421                                                        for (int i = 0; i < n; i++) {
422                                                                t = cs*V[i][j] + sn*V[i][j+1];
423                                                                V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1];
424                                                                V[i][j] = t;
425                                                        }
426                                                }
427                                                t = Maths.hypot(f,g);
428                                                cs = f/t;
429                                                sn = g/t;
430                                                s[j] = t;
431                                                f = cs*e[j] + sn*s[j+1];
432                                                s[j+1] = -sn*e[j] + cs*s[j+1];
433                                                g = sn*e[j+1];
434                                                e[j+1] = cs*e[j+1];
435                                                if (wantu && (j < m-1)) {
436                                                        for (int i = 0; i < m; i++) {
437                                                                t = cs*U[i][j] + sn*U[i][j+1];
438                                                                U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1];
439                                                                U[i][j] = t;
440                                                        }
441                                                }
442                                        }
443                                        e[p-2] = f;
444                                        iter = iter + 1;
445                                }
446                                break;
447
448                                // Convergence.
449
450                                case 4: {
451
452                                        // Make the singular values positive.
453
454                                        if (s[k] <= 0.0) {
455                                                s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
456                                                if (wantv) {
457                                                        for (int i = 0; i <= pp; i++) {
458                                                                V[i][k] = -V[i][k];
459                                                        }
460                                                }
461                                        }
462
463                                        // Order the singular values.
464
465                                        while (k < pp) {
466                                                if (s[k] >= s[k+1]) {
467                                                        break;
468                                                }
469                                                double t = s[k];
470                                                s[k] = s[k+1];
471                                                s[k+1] = t;
472                                                if (wantv && (k < n-1)) {
473                                                        for (int i = 0; i < n; i++) {
474                                                                t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t;
475                                                        }
476                                                }
477                                                if (wantu && (k < m-1)) {
478                                                        for (int i = 0; i < m; i++) {
479                                                                t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t;
480                                                        }
481                                                }
482                                                k++;
483                                        }
484                                        iter = 0;
485                                        p--;
486                                }
487                                break;
488                        }
489                }
490        }
491
492/* ------------------------
493        Public Methods
494 * ------------------------ */
495
496        /** Return the left singular vectors
497        @return     U
498        */
499
500        public Matrix getU () {
501                return new Matrix(U,m,Math.min(m+1,n));
502        }
503
504        /** Return the right singular vectors
505        @return     V
506        */
507
508        public Matrix getV () {
509                return new Matrix(V,n,n);
510        }
511
512        /** Return the one-dimensional array of singular values
513        @return     diagonal of S.
514        */
515
516        public double[] getSingularValues () {
517                return s;
518        }
519
520        /** Return the diagonal matrix of singular values
521        @return     S
522        */
523
524        public Matrix getS () {
525                Matrix X = new Matrix(n,n);
526                double[][] S = X.getArray();
527                for (int i = 0; i < n; i++) {
528                        for (int j = 0; j < n; j++) {
529                                S[i][j] = 0.0;
530                        }
531                        S[i][i] = this.s[i];
532                }
533                return X;
534        }
535
536        /** Two norm
537        @return     max(S)
538        */
539
540        public double norm2 () {
541                return s[0];
542        }
543
544        /** Two norm condition number
545        @return     max(S)/min(S)
546        */
547
548        public double cond () {
549                return s[0]/s[Math.min(m,n)-1];
550        }
551
552        /** Effective numerical matrix rank
553        @return     Number of nonnegligible singular values.
554        */
555
556        public int rank () {
557                double eps = Math.pow(2.0,-52.0);
558                double tol = Math.max(m,n)*s[0]*eps;
559                int r = 0;
560                for (int i = 0; i < s.length; i++) {
561                        if (s[i] > tol) {
562                                r++;
563                        }
564                }
565                return r;
566        }
567}