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.MathRuntimeException;
021    
022    /**
023     * Calculates the compact Singular Value Decomposition of a matrix.
024     * <p>
025     * The Singular Value Decomposition of matrix A is a set of three matrices: U,
026     * &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>. Let A be
027     * a m &times; n matrix, then U is a m &times; p orthogonal matrix, &Sigma; is a
028     * p &times; p diagonal matrix with positive or null elements, V is a p &times;
029     * n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
030     * p=min(m,n).
031     * </p>
032     * @version $Revision: 912413 $ $Date: 2010-02-21 16:46:12 -0500 (Sun, 21 Feb 2010) $
033     * @since 2.0
034     */
035    public class SingularValueDecompositionImpl implements
036            SingularValueDecomposition {
037    
038        /** Number of rows of the initial matrix. */
039        private int m;
040    
041        /** Number of columns of the initial matrix. */
042        private int n;
043    
044        /** Eigen decomposition of the tridiagonal matrix. */
045        private EigenDecomposition eigenDecomposition;
046    
047        /** Singular values. */
048        private double[] singularValues;
049    
050        /** Cached value of U. */
051        private RealMatrix cachedU;
052    
053        /** Cached value of U<sup>T</sup>. */
054        private RealMatrix cachedUt;
055    
056        /** Cached value of S. */
057        private RealMatrix cachedS;
058    
059        /** Cached value of V. */
060        private RealMatrix cachedV;
061    
062        /** Cached value of V<sup>T</sup>. */
063        private RealMatrix cachedVt;
064    
065        /**
066         * Calculates the compact Singular Value Decomposition of the given matrix.
067         * @param matrix
068         *            The matrix to decompose.
069         * @exception InvalidMatrixException
070         *                (wrapping a
071         *                {@link org.apache.commons.math.ConvergenceException} if
072         *                algorithm fails to converge
073         */
074        public SingularValueDecompositionImpl(final RealMatrix matrix)
075                throws InvalidMatrixException {
076    
077            m = matrix.getRowDimension();
078            n = matrix.getColumnDimension();
079    
080            cachedU = null;
081            cachedS = null;
082            cachedV = null;
083            cachedVt = null;
084    
085            double[][] localcopy = matrix.getData();
086            double[][] matATA = new double[n][n];
087            //
088            // create A^T*A
089            //
090            for (int i = 0; i < n; i++) {
091                for (int j = i; j < n; j++) {
092                    matATA[i][j] = 0.0;
093                    for (int k = 0; k < m; k++) {
094                        matATA[i][j] += localcopy[k][i] * localcopy[k][j];
095                    }
096                    matATA[j][i]=matATA[i][j];
097                }
098            }
099    
100            double[][] matAAT = new double[m][m];
101            //
102            // create A*A^T
103            //
104            for (int i = 0; i < m; i++) {
105                for (int j = i; j < m; j++) {
106                    matAAT[i][j] = 0.0;
107                    for (int k = 0; k < n; k++) {
108                        matAAT[i][j] += localcopy[i][k] * localcopy[j][k];
109                    }
110                    matAAT[j][i]=matAAT[i][j];
111                }
112            }
113            int p;
114            if (m>=n) {
115                p=n;
116                // compute eigen decomposition of A^T*A
117                eigenDecomposition = new EigenDecompositionImpl(
118                        new Array2DRowRealMatrix(matATA),1.0);
119                singularValues = eigenDecomposition.getRealEigenvalues();
120                cachedV = eigenDecomposition.getV();
121    
122                // compute eigen decomposition of A*A^T
123                eigenDecomposition = new EigenDecompositionImpl(
124                        new Array2DRowRealMatrix(matAAT),1.0);
125                cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);
126            } else {
127                p=m;
128                // compute eigen decomposition of A*A^T
129                eigenDecomposition = new EigenDecompositionImpl(
130                        new Array2DRowRealMatrix(matAAT),1.0);
131                singularValues = eigenDecomposition.getRealEigenvalues();
132                cachedU = eigenDecomposition.getV();
133    
134                // compute eigen decomposition of A^T*A
135                eigenDecomposition = new EigenDecompositionImpl(
136                        new Array2DRowRealMatrix(matATA),1.0);
137                cachedV = eigenDecomposition.getV().getSubMatrix(0,n-1,0,p-1);
138            }
139            for (int i = 0; i < p; i++) {
140                singularValues[i] = Math.sqrt(Math.abs(singularValues[i]));
141            }
142            // Up to this point, U and V are computed independently of each other.
143            // There still an sign indetermination of each column of, say, U.
144            // The sign is set such that A.V_i=sigma_i.U_i (i<=p)
145            // The right sign corresponds to a positive dot product of A.V_i and U_i
146            for (int i = 0; i < p; i++) {
147              RealVector tmp = cachedU.getColumnVector(i);
148              double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp);
149              if (product<0) {
150                cachedU.setColumnVector(i, tmp.mapMultiply(-1.0));
151              }
152            }
153        }
154    
155        /** {@inheritDoc} */
156        public RealMatrix getU() throws InvalidMatrixException {
157            // return the cached matrix
158            return cachedU;
159    
160        }
161    
162        /** {@inheritDoc} */
163        public RealMatrix getUT() throws InvalidMatrixException {
164    
165            if (cachedUt == null) {
166                cachedUt = getU().transpose();
167            }
168    
169            // return the cached matrix
170            return cachedUt;
171    
172        }
173    
174        /** {@inheritDoc} */
175        public RealMatrix getS() throws InvalidMatrixException {
176    
177            if (cachedS == null) {
178    
179                // cache the matrix for subsequent calls
180                cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
181    
182            }
183            return cachedS;
184        }
185    
186        /** {@inheritDoc} */
187        public double[] getSingularValues() throws InvalidMatrixException {
188            return singularValues.clone();
189        }
190    
191        /** {@inheritDoc} */
192        public RealMatrix getV() throws InvalidMatrixException {
193            // return the cached matrix
194            return cachedV;
195    
196        }
197    
198        /** {@inheritDoc} */
199        public RealMatrix getVT() throws InvalidMatrixException {
200    
201            if (cachedVt == null) {
202                cachedVt = getV().transpose();
203            }
204    
205            // return the cached matrix
206            return cachedVt;
207    
208        }
209    
210        /** {@inheritDoc} */
211        public RealMatrix getCovariance(final double minSingularValue) {
212    
213            // get the number of singular values to consider
214            final int p = singularValues.length;
215            int dimension = 0;
216            while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) {
217                ++dimension;
218            }
219    
220            if (dimension == 0) {
221                throw MathRuntimeException.createIllegalArgumentException(
222                        "cutoff singular value is {0}, should be at most {1}",
223                        minSingularValue, singularValues[0]);
224            }
225    
226            final double[][] data = new double[dimension][p];
227            getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
228                /** {@inheritDoc} */
229                @Override
230                public void visit(final int row, final int column,
231                        final double value) {
232                    data[row][column] = value / singularValues[row];
233                }
234            }, 0, dimension - 1, 0, p - 1);
235    
236            RealMatrix jv = new Array2DRowRealMatrix(data, false);
237            return jv.transpose().multiply(jv);
238    
239        }
240    
241        /** {@inheritDoc} */
242        public double getNorm() throws InvalidMatrixException {
243            return singularValues[0];
244        }
245    
246        /** {@inheritDoc} */
247        public double getConditionNumber() throws InvalidMatrixException {
248            return singularValues[0] / singularValues[singularValues.length - 1];
249        }
250    
251        /** {@inheritDoc} */
252        public int getRank() throws IllegalStateException {
253    
254            final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]);
255    
256            for (int i = singularValues.length - 1; i >= 0; --i) {
257                if (singularValues[i] > threshold) {
258                    return i + 1;
259                }
260            }
261            return 0;
262    
263        }
264    
265        /** {@inheritDoc} */
266        public DecompositionSolver getSolver() {
267            return new Solver(singularValues, getUT(), getV(), getRank() == Math
268                    .max(m, n));
269        }
270    
271        /** Specialized solver. */
272        private static class Solver implements DecompositionSolver {
273    
274            /** Pseudo-inverse of the initial matrix. */
275            private final RealMatrix pseudoInverse;
276    
277            /** Singularity indicator. */
278            private boolean nonSingular;
279    
280            /**
281             * Build a solver from decomposed matrix.
282             * @param singularValues
283             *            singularValues
284             * @param uT
285             *            U<sup>T</sup> matrix of the decomposition
286             * @param v
287             *            V matrix of the decomposition
288             * @param nonSingular
289             *            singularity indicator
290             */
291            private Solver(final double[] singularValues, final RealMatrix uT,
292                    final RealMatrix v, final boolean nonSingular) {
293                double[][] suT = uT.getData();
294                for (int i = 0; i < singularValues.length; ++i) {
295                    final double a;
296                    if (singularValues[i]>0) {
297                     a=1.0 / singularValues[i];
298                    } else {
299                     a=0.0;
300                    }
301                    final double[] suTi = suT[i];
302                    for (int j = 0; j < suTi.length; ++j) {
303                        suTi[j] *= a;
304                    }
305                }
306                pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));
307                this.nonSingular = nonSingular;
308            }
309    
310            /**
311             * Solve the linear equation A &times; X = B in least square sense.
312             * <p>
313             * The m&times;n matrix A may not be square, the solution X is such that
314             * ||A &times; X - B|| is minimal.
315             * </p>
316             * @param b
317             *            right-hand side of the equation A &times; X = B
318             * @return a vector X that minimizes the two norm of A &times; X - B
319             * @exception IllegalArgumentException
320             *                if matrices dimensions don't match
321             */
322            public double[] solve(final double[] b) throws IllegalArgumentException {
323                return pseudoInverse.operate(b);
324            }
325    
326            /**
327             * Solve the linear equation A &times; X = B in least square sense.
328             * <p>
329             * The m&times;n matrix A may not be square, the solution X is such that
330             * ||A &times; X - B|| is minimal.
331             * </p>
332             * @param b
333             *            right-hand side of the equation A &times; X = B
334             * @return a vector X that minimizes the two norm of A &times; X - B
335             * @exception IllegalArgumentException
336             *                if matrices dimensions don't match
337             */
338            public RealVector solve(final RealVector b)
339                    throws IllegalArgumentException {
340                return pseudoInverse.operate(b);
341            }
342    
343            /**
344             * Solve the linear equation A &times; X = B in least square sense.
345             * <p>
346             * The m&times;n matrix A may not be square, the solution X is such that
347             * ||A &times; X - B|| is minimal.
348             * </p>
349             * @param b
350             *            right-hand side of the equation A &times; X = B
351             * @return a matrix X that minimizes the two norm of A &times; X - B
352             * @exception IllegalArgumentException
353             *                if matrices dimensions don't match
354             */
355            public RealMatrix solve(final RealMatrix b)
356                    throws IllegalArgumentException {
357                return pseudoInverse.multiply(b);
358            }
359    
360            /**
361             * Check if the decomposed matrix is non-singular.
362             * @return true if the decomposed matrix is non-singular
363             */
364            public boolean isNonSingular() {
365                return nonSingular;
366            }
367    
368            /**
369             * Get the pseudo-inverse of the decomposed matrix.
370             * @return inverse matrix
371             */
372            public RealMatrix getInverse() {
373                return pseudoInverse;
374            }
375    
376        }
377    
378    }