    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    }