001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.linear;
019    
020    import java.util.Arrays;
021    import java.util.Random;
022    
023    import org.apache.commons.math.linear.EigenDecomposition;
024    import org.apache.commons.math.linear.EigenDecompositionImpl;
025    import org.apache.commons.math.linear.MatrixUtils;
026    import org.apache.commons.math.linear.RealMatrix;
027    import org.apache.commons.math.linear.RealVector;
028    import org.apache.commons.math.linear.TriDiagonalTransformer;
029    import org.apache.commons.math.util.MathUtils;
030    
031    import junit.framework.Test;
032    import junit.framework.TestCase;
033    import junit.framework.TestSuite;
034    
035    public class EigenDecompositionImplTest extends TestCase {
036    
037        private double[] refValues;
038        private RealMatrix matrix;
039    
040        public EigenDecompositionImplTest(String name) {
041            super(name);
042        }
043    
044        public static Test suite() {
045            TestSuite suite = new TestSuite(EigenDecompositionImplTest.class);
046            suite.setName("EigenDecompositionImpl Tests");
047            return suite;
048        }
049    
050        public void testDimension1() {
051            RealMatrix matrix =
052                MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
053            EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
054            assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
055        }
056    
057        public void testDimension2() {
058            RealMatrix matrix =
059                MatrixUtils.createRealMatrix(new double[][] {
060                        { 59.0, 12.0 },
061                        { 12.0, 66.0 }
062                });
063            EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
064            assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
065            assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
066        }
067    
068        public void testDimension3() {
069            RealMatrix matrix =
070                MatrixUtils.createRealMatrix(new double[][] {
071                                       {  39632.0, -4824.0, -16560.0 },
072                                       {  -4824.0,  8693.0,   7920.0 },
073                                       { -16560.0,  7920.0,  17300.0 }
074                                   });
075            EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
076            assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
077            assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
078            assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
079        }
080    
081        public void testDimension4WithSplit() {
082            RealMatrix matrix =
083                MatrixUtils.createRealMatrix(new double[][] {
084                                       {  0.784, -0.288,  0.000,  0.000 },
085                                       { -0.288,  0.616,  0.000,  0.000 },
086                                       {  0.000,  0.000,  0.164, -0.048 },
087                                       {  0.000,  0.000, -0.048,  0.136 }
088                                   });
089            EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
090            assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
091            assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
092            assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
093            assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
094        }
095    
096        public void testDimension4WithoutSplit() {
097            RealMatrix matrix =
098                MatrixUtils.createRealMatrix(new double[][] {
099                                       {  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    }