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  
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     /** test a matrix already in tridiagonal form. */
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     /** test dimensions */
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     /** test eigenvalues */
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     /** test eigenvalues for a big matrix. */
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     /** test eigenvectors */
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     /** test A = VDVt */
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     /** test that V is orthogonal */
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     /** test diagonal matrix */
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      * Matrix with eigenvalues {8, -1, -1}
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      * Matrix with eigenvalues {2, 0, 12}
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      * Verifies that the given EigenDecomposition has eigenvalues equivalent to
248      * the targetValues, ignoring the order of the values and allowing
249      * values to differ by tolerance.
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      * Returns true iff there is an entry within tolerance of value in
262      * searchArray.
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      * Returns true iff eigenVector is a scalar multiple of one of the columns
279      * of ed.getV().  Does not try linear combinations - i.e., should only be
280      * used to find vectors in one-dimensional eigenspaces.
281      */
282     protected void checkEigenVector(double[] eigenVector,
283             EigenDecomposition ed, double tolerance) {
284         assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
285     }
286     
287     /**
288      * Returns true iff there is a column that is a scalar multiple of column
289      * in searchMatrix (modulo tolerance)
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                 // Use the first entry where both are non-zero as scalar
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                 // generate randomly row I
348                 for (int j = 0; j < size; ++j) {
349                     dataI[j] = 2 * r.nextDouble() - 1;
350                 }
351 
352                 // project the row in the subspace orthogonal to previous rows
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                 // normalize the row
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 }