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