Home | History | Annotate | Download | only in linear 
      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.MathRuntimeException;
     21 import org.apache.commons.math.exception.util.LocalizedFormats;
     22 import org.apache.commons.math.util.FastMath;
     23 
     24 /**
     25  * Calculates the compact Singular Value Decomposition of a matrix.
     26  * <p>
     27  * The Singular Value Decomposition of matrix A is a set of three matrices: U,
     28  * &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>. Let A be
     29  * a m &times; n matrix, then U is a m &times; p orthogonal matrix, &Sigma; is a
     30  * p &times; p diagonal matrix with positive or null elements, V is a p &times;
     31  * n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
     32  * p=min(m,n).
     33  * </p>
     34  * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $
     35  * @since 2.0
     36  */
     37 public class SingularValueDecompositionImpl implements
     38         SingularValueDecomposition {
     39 
     40     /** Number of rows of the initial matrix. */
     41     private int m;
     42 
     43     /** Number of columns of the initial matrix. */
     44     private int n;
     45 
     46     /** Eigen decomposition of the tridiagonal matrix. */
     47     private EigenDecomposition eigenDecomposition;
     48 
     49     /** Singular values. */
     50     private double[] singularValues;
     51 
     52     /** Cached value of U. */
     53     private RealMatrix cachedU;
     54 
     55     /** Cached value of U<sup>T</sup>. */
     56     private RealMatrix cachedUt;
     57 
     58     /** Cached value of S. */
     59     private RealMatrix cachedS;
     60 
     61     /** Cached value of V. */
     62     private RealMatrix cachedV;
     63 
     64     /** Cached value of V<sup>T</sup>. */
     65     private RealMatrix cachedVt;
     66 
     67     /**
     68      * Calculates the compact Singular Value Decomposition of the given matrix.
     69      * @param matrix
     70      *            The matrix to decompose.
     71      * @exception InvalidMatrixException
     72      *                (wrapping a
     73      *                {@link org.apache.commons.math.ConvergenceException} if
     74      *                algorithm fails to converge
     75      */
     76     public SingularValueDecompositionImpl(final RealMatrix matrix)
     77             throws InvalidMatrixException {
     78 
     79         m = matrix.getRowDimension();
     80         n = matrix.getColumnDimension();
     81 
     82         cachedU = null;
     83         cachedS = null;
     84         cachedV = null;
     85         cachedVt = null;
     86 
     87         double[][] localcopy = matrix.getData();
     88         double[][] matATA = new double[n][n];
     89         //
     90         // create A^T*A
     91         //
     92         for (int i = 0; i < n; i++) {
     93             for (int j = i; j < n; j++) {
     94                 matATA[i][j] = 0.0;
     95                 for (int k = 0; k < m; k++) {
     96                     matATA[i][j] += localcopy[k][i] * localcopy[k][j];
     97                 }
     98                 matATA[j][i]=matATA[i][j];
     99             }
    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 }
    380