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 }