View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math.linear;
19  
20  import org.apache.commons.math.ConvergenceException;
21  import org.apache.commons.math.MathRuntimeException;
22  import org.apache.commons.math.util.MathUtils;
23  
24  /**
25   * Calculates the Singular Value Decomposition of a matrix.
26   * <p>The Singular Value Decomposition of matrix A is a set of three matrices:
27   * U, &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>.
28   * Let A be an m &times; n matrix, then U is an m &times; m orthogonal matrix,
29   * &Sigma; is a m &times; n diagonal matrix with positive diagonal elements,
30   * and V is an n &times; n orthogonal matrix.</p>
31   *
32   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
33   * @since 2.0
34   */
35  public class SingularValueDecompositionImpl implements SingularValueDecomposition {
36  
37      /** Number of rows of the initial matrix. */
38      private int m;
39  
40      /** Number of columns of the initial matrix. */
41      private int n;
42  
43      /** Transformer to bidiagonal. */
44      private BiDiagonalTransformer transformer;
45  
46      /** Main diagonal of the bidiagonal matrix. */
47      private double[] mainBidiagonal;
48  
49      /** Secondary diagonal of the bidiagonal matrix. */
50      private double[] secondaryBidiagonal;
51  
52      /** Main diagonal of the tridiagonal matrix. */
53      private double[] mainTridiagonal;
54  
55      /** Secondary diagonal of the tridiagonal matrix. */
56      private double[] secondaryTridiagonal;
57  
58      /** Eigen decomposition of the tridiagonal matrix. */
59      private EigenDecomposition eigenDecomposition;
60  
61      /** Singular values. */
62      private double[] singularValues;
63  
64      /** Cached value of U. */
65      private RealMatrix cachedU;
66  
67      /** Cached value of U<sup>T</sup>. */
68      private RealMatrix cachedUt;
69  
70      /** Cached value of S. */
71      private RealMatrix cachedS;
72  
73      /** Cached value of V. */
74      private RealMatrix cachedV;
75  
76      /** Cached value of V<sup>T</sup>. */
77      private RealMatrix cachedVt;
78  
79      /**
80       * Calculates the Singular Value Decomposition of the given matrix. 
81       * @param matrix The matrix to decompose.
82       * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
83       * if algorithm fails to converge
84       */
85      public SingularValueDecompositionImpl(RealMatrix matrix)
86          throws InvalidMatrixException {
87  
88          m = matrix.getRowDimension();
89          n = matrix.getColumnDimension();
90  
91          cachedU  = null;
92          cachedS  = null;
93          cachedV  = null;
94          cachedVt = null;
95  
96          // transform the matrix to bidiagonal
97          transformer         = new BiDiagonalTransformer(matrix);
98          mainBidiagonal      = transformer.getMainDiagonalRef();
99          secondaryBidiagonal = transformer.getSecondaryDiagonalRef();
100 
101         // compute Bt.B (if upper diagonal) or B.Bt (if lower diagonal)
102         mainTridiagonal      = new double[mainBidiagonal.length];
103         secondaryTridiagonal = new double[mainBidiagonal.length - 1];
104         double a = mainBidiagonal[0];
105         mainTridiagonal[0] = a * a;
106         for (int i = 1; i < mainBidiagonal.length; ++i) {
107             final double b  = secondaryBidiagonal[i - 1];
108             secondaryTridiagonal[i - 1] = a * b;
109             a = mainBidiagonal[i];
110             mainTridiagonal[i] = a * a + b * b;
111         }
112 
113         // compute singular values
114         eigenDecomposition =
115             new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal,
116                                        MathUtils.SAFE_MIN);
117         singularValues = eigenDecomposition.getRealEigenvalues();
118         for (int i = 0; i < singularValues.length; ++i) {
119             singularValues[i] = Math.sqrt(singularValues[i]);
120         }
121 
122     }
123 
124     /** {@inheritDoc} */
125     public RealMatrix getU()
126         throws InvalidMatrixException {
127 
128         if (cachedU == null) {
129 
130             if (m >= n) {
131                 // the tridiagonal matrix is Bt.B, where B is upper bidiagonal
132                 final double[][] eData = eigenDecomposition.getV().getData();
133                 final double[][] iData = new double[m][];
134                 double[] ei1 = eData[0];
135                 iData[0] = ei1;
136                 for (int i = 0; i < n - 1; ++i) {
137                     // compute Bt.E.S^(-1) where E is the eigenvectors matrix
138                     // we reuse the array from matrix E to store the result 
139                     final double[] ei0 = ei1;
140                     ei1 = eData[i + 1];
141                     iData[i + 1] = ei1;
142                     for (int j = 0; j < n; ++j) {
143                         ei0[j] = (mainBidiagonal[i] * ei0[j] +
144                                   secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
145                     }
146                 }
147                 // last row
148                 final double lastMain = mainBidiagonal[n - 1];
149                 for (int j = 0; j < n; ++j) {
150                     ei1[j] *= lastMain / singularValues[j];
151                 }
152                 for (int i = n; i < m; ++i) {
153                     iData[i] = new double[n];
154                 }
155                 cachedU =
156                     transformer.getU().multiply(MatrixUtils.createRealMatrix(iData));
157             } else {
158                 // the tridiagonal matrix is B.Bt, where B is lower bidiagonal
159                 cachedU = transformer.getU().multiply(eigenDecomposition.getV());
160             }
161 
162         }
163 
164         // return the cached matrix
165         return cachedU;
166 
167     }
168 
169     /** {@inheritDoc} */
170     public RealMatrix getUT()
171         throws InvalidMatrixException {
172 
173         if (cachedUt == null) {
174             cachedUt = getU().transpose();
175         }
176 
177         // return the cached matrix
178         return cachedUt;
179 
180     }
181 
182     /** {@inheritDoc} */
183     public RealMatrix getS()
184         throws InvalidMatrixException {
185 
186         if (cachedS == null) {
187 
188             // cache the matrix for subsequent calls
189             cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
190 
191         }
192         return cachedS;
193     }
194 
195     /** {@inheritDoc} */
196     public double[] getSingularValues()
197         throws InvalidMatrixException {
198         return singularValues.clone();
199     }
200 
201     /** {@inheritDoc} */
202     public RealMatrix getV()
203         throws InvalidMatrixException {
204 
205         if (cachedV == null) {
206 
207             if (m >= n) {
208                 // the tridiagonal matrix is Bt.B, where B is upper bidiagonal
209                 cachedV = transformer.getV().multiply(eigenDecomposition.getV());
210             } else {
211                 // the tridiagonal matrix is B.Bt, where B is lower bidiagonal
212                 final double[][] eData = eigenDecomposition.getV().getData();
213                 final double[][] iData = new double[n][];
214                 double[] ei1 = eData[0];
215                 iData[0] = ei1;
216                 for (int i = 0; i < m - 1; ++i) {
217                     // compute Bt.E.S^(-1) where E is the eigenvectors matrix
218                     // we reuse the array from matrix E to store the result 
219                     final double[] ei0 = ei1;
220                     ei1 = eData[i + 1];
221                     iData[i + 1] = ei1;
222                     for (int j = 0; j < m; ++j) {
223                         ei0[j] = (mainBidiagonal[i] * ei0[j] +
224                                   secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
225                     }
226                 }
227                 // last row
228                 final double lastMain = mainBidiagonal[m - 1];
229                 for (int j = 0; j < m; ++j) {
230                     ei1[j] *= lastMain / singularValues[j];
231                 }
232                 for (int i = m; i < n; ++i) {
233                     iData[i] = new double[m];
234                 }
235                 cachedV =
236                     transformer.getV().multiply(MatrixUtils.createRealMatrix(iData));
237             }
238 
239         }
240 
241         // return the cached matrix
242         return cachedV;
243 
244     }
245 
246     /** {@inheritDoc} */
247     public RealMatrix getVT()
248         throws InvalidMatrixException {
249 
250         if (cachedVt == null) {
251             cachedVt = getV().transpose();
252         }
253 
254         // return the cached matrix
255         return cachedVt;
256 
257     }
258 
259     /** {@inheritDoc} */
260     public RealMatrix getCovariance(final double minSingularValue) {
261 
262         // get the number of singular values to consider
263         int dimension = 0;
264         while ((dimension < n) && (singularValues[dimension] >= minSingularValue)) {
265             ++dimension;
266         }
267 
268         if (dimension == 0) {
269             throw MathRuntimeException.createIllegalArgumentException(
270                   "cutoff singular value is {0}, should be at most {1}",
271                   minSingularValue, singularValues[0]);
272         }
273 
274         final double[][] data = new double[dimension][n];
275         getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
276             /** {@inheritDoc} */
277             @Override
278             public void visit(final int row, final int column, final double value) {
279                 data[row][column] = value / singularValues[row];
280             }
281         }, 0, dimension - 1, 0, n - 1);
282 
283         RealMatrix jv = new Array2DRowRealMatrix(data, false);
284         return jv.transpose().multiply(jv);
285 
286     }
287 
288     /** {@inheritDoc} */
289     public double getNorm()
290         throws InvalidMatrixException {
291         return singularValues[0];
292     }
293 
294     /** {@inheritDoc} */
295     public double getConditionNumber()
296         throws InvalidMatrixException {
297         return singularValues[0] / singularValues[singularValues.length - 1];
298     }
299 
300     /** {@inheritDoc} */
301     public int getRank()
302         throws IllegalStateException {
303 
304         final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]);
305 
306         for (int i = singularValues.length - 1; i >= 0; --i) {
307            if (singularValues[i] > threshold) {
308               return i + 1;
309            }
310         }
311         return 0;
312 
313     }
314 
315     /** {@inheritDoc} */
316     public DecompositionSolver getSolver() {
317         return new Solver(singularValues, getUT(), getV(),
318                           getRank() == singularValues.length);
319     }
320 
321     /** Specialized solver. */
322     private static class Solver implements DecompositionSolver {
323         
324         /** Singular values. */
325         private final double[] singularValues;
326 
327         /** U<sup>T</sup> matrix of the decomposition. */
328         private final RealMatrix uT;
329 
330         /** V matrix of the decomposition. */
331         private final RealMatrix v;
332 
333         /** Singularity indicator. */
334         private boolean nonSingular;
335 
336         /**
337          * Build a solver from decomposed matrix.
338          * @param singularValues singularValues
339          * @param uT U<sup>T</sup> matrix of the decomposition
340          * @param v V matrix of the decomposition
341          * @param nonSingular singularity indicator
342          */
343         private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v,
344                        final boolean nonSingular) {
345             this.singularValues = singularValues;
346             this.uT             = uT;
347             this.v              = v;
348             this.nonSingular    = nonSingular;
349         }
350 
351         /** Solve the linear equation A &times; X = B in least square sense.
352          * <p>The m&times;n matrix A may not be square, the solution X is
353          * such that ||A &times; X - B|| is minimal.</p>
354          * @param b right-hand side of the equation A &times; X = B
355          * @return a vector X that minimizes the two norm of A &times; X - B
356          * @exception IllegalArgumentException if matrices dimensions don't match
357          * @exception InvalidMatrixException if decomposed matrix is singular
358          */
359         public double[] solve(final double[] b)
360             throws IllegalArgumentException, InvalidMatrixException {
361 
362             if (b.length != singularValues.length) {
363                 throw MathRuntimeException.createIllegalArgumentException(
364                         "vector length mismatch: got {0} but expected {1}",
365                         b.length, singularValues.length);
366             }
367 
368             final double[] w = uT.operate(b);
369             for (int i = 0; i < singularValues.length; ++i) {
370                 final double si = singularValues[i];
371                 if (si == 0) {
372                     throw new SingularMatrixException();
373                 }
374                 w[i] /= si;
375             }
376             return v.operate(w);
377 
378         }
379 
380         /** Solve the linear equation A &times; X = B in least square sense.
381          * <p>The m&times;n matrix A may not be square, the solution X is
382          * such that ||A &times; X - B|| is minimal.</p>
383          * @param b right-hand side of the equation A &times; X = B
384          * @return a vector X that minimizes the two norm of A &times; X - B
385          * @exception IllegalArgumentException if matrices dimensions don't match
386          * @exception InvalidMatrixException if decomposed matrix is singular
387          */
388         public RealVector solve(final RealVector b)
389             throws IllegalArgumentException, InvalidMatrixException {
390 
391             if (b.getDimension() != singularValues.length) {
392                 throw MathRuntimeException.createIllegalArgumentException(
393                         "vector length mismatch: got {0} but expected {1}",
394                          b.getDimension(), singularValues.length);
395             }
396 
397             final RealVector w = uT.operate(b);
398             for (int i = 0; i < singularValues.length; ++i) {
399                 final double si = singularValues[i];
400                 if (si == 0) {
401                     throw new SingularMatrixException();
402                 }
403                 w.setEntry(i, w.getEntry(i) / si);
404             }
405             return v.operate(w);
406 
407         }
408 
409         /** Solve the linear equation A &times; X = B in least square sense.
410          * <p>The m&times;n matrix A may not be square, the solution X is
411          * such that ||A &times; X - B|| is minimal.</p>
412          * @param b right-hand side of the equation A &times; X = B
413          * @return a matrix X that minimizes the two norm of A &times; X - B
414          * @exception IllegalArgumentException if matrices dimensions don't match
415          * @exception InvalidMatrixException if decomposed matrix is singular
416          */
417         public RealMatrix solve(final RealMatrix b)
418             throws IllegalArgumentException, InvalidMatrixException {
419 
420             if (b.getRowDimension() != singularValues.length) {
421                 throw MathRuntimeException.createIllegalArgumentException(
422                         "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
423                         b.getRowDimension(), b.getColumnDimension(),
424                         singularValues.length, "n");
425             }
426 
427             final RealMatrix w = uT.multiply(b);
428             for (int i = 0; i < singularValues.length; ++i) {
429                 final double si  = singularValues[i];
430                 if (si == 0) {
431                     throw new SingularMatrixException();
432                 }
433                 final double inv = 1.0 / si;
434                 for (int j = 0; j < b.getColumnDimension(); ++j) {
435                     w.multiplyEntry(i, j, inv);
436                 }
437             }
438             return v.multiply(w);
439 
440         }
441 
442         /**
443          * Check if the decomposed matrix is non-singular.
444          * @return true if the decomposed matrix is non-singular
445          */
446         public boolean isNonSingular() {
447             return nonSingular;
448         }
449 
450         /** Get the pseudo-inverse of the decomposed matrix.
451          * @return inverse matrix
452          * @throws InvalidMatrixException if decomposed matrix is singular
453          */
454         public RealMatrix getInverse()
455             throws InvalidMatrixException {
456 
457             if (!isNonSingular()) {
458                 throw new SingularMatrixException();
459             }
460 
461             return solve(MatrixUtils.createRealIdentityMatrix(singularValues.length));
462 
463         }
464 
465     }
466 
467 }