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 * Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be
027 * a m × n matrix, then U is a m × p orthogonal matrix, Σ is a
028 * p × p diagonal matrix with positive or null elements, V is a p ×
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 × X = B in least square sense.
312 * <p>
313 * The m×n matrix A may not be square, the solution X is such that
314 * ||A × X - B|| is minimal.
315 * </p>
316 * @param b
317 * right-hand side of the equation A × X = B
318 * @return a vector X that minimizes the two norm of A × 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 × X = B in least square sense.
328 * <p>
329 * The m×n matrix A may not be square, the solution X is such that
330 * ||A × X - B|| is minimal.
331 * </p>
332 * @param b
333 * right-hand side of the equation A × X = B
334 * @return a vector X that minimizes the two norm of A × 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 × X = B in least square sense.
345 * <p>
346 * The m×n matrix A may not be square, the solution X is such that
347 * ||A × X - B|| is minimal.
348 * </p>
349 * @param b
350 * right-hand side of the equation A × X = B
351 * @return a matrix X that minimizes the two norm of A × 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 }