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 }