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