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