1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.linear;
19
20 import java.util.Arrays;
21 import java.util.Random;
22
23 import org.apache.commons.statistics.distribution.ContinuousDistribution;
24 import org.apache.commons.statistics.distribution.NormalDistribution;
25 import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
26 import org.apache.commons.math4.core.jdkmath.JdkMath;
27 import org.apache.commons.math4.legacy.core.MathArrays;
28 import org.apache.commons.numbers.core.Precision;
29 import org.apache.commons.rng.simple.RandomSource;
30 import org.junit.After;
31 import org.junit.Assert;
32 import org.junit.Before;
33 import org.junit.Ignore;
34 import org.junit.Test;
35
36 public class EigenDecompositionTest {
37
38 private double[] refValues;
39 private RealMatrix matrix;
40
41 @Test
42 public void testDimension1() {
43 RealMatrix matrix =
44 MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
45 EigenDecomposition ed;
46 ed = new EigenDecomposition(matrix);
47 Assert.assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
48 }
49
50 @Test
51 public void testDimension2() {
52 RealMatrix matrix =
53 MatrixUtils.createRealMatrix(new double[][] {
54 { 59.0, 12.0 },
55 { 12.0, 66.0 }
56 });
57 EigenDecomposition ed;
58 ed = new EigenDecomposition(matrix);
59 Assert.assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
60 Assert.assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
61 }
62
63 @Test
64 public void testDimension3() {
65 RealMatrix matrix =
66 MatrixUtils.createRealMatrix(new double[][] {
67 { 39632.0, -4824.0, -16560.0 },
68 { -4824.0, 8693.0, 7920.0 },
69 { -16560.0, 7920.0, 17300.0 }
70 });
71 EigenDecomposition ed;
72 ed = new EigenDecomposition(matrix);
73 Assert.assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
74 Assert.assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
75 Assert.assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
76 }
77
78 @Test
79 public void testDimension3MultipleRoot() {
80 RealMatrix matrix =
81 MatrixUtils.createRealMatrix(new double[][] {
82 { 5, 10, 15 },
83 { 10, 20, 30 },
84 { 15, 30, 45 }
85 });
86 EigenDecomposition ed;
87 ed = new EigenDecomposition(matrix);
88 Assert.assertEquals(70.0, ed.getRealEigenvalue(0), 3.0e-11);
89 Assert.assertEquals(0.0, ed.getRealEigenvalue(1), 3.0e-11);
90 Assert.assertEquals(0.0, ed.getRealEigenvalue(2), 3.0e-11);
91 }
92
93 @Test
94 public void testDimension4WithSplit() {
95 RealMatrix matrix =
96 MatrixUtils.createRealMatrix(new double[][] {
97 { 0.784, -0.288, 0.000, 0.000 },
98 { -0.288, 0.616, 0.000, 0.000 },
99 { 0.000, 0.000, 0.164, -0.048 },
100 { 0.000, 0.000, -0.048, 0.136 }
101 });
102 EigenDecomposition ed;
103 ed = new EigenDecomposition(matrix);
104 Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
105 Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
106 Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
107 Assert.assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
108 }
109
110 @Test
111 public void testDimension4WithoutSplit() {
112 RealMatrix matrix =
113 MatrixUtils.createRealMatrix(new double[][] {
114 { 0.5608, -0.2016, 0.1152, -0.2976 },
115 { -0.2016, 0.4432, -0.2304, 0.1152 },
116 { 0.1152, -0.2304, 0.3088, -0.1344 },
117 { -0.2976, 0.1152, -0.1344, 0.3872 }
118 });
119 EigenDecomposition ed;
120 ed = new EigenDecomposition(matrix);
121 Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
122 Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
123 Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
124 Assert.assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
125 }
126
127
128 @Test
129 public void testMath308() {
130
131 double[] mainTridiagonal = {
132 22.330154644539597, 46.65485522478641, 17.393672330044705, 54.46687435351116, 80.17800767709437
133 };
134 double[] secondaryTridiagonal = {
135 13.04450406501361, -5.977590941539671, 2.9040909856707517, 7.1570352792841225
136 };
137
138
139
140 double[] refEigenValues = {
141 82.044413207204002, 53.456697699894512, 52.536278520113882, 18.847969733754262, 14.138204224043099
142 };
143 RealVector[] refEigenVectors = {
144 new ArrayRealVector(new double[] { -0.000462690386766, -0.002118073109055, 0.011530080757413, 0.252322434584915, 0.967572088232592 }),
145 new ArrayRealVector(new double[] { 0.314647769490148, 0.750806415553905, -0.167700312025760, -0.537092972407375, 0.143854968127780 }),
146 new ArrayRealVector(new double[] { 0.222368839324646, 0.514921891363332, -0.021377019336614, 0.801196801016305, -0.207446991247740 }),
147 new ArrayRealVector(new double[] { -0.713933751051495, 0.190582113553930, -0.671410443368332, 0.056056055955050, -0.006541576993581 }),
148 new ArrayRealVector(new double[] { -0.584677060845929, 0.367177264979103, 0.721453187784497, -0.052971054621812, 0.005740715188257 })
149 };
150
151 EigenDecomposition decomposition;
152 decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal);
153
154 double[] eigenValues = decomposition.getRealEigenvalues();
155 for (int i = 0; i < refEigenValues.length; ++i) {
156 Assert.assertEquals(refEigenValues[i], eigenValues[i], 1.0e-5);
157 Assert.assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 2.0e-7);
158 }
159 }
160
161 @Test
162 public void testMathpbx02() {
163
164 double[] mainTridiagonal = {
165 7484.860960227216, 18405.28129035345, 13855.225609560746,
166 10016.708722343366, 559.8117399576674, 6750.190788301587,
167 71.21428769782159
168 };
169 double[] secondaryTridiagonal = {
170 -4175.088570476366,1975.7955858241994,5193.178422374075,
171 1995.286659169179,75.34535882933804,-234.0808002076056
172 };
173
174
175
176 double[] refEigenValues = {
177 20654.744890306974412,16828.208208485466457,
178 6893.155912634994820,6757.083016675340332,
179 5887.799885688558788,64.309089923240379,
180 57.992628792736340
181 };
182 RealVector[] refEigenVectors = {
183 new ArrayRealVector(new double[] {-0.270356342026904, 0.852811091326997, 0.399639490702077, 0.198794657813990, 0.019739323307666, 0.000106983022327, -0.000001216636321}),
184 new ArrayRealVector(new double[] {0.179995273578326,-0.402807848153042,0.701870993525734,0.555058211014888,0.068079148898236,0.000509139115227,-0.000007112235617}),
185 new ArrayRealVector(new double[] {-0.399582721284727,-0.056629954519333,-0.514406488522827,0.711168164518580,0.225548081276367,0.125943999652923,-0.004321507456014}),
186 new ArrayRealVector(new double[] {0.058515721572821,0.010200130057739,0.063516274916536,-0.090696087449378,-0.017148420432597,0.991318870265707,-0.034707338554096}),
187 new ArrayRealVector(new double[] {0.855205995537564,0.327134656629775,-0.265382397060548,0.282690729026706,0.105736068025572,-0.009138126622039,0.000367751821196}),
188 new ArrayRealVector(new double[] {-0.002913069901144,-0.005177515777101,0.041906334478672,-0.109315918416258,0.436192305456741,0.026307315639535,0.891797507436344}),
189 new ArrayRealVector(new double[] {-0.005738311176435,-0.010207611670378,0.082662420517928,-0.215733886094368,0.861606487840411,-0.025478530652759,-0.451080697503958})
190 };
191
192
193 EigenDecomposition decomposition;
194 decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal);
195
196 double[] eigenValues = decomposition.getRealEigenvalues();
197 for (int i = 0; i < refEigenValues.length; ++i) {
198 Assert.assertEquals(refEigenValues[i], eigenValues[i], 1.0e-3);
199 if (refEigenVectors[i].dotProduct(decomposition.getEigenvector(i)) < 0) {
200 Assert.assertEquals(0, refEigenVectors[i].add(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
201 } else {
202 Assert.assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
203 }
204 }
205 }
206
207 @Test
208 public void testMathpbx03() {
209
210 double[] mainTridiagonal = {
211 1809.0978259647177,3395.4763425956166,1832.1894584712693,3804.364873592377,
212 806.0482458637571,2403.656427234185,28.48691431556015
213 };
214 double[] secondaryTridiagonal = {
215 -656.8932064545833,-469.30804108920734,-1021.7714889369421,
216 -1152.540497328983,-939.9765163817368,-12.885877015422391
217 };
218
219
220
221 double[] refEigenValues = {
222 4603.121913685183245,3691.195818048970978,2743.442955402465032,1657.596442107321764,
223 1336.797819095331306,30.129865209677519,17.035352085224986
224 };
225
226 RealVector[] refEigenVectors = {
227 new ArrayRealVector(new double[] {-0.036249830202337,0.154184732411519,-0.346016328392363,0.867540105133093,-0.294483395433451,0.125854235969548,-0.000354507444044}),
228 new ArrayRealVector(new double[] {-0.318654191697157,0.912992309960507,-0.129270874079777,-0.184150038178035,0.096521712579439,-0.070468788536461,0.000247918177736}),
229 new ArrayRealVector(new double[] {-0.051394668681147,0.073102235876933,0.173502042943743,-0.188311980310942,-0.327158794289386,0.905206581432676,-0.004296342252659}),
230 new ArrayRealVector(new double[] {0.838150199198361,0.193305209055716,-0.457341242126146,-0.166933875895419,0.094512811358535,0.119062381338757,-0.000941755685226}),
231 new ArrayRealVector(new double[] {0.438071395458547,0.314969169786246,0.768480630802146,0.227919171600705,-0.193317045298647,-0.170305467485594,0.001677380536009}),
232 new ArrayRealVector(new double[] {-0.003726503878741,-0.010091946369146,-0.067152015137611,-0.113798146542187,-0.313123000097908,-0.118940107954918,0.932862311396062}),
233 new ArrayRealVector(new double[] {0.009373003194332,0.025570377559400,0.170955836081348,0.291954519805750,0.807824267665706,0.320108347088646,0.360202112392266}),
234 };
235
236
237 EigenDecomposition decomposition;
238 decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal);
239
240 double[] eigenValues = decomposition.getRealEigenvalues();
241 for (int i = 0; i < refEigenValues.length; ++i) {
242 Assert.assertEquals(refEigenValues[i], eigenValues[i], 1.0e-4);
243 if (refEigenVectors[i].dotProduct(decomposition.getEigenvector(i)) < 0) {
244 Assert.assertEquals(0, refEigenVectors[i].add(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
245 } else {
246 Assert.assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 1.0e-5);
247 }
248 }
249 }
250
251
252 @Test
253 public void testTridiagonal() {
254 Random r = new Random(4366663527842L);
255 double[] ref = new double[30];
256 for (int i = 0; i < ref.length; ++i) {
257 if (i < 5) {
258 ref[i] = 2 * r.nextDouble() - 1;
259 } else {
260 ref[i] = 0.0001 * r.nextDouble() + 6;
261 }
262 }
263 Arrays.sort(ref);
264 TriDiagonalTransformer t =
265 new TriDiagonalTransformer(createTestMatrix(r, ref));
266 EigenDecomposition ed;
267 ed = new EigenDecomposition(t.getMainDiagonalRef(), t.getSecondaryDiagonalRef());
268 double[] eigenValues = ed.getRealEigenvalues();
269 Assert.assertEquals(ref.length, eigenValues.length);
270 for (int i = 0; i < ref.length; ++i) {
271 Assert.assertEquals(ref[ref.length - i - 1], eigenValues[i], 2.0e-14);
272 }
273 }
274
275
276 @Test
277 public void testDimensions() {
278 final int m = matrix.getRowDimension();
279 EigenDecomposition ed;
280 ed = new EigenDecomposition(matrix);
281 Assert.assertEquals(m, ed.getV().getRowDimension());
282 Assert.assertEquals(m, ed.getV().getColumnDimension());
283 Assert.assertEquals(m, ed.getD().getColumnDimension());
284 Assert.assertEquals(m, ed.getD().getColumnDimension());
285 Assert.assertEquals(m, ed.getVT().getRowDimension());
286 Assert.assertEquals(m, ed.getVT().getColumnDimension());
287 }
288
289
290 @Test
291 public void testEigenvalues() {
292 EigenDecomposition ed;
293 ed = new EigenDecomposition(matrix);
294 double[] eigenValues = ed.getRealEigenvalues();
295 Assert.assertEquals(refValues.length, eigenValues.length);
296 for (int i = 0; i < refValues.length; ++i) {
297 Assert.assertEquals(refValues[i], eigenValues[i], 3.0e-15);
298 }
299 }
300
301
302 @Test
303 public void testBigMatrix() {
304 Random r = new Random(17748333525117L);
305 double[] bigValues = new double[200];
306 for (int i = 0; i < bigValues.length; ++i) {
307 bigValues[i] = 2 * r.nextDouble() - 1;
308 }
309 Arrays.sort(bigValues);
310 EigenDecomposition ed;
311 ed = new EigenDecomposition(createTestMatrix(r, bigValues));
312 double[] eigenValues = ed.getRealEigenvalues();
313 Assert.assertEquals(bigValues.length, eigenValues.length);
314 for (int i = 0; i < bigValues.length; ++i) {
315 Assert.assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14);
316 }
317 }
318
319 @Test
320 public void testSymmetric() {
321 RealMatrix symmetric = MatrixUtils.createRealMatrix(new double[][] {
322 {4, 1, 1},
323 {1, 2, 3},
324 {1, 3, 6}
325 });
326
327 EigenDecomposition ed;
328 ed = new EigenDecomposition(symmetric);
329
330 RealMatrix d = ed.getD();
331 RealMatrix v = ed.getV();
332 RealMatrix vT = ed.getVT();
333
334 double norm = v.multiply(d).multiply(vT).subtract(symmetric).getNorm();
335 Assert.assertEquals(0, norm, 6.0e-13);
336 }
337
338 @Test
339 public void testSquareRoot() {
340 final double[][] data = {
341 { 33, 24, 7 },
342 { 24, 57, 11 },
343 { 7, 11, 9 }
344 };
345
346 final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data));
347 final RealMatrix sqrtM = dec.getSquareRoot();
348
349
350 final RealMatrix m = sqrtM.multiply(sqrtM);
351
352 final int dim = data.length;
353 for (int r = 0; r < dim; r++) {
354 for (int c = 0; c < dim; c++) {
355 Assert.assertEquals("m[" + r + "][" + c + "]",
356 data[r][c], m.getEntry(r, c), 1e-13);
357 }
358 }
359 }
360
361 @Test(expected=MathUnsupportedOperationException.class)
362 public void testSquareRootNonSymmetric() {
363 final double[][] data = {
364 { 1, 2, 4 },
365 { 2, 3, 5 },
366 { 11, 5, 9 }
367 };
368
369 final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data));
370 @SuppressWarnings("unused")
371 final RealMatrix sqrtM = dec.getSquareRoot();
372 }
373
374 @Test(expected=MathUnsupportedOperationException.class)
375 public void testSquareRootNonPositiveDefinite() {
376 final double[][] data = {
377 { 1, 2, 4 },
378 { 2, 3, 5 },
379 { 4, 5, -9 }
380 };
381
382 final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data));
383 @SuppressWarnings("unused")
384 final RealMatrix sqrtM = dec.getSquareRoot();
385 }
386
387 @Test
388 public void testUnsymmetric() {
389
390 double[][] vData = { { -1.0, 1.0, -1.0, 1.0 },
391 { -8.0, 4.0, -2.0, 1.0 },
392 { 27.0, 9.0, 3.0, 1.0 },
393 { 64.0, 16.0, 4.0, 1.0 } };
394 checkUnsymmetricMatrix(MatrixUtils.createRealMatrix(vData));
395
396 RealMatrix randMatrix = MatrixUtils.createRealMatrix(new double[][] {
397 {0, 1, 0, 0},
398 {1, 0, 2.e-7, 0},
399 {0, -2.e-7, 0, 1},
400 {0, 0, 1, 0}
401 });
402 checkUnsymmetricMatrix(randMatrix);
403
404
405 double[][] randData2 = {
406 { 0.680, -0.3300, -0.2700, -0.717, -0.687, 0.0259 },
407 { -0.211, 0.5360, 0.0268, 0.214, -0.198, 0.6780 },
408 { 0.566, -0.4440, 0.9040, -0.967, -0.740, 0.2250 },
409 { 0.597, 0.1080, 0.8320, -0.514, -0.782, -0.4080 },
410 { 0.823, -0.0452, 0.2710, -0.726, 0.998, 0.2750 },
411 { -0.605, 0.2580, 0.4350, 0.608, -0.563, 0.0486 }
412 };
413 checkUnsymmetricMatrix(MatrixUtils.createRealMatrix(randData2));
414 }
415
416 @Test
417 @Ignore
418 public void testRandomUnsymmetricMatrix() {
419 for (int run = 0; run < 100; run++) {
420 Random r = new Random(System.currentTimeMillis());
421
422
423 int size = r.nextInt(20) + 4;
424
425 double[][] data = new double[size][size];
426 for (int i = 0; i < size; i++) {
427 for (int j = 0; j < size; j++) {
428 data[i][j] = r.nextInt(100);
429 }
430 }
431
432 RealMatrix m = MatrixUtils.createRealMatrix(data);
433 checkUnsymmetricMatrix(m);
434 }
435 }
436
437
438
439
440
441
442
443 @Test
444 public void testMath1051() {
445 double[][] data = {
446 {0,0,0,0,0},
447 {0,0,0,0,1},
448 {0,0,0,1,0},
449 {1,1,0,0,1},
450 {1,0,1,0,1}
451 };
452
453 RealMatrix m = MatrixUtils.createRealMatrix(data);
454 checkUnsymmetricMatrix(m);
455 }
456
457 @Test
458 @Ignore
459 public void testNormalDistributionUnsymmetricMatrix() {
460 for (int run = 0; run < 100; run++) {
461 Random r = new Random(System.currentTimeMillis());
462 ContinuousDistribution.Sampler dist
463 = NormalDistribution.of(0.0, r.nextDouble() * 5).createSampler(RandomSource.WELL_512_A.create(64925784252L));
464
465
466 int size = r.nextInt(20) + 4;
467
468 double[][] data = new double[size][size];
469 for (int i = 0; i < size; i++) {
470 for (int j = 0; j < size; j++) {
471 data[i][j] = dist.sample();
472 }
473 }
474
475 RealMatrix m = MatrixUtils.createRealMatrix(data);
476 checkUnsymmetricMatrix(m);
477 }
478 }
479
480 @Test
481 public void testMath848() {
482 double[][] data = {
483 { 0.1849449280, -0.0646971046, 0.0774755812, -0.0969651755, -0.0692648806, 0.3282344352, -0.0177423074, 0.2063136340},
484 {-0.0742700134, -0.0289063030, -0.0017269460, -0.0375550146, -0.0487737922, -0.2616837868, -0.0821201295, -0.2530000167},
485 { 0.2549910127, 0.0995733692, -0.0009718388, 0.0149282808, 0.1791878897, -0.0823182816, 0.0582629256, 0.3219545182},
486 {-0.0694747557, -0.1880649148, -0.2740630911, 0.0720096468, -0.1800836914, -0.3518996425, 0.2486747833, 0.6257938167},
487 { 0.0536360918, -0.1339297778, 0.2241579764, -0.0195327484, -0.0054103808, 0.0347564518, 0.5120802482, -0.0329902864},
488 {-0.5933332356, -0.2488721082, 0.2357173629, 0.0177285473, 0.0856630593, -0.3567126300, -0.1600668126, -0.1010899621},
489 {-0.0514349819, -0.0854319435, 0.1125050061, 0.0063453560, -0.2250000688, -0.2209343090, 0.1964623477, -0.1512329924},
490 { 0.0197395947, -0.1997170581, -0.1425959019, -0.2749477910, -0.0969467073, 0.0603688520, -0.2826905192, 0.1794315473}};
491 RealMatrix m = MatrixUtils.createRealMatrix(data);
492 checkUnsymmetricMatrix(m);
493 }
494
495
496
497
498
499 private void checkUnsymmetricMatrix(final RealMatrix m) {
500 try {
501 EigenDecomposition ed = new EigenDecomposition(m);
502
503 RealMatrix d = ed.getD();
504 RealMatrix v = ed.getV();
505
506
507 RealMatrix x = m.multiply(v);
508 RealMatrix y = v.multiply(d);
509
510 double diffNorm = x.subtract(y).getNorm();
511 Assert.assertTrue("The norm of (X-Y) is too large: " + diffNorm + ", matrix=" + m.toString(),
512 x.subtract(y).getNorm() < 1000 * Precision.EPSILON * JdkMath.max(x.getNorm(), y.getNorm()));
513
514 RealMatrix invV = new LUDecomposition(v).getSolver().getInverse();
515 double norm = v.multiply(d).multiply(invV).subtract(m).getNorm();
516 Assert.assertEquals(0.0, norm, 1.0e-10);
517 } catch (Exception e) {
518 Assert.fail("Failed to create EigenDecomposition for matrix " + m.toString() + ", ex=" + e.toString());
519 }
520 }
521
522
523 @Test
524 public void testEigenvectors() {
525 EigenDecomposition ed;
526 ed = new EigenDecomposition(matrix);
527 for (int i = 0; i < matrix.getRowDimension(); ++i) {
528 double lambda = ed.getRealEigenvalue(i);
529 RealVector v = ed.getEigenvector(i);
530 RealVector mV = matrix.operate(v);
531 Assert.assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13);
532 }
533 }
534
535
536 @Test
537 public void testAEqualVDVt() {
538 EigenDecomposition ed;
539 ed = new EigenDecomposition(matrix);
540 RealMatrix v = ed.getV();
541 RealMatrix d = ed.getD();
542 RealMatrix vT = ed.getVT();
543 double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm();
544 Assert.assertEquals(0, norm, 6.0e-13);
545 }
546
547
548 @Test
549 public void testVOrthogonal() {
550 RealMatrix v = new EigenDecomposition(matrix).getV();
551 RealMatrix vTv = v.transpose().multiply(v);
552 RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
553 Assert.assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
554 }
555
556
557 @Test
558 public void testDiagonal() {
559 double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
560 RealMatrix m = MatrixUtils.createRealDiagonalMatrix(diagonal);
561 EigenDecomposition ed;
562 ed = new EigenDecomposition(m);
563 Assert.assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
564 Assert.assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
565 Assert.assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
566 Assert.assertEquals(diagonal[3], ed.getRealEigenvalue(0), 2.0e-15);
567 }
568
569
570
571
572 @Test
573 public void testRepeatedEigenvalue() {
574 RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] {
575 {3, 2, 4},
576 {2, 0, 2},
577 {4, 2, 3}
578 });
579 EigenDecomposition ed;
580 ed = new EigenDecomposition(repeated);
581 checkEigenValues(new double[] {8, -1, -1}, ed, 1E-12);
582 checkEigenVector(new double[] {2, 1, 2}, ed, 1E-12);
583 }
584
585
586
587
588 @Test
589 public void testDistinctEigenvalues() {
590 RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
591 {3, 1, -4},
592 {1, 3, -4},
593 {-4, -4, 8}
594 });
595 EigenDecomposition ed;
596 ed = new EigenDecomposition(distinct);
597 checkEigenValues(new double[] {2, 0, 12}, ed, 1E-12);
598 checkEigenVector(new double[] {1, -1, 0}, ed, 1E-12);
599 checkEigenVector(new double[] {1, 1, 1}, ed, 1E-12);
600 checkEigenVector(new double[] {-1, -1, 2}, ed, 1E-12);
601 }
602
603
604
605
606 @Test
607 public void testZeroDivide() {
608 RealMatrix indefinite = MatrixUtils.createRealMatrix(new double [][] {
609 { 0.0, 1.0, -1.0 },
610 { 1.0, 1.0, 0.0 },
611 { -1.0,0.0, 1.0 }
612 });
613 EigenDecomposition ed;
614 ed = new EigenDecomposition(indefinite);
615 checkEigenValues(new double[] {2, 1, -1}, ed, 1E-12);
616 double isqrt3 = 1/JdkMath.sqrt(3.0);
617 checkEigenVector(new double[] {isqrt3,isqrt3,-isqrt3}, ed, 1E-12);
618 double isqrt2 = 1/JdkMath.sqrt(2.0);
619 checkEigenVector(new double[] {0.0,-isqrt2,-isqrt2}, ed, 1E-12);
620 double isqrt6 = 1/JdkMath.sqrt(6.0);
621 checkEigenVector(new double[] {2*isqrt6,-isqrt6,isqrt6}, ed, 1E-12);
622 }
623
624
625
626
627
628 @Test
629 public void testTinyValues() {
630 final double tiny = 1e-100;
631 RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
632 {3, 1, -4},
633 {1, 3, -4},
634 {-4, -4, 8}
635 });
636 distinct = distinct.scalarMultiply(tiny);
637
638 final EigenDecomposition ed = new EigenDecomposition(distinct);
639 checkEigenValues(MathArrays.scale(tiny, new double[] {2, 0, 12}), ed, 1e-12 * tiny);
640 checkEigenVector(new double[] {1, -1, 0}, ed, 1e-12);
641 checkEigenVector(new double[] {1, 1, 1}, ed, 1e-12);
642 checkEigenVector(new double[] {-1, -1, 2}, ed, 1e-12);
643 }
644
645
646
647
648
649
650 protected void checkEigenValues(double[] targetValues,
651 EigenDecomposition ed, double tolerance) {
652 double[] observed = ed.getRealEigenvalues();
653 for (int i = 0; i < observed.length; i++) {
654 Assert.assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
655 Assert.assertTrue(isIncludedValue(targetValues[i], observed, tolerance));
656 }
657 }
658
659
660
661
662
663
664 private boolean isIncludedValue(double value, double[] searchArray,
665 double tolerance) {
666 boolean found = false;
667 int i = 0;
668 while (!found && i < searchArray.length) {
669 if (JdkMath.abs(value - searchArray[i]) < tolerance) {
670 found = true;
671 }
672 i++;
673 }
674 return found;
675 }
676
677
678
679
680
681
682 protected void checkEigenVector(double[] eigenVector,
683 EigenDecomposition ed, double tolerance) {
684 Assert.assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
685 }
686
687
688
689
690
691 private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix,
692 double tolerance) {
693 boolean found = false;
694 int i = 0;
695 while (!found && i < searchMatrix.getColumnDimension()) {
696 double multiplier = 1.0;
697 boolean matching = true;
698 int j = 0;
699 while (matching && j < searchMatrix.getRowDimension()) {
700 double colEntry = searchMatrix.getEntry(j, i);
701
702 if (JdkMath.abs(multiplier - 1.0) <= JdkMath.ulp(1.0) && JdkMath.abs(colEntry) > 1E-14
703 && JdkMath.abs(column[j]) > 1e-14) {
704 multiplier = colEntry / column[j];
705 }
706 if (JdkMath.abs(column[j] * multiplier - colEntry) > tolerance) {
707 matching = false;
708 }
709 j++;
710 }
711 found = matching;
712 i++;
713 }
714 return found;
715 }
716
717 @Before
718 public void setUp() {
719 refValues = new double[] {
720 2.003, 2.002, 2.001, 1.001, 1.000, 0.001
721 };
722 matrix = createTestMatrix(new Random(35992629946426L), refValues);
723 }
724
725 @After
726 public void tearDown() {
727 refValues = null;
728 matrix = null;
729 }
730
731 static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
732 final int n = eigenValues.length;
733 final RealMatrix v = createOrthogonalMatrix(r, n);
734 final RealMatrix d = MatrixUtils.createRealDiagonalMatrix(eigenValues);
735 return v.multiply(d).multiply(v.transpose());
736 }
737
738 public static RealMatrix createOrthogonalMatrix(final Random r, final int size) {
739
740 final double[][] data = new double[size][size];
741
742 for (int i = 0; i < size; ++i) {
743 final double[] dataI = data[i];
744 double norm2 = 0;
745 do {
746
747
748 for (int j = 0; j < size; ++j) {
749 dataI[j] = 2 * r.nextDouble() - 1;
750 }
751
752
753 for (int k = 0; k < i; ++k) {
754 final double[] dataK = data[k];
755 double dotProduct = 0;
756 for (int j = 0; j < size; ++j) {
757 dotProduct += dataI[j] * dataK[j];
758 }
759 for (int j = 0; j < size; ++j) {
760 dataI[j] -= dotProduct * dataK[j];
761 }
762 }
763
764
765 norm2 = 0;
766 for (final double dataIJ : dataI) {
767 norm2 += dataIJ * dataIJ;
768 }
769 final double inv = 1.0 / JdkMath.sqrt(norm2);
770 for (int j = 0; j < size; ++j) {
771 dataI[j] *= inv;
772 }
773 } while (norm2 * size < 0.01);
774 }
775
776 return MatrixUtils.createRealMatrix(data);
777 }
778 }