source: src/main/java/agents/org/apache/commons/math/linear/SingularValueDecompositionImpl.java

Last change on this file was 1, checked in by Wouter Pasman, 7 years ago

Initial import : Genius 9.0.0
File size: 12.6 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package agents.org.apache.commons.math.linear;
19
20import agents.org.apache.commons.math.MathRuntimeException;
21import agents.org.apache.commons.math.exception.util.LocalizedFormats;
22import agents.org.apache.commons.math.util.FastMath;
23
24/**
25 * Calculates the compact Singular Value Decomposition of a matrix.
26 * <p>
27 * The Singular Value Decomposition of matrix A is a set of three matrices: U,
28 * &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>. Let A be
29 * a m &times; n matrix, then U is a m &times; p orthogonal matrix, &Sigma; is a
30 * p &times; p diagonal matrix with positive or null elements, V is a p &times;
31 * n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
32 * p=min(m,n).
33 * </p>
34 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
35 * @since 2.0
36 */
37public class SingularValueDecompositionImpl implements
38 SingularValueDecomposition {
39
40 /** Number of rows of the initial matrix. */
41 private int m;
42
43 /** Number of columns of the initial matrix. */
44 private int n;
45
46 /** Eigen decomposition of the tridiagonal matrix. */
47 private EigenDecomposition eigenDecomposition;
48
49 /** Singular values. */
50 private double[] singularValues;
51
52 /** Cached value of U. */
53 private RealMatrix cachedU;
54
55 /** Cached value of U<sup>T</sup>. */
56 private RealMatrix cachedUt;
57
58 /** Cached value of S. */
59 private RealMatrix cachedS;
60
61 /** Cached value of V. */
62 private RealMatrix cachedV;
63
64 /** Cached value of V<sup>T</sup>. */
65 private RealMatrix cachedVt;
66
67 /**
68 * Calculates the compact Singular Value Decomposition of the given matrix.
69 * @param matrix
70 * The matrix to decompose.
71 * @exception InvalidMatrixException
72 * (wrapping a
73 * {@link agents.org.apache.commons.math.ConvergenceException} if
74 * algorithm fails to converge
75 */
76 public SingularValueDecompositionImpl(final RealMatrix matrix)
77 throws InvalidMatrixException {
78
79 m = matrix.getRowDimension();
80 n = matrix.getColumnDimension();
81
82 cachedU = null;
83 cachedS = null;
84 cachedV = null;
85 cachedVt = null;
86
87 double[][] localcopy = matrix.getData();
88 double[][] matATA = new double[n][n];
89 //
90 // create A^T*A
91 //
92 for (int i = 0; i < n; i++) {
93 for (int j = i; j < n; j++) {
94 matATA[i][j] = 0.0;
95 for (int k = 0; k < m; k++) {
96 matATA[i][j] += localcopy[k][i] * localcopy[k][j];
97 }
98 matATA[j][i]=matATA[i][j];
99 }
100 }
101
102 double[][] matAAT = new double[m][m];
103 //
104 // create A*A^T
105 //
106 for (int i = 0; i < m; i++) {
107 for (int j = i; j < m; j++) {
108 matAAT[i][j] = 0.0;
109 for (int k = 0; k < n; k++) {
110 matAAT[i][j] += localcopy[i][k] * localcopy[j][k];
111 }
112 matAAT[j][i]=matAAT[i][j];
113 }
114 }
115 int p;
116 if (m>=n) {
117 p=n;
118 // compute eigen decomposition of A^T*A
119 eigenDecomposition = new EigenDecompositionImpl(
120 new Array2DRowRealMatrix(matATA),1.0);
121 singularValues = eigenDecomposition.getRealEigenvalues();
122 cachedV = eigenDecomposition.getV();
123 // compute eigen decomposition of A*A^T
124 eigenDecomposition = new EigenDecompositionImpl(
125 new Array2DRowRealMatrix(matAAT),1.0);
126 cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);
127 } else {
128 p=m;
129 // compute eigen decomposition of A*A^T
130 eigenDecomposition = new EigenDecompositionImpl(
131 new Array2DRowRealMatrix(matAAT),1.0);
132 singularValues = eigenDecomposition.getRealEigenvalues();
133 cachedU = eigenDecomposition.getV();
134
135 // compute eigen decomposition of A^T*A
136 eigenDecomposition = new EigenDecompositionImpl(
137 new Array2DRowRealMatrix(matATA),1.0);
138 cachedV = eigenDecomposition.getV().getSubMatrix(0,n-1,0,p-1);
139 }
140 for (int i = 0; i < p; i++) {
141 singularValues[i] = FastMath.sqrt(FastMath.abs(singularValues[i]));
142 }
143 // Up to this point, U and V are computed independently of each other.
144 // There still a sign indetermination of each column of, say, U.
145 // The sign is set such that A.V_i=sigma_i.U_i (i<=p)
146 // The right sign corresponds to a positive dot product of A.V_i and U_i
147 for (int i = 0; i < p; i++) {
148 RealVector tmp = cachedU.getColumnVector(i);
149 double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp);
150 if (product<0) {
151 cachedU.setColumnVector(i, tmp.mapMultiply(-1.0));
152 }
153 }
154 }
155
156 /** {@inheritDoc} */
157 public RealMatrix getU() throws InvalidMatrixException {
158 // return the cached matrix
159 return cachedU;
160
161 }
162
163 /** {@inheritDoc} */
164 public RealMatrix getUT() throws InvalidMatrixException {
165
166 if (cachedUt == null) {
167 cachedUt = getU().transpose();
168 }
169
170 // return the cached matrix
171 return cachedUt;
172
173 }
174
175 /** {@inheritDoc} */
176 public RealMatrix getS() throws InvalidMatrixException {
177
178 if (cachedS == null) {
179
180 // cache the matrix for subsequent calls
181 cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
182
183 }
184 return cachedS;
185 }
186
187 /** {@inheritDoc} */
188 public double[] getSingularValues() throws InvalidMatrixException {
189 return singularValues.clone();
190 }
191
192 /** {@inheritDoc} */
193 public RealMatrix getV() throws InvalidMatrixException {
194 // return the cached matrix
195 return cachedV;
196
197 }
198
199 /** {@inheritDoc} */
200 public RealMatrix getVT() throws InvalidMatrixException {
201
202 if (cachedVt == null) {
203 cachedVt = getV().transpose();
204 }
205
206 // return the cached matrix
207 return cachedVt;
208
209 }
210
211 /** {@inheritDoc} */
212 public RealMatrix getCovariance(final double minSingularValue) {
213
214 // get the number of singular values to consider
215 final int p = singularValues.length;
216 int dimension = 0;
217 while ((dimension < p) && (singularValues[dimension] >= minSingularValue)) {
218 ++dimension;
219 }
220
221 if (dimension == 0) {
222 throw MathRuntimeException.createIllegalArgumentException(
223 LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE,
224 minSingularValue, singularValues[0]);
225 }
226
227 final double[][] data = new double[dimension][p];
228 getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
229 /** {@inheritDoc} */
230 @Override
231 public void visit(final int row, final int column,
232 final double value) {
233 data[row][column] = value / singularValues[row];
234 }
235 }, 0, dimension - 1, 0, p - 1);
236
237 RealMatrix jv = new Array2DRowRealMatrix(data, false);
238 return jv.transpose().multiply(jv);
239
240 }
241
242 /** {@inheritDoc} */
243 public double getNorm() throws InvalidMatrixException {
244 return singularValues[0];
245 }
246
247 /** {@inheritDoc} */
248 public double getConditionNumber() throws InvalidMatrixException {
249 return singularValues[0] / singularValues[singularValues.length - 1];
250 }
251
252 /** {@inheritDoc} */
253 public int getRank() throws IllegalStateException {
254
255 final double threshold = FastMath.max(m, n) * FastMath.ulp(singularValues[0]);
256
257 for (int i = singularValues.length - 1; i >= 0; --i) {
258 if (singularValues[i] > threshold) {
259 return i + 1;
260 }
261 }
262 return 0;
263
264 }
265
266 /** {@inheritDoc} */
267 public DecompositionSolver getSolver() {
268 return new Solver(singularValues, getUT(), getV(), getRank() == Math
269 .max(m, n));
270 }
271
272 /** Specialized solver. */
273 private static class Solver implements DecompositionSolver {
274
275 /** Pseudo-inverse of the initial matrix. */
276 private final RealMatrix pseudoInverse;
277
278 /** Singularity indicator. */
279 private boolean nonSingular;
280
281 /**
282 * Build a solver from decomposed matrix.
283 * @param singularValues
284 * singularValues
285 * @param uT
286 * U<sup>T</sup> matrix of the decomposition
287 * @param v
288 * V matrix of the decomposition
289 * @param nonSingular
290 * singularity indicator
291 */
292 private Solver(final double[] singularValues, final RealMatrix uT,
293 final RealMatrix v, final boolean nonSingular) {
294 double[][] suT = uT.getData();
295 for (int i = 0; i < singularValues.length; ++i) {
296 final double a;
297 if (singularValues[i]>0) {
298 a=1.0 / singularValues[i];
299 } else {
300 a=0.0;
301 }
302 final double[] suTi = suT[i];
303 for (int j = 0; j < suTi.length; ++j) {
304 suTi[j] *= a;
305 }
306 }
307 pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));
308 this.nonSingular = nonSingular;
309 }
310
311 /**
312 * Solve the linear equation A &times; X = B in least square sense.
313 * <p>
314 * The m&times;n matrix A may not be square, the solution X is such that
315 * ||A &times; X - B|| is minimal.
316 * </p>
317 * @param b
318 * right-hand side of the equation A &times; X = B
319 * @return a vector X that minimizes the two norm of A &times; X - B
320 * @exception IllegalArgumentException
321 * if matrices dimensions don't match
322 */
323 public double[] solve(final double[] b) throws IllegalArgumentException {
324 return pseudoInverse.operate(b);
325 }
326
327 /**
328 * Solve the linear equation A &times; X = B in least square sense.
329 * <p>
330 * The m&times;n matrix A may not be square, the solution X is such that
331 * ||A &times; X - B|| is minimal.
332 * </p>
333 * @param b
334 * right-hand side of the equation A &times; X = B
335 * @return a vector X that minimizes the two norm of A &times; X - B
336 * @exception IllegalArgumentException
337 * if matrices dimensions don't match
338 */
339 public RealVector solve(final RealVector b)
340 throws IllegalArgumentException {
341 return pseudoInverse.operate(b);
342 }
343
344 /**
345 * Solve the linear equation A &times; X = B in least square sense.
346 * <p>
347 * The m&times;n matrix A may not be square, the solution X is such that
348 * ||A &times; X - B|| is minimal.
349 * </p>
350 * @param b
351 * right-hand side of the equation A &times; X = B
352 * @return a matrix X that minimizes the two norm of A &times; X - B
353 * @exception IllegalArgumentException
354 * if matrices dimensions don't match
355 */
356 public RealMatrix solve(final RealMatrix b)
357 throws IllegalArgumentException {
358 return pseudoInverse.multiply(b);
359 }
360
361 /**
362 * Check if the decomposed matrix is non-singular.
363 * @return true if the decomposed matrix is non-singular
364 */
365 public boolean isNonSingular() {
366 return nonSingular;
367 }
368
369 /**
370 * Get the pseudo-inverse of the decomposed matrix.
371 * @return inverse matrix
372 */
373 public RealMatrix getInverse() {
374 return pseudoInverse;
375 }
376
377 }
378
379}
Note: See TracBrowser for help on using the repository browser.