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 | package agents.anac.y2019.harddealer.math3.linear;
|
---|
18 |
|
---|
19 | import agents.anac.y2019.harddealer.math3.exception.NumberIsTooLargeException;
|
---|
20 | import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
|
---|
21 | import agents.anac.y2019.harddealer.math3.util.FastMath;
|
---|
22 | import agents.anac.y2019.harddealer.math3.util.Precision;
|
---|
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 | * Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be
|
---|
29 | * a m × n matrix, then U is a m × p orthogonal matrix, Σ is a
|
---|
30 | * p × p diagonal matrix with positive or null elements, V is a p ×
|
---|
31 | * n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
|
---|
32 | * p=min(m,n).
|
---|
33 | * </p>
|
---|
34 | * <p>This class is similar to the class with similar name from the
|
---|
35 | * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
|
---|
36 | * following changes:</p>
|
---|
37 | * <ul>
|
---|
38 | * <li>the {@code norm2} method which has been renamed as {@link #getNorm()
|
---|
39 | * getNorm},</li>
|
---|
40 | * <li>the {@code cond} method which has been renamed as {@link
|
---|
41 | * #getConditionNumber() getConditionNumber},</li>
|
---|
42 | * <li>the {@code rank} method which has been renamed as {@link #getRank()
|
---|
43 | * getRank},</li>
|
---|
44 | * <li>a {@link #getUT() getUT} method has been added,</li>
|
---|
45 | * <li>a {@link #getVT() getVT} method has been added,</li>
|
---|
46 | * <li>a {@link #getSolver() getSolver} method has been added,</li>
|
---|
47 | * <li>a {@link #getCovariance(double) getCovariance} method has been added.</li>
|
---|
48 | * </ul>
|
---|
49 | * @see <a href="http://mathworld.wolfram.com/SingularValueDecomposition.html">MathWorld</a>
|
---|
50 | * @see <a href="http://en.wikipedia.org/wiki/Singular_value_decomposition">Wikipedia</a>
|
---|
51 | * @since 2.0 (changed to concrete class in 3.0)
|
---|
52 | */
|
---|
53 | public class SingularValueDecomposition {
|
---|
54 | /** Relative threshold for small singular values. */
|
---|
55 | private static final double EPS = 0x1.0p-52;
|
---|
56 | /** Absolute threshold for small singular values. */
|
---|
57 | private static final double TINY = 0x1.0p-966;
|
---|
58 | /** Computed singular values. */
|
---|
59 | private final double[] singularValues;
|
---|
60 | /** max(row dimension, column dimension). */
|
---|
61 | private final int m;
|
---|
62 | /** min(row dimension, column dimension). */
|
---|
63 | private final int n;
|
---|
64 | /** Indicator for transposed matrix. */
|
---|
65 | private final boolean transposed;
|
---|
66 | /** Cached value of U matrix. */
|
---|
67 | private final RealMatrix cachedU;
|
---|
68 | /** Cached value of transposed U matrix. */
|
---|
69 | private RealMatrix cachedUt;
|
---|
70 | /** Cached value of S (diagonal) matrix. */
|
---|
71 | private RealMatrix cachedS;
|
---|
72 | /** Cached value of V matrix. */
|
---|
73 | private final RealMatrix cachedV;
|
---|
74 | /** Cached value of transposed V matrix. */
|
---|
75 | private RealMatrix cachedVt;
|
---|
76 | /**
|
---|
77 | * Tolerance value for small singular values, calculated once we have
|
---|
78 | * populated "singularValues".
|
---|
79 | **/
|
---|
80 | private final double tol;
|
---|
81 |
|
---|
82 | /**
|
---|
83 | * Calculates the compact Singular Value Decomposition of the given matrix.
|
---|
84 | *
|
---|
85 | * @param matrix Matrix to decompose.
|
---|
86 | */
|
---|
87 | public SingularValueDecomposition(final RealMatrix matrix) {
|
---|
88 | final double[][] A;
|
---|
89 |
|
---|
90 | // "m" is always the largest dimension.
|
---|
91 | if (matrix.getRowDimension() < matrix.getColumnDimension()) {
|
---|
92 | transposed = true;
|
---|
93 | A = matrix.transpose().getData();
|
---|
94 | m = matrix.getColumnDimension();
|
---|
95 | n = matrix.getRowDimension();
|
---|
96 | } else {
|
---|
97 | transposed = false;
|
---|
98 | A = matrix.getData();
|
---|
99 | m = matrix.getRowDimension();
|
---|
100 | n = matrix.getColumnDimension();
|
---|
101 | }
|
---|
102 |
|
---|
103 | singularValues = new double[n];
|
---|
104 | final double[][] U = new double[m][n];
|
---|
105 | final double[][] V = new double[n][n];
|
---|
106 | final double[] e = new double[n];
|
---|
107 | final double[] work = new double[m];
|
---|
108 | // Reduce A to bidiagonal form, storing the diagonal elements
|
---|
109 | // in s and the super-diagonal elements in e.
|
---|
110 | final int nct = FastMath.min(m - 1, n);
|
---|
111 | final int nrt = FastMath.max(0, n - 2);
|
---|
112 | for (int k = 0; k < FastMath.max(nct, nrt); k++) {
|
---|
113 | if (k < nct) {
|
---|
114 | // Compute the transformation for the k-th column and
|
---|
115 | // place the k-th diagonal in s[k].
|
---|
116 | // Compute 2-norm of k-th column without under/overflow.
|
---|
117 | singularValues[k] = 0;
|
---|
118 | for (int i = k; i < m; i++) {
|
---|
119 | singularValues[k] = FastMath.hypot(singularValues[k], A[i][k]);
|
---|
120 | }
|
---|
121 | if (singularValues[k] != 0) {
|
---|
122 | if (A[k][k] < 0) {
|
---|
123 | singularValues[k] = -singularValues[k];
|
---|
124 | }
|
---|
125 | for (int i = k; i < m; i++) {
|
---|
126 | A[i][k] /= singularValues[k];
|
---|
127 | }
|
---|
128 | A[k][k] += 1;
|
---|
129 | }
|
---|
130 | singularValues[k] = -singularValues[k];
|
---|
131 | }
|
---|
132 | for (int j = k + 1; j < n; j++) {
|
---|
133 | if (k < nct &&
|
---|
134 | singularValues[k] != 0) {
|
---|
135 | // Apply the transformation.
|
---|
136 | double t = 0;
|
---|
137 | for (int i = k; i < m; i++) {
|
---|
138 | t += A[i][k] * A[i][j];
|
---|
139 | }
|
---|
140 | t = -t / A[k][k];
|
---|
141 | for (int i = k; i < m; i++) {
|
---|
142 | A[i][j] += t * A[i][k];
|
---|
143 | }
|
---|
144 | }
|
---|
145 | // Place the k-th row of A into e for the
|
---|
146 | // subsequent calculation of the row transformation.
|
---|
147 | e[j] = A[k][j];
|
---|
148 | }
|
---|
149 | if (k < nct) {
|
---|
150 | // Place the transformation in U for subsequent back
|
---|
151 | // multiplication.
|
---|
152 | for (int i = k; i < m; i++) {
|
---|
153 | U[i][k] = A[i][k];
|
---|
154 | }
|
---|
155 | }
|
---|
156 | if (k < nrt) {
|
---|
157 | // Compute the k-th row transformation and place the
|
---|
158 | // k-th super-diagonal in e[k].
|
---|
159 | // Compute 2-norm without under/overflow.
|
---|
160 | e[k] = 0;
|
---|
161 | for (int i = k + 1; i < n; i++) {
|
---|
162 | e[k] = FastMath.hypot(e[k], e[i]);
|
---|
163 | }
|
---|
164 | if (e[k] != 0) {
|
---|
165 | if (e[k + 1] < 0) {
|
---|
166 | e[k] = -e[k];
|
---|
167 | }
|
---|
168 | for (int i = k + 1; i < n; i++) {
|
---|
169 | e[i] /= e[k];
|
---|
170 | }
|
---|
171 | e[k + 1] += 1;
|
---|
172 | }
|
---|
173 | e[k] = -e[k];
|
---|
174 | if (k + 1 < m &&
|
---|
175 | e[k] != 0) {
|
---|
176 | // Apply the transformation.
|
---|
177 | for (int i = k + 1; i < m; i++) {
|
---|
178 | work[i] = 0;
|
---|
179 | }
|
---|
180 | for (int j = k + 1; j < n; j++) {
|
---|
181 | for (int i = k + 1; i < m; i++) {
|
---|
182 | work[i] += e[j] * A[i][j];
|
---|
183 | }
|
---|
184 | }
|
---|
185 | for (int j = k + 1; j < n; j++) {
|
---|
186 | final double t = -e[j] / e[k + 1];
|
---|
187 | for (int i = k + 1; i < m; i++) {
|
---|
188 | A[i][j] += t * work[i];
|
---|
189 | }
|
---|
190 | }
|
---|
191 | }
|
---|
192 |
|
---|
193 | // Place the transformation in V for subsequent
|
---|
194 | // back multiplication.
|
---|
195 | for (int i = k + 1; i < n; i++) {
|
---|
196 | V[i][k] = e[i];
|
---|
197 | }
|
---|
198 | }
|
---|
199 | }
|
---|
200 | // Set up the final bidiagonal matrix or order p.
|
---|
201 | int p = n;
|
---|
202 | if (nct < n) {
|
---|
203 | singularValues[nct] = A[nct][nct];
|
---|
204 | }
|
---|
205 | if (m < p) {
|
---|
206 | singularValues[p - 1] = 0;
|
---|
207 | }
|
---|
208 | if (nrt + 1 < p) {
|
---|
209 | e[nrt] = A[nrt][p - 1];
|
---|
210 | }
|
---|
211 | e[p - 1] = 0;
|
---|
212 |
|
---|
213 | // Generate U.
|
---|
214 | for (int j = nct; j < n; j++) {
|
---|
215 | for (int i = 0; i < m; i++) {
|
---|
216 | U[i][j] = 0;
|
---|
217 | }
|
---|
218 | U[j][j] = 1;
|
---|
219 | }
|
---|
220 | for (int k = nct - 1; k >= 0; k--) {
|
---|
221 | if (singularValues[k] != 0) {
|
---|
222 | for (int j = k + 1; j < n; j++) {
|
---|
223 | double t = 0;
|
---|
224 | for (int i = k; i < m; i++) {
|
---|
225 | t += U[i][k] * U[i][j];
|
---|
226 | }
|
---|
227 | t = -t / U[k][k];
|
---|
228 | for (int i = k; i < m; i++) {
|
---|
229 | U[i][j] += t * U[i][k];
|
---|
230 | }
|
---|
231 | }
|
---|
232 | for (int i = k; i < m; i++) {
|
---|
233 | U[i][k] = -U[i][k];
|
---|
234 | }
|
---|
235 | U[k][k] = 1 + U[k][k];
|
---|
236 | for (int i = 0; i < k - 1; i++) {
|
---|
237 | U[i][k] = 0;
|
---|
238 | }
|
---|
239 | } else {
|
---|
240 | for (int i = 0; i < m; i++) {
|
---|
241 | U[i][k] = 0;
|
---|
242 | }
|
---|
243 | U[k][k] = 1;
|
---|
244 | }
|
---|
245 | }
|
---|
246 |
|
---|
247 | // Generate V.
|
---|
248 | for (int k = n - 1; k >= 0; k--) {
|
---|
249 | if (k < nrt &&
|
---|
250 | e[k] != 0) {
|
---|
251 | for (int j = k + 1; j < n; j++) {
|
---|
252 | double t = 0;
|
---|
253 | for (int i = k + 1; i < n; i++) {
|
---|
254 | t += V[i][k] * V[i][j];
|
---|
255 | }
|
---|
256 | t = -t / V[k + 1][k];
|
---|
257 | for (int i = k + 1; i < n; i++) {
|
---|
258 | V[i][j] += t * V[i][k];
|
---|
259 | }
|
---|
260 | }
|
---|
261 | }
|
---|
262 | for (int i = 0; i < n; i++) {
|
---|
263 | V[i][k] = 0;
|
---|
264 | }
|
---|
265 | V[k][k] = 1;
|
---|
266 | }
|
---|
267 |
|
---|
268 | // Main iteration loop for the singular values.
|
---|
269 | final int pp = p - 1;
|
---|
270 | while (p > 0) {
|
---|
271 | int k;
|
---|
272 | int kase;
|
---|
273 | // Here is where a test for too many iterations would go.
|
---|
274 | // This section of the program inspects for
|
---|
275 | // negligible elements in the s and e arrays. On
|
---|
276 | // completion the variables kase and k are set as follows.
|
---|
277 | // kase = 1 if s(p) and e[k-1] are negligible and k<p
|
---|
278 | // kase = 2 if s(k) is negligible and k<p
|
---|
279 | // kase = 3 if e[k-1] is negligible, k<p, and
|
---|
280 | // s(k), ..., s(p) are not negligible (qr step).
|
---|
281 | // kase = 4 if e(p-1) is negligible (convergence).
|
---|
282 | for (k = p - 2; k >= 0; k--) {
|
---|
283 | final double threshold
|
---|
284 | = TINY + EPS * (FastMath.abs(singularValues[k]) +
|
---|
285 | FastMath.abs(singularValues[k + 1]));
|
---|
286 |
|
---|
287 | // the following condition is written this way in order
|
---|
288 | // to break out of the loop when NaN occurs, writing it
|
---|
289 | // as "if (FastMath.abs(e[k]) <= threshold)" would loop
|
---|
290 | // indefinitely in case of NaNs because comparison on NaNs
|
---|
291 | // always return false, regardless of what is checked
|
---|
292 | // see issue MATH-947
|
---|
293 | if (!(FastMath.abs(e[k]) > threshold)) {
|
---|
294 | e[k] = 0;
|
---|
295 | break;
|
---|
296 | }
|
---|
297 |
|
---|
298 | }
|
---|
299 |
|
---|
300 | if (k == p - 2) {
|
---|
301 | kase = 4;
|
---|
302 | } else {
|
---|
303 | int ks;
|
---|
304 | for (ks = p - 1; ks >= k; ks--) {
|
---|
305 | if (ks == k) {
|
---|
306 | break;
|
---|
307 | }
|
---|
308 | final double t = (ks != p ? FastMath.abs(e[ks]) : 0) +
|
---|
309 | (ks != k + 1 ? FastMath.abs(e[ks - 1]) : 0);
|
---|
310 | if (FastMath.abs(singularValues[ks]) <= TINY + EPS * t) {
|
---|
311 | singularValues[ks] = 0;
|
---|
312 | break;
|
---|
313 | }
|
---|
314 | }
|
---|
315 | if (ks == k) {
|
---|
316 | kase = 3;
|
---|
317 | } else if (ks == p - 1) {
|
---|
318 | kase = 1;
|
---|
319 | } else {
|
---|
320 | kase = 2;
|
---|
321 | k = ks;
|
---|
322 | }
|
---|
323 | }
|
---|
324 | k++;
|
---|
325 | // Perform the task indicated by kase.
|
---|
326 | switch (kase) {
|
---|
327 | // Deflate negligible s(p).
|
---|
328 | case 1: {
|
---|
329 | double f = e[p - 2];
|
---|
330 | e[p - 2] = 0;
|
---|
331 | for (int j = p - 2; j >= k; j--) {
|
---|
332 | double t = FastMath.hypot(singularValues[j], f);
|
---|
333 | final double cs = singularValues[j] / t;
|
---|
334 | final double sn = f / t;
|
---|
335 | singularValues[j] = t;
|
---|
336 | if (j != k) {
|
---|
337 | f = -sn * e[j - 1];
|
---|
338 | e[j - 1] = cs * e[j - 1];
|
---|
339 | }
|
---|
340 |
|
---|
341 | for (int i = 0; i < n; i++) {
|
---|
342 | t = cs * V[i][j] + sn * V[i][p - 1];
|
---|
343 | V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
|
---|
344 | V[i][j] = t;
|
---|
345 | }
|
---|
346 | }
|
---|
347 | }
|
---|
348 | break;
|
---|
349 | // Split at negligible s(k).
|
---|
350 | case 2: {
|
---|
351 | double f = e[k - 1];
|
---|
352 | e[k - 1] = 0;
|
---|
353 | for (int j = k; j < p; j++) {
|
---|
354 | double t = FastMath.hypot(singularValues[j], f);
|
---|
355 | final double cs = singularValues[j] / t;
|
---|
356 | final double sn = f / t;
|
---|
357 | singularValues[j] = t;
|
---|
358 | f = -sn * e[j];
|
---|
359 | e[j] = cs * e[j];
|
---|
360 |
|
---|
361 | for (int i = 0; i < m; i++) {
|
---|
362 | t = cs * U[i][j] + sn * U[i][k - 1];
|
---|
363 | U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
|
---|
364 | U[i][j] = t;
|
---|
365 | }
|
---|
366 | }
|
---|
367 | }
|
---|
368 | break;
|
---|
369 | // Perform one qr step.
|
---|
370 | case 3: {
|
---|
371 | // Calculate the shift.
|
---|
372 | final double maxPm1Pm2 = FastMath.max(FastMath.abs(singularValues[p - 1]),
|
---|
373 | FastMath.abs(singularValues[p - 2]));
|
---|
374 | final double scale = FastMath.max(FastMath.max(FastMath.max(maxPm1Pm2,
|
---|
375 | FastMath.abs(e[p - 2])),
|
---|
376 | FastMath.abs(singularValues[k])),
|
---|
377 | FastMath.abs(e[k]));
|
---|
378 | final double sp = singularValues[p - 1] / scale;
|
---|
379 | final double spm1 = singularValues[p - 2] / scale;
|
---|
380 | final double epm1 = e[p - 2] / scale;
|
---|
381 | final double sk = singularValues[k] / scale;
|
---|
382 | final double ek = e[k] / scale;
|
---|
383 | final double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
|
---|
384 | final double c = (sp * epm1) * (sp * epm1);
|
---|
385 | double shift = 0;
|
---|
386 | if (b != 0 ||
|
---|
387 | c != 0) {
|
---|
388 | shift = FastMath.sqrt(b * b + c);
|
---|
389 | if (b < 0) {
|
---|
390 | shift = -shift;
|
---|
391 | }
|
---|
392 | shift = c / (b + shift);
|
---|
393 | }
|
---|
394 | double f = (sk + sp) * (sk - sp) + shift;
|
---|
395 | double g = sk * ek;
|
---|
396 | // Chase zeros.
|
---|
397 | for (int j = k; j < p - 1; j++) {
|
---|
398 | double t = FastMath.hypot(f, g);
|
---|
399 | double cs = f / t;
|
---|
400 | double sn = g / t;
|
---|
401 | if (j != k) {
|
---|
402 | e[j - 1] = t;
|
---|
403 | }
|
---|
404 | f = cs * singularValues[j] + sn * e[j];
|
---|
405 | e[j] = cs * e[j] - sn * singularValues[j];
|
---|
406 | g = sn * singularValues[j + 1];
|
---|
407 | singularValues[j + 1] = cs * singularValues[j + 1];
|
---|
408 |
|
---|
409 | for (int i = 0; i < n; i++) {
|
---|
410 | t = cs * V[i][j] + sn * V[i][j + 1];
|
---|
411 | V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
|
---|
412 | V[i][j] = t;
|
---|
413 | }
|
---|
414 | t = FastMath.hypot(f, g);
|
---|
415 | cs = f / t;
|
---|
416 | sn = g / t;
|
---|
417 | singularValues[j] = t;
|
---|
418 | f = cs * e[j] + sn * singularValues[j + 1];
|
---|
419 | singularValues[j + 1] = -sn * e[j] + cs * singularValues[j + 1];
|
---|
420 | g = sn * e[j + 1];
|
---|
421 | e[j + 1] = cs * e[j + 1];
|
---|
422 | if (j < m - 1) {
|
---|
423 | for (int i = 0; i < m; i++) {
|
---|
424 | t = cs * U[i][j] + sn * U[i][j + 1];
|
---|
425 | U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
|
---|
426 | U[i][j] = t;
|
---|
427 | }
|
---|
428 | }
|
---|
429 | }
|
---|
430 | e[p - 2] = f;
|
---|
431 | }
|
---|
432 | break;
|
---|
433 | // Convergence.
|
---|
434 | default: {
|
---|
435 | // Make the singular values positive.
|
---|
436 | if (singularValues[k] <= 0) {
|
---|
437 | singularValues[k] = singularValues[k] < 0 ? -singularValues[k] : 0;
|
---|
438 |
|
---|
439 | for (int i = 0; i <= pp; i++) {
|
---|
440 | V[i][k] = -V[i][k];
|
---|
441 | }
|
---|
442 | }
|
---|
443 | // Order the singular values.
|
---|
444 | while (k < pp) {
|
---|
445 | if (singularValues[k] >= singularValues[k + 1]) {
|
---|
446 | break;
|
---|
447 | }
|
---|
448 | double t = singularValues[k];
|
---|
449 | singularValues[k] = singularValues[k + 1];
|
---|
450 | singularValues[k + 1] = t;
|
---|
451 | if (k < n - 1) {
|
---|
452 | for (int i = 0; i < n; i++) {
|
---|
453 | t = V[i][k + 1];
|
---|
454 | V[i][k + 1] = V[i][k];
|
---|
455 | V[i][k] = t;
|
---|
456 | }
|
---|
457 | }
|
---|
458 | if (k < m - 1) {
|
---|
459 | for (int i = 0; i < m; i++) {
|
---|
460 | t = U[i][k + 1];
|
---|
461 | U[i][k + 1] = U[i][k];
|
---|
462 | U[i][k] = t;
|
---|
463 | }
|
---|
464 | }
|
---|
465 | k++;
|
---|
466 | }
|
---|
467 | p--;
|
---|
468 | }
|
---|
469 | break;
|
---|
470 | }
|
---|
471 | }
|
---|
472 |
|
---|
473 | // Set the small value tolerance used to calculate rank and pseudo-inverse
|
---|
474 | tol = FastMath.max(m * singularValues[0] * EPS,
|
---|
475 | FastMath.sqrt(Precision.SAFE_MIN));
|
---|
476 |
|
---|
477 | if (!transposed) {
|
---|
478 | cachedU = MatrixUtils.createRealMatrix(U);
|
---|
479 | cachedV = MatrixUtils.createRealMatrix(V);
|
---|
480 | } else {
|
---|
481 | cachedU = MatrixUtils.createRealMatrix(V);
|
---|
482 | cachedV = MatrixUtils.createRealMatrix(U);
|
---|
483 | }
|
---|
484 | }
|
---|
485 |
|
---|
486 | /**
|
---|
487 | * Returns the matrix U of the decomposition.
|
---|
488 | * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
|
---|
489 | * @return the U matrix
|
---|
490 | * @see #getUT()
|
---|
491 | */
|
---|
492 | public RealMatrix getU() {
|
---|
493 | // return the cached matrix
|
---|
494 | return cachedU;
|
---|
495 |
|
---|
496 | }
|
---|
497 |
|
---|
498 | /**
|
---|
499 | * Returns the transpose of the matrix U of the decomposition.
|
---|
500 | * <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
|
---|
501 | * @return the U matrix (or null if decomposed matrix is singular)
|
---|
502 | * @see #getU()
|
---|
503 | */
|
---|
504 | public RealMatrix getUT() {
|
---|
505 | if (cachedUt == null) {
|
---|
506 | cachedUt = getU().transpose();
|
---|
507 | }
|
---|
508 | // return the cached matrix
|
---|
509 | return cachedUt;
|
---|
510 | }
|
---|
511 |
|
---|
512 | /**
|
---|
513 | * Returns the diagonal matrix Σ of the decomposition.
|
---|
514 | * <p>Σ is a diagonal matrix. The singular values are provided in
|
---|
515 | * non-increasing order, for compatibility with Jama.</p>
|
---|
516 | * @return the Σ matrix
|
---|
517 | */
|
---|
518 | public RealMatrix getS() {
|
---|
519 | if (cachedS == null) {
|
---|
520 | // cache the matrix for subsequent calls
|
---|
521 | cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
|
---|
522 | }
|
---|
523 | return cachedS;
|
---|
524 | }
|
---|
525 |
|
---|
526 | /**
|
---|
527 | * Returns the diagonal elements of the matrix Σ of the decomposition.
|
---|
528 | * <p>The singular values are provided in non-increasing order, for
|
---|
529 | * compatibility with Jama.</p>
|
---|
530 | * @return the diagonal elements of the Σ matrix
|
---|
531 | */
|
---|
532 | public double[] getSingularValues() {
|
---|
533 | return singularValues.clone();
|
---|
534 | }
|
---|
535 |
|
---|
536 | /**
|
---|
537 | * Returns the matrix V of the decomposition.
|
---|
538 | * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
|
---|
539 | * @return the V matrix (or null if decomposed matrix is singular)
|
---|
540 | * @see #getVT()
|
---|
541 | */
|
---|
542 | public RealMatrix getV() {
|
---|
543 | // return the cached matrix
|
---|
544 | return cachedV;
|
---|
545 | }
|
---|
546 |
|
---|
547 | /**
|
---|
548 | * Returns the transpose of the matrix V of the decomposition.
|
---|
549 | * <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
|
---|
550 | * @return the V matrix (or null if decomposed matrix is singular)
|
---|
551 | * @see #getV()
|
---|
552 | */
|
---|
553 | public RealMatrix getVT() {
|
---|
554 | if (cachedVt == null) {
|
---|
555 | cachedVt = getV().transpose();
|
---|
556 | }
|
---|
557 | // return the cached matrix
|
---|
558 | return cachedVt;
|
---|
559 | }
|
---|
560 |
|
---|
561 | /**
|
---|
562 | * Returns the n × n covariance matrix.
|
---|
563 | * <p>The covariance matrix is V × J × V<sup>T</sup>
|
---|
564 | * where J is the diagonal matrix of the inverse of the squares of
|
---|
565 | * the singular values.</p>
|
---|
566 | * @param minSingularValue value below which singular values are ignored
|
---|
567 | * (a 0 or negative value implies all singular value will be used)
|
---|
568 | * @return covariance matrix
|
---|
569 | * @exception IllegalArgumentException if minSingularValue is larger than
|
---|
570 | * the largest singular value, meaning all singular values are ignored
|
---|
571 | */
|
---|
572 | public RealMatrix getCovariance(final double minSingularValue) {
|
---|
573 | // get the number of singular values to consider
|
---|
574 | final int p = singularValues.length;
|
---|
575 | int dimension = 0;
|
---|
576 | while (dimension < p &&
|
---|
577 | singularValues[dimension] >= minSingularValue) {
|
---|
578 | ++dimension;
|
---|
579 | }
|
---|
580 |
|
---|
581 | if (dimension == 0) {
|
---|
582 | throw new NumberIsTooLargeException(LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE,
|
---|
583 | minSingularValue, singularValues[0], true);
|
---|
584 | }
|
---|
585 |
|
---|
586 | final double[][] data = new double[dimension][p];
|
---|
587 | getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
|
---|
588 | /** {@inheritDoc} */
|
---|
589 | @Override
|
---|
590 | public void visit(final int row, final int column,
|
---|
591 | final double value) {
|
---|
592 | data[row][column] = value / singularValues[row];
|
---|
593 | }
|
---|
594 | }, 0, dimension - 1, 0, p - 1);
|
---|
595 |
|
---|
596 | RealMatrix jv = new Array2DRowRealMatrix(data, false);
|
---|
597 | return jv.transpose().multiply(jv);
|
---|
598 | }
|
---|
599 |
|
---|
600 | /**
|
---|
601 | * Returns the L<sub>2</sub> norm of the matrix.
|
---|
602 | * <p>The L<sub>2</sub> norm is max(|A × u|<sub>2</sub> /
|
---|
603 | * |u|<sub>2</sub>), where |.|<sub>2</sub> denotes the vectorial 2-norm
|
---|
604 | * (i.e. the traditional euclidian norm).</p>
|
---|
605 | * @return norm
|
---|
606 | */
|
---|
607 | public double getNorm() {
|
---|
608 | return singularValues[0];
|
---|
609 | }
|
---|
610 |
|
---|
611 | /**
|
---|
612 | * Return the condition number of the matrix.
|
---|
613 | * @return condition number of the matrix
|
---|
614 | */
|
---|
615 | public double getConditionNumber() {
|
---|
616 | return singularValues[0] / singularValues[n - 1];
|
---|
617 | }
|
---|
618 |
|
---|
619 | /**
|
---|
620 | * Computes the inverse of the condition number.
|
---|
621 | * In cases of rank deficiency, the {@link #getConditionNumber() condition
|
---|
622 | * number} will become undefined.
|
---|
623 | *
|
---|
624 | * @return the inverse of the condition number.
|
---|
625 | */
|
---|
626 | public double getInverseConditionNumber() {
|
---|
627 | return singularValues[n - 1] / singularValues[0];
|
---|
628 | }
|
---|
629 |
|
---|
630 | /**
|
---|
631 | * Return the effective numerical matrix rank.
|
---|
632 | * <p>The effective numerical rank is the number of non-negligible
|
---|
633 | * singular values. The threshold used to identify non-negligible
|
---|
634 | * terms is max(m,n) × ulp(s<sub>1</sub>) where ulp(s<sub>1</sub>)
|
---|
635 | * is the least significant bit of the largest singular value.</p>
|
---|
636 | * @return effective numerical matrix rank
|
---|
637 | */
|
---|
638 | public int getRank() {
|
---|
639 | int r = 0;
|
---|
640 | for (int i = 0; i < singularValues.length; i++) {
|
---|
641 | if (singularValues[i] > tol) {
|
---|
642 | r++;
|
---|
643 | }
|
---|
644 | }
|
---|
645 | return r;
|
---|
646 | }
|
---|
647 |
|
---|
648 | /**
|
---|
649 | * Get a solver for finding the A × X = B solution in least square sense.
|
---|
650 | * @return a solver
|
---|
651 | */
|
---|
652 | public DecompositionSolver getSolver() {
|
---|
653 | return new Solver(singularValues, getUT(), getV(), getRank() == m, tol);
|
---|
654 | }
|
---|
655 |
|
---|
656 | /** Specialized solver. */
|
---|
657 | private static class Solver implements DecompositionSolver {
|
---|
658 | /** Pseudo-inverse of the initial matrix. */
|
---|
659 | private final RealMatrix pseudoInverse;
|
---|
660 | /** Singularity indicator. */
|
---|
661 | private boolean nonSingular;
|
---|
662 |
|
---|
663 | /**
|
---|
664 | * Build a solver from decomposed matrix.
|
---|
665 | *
|
---|
666 | * @param singularValues Singular values.
|
---|
667 | * @param uT U<sup>T</sup> matrix of the decomposition.
|
---|
668 | * @param v V matrix of the decomposition.
|
---|
669 | * @param nonSingular Singularity indicator.
|
---|
670 | * @param tol tolerance for singular values
|
---|
671 | */
|
---|
672 | private Solver(final double[] singularValues, final RealMatrix uT,
|
---|
673 | final RealMatrix v, final boolean nonSingular, final double tol) {
|
---|
674 | final double[][] suT = uT.getData();
|
---|
675 | for (int i = 0; i < singularValues.length; ++i) {
|
---|
676 | final double a;
|
---|
677 | if (singularValues[i] > tol) {
|
---|
678 | a = 1 / singularValues[i];
|
---|
679 | } else {
|
---|
680 | a = 0;
|
---|
681 | }
|
---|
682 | final double[] suTi = suT[i];
|
---|
683 | for (int j = 0; j < suTi.length; ++j) {
|
---|
684 | suTi[j] *= a;
|
---|
685 | }
|
---|
686 | }
|
---|
687 | pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));
|
---|
688 | this.nonSingular = nonSingular;
|
---|
689 | }
|
---|
690 |
|
---|
691 | /**
|
---|
692 | * Solve the linear equation A × X = B in least square sense.
|
---|
693 | * <p>
|
---|
694 | * The m×n matrix A may not be square, the solution X is such that
|
---|
695 | * ||A × X - B|| is minimal.
|
---|
696 | * </p>
|
---|
697 | * @param b Right-hand side of the equation A × X = B
|
---|
698 | * @return a vector X that minimizes the two norm of A × X - B
|
---|
699 | * @throws agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException
|
---|
700 | * if the matrices dimensions do not match.
|
---|
701 | */
|
---|
702 | public RealVector solve(final RealVector b) {
|
---|
703 | return pseudoInverse.operate(b);
|
---|
704 | }
|
---|
705 |
|
---|
706 | /**
|
---|
707 | * Solve the linear equation A × X = B in least square sense.
|
---|
708 | * <p>
|
---|
709 | * The m×n matrix A may not be square, the solution X is such that
|
---|
710 | * ||A × X - B|| is minimal.
|
---|
711 | * </p>
|
---|
712 | *
|
---|
713 | * @param b Right-hand side of the equation A × X = B
|
---|
714 | * @return a matrix X that minimizes the two norm of A × X - B
|
---|
715 | * @throws agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException
|
---|
716 | * if the matrices dimensions do not match.
|
---|
717 | */
|
---|
718 | public RealMatrix solve(final RealMatrix b) {
|
---|
719 | return pseudoInverse.multiply(b);
|
---|
720 | }
|
---|
721 |
|
---|
722 | /**
|
---|
723 | * Check if the decomposed matrix is non-singular.
|
---|
724 | *
|
---|
725 | * @return {@code true} if the decomposed matrix is non-singular.
|
---|
726 | */
|
---|
727 | public boolean isNonSingular() {
|
---|
728 | return nonSingular;
|
---|
729 | }
|
---|
730 |
|
---|
731 | /**
|
---|
732 | * Get the pseudo-inverse of the decomposed matrix.
|
---|
733 | *
|
---|
734 | * @return the inverse matrix.
|
---|
735 | */
|
---|
736 | public RealMatrix getInverse() {
|
---|
737 | return pseudoInverse;
|
---|
738 | }
|
---|
739 | }
|
---|
740 | }
|
---|