18
18
package org.apache.commons.math.linear;
20
import org.apache.commons.math.ConvergenceException;
21
20
import org.apache.commons.math.MathRuntimeException;
22
import org.apache.commons.math.util.MathUtils;
25
* Calculates the Singular Value Decomposition of a matrix.
26
* <p>The Singular Value Decomposition of matrix A is a set of three matrices:
27
* U, Σ and V such that A = U × Σ × V<sup>T</sup>.
28
* Let A be an m × n matrix, then U is an m × m orthogonal matrix,
29
* Σ is a m × n diagonal matrix with positive diagonal elements,
30
* and V is an n × n orthogonal matrix.</p>
32
* @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
23
* Calculates the compact Singular Value Decomposition of a matrix.
25
* The Singular Value Decomposition of matrix A is a set of three matrices: U,
26
* Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be
27
* a m × n matrix, then U is a m × p orthogonal matrix, Σ is a
28
* p × p diagonal matrix with positive or null elements, V is a p ×
29
* n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
32
* @version $Revision: 912413 $ $Date: 2010-02-21 16:46:12 -0500 (Sun, 21 Feb 2010) $
35
public class SingularValueDecompositionImpl implements SingularValueDecomposition {
35
public class SingularValueDecompositionImpl implements
36
SingularValueDecomposition {
37
38
/** Number of rows of the initial matrix. */
77
63
private RealMatrix cachedVt;
80
* Calculates the Singular Value Decomposition of the given matrix.
81
* @param matrix The matrix to decompose.
82
* @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
83
* if algorithm fails to converge
66
* Calculates the compact Singular Value Decomposition of the given matrix.
68
* The matrix to decompose.
69
* @exception InvalidMatrixException
71
* {@link org.apache.commons.math.ConvergenceException} if
72
* algorithm fails to converge
85
public SingularValueDecompositionImpl(RealMatrix matrix)
86
throws InvalidMatrixException {
74
public SingularValueDecompositionImpl(final RealMatrix matrix)
75
throws InvalidMatrixException {
88
77
m = matrix.getRowDimension();
89
78
n = matrix.getColumnDimension();
96
// transform the matrix to bidiagonal
97
transformer = new BiDiagonalTransformer(matrix);
98
mainBidiagonal = transformer.getMainDiagonalRef();
99
secondaryBidiagonal = transformer.getSecondaryDiagonalRef();
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;
113
// compute singular values
115
new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal,
117
singularValues = eigenDecomposition.getRealEigenvalues();
118
for (int i = 0; i < singularValues.length; ++i) {
119
singularValues[i] = Math.sqrt(singularValues[i]);
85
double[][] localcopy = matrix.getData();
86
double[][] matATA = new double[n][n];
90
for (int i = 0; i < n; i++) {
91
for (int j = i; j < n; j++) {
93
for (int k = 0; k < m; k++) {
94
matATA[i][j] += localcopy[k][i] * localcopy[k][j];
96
matATA[j][i]=matATA[i][j];
100
double[][] matAAT = new double[m][m];
104
for (int i = 0; i < m; i++) {
105
for (int j = i; j < m; j++) {
107
for (int k = 0; k < n; k++) {
108
matAAT[i][j] += localcopy[i][k] * localcopy[j][k];
110
matAAT[j][i]=matAAT[i][j];
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();
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);
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();
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);
139
for (int i = 0; i < p; i++) {
140
singularValues[i] = Math.sqrt(Math.abs(singularValues[i]));
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);
150
cachedU.setColumnVector(i, tmp.mapMultiply(-1.0));
124
155
/** {@inheritDoc} */
125
public RealMatrix getU()
126
throws InvalidMatrixException {
128
if (cachedU == null) {
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];
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;
142
for (int j = 0; j < n; ++j) {
143
ei0[j] = (mainBidiagonal[i] * ei0[j] +
144
secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
148
final double lastMain = mainBidiagonal[n - 1];
149
for (int j = 0; j < n; ++j) {
150
ei1[j] *= lastMain / singularValues[j];
152
for (int i = n; i < m; ++i) {
153
iData[i] = new double[n];
156
transformer.getU().multiply(MatrixUtils.createRealMatrix(iData));
158
// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
159
cachedU = transformer.getU().multiply(eigenDecomposition.getV());
156
public RealMatrix getU() throws InvalidMatrixException {
164
157
// return the cached matrix
169
162
/** {@inheritDoc} */
170
public RealMatrix getUT()
171
throws InvalidMatrixException {
163
public RealMatrix getUT() throws InvalidMatrixException {
173
165
if (cachedUt == null) {
174
166
cachedUt = getU().transpose();
195
186
/** {@inheritDoc} */
196
public double[] getSingularValues()
197
throws InvalidMatrixException {
187
public double[] getSingularValues() throws InvalidMatrixException {
198
188
return singularValues.clone();
201
191
/** {@inheritDoc} */
202
public RealMatrix getV()
203
throws InvalidMatrixException {
205
if (cachedV == null) {
208
// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
209
cachedV = transformer.getV().multiply(eigenDecomposition.getV());
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];
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;
222
for (int j = 0; j < m; ++j) {
223
ei0[j] = (mainBidiagonal[i] * ei0[j] +
224
secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
228
final double lastMain = mainBidiagonal[m - 1];
229
for (int j = 0; j < m; ++j) {
230
ei1[j] *= lastMain / singularValues[j];
232
for (int i = m; i < n; ++i) {
233
iData[i] = new double[m];
236
transformer.getV().multiply(MatrixUtils.createRealMatrix(iData));
192
public RealMatrix getV() throws InvalidMatrixException {
241
193
// return the cached matrix
246
198
/** {@inheritDoc} */
247
public RealMatrix getVT()
248
throws InvalidMatrixException {
199
public RealMatrix getVT() throws InvalidMatrixException {
250
201
if (cachedVt == null) {
251
202
cachedVt = getV().transpose();
260
211
public RealMatrix getCovariance(final double minSingularValue) {
262
213
// get the number of singular values to consider
214
final int p = singularValues.length;
263
215
int dimension = 0;
264
while ((dimension < n) && (singularValues[dimension] >= minSingularValue)) {
216
while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) {
268
220
if (dimension == 0) {
269
221
throw MathRuntimeException.createIllegalArgumentException(
270
"cutoff singular value is {0}, should be at most {1}",
271
minSingularValue, singularValues[0]);
222
"cutoff singular value is {0}, should be at most {1}",
223
minSingularValue, singularValues[0]);
274
final double[][] data = new double[dimension][n];
226
final double[][] data = new double[dimension][p];
275
227
getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
276
228
/** {@inheritDoc} */
278
public void visit(final int row, final int column, final double value) {
230
public void visit(final int row, final int column,
231
final double value) {
279
232
data[row][column] = value / singularValues[row];
281
}, 0, dimension - 1, 0, n - 1);
234
}, 0, dimension - 1, 0, p - 1);
283
236
RealMatrix jv = new Array2DRowRealMatrix(data, false);
284
237
return jv.transpose().multiply(jv);
315
265
/** {@inheritDoc} */
316
266
public DecompositionSolver getSolver() {
317
return new Solver(singularValues, getUT(), getV(),
318
getRank() == singularValues.length);
267
return new Solver(singularValues, getUT(), getV(), getRank() == Math
321
271
/** Specialized solver. */
322
272
private static class Solver implements DecompositionSolver {
324
/** Singular values. */
325
private final double[] singularValues;
327
/** U<sup>T</sup> matrix of the decomposition. */
328
private final RealMatrix uT;
330
/** V matrix of the decomposition. */
331
private final RealMatrix v;
274
/** Pseudo-inverse of the initial matrix. */
275
private final RealMatrix pseudoInverse;
333
277
/** Singularity indicator. */
334
278
private boolean nonSingular;
337
281
* 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
343
private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v,
344
final boolean nonSingular) {
345
this.singularValues = singularValues;
348
this.nonSingular = nonSingular;
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
359
public double[] solve(final double[] b)
360
throws IllegalArgumentException, InvalidMatrixException {
362
if (b.length != singularValues.length) {
363
throw MathRuntimeException.createIllegalArgumentException(
364
"vector length mismatch: got {0} but expected {1}",
365
b.length, singularValues.length);
368
final double[] w = uT.operate(b);
282
* @param singularValues
285
* U<sup>T</sup> matrix of the decomposition
287
* V matrix of the decomposition
289
* singularity indicator
291
private Solver(final double[] singularValues, final RealMatrix uT,
292
final RealMatrix v, final boolean nonSingular) {
293
double[][] suT = uT.getData();
369
294
for (int i = 0; i < singularValues.length; ++i) {
370
final double si = singularValues[i];
372
throw new SingularMatrixException();
296
if (singularValues[i]>0) {
297
a=1.0 / singularValues[i];
301
final double[] suTi = suT[i];
302
for (int j = 0; j < suTi.length; ++j) {
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
306
pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));
307
this.nonSingular = nonSingular;
311
* Solve the linear equation A × X = B in least square sense.
313
* The m×n matrix A may not be square, the solution X is such that
314
* ||A × X - B|| is minimal.
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
322
public double[] solve(final double[] b) throws IllegalArgumentException {
323
return pseudoInverse.operate(b);
327
* Solve the linear equation A × X = B in least square sense.
329
* The m×n matrix A may not be square, the solution X is such that
330
* ||A × X - B|| is minimal.
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
388
338
public RealVector solve(final RealVector b)
389
throws IllegalArgumentException, InvalidMatrixException {
391
if (b.getDimension() != singularValues.length) {
392
throw MathRuntimeException.createIllegalArgumentException(
393
"vector length mismatch: got {0} but expected {1}",
394
b.getDimension(), singularValues.length);
397
final RealVector w = uT.operate(b);
398
for (int i = 0; i < singularValues.length; ++i) {
399
final double si = singularValues[i];
401
throw new SingularMatrixException();
403
w.setEntry(i, w.getEntry(i) / si);
339
throws IllegalArgumentException {
340
return pseudoInverse.operate(b);
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
344
* Solve the linear equation A × X = B in least square sense.
346
* The m×n matrix A may not be square, the solution X is such that
347
* ||A × X - B|| is minimal.
350
* right-hand side of the equation A × X = B
413
351
* @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
352
* @exception IllegalArgumentException
353
* if matrices dimensions don't match
417
355
public RealMatrix solve(final RealMatrix b)
418
throws IllegalArgumentException, InvalidMatrixException {
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");
427
final RealMatrix w = uT.multiply(b);
428
for (int i = 0; i < singularValues.length; ++i) {
429
final double si = singularValues[i];
431
throw new SingularMatrixException();
433
final double inv = 1.0 / si;
434
for (int j = 0; j < b.getColumnDimension(); ++j) {
435
w.multiplyEntry(i, j, inv);
438
return v.multiply(w);
356
throws IllegalArgumentException {
357
return pseudoInverse.multiply(b);
added
removed
src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
