1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.linear;
19
20 import java.util.Arrays;
21 import java.util.Random;
22
23 import org.apache.commons.math.linear.EigenDecomposition;
24 import org.apache.commons.math.linear.EigenDecompositionImpl;
25 import org.apache.commons.math.linear.MatrixUtils;
26 import org.apache.commons.math.linear.RealMatrix;
27 import org.apache.commons.math.linear.RealVector;
28 import org.apache.commons.math.linear.TriDiagonalTransformer;
29 import org.apache.commons.math.util.MathUtils;
30
31 import junit.framework.Test;
32 import junit.framework.TestCase;
33 import junit.framework.TestSuite;
34
35 public class EigenDecompositionImplTest extends TestCase {
36
37 private double[] refValues;
38 private RealMatrix matrix;
39
40 public EigenDecompositionImplTest(String name) {
41 super(name);
42 }
43
44 public static Test suite() {
45 TestSuite suite = new TestSuite(EigenDecompositionImplTest.class);
46 suite.setName("EigenDecompositionImpl Tests");
47 return suite;
48 }
49
50 public void testDimension1() {
51 RealMatrix matrix =
52 MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
53 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
54 assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
55 }
56
57 public void testDimension2() {
58 RealMatrix matrix =
59 MatrixUtils.createRealMatrix(new double[][] {
60 { 59.0, 12.0 },
61 { 12.0, 66.0 }
62 });
63 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
64 assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
65 assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
66 }
67
68 public void testDimension3() {
69 RealMatrix matrix =
70 MatrixUtils.createRealMatrix(new double[][] {
71 { 39632.0, -4824.0, -16560.0 },
72 { -4824.0, 8693.0, 7920.0 },
73 { -16560.0, 7920.0, 17300.0 }
74 });
75 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
76 assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
77 assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
78 assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
79 }
80
81 public void testDimension4WithSplit() {
82 RealMatrix matrix =
83 MatrixUtils.createRealMatrix(new double[][] {
84 { 0.784, -0.288, 0.000, 0.000 },
85 { -0.288, 0.616, 0.000, 0.000 },
86 { 0.000, 0.000, 0.164, -0.048 },
87 { 0.000, 0.000, -0.048, 0.136 }
88 });
89 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
90 assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
91 assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
92 assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
93 assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
94 }
95
96 public void testDimension4WithoutSplit() {
97 RealMatrix matrix =
98 MatrixUtils.createRealMatrix(new double[][] {
99 { 0.5608, -0.2016, 0.1152, -0.2976 },
100 { -0.2016, 0.4432, -0.2304, 0.1152 },
101 { 0.1152, -0.2304, 0.3088, -0.1344 },
102 { -0.2976, 0.1152, -0.1344, 0.3872 }
103 });
104 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
105 assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
106 assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
107 assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
108 assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
109 }
110
111
112 public void testTridiagonal() {
113 Random r = new Random(4366663527842l);
114 double[] ref = new double[30];
115 for (int i = 0; i < ref.length; ++i) {
116 if (i < 5) {
117 ref[i] = 2 * r.nextDouble() - 1;
118 } else {
119 ref[i] = 0.0001 * r.nextDouble() + 6;
120 }
121 }
122 Arrays.sort(ref);
123 TriDiagonalTransformer t =
124 new TriDiagonalTransformer(createTestMatrix(r, ref));
125 EigenDecomposition ed =
126 new EigenDecompositionImpl(t.getMainDiagonalRef(),
127 t.getSecondaryDiagonalRef(),
128 MathUtils.SAFE_MIN);
129 double[] eigenValues = ed.getRealEigenvalues();
130 assertEquals(ref.length, eigenValues.length);
131 for (int i = 0; i < ref.length; ++i) {
132 assertEquals(ref[ref.length - i - 1], eigenValues[i], 2.0e-14);
133 }
134
135 }
136
137
138 public void testDimensions() {
139 final int m = matrix.getRowDimension();
140 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
141 assertEquals(m, ed.getV().getRowDimension());
142 assertEquals(m, ed.getV().getColumnDimension());
143 assertEquals(m, ed.getD().getColumnDimension());
144 assertEquals(m, ed.getD().getColumnDimension());
145 assertEquals(m, ed.getVT().getRowDimension());
146 assertEquals(m, ed.getVT().getColumnDimension());
147 }
148
149
150 public void testEigenvalues() {
151 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
152 double[] eigenValues = ed.getRealEigenvalues();
153 assertEquals(refValues.length, eigenValues.length);
154 for (int i = 0; i < refValues.length; ++i) {
155 assertEquals(refValues[i], eigenValues[i], 3.0e-15);
156 }
157 }
158
159
160 public void testBigMatrix() {
161 Random r = new Random(17748333525117l);
162 double[] bigValues = new double[200];
163 for (int i = 0; i < bigValues.length; ++i) {
164 bigValues[i] = 2 * r.nextDouble() - 1;
165 }
166 Arrays.sort(bigValues);
167 EigenDecomposition ed =
168 new EigenDecompositionImpl(createTestMatrix(r, bigValues), MathUtils.SAFE_MIN);
169 double[] eigenValues = ed.getRealEigenvalues();
170 assertEquals(bigValues.length, eigenValues.length);
171 for (int i = 0; i < bigValues.length; ++i) {
172 assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14);
173 }
174 }
175
176
177 public void testEigenvectors() {
178 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
179 for (int i = 0; i < matrix.getRowDimension(); ++i) {
180 double lambda = ed.getRealEigenvalue(i);
181 RealVector v = ed.getEigenvector(i);
182 RealVector mV = matrix.operate(v);
183 assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13);
184 }
185 }
186
187
188 public void testAEqualVDVt() {
189 EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
190 RealMatrix v = ed.getV();
191 RealMatrix d = ed.getD();
192 RealMatrix vT = ed.getVT();
193 double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm();
194 assertEquals(0, norm, 6.0e-13);
195 }
196
197
198 public void testVOrthogonal() {
199 RealMatrix v = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN).getV();
200 RealMatrix vTv = v.transpose().multiply(v);
201 RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
202 assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
203 }
204
205
206 public void testDiagonal() {
207 double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
208 RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
209 EigenDecomposition ed = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN);
210 assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
211 assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
212 assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
213 assertEquals(diagonal[3], ed.getRealEigenvalue(0), 2.0e-15);
214 }
215
216
217
218
219 public void testRepeatedEigenvalue() {
220 RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] {
221 {3, 2, 4},
222 {2, 0, 2},
223 {4, 2, 3}
224 });
225 EigenDecomposition ed = new EigenDecompositionImpl(repeated, MathUtils.SAFE_MIN);
226 checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
227 checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
228 }
229
230
231
232
233 public void testDistinctEigenvalues() {
234 RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
235 {3, 1, -4},
236 {1, 3, -4},
237 {-4, -4, 8}
238 });
239 EigenDecomposition ed = new EigenDecompositionImpl(distinct, MathUtils.SAFE_MIN);
240 checkEigenValues((new double[] {2, 0, 12}), ed, 1E-12);
241 checkEigenVector((new double[] {1, -1, 0}), ed, 1E-12);
242 checkEigenVector((new double[] {1, 1, 1}), ed, 1E-12);
243 checkEigenVector((new double[] {-1, -1, 2}), ed, 1E-12);
244 }
245
246
247
248
249
250
251 protected void checkEigenValues(double[] targetValues,
252 EigenDecomposition ed, double tolerance) {
253 double[] observed = ed.getRealEigenvalues();
254 for (int i = 0; i < observed.length; i++) {
255 assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
256 assertTrue(isIncludedValue(targetValues[i], observed, tolerance));
257 }
258 }
259
260
261
262
263
264 private boolean isIncludedValue(double value, double[] searchArray,
265 double tolerance) {
266 boolean found = false;
267 int i = 0;
268 while (!found && i < searchArray.length) {
269 if (Math.abs(value - searchArray[i]) < tolerance) {
270 found = true;
271 }
272 i++;
273 }
274 return found;
275 }
276
277
278
279
280
281
282 protected void checkEigenVector(double[] eigenVector,
283 EigenDecomposition ed, double tolerance) {
284 assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
285 }
286
287
288
289
290
291 private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix,
292 double tolerance) {
293 boolean found = false;
294 int i = 0;
295 while (!found && i < searchMatrix.getColumnDimension()) {
296 double multiplier = 1.0;
297 boolean matching = true;
298 int j = 0;
299 while (matching && j < searchMatrix.getRowDimension()) {
300 double colEntry = searchMatrix.getEntry(j, i);
301
302 if (Math.abs(multiplier - 1.0) <= Math.ulp(1.0) && Math.abs(colEntry) > 1E-14
303 && Math.abs(column[j]) > 1e-14) {
304 multiplier = colEntry / column[j];
305 }
306 if (Math.abs(column[j] * multiplier - colEntry) > tolerance) {
307 matching = false;
308 }
309 j++;
310 }
311 found = matching;
312 i++;
313 }
314 return found;
315 }
316
317 @Override
318 public void setUp() {
319 refValues = new double[] {
320 2.003, 2.002, 2.001, 1.001, 1.000, 0.001
321 };
322 matrix = createTestMatrix(new Random(35992629946426l), refValues);
323 }
324
325 @Override
326 public void tearDown() {
327 refValues = null;
328 matrix = null;
329 }
330
331 static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
332 final int n = eigenValues.length;
333 final RealMatrix v = createOrthogonalMatrix(r, n);
334 final RealMatrix d = createDiagonalMatrix(eigenValues, n, n);
335 return v.multiply(d).multiply(v.transpose());
336 }
337
338 public static RealMatrix createOrthogonalMatrix(final Random r, final int size) {
339
340 final double[][] data = new double[size][size];
341
342 for (int i = 0; i < size; ++i) {
343 final double[] dataI = data[i];
344 double norm2 = 0;
345 do {
346
347
348 for (int j = 0; j < size; ++j) {
349 dataI[j] = 2 * r.nextDouble() - 1;
350 }
351
352
353 for (int k = 0; k < i; ++k) {
354 final double[] dataK = data[k];
355 double dotProduct = 0;
356 for (int j = 0; j < size; ++j) {
357 dotProduct += dataI[j] * dataK[j];
358 }
359 for (int j = 0; j < size; ++j) {
360 dataI[j] -= dotProduct * dataK[j];
361 }
362 }
363
364
365 norm2 = 0;
366 for (final double dataIJ : dataI) {
367 norm2 += dataIJ * dataIJ;
368 }
369 final double inv = 1.0 / Math.sqrt(norm2);
370 for (int j = 0; j < size; ++j) {
371 dataI[j] *= inv;
372 }
373
374 } while (norm2 * size < 0.01);
375 }
376
377 return MatrixUtils.createRealMatrix(data);
378
379 }
380
381 public static RealMatrix createDiagonalMatrix(final double[] diagonal,
382 final int rows, final int columns) {
383 final double[][] dData = new double[rows][columns];
384 for (int i = 0; i < Math.min(rows, columns); ++i) {
385 dData[i][i] = diagonal[i];
386 }
387 return MatrixUtils.createRealMatrix(dData);
388 }
389
390 }