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.ConvergenceException;
021 import org.apache.commons.math.MathRuntimeException;
022 import org.apache.commons.math.util.MathUtils;
023
024 /**
025 * Calculates the Singular Value Decomposition of a matrix.
026 * <p>The Singular Value Decomposition of matrix A is a set of three matrices:
027 * U, Σ and V such that A = U × Σ × V<sup>T</sup>.
028 * Let A be an m × n matrix, then U is an m × m orthogonal matrix,
029 * Σ is a m × n diagonal matrix with positive diagonal elements,
030 * and V is an n × n orthogonal matrix.</p>
031 *
032 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
033 * @since 2.0
034 */
035 public class SingularValueDecompositionImpl implements SingularValueDecomposition {
036
037 /** Number of rows of the initial matrix. */
038 private int m;
039
040 /** Number of columns of the initial matrix. */
041 private int n;
042
043 /** Transformer to bidiagonal. */
044 private BiDiagonalTransformer transformer;
045
046 /** Main diagonal of the bidiagonal matrix. */
047 private double[] mainBidiagonal;
048
049 /** Secondary diagonal of the bidiagonal matrix. */
050 private double[] secondaryBidiagonal;
051
052 /** Main diagonal of the tridiagonal matrix. */
053 private double[] mainTridiagonal;
054
055 /** Secondary diagonal of the tridiagonal matrix. */
056 private double[] secondaryTridiagonal;
057
058 /** Eigen decomposition of the tridiagonal matrix. */
059 private EigenDecomposition eigenDecomposition;
060
061 /** Singular values. */
062 private double[] singularValues;
063
064 /** Cached value of U. */
065 private RealMatrix cachedU;
066
067 /** Cached value of U<sup>T</sup>. */
068 private RealMatrix cachedUt;
069
070 /** Cached value of S. */
071 private RealMatrix cachedS;
072
073 /** Cached value of V. */
074 private RealMatrix cachedV;
075
076 /** Cached value of V<sup>T</sup>. */
077 private RealMatrix cachedVt;
078
079 /**
080 * Calculates the Singular Value Decomposition of the given matrix.
081 * @param matrix The matrix to decompose.
082 * @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
083 * if algorithm fails to converge
084 */
085 public SingularValueDecompositionImpl(RealMatrix matrix)
086 throws InvalidMatrixException {
087
088 m = matrix.getRowDimension();
089 n = matrix.getColumnDimension();
090
091 cachedU = null;
092 cachedS = null;
093 cachedV = null;
094 cachedVt = null;
095
096 // transform the matrix to bidiagonal
097 transformer = new BiDiagonalTransformer(matrix);
098 mainBidiagonal = transformer.getMainDiagonalRef();
099 secondaryBidiagonal = transformer.getSecondaryDiagonalRef();
100
101 // compute Bt.B (if upper diagonal) or B.Bt (if lower diagonal)
102 mainTridiagonal = new double[mainBidiagonal.length];
103 secondaryTridiagonal = new double[mainBidiagonal.length - 1];
104 double a = mainBidiagonal[0];
105 mainTridiagonal[0] = a * a;
106 for (int i = 1; i < mainBidiagonal.length; ++i) {
107 final double b = secondaryBidiagonal[i - 1];
108 secondaryTridiagonal[i - 1] = a * b;
109 a = mainBidiagonal[i];
110 mainTridiagonal[i] = a * a + b * b;
111 }
112
113 // compute singular values
114 eigenDecomposition =
115 new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal,
116 MathUtils.SAFE_MIN);
117 singularValues = eigenDecomposition.getRealEigenvalues();
118 for (int i = 0; i < singularValues.length; ++i) {
119 singularValues[i] = Math.sqrt(singularValues[i]);
120 }
121
122 }
123
124 /** {@inheritDoc} */
125 public RealMatrix getU()
126 throws InvalidMatrixException {
127
128 if (cachedU == null) {
129
130 if (m >= n) {
131 // the tridiagonal matrix is Bt.B, where B is upper bidiagonal
132 final double[][] eData = eigenDecomposition.getV().getData();
133 final double[][] iData = new double[m][];
134 double[] ei1 = eData[0];
135 iData[0] = ei1;
136 for (int i = 0; i < n - 1; ++i) {
137 // compute Bt.E.S^(-1) where E is the eigenvectors matrix
138 // we reuse the array from matrix E to store the result
139 final double[] ei0 = ei1;
140 ei1 = eData[i + 1];
141 iData[i + 1] = ei1;
142 for (int j = 0; j < n; ++j) {
143 ei0[j] = (mainBidiagonal[i] * ei0[j] +
144 secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
145 }
146 }
147 // last row
148 final double lastMain = mainBidiagonal[n - 1];
149 for (int j = 0; j < n; ++j) {
150 ei1[j] *= lastMain / singularValues[j];
151 }
152 for (int i = n; i < m; ++i) {
153 iData[i] = new double[n];
154 }
155 cachedU =
156 transformer.getU().multiply(MatrixUtils.createRealMatrix(iData));
157 } else {
158 // the tridiagonal matrix is B.Bt, where B is lower bidiagonal
159 cachedU = transformer.getU().multiply(eigenDecomposition.getV());
160 }
161
162 }
163
164 // return the cached matrix
165 return cachedU;
166
167 }
168
169 /** {@inheritDoc} */
170 public RealMatrix getUT()
171 throws InvalidMatrixException {
172
173 if (cachedUt == null) {
174 cachedUt = getU().transpose();
175 }
176
177 // return the cached matrix
178 return cachedUt;
179
180 }
181
182 /** {@inheritDoc} */
183 public RealMatrix getS()
184 throws InvalidMatrixException {
185
186 if (cachedS == null) {
187
188 // cache the matrix for subsequent calls
189 cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
190
191 }
192 return cachedS;
193 }
194
195 /** {@inheritDoc} */
196 public double[] getSingularValues()
197 throws InvalidMatrixException {
198 return singularValues.clone();
199 }
200
201 /** {@inheritDoc} */
202 public RealMatrix getV()
203 throws InvalidMatrixException {
204
205 if (cachedV == null) {
206
207 if (m >= n) {
208 // the tridiagonal matrix is Bt.B, where B is upper bidiagonal
209 cachedV = transformer.getV().multiply(eigenDecomposition.getV());
210 } else {
211 // the tridiagonal matrix is B.Bt, where B is lower bidiagonal
212 final double[][] eData = eigenDecomposition.getV().getData();
213 final double[][] iData = new double[n][];
214 double[] ei1 = eData[0];
215 iData[0] = ei1;
216 for (int i = 0; i < m - 1; ++i) {
217 // compute Bt.E.S^(-1) where E is the eigenvectors matrix
218 // we reuse the array from matrix E to store the result
219 final double[] ei0 = ei1;
220 ei1 = eData[i + 1];
221 iData[i + 1] = ei1;
222 for (int j = 0; j < m; ++j) {
223 ei0[j] = (mainBidiagonal[i] * ei0[j] +
224 secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
225 }
226 }
227 // last row
228 final double lastMain = mainBidiagonal[m - 1];
229 for (int j = 0; j < m; ++j) {
230 ei1[j] *= lastMain / singularValues[j];
231 }
232 for (int i = m; i < n; ++i) {
233 iData[i] = new double[m];
234 }
235 cachedV =
236 transformer.getV().multiply(MatrixUtils.createRealMatrix(iData));
237 }
238
239 }
240
241 // return the cached matrix
242 return cachedV;
243
244 }
245
246 /** {@inheritDoc} */
247 public RealMatrix getVT()
248 throws InvalidMatrixException {
249
250 if (cachedVt == null) {
251 cachedVt = getV().transpose();
252 }
253
254 // return the cached matrix
255 return cachedVt;
256
257 }
258
259 /** {@inheritDoc} */
260 public RealMatrix getCovariance(final double minSingularValue) {
261
262 // get the number of singular values to consider
263 int dimension = 0;
264 while ((dimension < n) && (singularValues[dimension] >= minSingularValue)) {
265 ++dimension;
266 }
267
268 if (dimension == 0) {
269 throw MathRuntimeException.createIllegalArgumentException(
270 "cutoff singular value is {0}, should be at most {1}",
271 minSingularValue, singularValues[0]);
272 }
273
274 final double[][] data = new double[dimension][n];
275 getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
276 /** {@inheritDoc} */
277 @Override
278 public void visit(final int row, final int column, final double value) {
279 data[row][column] = value / singularValues[row];
280 }
281 }, 0, dimension - 1, 0, n - 1);
282
283 RealMatrix jv = new Array2DRowRealMatrix(data, false);
284 return jv.transpose().multiply(jv);
285
286 }
287
288 /** {@inheritDoc} */
289 public double getNorm()
290 throws InvalidMatrixException {
291 return singularValues[0];
292 }
293
294 /** {@inheritDoc} */
295 public double getConditionNumber()
296 throws InvalidMatrixException {
297 return singularValues[0] / singularValues[singularValues.length - 1];
298 }
299
300 /** {@inheritDoc} */
301 public int getRank()
302 throws IllegalStateException {
303
304 final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]);
305
306 for (int i = singularValues.length - 1; i >= 0; --i) {
307 if (singularValues[i] > threshold) {
308 return i + 1;
309 }
310 }
311 return 0;
312
313 }
314
315 /** {@inheritDoc} */
316 public DecompositionSolver getSolver() {
317 return new Solver(singularValues, getUT(), getV(),
318 getRank() == singularValues.length);
319 }
320
321 /** Specialized solver. */
322 private static class Solver implements DecompositionSolver {
323
324 /** Singular values. */
325 private final double[] singularValues;
326
327 /** U<sup>T</sup> matrix of the decomposition. */
328 private final RealMatrix uT;
329
330 /** V matrix of the decomposition. */
331 private final RealMatrix v;
332
333 /** Singularity indicator. */
334 private boolean nonSingular;
335
336 /**
337 * Build a solver from decomposed matrix.
338 * @param singularValues singularValues
339 * @param uT U<sup>T</sup> matrix of the decomposition
340 * @param v V matrix of the decomposition
341 * @param nonSingular singularity indicator
342 */
343 private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v,
344 final boolean nonSingular) {
345 this.singularValues = singularValues;
346 this.uT = uT;
347 this.v = v;
348 this.nonSingular = nonSingular;
349 }
350
351 /** Solve the linear equation A × X = B in least square sense.
352 * <p>The m×n matrix A may not be square, the solution X is
353 * such that ||A × X - B|| is minimal.</p>
354 * @param b right-hand side of the equation A × X = B
355 * @return a vector X that minimizes the two norm of A × X - B
356 * @exception IllegalArgumentException if matrices dimensions don't match
357 * @exception InvalidMatrixException if decomposed matrix is singular
358 */
359 public double[] solve(final double[] b)
360 throws IllegalArgumentException, InvalidMatrixException {
361
362 if (b.length != singularValues.length) {
363 throw MathRuntimeException.createIllegalArgumentException(
364 "vector length mismatch: got {0} but expected {1}",
365 b.length, singularValues.length);
366 }
367
368 final double[] w = uT.operate(b);
369 for (int i = 0; i < singularValues.length; ++i) {
370 final double si = singularValues[i];
371 if (si == 0) {
372 throw new SingularMatrixException();
373 }
374 w[i] /= si;
375 }
376 return v.operate(w);
377
378 }
379
380 /** Solve the linear equation A × X = B in least square sense.
381 * <p>The m×n matrix A may not be square, the solution X is
382 * such that ||A × X - B|| is minimal.</p>
383 * @param b right-hand side of the equation A × X = B
384 * @return a vector X that minimizes the two norm of A × X - B
385 * @exception IllegalArgumentException if matrices dimensions don't match
386 * @exception InvalidMatrixException if decomposed matrix is singular
387 */
388 public RealVector solve(final RealVector b)
389 throws IllegalArgumentException, InvalidMatrixException {
390
391 if (b.getDimension() != singularValues.length) {
392 throw MathRuntimeException.createIllegalArgumentException(
393 "vector length mismatch: got {0} but expected {1}",
394 b.getDimension(), singularValues.length);
395 }
396
397 final RealVector w = uT.operate(b);
398 for (int i = 0; i < singularValues.length; ++i) {
399 final double si = singularValues[i];
400 if (si == 0) {
401 throw new SingularMatrixException();
402 }
403 w.setEntry(i, w.getEntry(i) / si);
404 }
405 return v.operate(w);
406
407 }
408
409 /** Solve the linear equation A × X = B in least square sense.
410 * <p>The m×n matrix A may not be square, the solution X is
411 * such that ||A × X - B|| is minimal.</p>
412 * @param b right-hand side of the equation A × X = B
413 * @return a matrix X that minimizes the two norm of A × X - B
414 * @exception IllegalArgumentException if matrices dimensions don't match
415 * @exception InvalidMatrixException if decomposed matrix is singular
416 */
417 public RealMatrix solve(final RealMatrix b)
418 throws IllegalArgumentException, InvalidMatrixException {
419
420 if (b.getRowDimension() != singularValues.length) {
421 throw MathRuntimeException.createIllegalArgumentException(
422 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
423 b.getRowDimension(), b.getColumnDimension(),
424 singularValues.length, "n");
425 }
426
427 final RealMatrix w = uT.multiply(b);
428 for (int i = 0; i < singularValues.length; ++i) {
429 final double si = singularValues[i];
430 if (si == 0) {
431 throw new SingularMatrixException();
432 }
433 final double inv = 1.0 / si;
434 for (int j = 0; j < b.getColumnDimension(); ++j) {
435 w.multiplyEntry(i, j, inv);
436 }
437 }
438 return v.multiply(w);
439
440 }
441
442 /**
443 * Check if the decomposed matrix is non-singular.
444 * @return true if the decomposed matrix is non-singular
445 */
446 public boolean isNonSingular() {
447 return nonSingular;
448 }
449
450 /** Get the pseudo-inverse of the decomposed matrix.
451 * @return inverse matrix
452 * @throws InvalidMatrixException if decomposed matrix is singular
453 */
454 public RealMatrix getInverse()
455 throws InvalidMatrixException {
456
457 if (!isNonSingular()) {
458 throw new SingularMatrixException();
459 }
460
461 return solve(MatrixUtils.createRealIdentityMatrix(singularValues.length));
462
463 }
464
465 }
466
467 }