1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.linear;
19
20 import org.apache.commons.math.ConvergenceException;
21 import org.apache.commons.math.MathRuntimeException;
22 import org.apache.commons.math.util.MathUtils;
23
24
25
26
27
28
29
30
31
32
33
34
35 public class SingularValueDecompositionImpl implements SingularValueDecomposition {
36
37
38 private int m;
39
40
41 private int n;
42
43
44 private BiDiagonalTransformer transformer;
45
46
47 private double[] mainBidiagonal;
48
49
50 private double[] secondaryBidiagonal;
51
52
53 private double[] mainTridiagonal;
54
55
56 private double[] secondaryTridiagonal;
57
58
59 private EigenDecomposition eigenDecomposition;
60
61
62 private double[] singularValues;
63
64
65 private RealMatrix cachedU;
66
67
68 private RealMatrix cachedUt;
69
70
71 private RealMatrix cachedS;
72
73
74 private RealMatrix cachedV;
75
76
77 private RealMatrix cachedVt;
78
79
80
81
82
83
84
85 public SingularValueDecompositionImpl(RealMatrix matrix)
86 throws InvalidMatrixException {
87
88 m = matrix.getRowDimension();
89 n = matrix.getColumnDimension();
90
91 cachedU = null;
92 cachedS = null;
93 cachedV = null;
94 cachedVt = null;
95
96
97 transformer = new BiDiagonalTransformer(matrix);
98 mainBidiagonal = transformer.getMainDiagonalRef();
99 secondaryBidiagonal = transformer.getSecondaryDiagonalRef();
100
101
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
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
125 public RealMatrix getU()
126 throws InvalidMatrixException {
127
128 if (cachedU == null) {
129
130 if (m >= n) {
131
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
138
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
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
159 cachedU = transformer.getU().multiply(eigenDecomposition.getV());
160 }
161
162 }
163
164
165 return cachedU;
166
167 }
168
169
170 public RealMatrix getUT()
171 throws InvalidMatrixException {
172
173 if (cachedUt == null) {
174 cachedUt = getU().transpose();
175 }
176
177
178 return cachedUt;
179
180 }
181
182
183 public RealMatrix getS()
184 throws InvalidMatrixException {
185
186 if (cachedS == null) {
187
188
189 cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
190
191 }
192 return cachedS;
193 }
194
195
196 public double[] getSingularValues()
197 throws InvalidMatrixException {
198 return singularValues.clone();
199 }
200
201
202 public RealMatrix getV()
203 throws InvalidMatrixException {
204
205 if (cachedV == null) {
206
207 if (m >= n) {
208
209 cachedV = transformer.getV().multiply(eigenDecomposition.getV());
210 } else {
211
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
218
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
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
242 return cachedV;
243
244 }
245
246
247 public RealMatrix getVT()
248 throws InvalidMatrixException {
249
250 if (cachedVt == null) {
251 cachedVt = getV().transpose();
252 }
253
254
255 return cachedVt;
256
257 }
258
259
260 public RealMatrix getCovariance(final double minSingularValue) {
261
262
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
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
289 public double getNorm()
290 throws InvalidMatrixException {
291 return singularValues[0];
292 }
293
294
295 public double getConditionNumber()
296 throws InvalidMatrixException {
297 return singularValues[0] / singularValues[singularValues.length - 1];
298 }
299
300
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
316 public DecompositionSolver getSolver() {
317 return new Solver(singularValues, getUT(), getV(),
318 getRank() == singularValues.length);
319 }
320
321
322 private static class Solver implements DecompositionSolver {
323
324
325 private final double[] singularValues;
326
327
328 private final RealMatrix uT;
329
330
331 private final RealMatrix v;
332
333
334 private boolean nonSingular;
335
336
337
338
339
340
341
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
352
353
354
355
356
357
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
381
382
383
384
385
386
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
410
411
412
413
414
415
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
444
445
446 public boolean isNonSingular() {
447 return nonSingular;
448 }
449
450
451
452
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 }