1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.linear;
24
25 import org.hipparchus.exception.LocalizedCoreFormats;
26 import org.hipparchus.exception.MathIllegalArgumentException;
27 import org.hipparchus.exception.MathRuntimeException;
28 import org.hipparchus.util.FastMath;
29 import org.hipparchus.util.MathArrays;
30 import org.junit.jupiter.api.AfterEach;
31 import org.junit.jupiter.api.BeforeEach;
32 import org.junit.jupiter.api.Test;
33
34 import java.util.Arrays;
35 import java.util.Random;
36
37 import static org.junit.jupiter.api.Assertions.assertEquals;
38 import static org.junit.jupiter.api.Assertions.assertFalse;
39 import static org.junit.jupiter.api.Assertions.assertThrows;
40 import static org.junit.jupiter.api.Assertions.assertTrue;
41 import static org.junit.jupiter.api.Assertions.fail;
42
43 public class EigenDecompositionSymmetricTest {
44
45 private double[] refValues;
46 private RealMatrix matrix;
47
48 @Test
49 void testDimension1() {
50 RealMatrix matrix =
51 MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
52 EigenDecompositionSymmetric ed;
53 ed = new EigenDecompositionSymmetric(matrix);
54 assertEquals(1.5, ed.getEigenvalue(0), 1.0e-15);
55 }
56
57 @Test
58 void testDimension2() {
59 RealMatrix matrix =
60 MatrixUtils.createRealMatrix(new double[][] {
61 { 59.0, 12.0 },
62 { 12.0, 66.0 }
63 });
64 EigenDecompositionSymmetric ed;
65 ed = new EigenDecompositionSymmetric(matrix);
66 assertEquals(75.0, ed.getEigenvalue(0), 1.0e-15);
67 assertEquals(50.0, ed.getEigenvalue(1), 1.0e-15);
68 }
69
70 @Test
71 void testDimension3() {
72 RealMatrix matrix =
73 MatrixUtils.createRealMatrix(new double[][] {
74 { 39632.0, -4824.0, -16560.0 },
75 { -4824.0, 8693.0, 7920.0 },
76 { -16560.0, 7920.0, 17300.0 }
77 });
78 EigenDecompositionSymmetric ed;
79 ed = new EigenDecompositionSymmetric(matrix);
80 assertEquals(50000.0, ed.getEigenvalue(0), 3.0e-11);
81 assertEquals(12500.0, ed.getEigenvalue(1), 3.0e-11);
82 assertEquals( 3125.0, ed.getEigenvalue(2), 3.0e-11);
83 }
84
85 @Test
86 void testDimension3MultipleRoot() {
87 RealMatrix matrix =
88 MatrixUtils.createRealMatrix(new double[][] {
89 { 5, 10, 15 },
90 { 10, 20, 30 },
91 { 15, 30, 45 }
92 });
93 EigenDecompositionSymmetric ed;
94 ed = new EigenDecompositionSymmetric(matrix);
95 assertEquals(70.0, ed.getEigenvalue(0), 3.0e-11);
96 assertEquals(0.0, ed.getEigenvalue(1), 3.0e-11);
97 assertEquals(0.0, ed.getEigenvalue(2), 3.0e-11);
98 assertEquals(matrix.getRowDimension(), ed.getSolver().getRowDimension());
99 assertEquals(matrix.getColumnDimension(), ed.getSolver().getColumnDimension());
100 }
101
102 @Test
103 void testDimension4WithSplit() {
104 RealMatrix matrix =
105 MatrixUtils.createRealMatrix(new double[][] {
106 { 0.784, -0.288, 0.000, 0.000 },
107 { -0.288, 0.616, 0.000, 0.000 },
108 { 0.000, 0.000, 0.164, -0.048 },
109 { 0.000, 0.000, -0.048, 0.136 }
110 });
111 EigenDecompositionSymmetric ed;
112 ed = new EigenDecompositionSymmetric(matrix);
113 assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
114 assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
115 assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
116 assertEquals(0.1, ed.getEigenvalue(3), 1.0e-15);
117 assertEquals(matrix.getRowDimension(), ed.getSolver().getRowDimension());
118 assertEquals(matrix.getColumnDimension(), ed.getSolver().getColumnDimension());
119 }
120
121 @Test
122 void testDimension4WithoutSplit() {
123 RealMatrix matrix =
124 MatrixUtils.createRealMatrix(new double[][] {
125 { 0.5608, -0.2016, 0.1152, -0.2976 },
126 { -0.2016, 0.4432, -0.2304, 0.1152 },
127 { 0.1152, -0.2304, 0.3088, -0.1344 },
128 { -0.2976, 0.1152, -0.1344, 0.3872 }
129 });
130 EigenDecompositionSymmetric ed;
131 ed = new EigenDecompositionSymmetric(matrix);
132 assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
133 assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
134 assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
135 assertEquals(0.1, ed.getEigenvalue(3), 1.0e-15);
136 assertEquals(matrix.getRowDimension(), ed.getSolver().getRowDimension());
137 assertEquals(matrix.getColumnDimension(), ed.getSolver().getColumnDimension());
138 }
139
140 @Test
141 void testMath308() {
142
143 double[] mainTridiagonal = {
144 22.330154644539597, 46.65485522478641, 17.393672330044705, 54.46687435351116, 80.17800767709437
145 };
146 double[] secondaryTridiagonal = {
147 13.04450406501361, -5.977590941539671, 2.9040909856707517, 7.1570352792841225
148 };
149
150
151
152 double[] refEigenValues = {
153 82.044413207204002, 53.456697699894512, 52.536278520113882, 18.847969733754262, 14.138204224043099
154 };
155 RealVector[] refEigenVectors = {
156 new ArrayRealVector(new double[] { -0.000462690386766, -0.002118073109055, 0.011530080757413, 0.252322434584915, 0.967572088232592 }),
157 new ArrayRealVector(new double[] { 0.314647769490148, 0.750806415553905, -0.167700312025760, -0.537092972407375, 0.143854968127780 }),
158 new ArrayRealVector(new double[] { 0.222368839324646, 0.514921891363332, -0.021377019336614, 0.801196801016305, -0.207446991247740 }),
159 new ArrayRealVector(new double[] { -0.713933751051495, 0.190582113553930, -0.671410443368332, 0.056056055955050, -0.006541576993581 }),
160 new ArrayRealVector(new double[] { -0.584677060845929, 0.367177264979103, 0.721453187784497, -0.052971054621812, 0.005740715188257 })
161 };
162
163 EigenDecompositionSymmetric decomposition;
164 decomposition = new EigenDecompositionSymmetric(mainTridiagonal, secondaryTridiagonal);
165
166 double[] eigenValues = decomposition.getEigenvalues();
167 for (int i = 0; i < refEigenValues.length; ++i) {
168 assertEquals(refEigenValues[i], eigenValues[i], 1.0e-5);
169 assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 2.0e-7);
170 }
171
172 assertEquals(mainTridiagonal.length, decomposition.getSolver().getRowDimension());
173 assertEquals(mainTridiagonal.length, decomposition.getSolver().getColumnDimension());
174 }
175
176 @Test
177 void testMathpbx02() {
178
179 double[] mainTridiagonal = {
180 7484.860960227216, 18405.28129035345, 13855.225609560746,
181 10016.708722343366, 559.8117399576674, 6750.190788301587,
182 71.21428769782159
183 };
184 double[] secondaryTridiagonal = {
185 -4175.088570476366,1975.7955858241994,5193.178422374075,
186 1995.286659169179,75.34535882933804,-234.0808002076056
187 };
188
189
190
191 double[] refEigenValues = {
192 20654.744890306974412,16828.208208485466457,
193 6893.155912634994820,6757.083016675340332,
194 5887.799885688558788,64.309089923240379,
195 57.992628792736340
196 };
197 RealVector[] refEigenVectors = {
198 new ArrayRealVector(new double[] {-0.270356342026904, 0.852811091326997, 0.399639490702077, 0.198794657813990, 0.019739323307666, 0.000106983022327, -0.000001216636321}),
199 new ArrayRealVector(new double[] {0.179995273578326,-0.402807848153042,0.701870993525734,0.555058211014888,0.068079148898236,0.000509139115227,-0.000007112235617}),
200 new ArrayRealVector(new double[] {-0.399582721284727,-0.056629954519333,-0.514406488522827,0.711168164518580,0.225548081276367,0.125943999652923,-0.004321507456014}),
201 new ArrayRealVector(new double[] {0.058515721572821,0.010200130057739,0.063516274916536,-0.090696087449378,-0.017148420432597,0.991318870265707,-0.034707338554096}),
202 new ArrayRealVector(new double[] {0.855205995537564,0.327134656629775,-0.265382397060548,0.282690729026706,0.105736068025572,-0.009138126622039,0.000367751821196}),
203 new ArrayRealVector(new double[] {-0.002913069901144,-0.005177515777101,0.041906334478672,-0.109315918416258,0.436192305456741,0.026307315639535,0.891797507436344}),
204 new ArrayRealVector(new double[] {-0.005738311176435,-0.010207611670378,0.082662420517928,-0.215733886094368,0.861606487840411,-0.025478530652759,-0.451080697503958})
205 };
206
207
208 EigenDecompositionSymmetric decomposition;
209 decomposition = new EigenDecompositionSymmetric(mainTridiagonal, secondaryTridiagonal);
210
211 double[] eigenValues = decomposition.getEigenvalues();
212 for (int i = 0; i < refEigenValues.length; ++i) {
213 assertEquals(refEigenValues[i], eigenValues[i], 1.0e-3);
214 if (refEigenVectors[i].dotProduct(decomposition.getEigenvector(i)) < 0) {
215 assertEquals(0, refEigenVectors[i].add(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
216 } else {
217 assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
218 }
219 }
220
221 }
222
223 @Test
224 void testMathpbx03() {
225
226 double[] mainTridiagonal = {
227 1809.0978259647177,3395.4763425956166,1832.1894584712693,3804.364873592377,
228 806.0482458637571,2403.656427234185,28.48691431556015
229 };
230 double[] secondaryTridiagonal = {
231 -656.8932064545833,-469.30804108920734,-1021.7714889369421,
232 -1152.540497328983,-939.9765163817368,-12.885877015422391
233 };
234
235
236
237 double[] refEigenValues = {
238 4603.121913685183245,3691.195818048970978,2743.442955402465032,1657.596442107321764,
239 1336.797819095331306,30.129865209677519,17.035352085224986
240 };
241
242 RealVector[] refEigenVectors = {
243 new ArrayRealVector(new double[] {-0.036249830202337,0.154184732411519,-0.346016328392363,0.867540105133093,-0.294483395433451,0.125854235969548,-0.000354507444044}),
244 new ArrayRealVector(new double[] {-0.318654191697157,0.912992309960507,-0.129270874079777,-0.184150038178035,0.096521712579439,-0.070468788536461,0.000247918177736}),
245 new ArrayRealVector(new double[] {-0.051394668681147,0.073102235876933,0.173502042943743,-0.188311980310942,-0.327158794289386,0.905206581432676,-0.004296342252659}),
246 new ArrayRealVector(new double[] {0.838150199198361,0.193305209055716,-0.457341242126146,-0.166933875895419,0.094512811358535,0.119062381338757,-0.000941755685226}),
247 new ArrayRealVector(new double[] {0.438071395458547,0.314969169786246,0.768480630802146,0.227919171600705,-0.193317045298647,-0.170305467485594,0.001677380536009}),
248 new ArrayRealVector(new double[] {-0.003726503878741,-0.010091946369146,-0.067152015137611,-0.113798146542187,-0.313123000097908,-0.118940107954918,0.932862311396062}),
249 new ArrayRealVector(new double[] {0.009373003194332,0.025570377559400,0.170955836081348,0.291954519805750,0.807824267665706,0.320108347088646,0.360202112392266}),
250 };
251
252
253 EigenDecompositionSymmetric decomposition;
254 decomposition = new EigenDecompositionSymmetric(mainTridiagonal, secondaryTridiagonal);
255
256 double[] eigenValues = decomposition.getEigenvalues();
257 for (int i = 0; i < refEigenValues.length; ++i) {
258 assertEquals(refEigenValues[i], eigenValues[i], 1.0e-4);
259 if (refEigenVectors[i].dotProduct(decomposition.getEigenvector(i)) < 0) {
260 assertEquals(0, refEigenVectors[i].add(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
261 } else {
262 assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
263 }
264 }
265
266 }
267
268
269 @Test
270 void testTridiagonal() {
271 Random r = new Random(4366663527842l);
272 double[] ref = new double[30];
273 for (int i = 0; i < ref.length; ++i) {
274 if (i < 5) {
275 ref[i] = 2 * r.nextDouble() - 1;
276 } else {
277 ref[i] = 0.0001 * r.nextDouble() + 6;
278 }
279 }
280 Arrays.sort(ref);
281 TriDiagonalTransformer t =
282 new TriDiagonalTransformer(createTestMatrix(r, ref));
283 EigenDecompositionSymmetric ed;
284 ed = new EigenDecompositionSymmetric(t.getMainDiagonalRef(), t.getSecondaryDiagonalRef());
285 double[] eigenValues = ed.getEigenvalues();
286 assertEquals(ref.length, eigenValues.length);
287 for (int i = 0; i < ref.length; ++i) {
288 assertEquals(ref[ref.length - i - 1], eigenValues[i], 2.0e-14);
289 }
290
291 }
292
293 @Test
294 void testGitHubIssue30() {
295 RealMatrix m = MatrixUtils.createRealMatrix(new double[][] {
296 { 23473.684554963584, 4273.093076392109 },
297 { 4273.093076392048, 4462.13956661408 }
298 });
299 EigenDecompositionSymmetric ed = new EigenDecompositionSymmetric(m, 1.0e-13, true);
300 assertEquals(0.0,
301 ed.getV().multiply(ed.getD()).multiply(ed.getVT()).subtract(m).getNorm1(),
302 1.0e-10);
303 }
304
305
306 @Test
307 void testDimensions() {
308 final int m = matrix.getRowDimension();
309 EigenDecompositionSymmetric ed;
310 ed = new EigenDecompositionSymmetric(matrix);
311 assertEquals(m, ed.getV().getRowDimension());
312 assertEquals(m, ed.getV().getColumnDimension());
313 assertEquals(m, ed.getD().getColumnDimension());
314 assertEquals(m, ed.getD().getColumnDimension());
315 assertEquals(m, ed.getVT().getRowDimension());
316 assertEquals(m, ed.getVT().getColumnDimension());
317 }
318
319
320 @Test
321 void testEigenvalues() {
322 EigenDecompositionSymmetric ed;
323 ed = new EigenDecompositionSymmetric(matrix);
324 double[] eigenValues = ed.getEigenvalues();
325 assertEquals(refValues.length, eigenValues.length);
326 for (int i = 0; i < refValues.length; ++i) {
327 assertEquals(refValues[i], eigenValues[i], 3.0e-15);
328 }
329 }
330
331 @Test
332 void testSymmetric() {
333 RealMatrix symmetric = MatrixUtils.createRealMatrix(new double[][] {
334 {4, 1, 1},
335 {1, 2, 3},
336 {1, 3, 6}
337 });
338
339 EigenDecompositionSymmetric ed;
340 ed = new EigenDecompositionSymmetric(symmetric);
341
342 DiagonalMatrix d = ed.getD();
343 RealMatrix v = ed.getV();
344 RealMatrix vT = ed.getVT();
345
346 double norm = v.multiply(d).multiply(vT).subtract(symmetric).getNorm1();
347 assertEquals(0, norm, 6.0e-13);
348 }
349
350 @Test
351 void testSquareRoot() {
352 final double[][] data = {
353 { 33, 24, 7 },
354 { 24, 57, 11 },
355 { 7, 11, 9 }
356 };
357
358 final EigenDecompositionSymmetric dec = new EigenDecompositionSymmetric(MatrixUtils.createRealMatrix(data));
359 final RealMatrix sqrtM = dec.getSquareRoot();
360
361
362 final RealMatrix m = sqrtM.multiply(sqrtM);
363
364 final int dim = data.length;
365 for (int r = 0; r < dim; r++) {
366 for (int c = 0; c < dim; c++) {
367 assertEquals(data[r][c], m.getEntry(r, c), 1e-13, "m[" + r + "][" + c + "]");
368 }
369 }
370 }
371
372 @Test
373 void testSquareRootNonSymmetric() {
374 assertThrows(MathRuntimeException.class, () -> {
375 final double[][] data = {
376 {1, 2, 4},
377 {2, 3, 5},
378 {11, 5, 9}
379 };
380
381 final EigenDecompositionSymmetric dec = new EigenDecompositionSymmetric(MatrixUtils.createRealMatrix(data));
382 @SuppressWarnings("unused")
383 final RealMatrix sqrtM = dec.getSquareRoot();
384 });
385 }
386
387 @Test
388 void testSquareRootNonPositiveDefinite() {
389 assertThrows(MathRuntimeException.class, () -> {
390 final double[][] data = {
391 {1, 2, 4},
392 {2, 3, 5},
393 {4, 5, -9}
394 };
395
396 final EigenDecompositionSymmetric dec = new EigenDecompositionSymmetric(MatrixUtils.createRealMatrix(data));
397 @SuppressWarnings("unused")
398 final RealMatrix sqrtM = dec.getSquareRoot();
399 });
400 }
401
402 @Test
403 void testIsNonSingularTinyOutOfOrderEigenvalue() {
404 try {
405 new EigenDecompositionSymmetric(MatrixUtils.createRealMatrix(new double[][] {
406 { 1e-13, 0 },
407 { 1, 1 },
408 }));
409 fail("an exception should have been thrown");
410 } catch (MathIllegalArgumentException miae) {
411 assertEquals(LocalizedCoreFormats.NON_SYMMETRIC_MATRIX, miae.getSpecifier());
412 }
413 }
414
415
416 @Test
417 void testEigenvectors() {
418 EigenDecompositionSymmetric ed;
419 ed = new EigenDecompositionSymmetric(matrix);
420 for (int i = 0; i < matrix.getRowDimension(); ++i) {
421 double lambda = ed.getEigenvalue(i);
422 RealVector v = ed.getEigenvector(i);
423 RealVector mV = matrix.operate(v);
424 assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13);
425 }
426 }
427
428
429 @Test
430 void testAEqualVDVt() {
431 EigenDecompositionSymmetric ed;
432 ed = new EigenDecompositionSymmetric(matrix);
433 RealMatrix v = ed.getV();
434 DiagonalMatrix d = ed.getD();
435 RealMatrix vT = ed.getVT();
436 double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm1();
437 assertEquals(0, norm, 6.0e-13);
438 }
439
440
441 @Test
442 void testVOrthogonal() {
443 RealMatrix v = new EigenDecompositionSymmetric(matrix).getV();
444 RealMatrix vTv = v.transposeMultiply(v);
445 RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
446 assertEquals(0, vTv.subtract(id).getNorm1(), 2.0e-13);
447 }
448
449
450 @Test
451 void testDiagonal() {
452 double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
453 RealMatrix m = MatrixUtils.createRealDiagonalMatrix(diagonal);
454 EigenDecompositionSymmetric ed;
455 ed = new EigenDecompositionSymmetric(m);
456 assertEquals(diagonal[0], ed.getEigenvalue(3), 2.0e-15);
457 assertEquals(diagonal[1], ed.getEigenvalue(2), 2.0e-15);
458 assertEquals(diagonal[2], ed.getEigenvalue(1), 2.0e-15);
459 assertEquals(diagonal[3], ed.getEigenvalue(0), 2.0e-15);
460 }
461
462
463
464
465 @Test
466 void testRepeatedEigenvalue() {
467 RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] {
468 {3, 2, 4},
469 {2, 0, 2},
470 {4, 2, 3}
471 });
472 EigenDecompositionSymmetric ed;
473 ed = new EigenDecompositionSymmetric(repeated);
474 checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
475 checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
476 }
477
478
479
480
481 @Test
482 void testDistinctEigenvalues() {
483 RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
484 {3, 1, -4},
485 {1, 3, -4},
486 {-4, -4, 8}
487 });
488 EigenDecompositionSymmetric ed;
489 ed = new EigenDecompositionSymmetric(distinct);
490 checkEigenValues((new double[] {2, 0, 12}), ed, 1E-12);
491 checkEigenVector((new double[] {1, -1, 0}), ed, 1E-12);
492 checkEigenVector((new double[] {1, 1, 1}), ed, 1E-12);
493 checkEigenVector((new double[] {-1, -1, 2}), ed, 1E-12);
494 }
495
496
497
498
499 @Test
500 void testZeroDivide() {
501 RealMatrix indefinite = MatrixUtils.createRealMatrix(new double [][] {
502 { 0.0, 1.0, -1.0 },
503 { 1.0, 1.0, 0.0 },
504 { -1.0,0.0, 1.0 }
505 });
506 EigenDecompositionSymmetric ed;
507 ed = new EigenDecompositionSymmetric(indefinite);
508 checkEigenValues((new double[] {2, 1, -1}), ed, 1E-12);
509 double isqrt3 = 1/FastMath.sqrt(3.0);
510 checkEigenVector((new double[] {isqrt3,isqrt3,-isqrt3}), ed, 1E-12);
511 double isqrt2 = 1/FastMath.sqrt(2.0);
512 checkEigenVector((new double[] {0.0,-isqrt2,-isqrt2}), ed, 1E-12);
513 double isqrt6 = 1/FastMath.sqrt(6.0);
514 checkEigenVector((new double[] {2*isqrt6,-isqrt6,isqrt6}), ed, 1E-12);
515 }
516
517
518
519
520
521 @Test
522 void testTinyValues() {
523 final double tiny = 1e-100;
524 RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
525 {3, 1, -4},
526 {1, 3, -4},
527 {-4, -4, 8}
528 });
529 distinct = distinct.scalarMultiply(tiny);
530
531 final EigenDecompositionSymmetric ed = new EigenDecompositionSymmetric(distinct);
532 checkEigenValues(MathArrays.scale(tiny, new double[] {2, 0, 12}), ed, 1e-12 * tiny);
533 checkEigenVector(new double[] {1, -1, 0}, ed, 1e-12);
534 checkEigenVector(new double[] {1, 1, 1}, ed, 1e-12);
535 checkEigenVector(new double[] {-1, -1, 2}, ed, 1e-12);
536 }
537
538
539
540
541 @Test
542 void testCustomEpsilon() {
543 RealMatrix matrix = MatrixUtils.createRealMatrix(new double[][] {
544 { 1.76208738E-13, -9.37625373E-13, -1.94760551E-12, -2.56572222E-11, -7.30093964E-11, -1.98340808E-09 },
545 { -9.37625373E-13, 5.00812620E-12, 1.06017205E-11, 1.40431472E-10, 3.62452521E-10, 1.05830167E-08 },
546 { -1.94760551E-12, 1.06017205E-11, 3.15658331E-11, 2.32155752E-09, -1.53067748E-09, 2.23110293E-08 },
547 { -2.56572222E-11, 1.40431472E-10, 2.32155752E-09, 8.81161492E-07, -8.70304198E-07, 2.93564832E-07 },
548 { -7.30093964E-11, 3.62452521E-10, -1.53067748E-09, -8.70304198E-07, 9.42413982E-07, 7.81029359E-07 },
549 { -1.98340808E-09, 1.05830167E-08, 2.23110293E-08, 2.93564832E-07, 7.81029359E-07, 2.23721205E-05 }
550 });
551
552 final EigenDecompositionSymmetric defaultEd = new EigenDecompositionSymmetric(matrix);
553 assertFalse(defaultEd.getSolver().isNonSingular());
554
555 final double customEpsilon = 1e-20;
556 final EigenDecompositionSymmetric customEd = new EigenDecompositionSymmetric(matrix, customEpsilon, false);
557 assertTrue(customEd.getSolver().isNonSingular());
558 assertEquals(customEpsilon, customEd.getEpsilon(), 1.0e-25);
559 }
560
561 @Test
562 void testSingular() {
563 final RealMatrix m = MatrixUtils.createRealMatrix(new double[][] { { 1, 0 }, { 0, 0 } });
564 DecompositionSolver solver = new EigenDecompositionSymmetric(m).getSolver();
565 try {
566 solver.solve(new ArrayRealVector(new double[] { 1.0, 1.0 }));
567 fail("an exception should have been thrown");
568 } catch (MathIllegalArgumentException miae) {
569 assertEquals(LocalizedCoreFormats.SINGULAR_MATRIX, miae.getSpecifier());
570 }
571 try {
572 solver.solve(new Array2DRowRealMatrix(new double[][] { { 1.0, 1.0 }, { 1.0, 2.0 } }));
573 fail("an exception should have been thrown");
574 } catch (MathIllegalArgumentException miae) {
575 assertEquals(LocalizedCoreFormats.SINGULAR_MATRIX, miae.getSpecifier());
576 }
577 }
578
579 @Test
580 void testIncreasingOrder() {
581
582 final RealMatrix m = MatrixUtils.createRealMatrix(new double[][] {
583 { 81.0, 63.0, 55.0, 49.0, 0.0, 0.0},
584 { 63.0, 82.0, 80.0, 69.0, 0.0, 0.0},
585 { 55.0, 80.0, 92.0, 75.0, 0.0, 0.0},
586 { 49.0, 69.0, 75.0, 73.0, 0.0, 0.0},
587 { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
588 { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}
589 });
590
591 EigenDecompositionSymmetric ed = new EigenDecompositionSymmetric(m,
592 EigenDecompositionSymmetric.DEFAULT_EPSILON,
593 false);
594
595 assertEquals( 0.0, ed.getD().getEntry(0, 0), 1.0e-14);
596 assertEquals( 0.0, ed.getD().getEntry(1, 1), 1.0e-14);
597 assertEquals( 4.793253233672134, ed.getD().getEntry(2, 2), 1.0e-14);
598 assertEquals( 7.429392849462756, ed.getD().getEntry(3, 3), 1.0e-14);
599 assertEquals( 36.43053556404571, ed.getD().getEntry(4, 4), 1.0e-14);
600 assertEquals(279.34681835281935, ed.getD().getEntry(5, 5), 1.0e-14);
601
602 }
603
604
605
606
607
608
609 protected void checkEigenValues(double[] targetValues,
610 EigenDecompositionSymmetric ed, double tolerance) {
611 double[] observed = ed.getEigenvalues();
612 for (int i = 0; i < observed.length; i++) {
613 assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
614 assertTrue(isIncludedValue(targetValues[i], observed, tolerance));
615 }
616 }
617
618
619
620
621
622
623 private boolean isIncludedValue(double value, double[] searchArray,
624 double tolerance) {
625 boolean found = false;
626 int i = 0;
627 while (!found && i < searchArray.length) {
628 if (FastMath.abs(value - searchArray[i]) < tolerance) {
629 found = true;
630 }
631 i++;
632 }
633 return found;
634 }
635
636
637
638
639
640
641 protected void checkEigenVector(double[] eigenVector,
642 EigenDecompositionSymmetric ed, double tolerance) {
643 assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
644 }
645
646
647
648
649
650 private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix,
651 double tolerance) {
652 boolean found = false;
653 int i = 0;
654 while (!found && i < searchMatrix.getColumnDimension()) {
655 double multiplier = 1.0;
656 boolean matching = true;
657 int j = 0;
658 while (matching && j < searchMatrix.getRowDimension()) {
659 double colEntry = searchMatrix.getEntry(j, i);
660
661 if (FastMath.abs(multiplier - 1.0) <= FastMath.ulp(1.0) && FastMath.abs(colEntry) > 1E-14
662 && FastMath.abs(column[j]) > 1e-14) {
663 multiplier = colEntry / column[j];
664 }
665 if (FastMath.abs(column[j] * multiplier - colEntry) > tolerance) {
666 matching = false;
667 }
668 j++;
669 }
670 found = matching;
671 i++;
672 }
673 return found;
674 }
675
676 @BeforeEach
677 void setUp() {
678 refValues = new double[] {
679 2.003, 2.002, 2.001, 1.001, 1.000, 0.001
680 };
681 matrix = createTestMatrix(new Random(35992629946426l), refValues);
682 }
683
684 @AfterEach
685 void tearDown() {
686 refValues = null;
687 matrix = null;
688 }
689
690 static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
691 final int n = eigenValues.length;
692 final RealMatrix v = createOrthogonalMatrix(r, n);
693 final RealMatrix d = MatrixUtils.createRealDiagonalMatrix(eigenValues);
694 return v.multiply(d).multiplyTransposed(v);
695 }
696
697 public static RealMatrix createOrthogonalMatrix(final Random r, final int size) {
698
699 final double[][] data = new double[size][size];
700
701 for (int i = 0; i < size; ++i) {
702 final double[] dataI = data[i];
703 double norm2 = 0;
704 do {
705
706
707 for (int j = 0; j < size; ++j) {
708 dataI[j] = 2 * r.nextDouble() - 1;
709 }
710
711
712 for (int k = 0; k < i; ++k) {
713 final double[] dataK = data[k];
714 double dotProduct = 0;
715 for (int j = 0; j < size; ++j) {
716 dotProduct += dataI[j] * dataK[j];
717 }
718 for (int j = 0; j < size; ++j) {
719 dataI[j] -= dotProduct * dataK[j];
720 }
721 }
722
723
724 norm2 = 0;
725 for (final double dataIJ : dataI) {
726 norm2 += dataIJ * dataIJ;
727 }
728 final double inv = 1.0 / FastMath.sqrt(norm2);
729 for (int j = 0; j < size; ++j) {
730 dataI[j] *= inv;
731 }
732
733 } while (norm2 * size < 0.01);
734 }
735
736 return MatrixUtils.createRealMatrix(data);
737
738 }
739 }