001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.linear;
019    
020    import org.apache.commons.math.ConvergenceException;
021    import org.apache.commons.math.MathRuntimeException;
022    import org.apache.commons.math.util.MathUtils;
023    
024    /**
025     * Calculates the Singular Value Decomposition of a matrix.
026     * <p>The Singular Value Decomposition of matrix A is a set of three matrices:
027     * U, &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>.
028     * Let A be an m &times; n matrix, then U is an m &times; m orthogonal matrix,
029     * &Sigma; is a m &times; n diagonal matrix with positive diagonal elements,
030     * and V is an n &times; n orthogonal matrix.</p>
031     *
032     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
033     * @since 2.0
034     */
035    public class SingularValueDecompositionImpl implements SingularValueDecomposition {
036    
037        /** Number of rows of the initial matrix. */
038        private int m;
039    
040        /** Number of columns of the initial matrix. */
041        private int n;
042    
043        /** Transformer to bidiagonal. */
044        private BiDiagonalTransformer transformer;
045    
046        /** Main diagonal of the bidiagonal matrix. */
047        private double[] mainBidiagonal;
048    
049        /** Secondary diagonal of the bidiagonal matrix. */
050        private double[] secondaryBidiagonal;
051    
052        /** Main diagonal of the tridiagonal matrix. */
053        private double[] mainTridiagonal;
054    
055        /** Secondary diagonal of the tridiagonal matrix. */
056        private double[] secondaryTridiagonal;
057    
058        /** Eigen decomposition of the tridiagonal matrix. */
059        private EigenDecomposition eigenDecomposition;
060    
061        /** Singular values. */
062        private double[] singularValues;
063    
064        /** Cached value of U. */
065        private RealMatrix cachedU;
066    
067        /** Cached value of U<sup>T</sup>. */
068        private RealMatrix cachedUt;
069    
070        /** Cached value of S. */
071        private RealMatrix cachedS;
072    
073        /** Cached value of V. */
074        private RealMatrix cachedV;
075    
076        /** Cached value of V<sup>T</sup>. */
077        private RealMatrix cachedVt;
078    
079        /**
080         * Calculates the Singular Value Decomposition of the given matrix. 
081         * @param matrix The matrix to decompose.
082         * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
083         * if algorithm fails to converge
084         */
085        public SingularValueDecompositionImpl(RealMatrix matrix)
086            throws InvalidMatrixException {
087    
088            m = matrix.getRowDimension();
089            n = matrix.getColumnDimension();
090    
091            cachedU  = null;
092            cachedS  = null;
093            cachedV  = null;
094            cachedVt = null;
095    
096            // transform the matrix to bidiagonal
097            transformer         = new BiDiagonalTransformer(matrix);
098            mainBidiagonal      = transformer.getMainDiagonalRef();
099            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    }