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.math3.linear; 19 20 import java.util.Arrays; 21 import java.util.Random; 22 23 24 import org.apache.commons.math3.distribution.NormalDistribution; 25 import org.apache.commons.math3.util.FastMath; 26 import org.apache.commons.math3.util.Precision; 27 import org.apache.commons.math3.util.MathArrays; 28 import org.apache.commons.math3.exception.MathUnsupportedOperationException; 29 import org.junit.After; 30 import org.junit.Assert; 31 import org.junit.Before; 32 import org.junit.Ignore; 33 import org.junit.Test; 34 35 public class EigenDecompositionTest { 36 37 private double[] refValues; 38 private RealMatrix matrix; 39 40 @Test testDimension1()41 public void testDimension1() { 42 RealMatrix matrix = 43 MatrixUtils.createRealMatrix(new double[][] { { 1.5 } }); 44 EigenDecomposition ed; 45 ed = new EigenDecomposition(matrix); 46 Assert.assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15); 47 } 48 49 @Test testDimension2()50 public void testDimension2() { 51 RealMatrix matrix = 52 MatrixUtils.createRealMatrix(new double[][] { 53 { 59.0, 12.0 }, 54 { 12.0, 66.0 } 55 }); 56 EigenDecomposition ed; 57 ed = new EigenDecomposition(matrix); 58 Assert.assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15); 59 Assert.assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15); 60 } 61 62 @Test testDimension3()63 public void testDimension3() { 64 RealMatrix matrix = 65 MatrixUtils.createRealMatrix(new double[][] { 66 { 39632.0, -4824.0, -16560.0 }, 67 { -4824.0, 8693.0, 7920.0 }, 68 { -16560.0, 7920.0, 17300.0 } 69 }); 70 EigenDecomposition ed; 71 ed = new EigenDecomposition(matrix); 72 Assert.assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11); 73 Assert.assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11); 74 Assert.assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11); 75 } 76 77 @Test testDimension3MultipleRoot()78 public void testDimension3MultipleRoot() { 79 RealMatrix matrix = 80 MatrixUtils.createRealMatrix(new double[][] { 81 { 5, 10, 15 }, 82 { 10, 20, 30 }, 83 { 15, 30, 45 } 84 }); 85 EigenDecomposition ed; 86 ed = new EigenDecomposition(matrix); 87 Assert.assertEquals(70.0, ed.getRealEigenvalue(0), 3.0e-11); 88 Assert.assertEquals(0.0, ed.getRealEigenvalue(1), 3.0e-11); 89 Assert.assertEquals(0.0, ed.getRealEigenvalue(2), 3.0e-11); 90 } 91 92 @Test testDimension4WithSplit()93 public void testDimension4WithSplit() { 94 RealMatrix matrix = 95 MatrixUtils.createRealMatrix(new double[][] { 96 { 0.784, -0.288, 0.000, 0.000 }, 97 { -0.288, 0.616, 0.000, 0.000 }, 98 { 0.000, 0.000, 0.164, -0.048 }, 99 { 0.000, 0.000, -0.048, 0.136 } 100 }); 101 EigenDecomposition ed; 102 ed = new EigenDecomposition(matrix); 103 Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15); 104 Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15); 105 Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15); 106 Assert.assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15); 107 } 108 109 @Test testDimension4WithoutSplit()110 public void testDimension4WithoutSplit() { 111 RealMatrix matrix = 112 MatrixUtils.createRealMatrix(new double[][] { 113 { 0.5608, -0.2016, 0.1152, -0.2976 }, 114 { -0.2016, 0.4432, -0.2304, 0.1152 }, 115 { 0.1152, -0.2304, 0.3088, -0.1344 }, 116 { -0.2976, 0.1152, -0.1344, 0.3872 } 117 }); 118 EigenDecomposition ed; 119 ed = new EigenDecomposition(matrix); 120 Assert.assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15); 121 Assert.assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15); 122 Assert.assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15); 123 Assert.assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15); 124 } 125 126 // the following test triggered an ArrayIndexOutOfBoundsException in commons-math 2.0 127 @Test testMath308()128 public void testMath308() { 129 130 double[] mainTridiagonal = { 131 22.330154644539597, 46.65485522478641, 17.393672330044705, 54.46687435351116, 80.17800767709437 132 }; 133 double[] secondaryTridiagonal = { 134 13.04450406501361, -5.977590941539671, 2.9040909856707517, 7.1570352792841225 135 }; 136 137 // the reference values have been computed using routine DSTEMR 138 // from the fortran library LAPACK version 3.2.1 139 double[] refEigenValues = { 140 82.044413207204002, 53.456697699894512, 52.536278520113882, 18.847969733754262, 14.138204224043099 141 }; 142 RealVector[] refEigenVectors = { 143 new ArrayRealVector(new double[] { -0.000462690386766, -0.002118073109055, 0.011530080757413, 0.252322434584915, 0.967572088232592 }), 144 new ArrayRealVector(new double[] { 0.314647769490148, 0.750806415553905, -0.167700312025760, -0.537092972407375, 0.143854968127780 }), 145 new ArrayRealVector(new double[] { 0.222368839324646, 0.514921891363332, -0.021377019336614, 0.801196801016305, -0.207446991247740 }), 146 new ArrayRealVector(new double[] { -0.713933751051495, 0.190582113553930, -0.671410443368332, 0.056056055955050, -0.006541576993581 }), 147 new ArrayRealVector(new double[] { -0.584677060845929, 0.367177264979103, 0.721453187784497, -0.052971054621812, 0.005740715188257 }) 148 }; 149 150 EigenDecomposition decomposition; 151 decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal); 152 153 double[] eigenValues = decomposition.getRealEigenvalues(); 154 for (int i = 0; i < refEigenValues.length; ++i) { 155 Assert.assertEquals(refEigenValues[i], eigenValues[i], 1.0e-5); 156 Assert.assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 2.0e-7); 157 } 158 159 } 160 161 @Test testMathpbx02()162 public void testMathpbx02() { 163 164 double[] mainTridiagonal = { 165 7484.860960227216, 18405.28129035345, 13855.225609560746, 166 10016.708722343366, 559.8117399576674, 6750.190788301587, 167 71.21428769782159 168 }; 169 double[] secondaryTridiagonal = { 170 -4175.088570476366,1975.7955858241994,5193.178422374075, 171 1995.286659169179,75.34535882933804,-234.0808002076056 172 }; 173 174 // the reference values have been computed using routine DSTEMR 175 // from the fortran library LAPACK version 3.2.1 176 double[] refEigenValues = { 177 20654.744890306974412,16828.208208485466457, 178 6893.155912634994820,6757.083016675340332, 179 5887.799885688558788,64.309089923240379, 180 57.992628792736340 181 }; 182 RealVector[] refEigenVectors = { 183 new ArrayRealVector(new double[] {-0.270356342026904, 0.852811091326997, 0.399639490702077, 0.198794657813990, 0.019739323307666, 0.000106983022327, -0.000001216636321}), 184 new ArrayRealVector(new double[] {0.179995273578326,-0.402807848153042,0.701870993525734,0.555058211014888,0.068079148898236,0.000509139115227,-0.000007112235617}), 185 new ArrayRealVector(new double[] {-0.399582721284727,-0.056629954519333,-0.514406488522827,0.711168164518580,0.225548081276367,0.125943999652923,-0.004321507456014}), 186 new ArrayRealVector(new double[] {0.058515721572821,0.010200130057739,0.063516274916536,-0.090696087449378,-0.017148420432597,0.991318870265707,-0.034707338554096}), 187 new ArrayRealVector(new double[] {0.855205995537564,0.327134656629775,-0.265382397060548,0.282690729026706,0.105736068025572,-0.009138126622039,0.000367751821196}), 188 new ArrayRealVector(new double[] {-0.002913069901144,-0.005177515777101,0.041906334478672,-0.109315918416258,0.436192305456741,0.026307315639535,0.891797507436344}), 189 new ArrayRealVector(new double[] {-0.005738311176435,-0.010207611670378,0.082662420517928,-0.215733886094368,0.861606487840411,-0.025478530652759,-0.451080697503958}) 190 }; 191 192 // the following line triggers the exception 193 EigenDecomposition decomposition; 194 decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal); 195 196 double[] eigenValues = decomposition.getRealEigenvalues(); 197 for (int i = 0; i < refEigenValues.length; ++i) { 198 Assert.assertEquals(refEigenValues[i], eigenValues[i], 1.0e-3); 199 if (refEigenVectors[i].dotProduct(decomposition.getEigenvector(i)) < 0) { 200 Assert.assertEquals(0, refEigenVectors[i].add(decomposition.getEigenvector(i)).getNorm(), 1.0e-5); 201 } else { 202 Assert.assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 1.0e-5); 203 } 204 } 205 206 } 207 208 @Test testMathpbx03()209 public void testMathpbx03() { 210 211 double[] mainTridiagonal = { 212 1809.0978259647177,3395.4763425956166,1832.1894584712693,3804.364873592377, 213 806.0482458637571,2403.656427234185,28.48691431556015 214 }; 215 double[] secondaryTridiagonal = { 216 -656.8932064545833,-469.30804108920734,-1021.7714889369421, 217 -1152.540497328983,-939.9765163817368,-12.885877015422391 218 }; 219 220 // the reference values have been computed using routine DSTEMR 221 // from the fortran library LAPACK version 3.2.1 222 double[] refEigenValues = { 223 4603.121913685183245,3691.195818048970978,2743.442955402465032,1657.596442107321764, 224 1336.797819095331306,30.129865209677519,17.035352085224986 225 }; 226 227 RealVector[] refEigenVectors = { 228 new ArrayRealVector(new double[] {-0.036249830202337,0.154184732411519,-0.346016328392363,0.867540105133093,-0.294483395433451,0.125854235969548,-0.000354507444044}), 229 new ArrayRealVector(new double[] {-0.318654191697157,0.912992309960507,-0.129270874079777,-0.184150038178035,0.096521712579439,-0.070468788536461,0.000247918177736}), 230 new ArrayRealVector(new double[] {-0.051394668681147,0.073102235876933,0.173502042943743,-0.188311980310942,-0.327158794289386,0.905206581432676,-0.004296342252659}), 231 new ArrayRealVector(new double[] {0.838150199198361,0.193305209055716,-0.457341242126146,-0.166933875895419,0.094512811358535,0.119062381338757,-0.000941755685226}), 232 new ArrayRealVector(new double[] {0.438071395458547,0.314969169786246,0.768480630802146,0.227919171600705,-0.193317045298647,-0.170305467485594,0.001677380536009}), 233 new ArrayRealVector(new double[] {-0.003726503878741,-0.010091946369146,-0.067152015137611,-0.113798146542187,-0.313123000097908,-0.118940107954918,0.932862311396062}), 234 new ArrayRealVector(new double[] {0.009373003194332,0.025570377559400,0.170955836081348,0.291954519805750,0.807824267665706,0.320108347088646,0.360202112392266}), 235 }; 236 237 // the following line triggers the exception 238 EigenDecomposition decomposition; 239 decomposition = new EigenDecomposition(mainTridiagonal, secondaryTridiagonal); 240 241 double[] eigenValues = decomposition.getRealEigenvalues(); 242 for (int i = 0; i < refEigenValues.length; ++i) { 243 Assert.assertEquals(refEigenValues[i], eigenValues[i], 1.0e-4); 244 if (refEigenVectors[i].dotProduct(decomposition.getEigenvector(i)) < 0) { 245 Assert.assertEquals(0, refEigenVectors[i].add(decomposition.getEigenvector(i)).getNorm(), 1.0e-5); 246 } else { 247 Assert.assertEquals(0, refEigenVectors[i].subtract(decomposition.getEigenvector(i)).getNorm(), 1.0e-5); 248 } 249 } 250 251 } 252 253 /** test a matrix already in tridiagonal form. */ 254 @Test testTridiagonal()255 public void testTridiagonal() { 256 Random r = new Random(4366663527842l); 257 double[] ref = new double[30]; 258 for (int i = 0; i < ref.length; ++i) { 259 if (i < 5) { 260 ref[i] = 2 * r.nextDouble() - 1; 261 } else { 262 ref[i] = 0.0001 * r.nextDouble() + 6; 263 } 264 } 265 Arrays.sort(ref); 266 TriDiagonalTransformer t = 267 new TriDiagonalTransformer(createTestMatrix(r, ref)); 268 EigenDecomposition ed; 269 ed = new EigenDecomposition(t.getMainDiagonalRef(), t.getSecondaryDiagonalRef()); 270 double[] eigenValues = ed.getRealEigenvalues(); 271 Assert.assertEquals(ref.length, eigenValues.length); 272 for (int i = 0; i < ref.length; ++i) { 273 Assert.assertEquals(ref[ref.length - i - 1], eigenValues[i], 2.0e-14); 274 } 275 276 } 277 278 /** test dimensions */ 279 @Test testDimensions()280 public void testDimensions() { 281 final int m = matrix.getRowDimension(); 282 EigenDecomposition ed; 283 ed = new EigenDecomposition(matrix); 284 Assert.assertEquals(m, ed.getV().getRowDimension()); 285 Assert.assertEquals(m, ed.getV().getColumnDimension()); 286 Assert.assertEquals(m, ed.getD().getColumnDimension()); 287 Assert.assertEquals(m, ed.getD().getColumnDimension()); 288 Assert.assertEquals(m, ed.getVT().getRowDimension()); 289 Assert.assertEquals(m, ed.getVT().getColumnDimension()); 290 } 291 292 /** test eigenvalues */ 293 @Test testEigenvalues()294 public void testEigenvalues() { 295 EigenDecomposition ed; 296 ed = new EigenDecomposition(matrix); 297 double[] eigenValues = ed.getRealEigenvalues(); 298 Assert.assertEquals(refValues.length, eigenValues.length); 299 for (int i = 0; i < refValues.length; ++i) { 300 Assert.assertEquals(refValues[i], eigenValues[i], 3.0e-15); 301 } 302 } 303 304 /** test eigenvalues for a big matrix. */ 305 @Test testBigMatrix()306 public void testBigMatrix() { 307 Random r = new Random(17748333525117l); 308 double[] bigValues = new double[200]; 309 for (int i = 0; i < bigValues.length; ++i) { 310 bigValues[i] = 2 * r.nextDouble() - 1; 311 } 312 Arrays.sort(bigValues); 313 EigenDecomposition ed; 314 ed = new EigenDecomposition(createTestMatrix(r, bigValues)); 315 double[] eigenValues = ed.getRealEigenvalues(); 316 Assert.assertEquals(bigValues.length, eigenValues.length); 317 for (int i = 0; i < bigValues.length; ++i) { 318 Assert.assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14); 319 } 320 } 321 322 @Test testSymmetric()323 public void testSymmetric() { 324 RealMatrix symmetric = MatrixUtils.createRealMatrix(new double[][] { 325 {4, 1, 1}, 326 {1, 2, 3}, 327 {1, 3, 6} 328 }); 329 330 EigenDecomposition ed; 331 ed = new EigenDecomposition(symmetric); 332 333 RealMatrix d = ed.getD(); 334 RealMatrix v = ed.getV(); 335 RealMatrix vT = ed.getVT(); 336 337 double norm = v.multiply(d).multiply(vT).subtract(symmetric).getNorm(); 338 Assert.assertEquals(0, norm, 6.0e-13); 339 } 340 341 @Test testSquareRoot()342 public void testSquareRoot() { 343 final double[][] data = { 344 { 33, 24, 7 }, 345 { 24, 57, 11 }, 346 { 7, 11, 9 } 347 }; 348 349 final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data)); 350 final RealMatrix sqrtM = dec.getSquareRoot(); 351 352 // Reconstruct initial matrix. 353 final RealMatrix m = sqrtM.multiply(sqrtM); 354 355 final int dim = data.length; 356 for (int r = 0; r < dim; r++) { 357 for (int c = 0; c < dim; c++) { 358 Assert.assertEquals("m[" + r + "][" + c + "]", 359 data[r][c], m.getEntry(r, c), 1e-13); 360 } 361 } 362 } 363 364 @Test(expected=MathUnsupportedOperationException.class) testSquareRootNonSymmetric()365 public void testSquareRootNonSymmetric() { 366 final double[][] data = { 367 { 1, 2, 4 }, 368 { 2, 3, 5 }, 369 { 11, 5, 9 } 370 }; 371 372 final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data)); 373 @SuppressWarnings("unused") 374 final RealMatrix sqrtM = dec.getSquareRoot(); 375 } 376 377 @Test(expected=MathUnsupportedOperationException.class) testSquareRootNonPositiveDefinite()378 public void testSquareRootNonPositiveDefinite() { 379 final double[][] data = { 380 { 1, 2, 4 }, 381 { 2, 3, 5 }, 382 { 4, 5, -9 } 383 }; 384 385 final EigenDecomposition dec = new EigenDecomposition(MatrixUtils.createRealMatrix(data)); 386 @SuppressWarnings("unused") 387 final RealMatrix sqrtM = dec.getSquareRoot(); 388 } 389 390 @Test testUnsymmetric()391 public void testUnsymmetric() { 392 // Vandermonde matrix V(x;i,j) = x_i^{n - j} with x = (-1,-2,3,4) 393 double[][] vData = { { -1.0, 1.0, -1.0, 1.0 }, 394 { -8.0, 4.0, -2.0, 1.0 }, 395 { 27.0, 9.0, 3.0, 1.0 }, 396 { 64.0, 16.0, 4.0, 1.0 } }; 397 checkUnsymmetricMatrix(MatrixUtils.createRealMatrix(vData)); 398 399 RealMatrix randMatrix = MatrixUtils.createRealMatrix(new double[][] { 400 {0, 1, 0, 0}, 401 {1, 0, 2.e-7, 0}, 402 {0, -2.e-7, 0, 1}, 403 {0, 0, 1, 0} 404 }); 405 checkUnsymmetricMatrix(randMatrix); 406 407 // from http://eigen.tuxfamily.org/dox/classEigen_1_1RealSchur.html 408 double[][] randData2 = { 409 { 0.680, -0.3300, -0.2700, -0.717, -0.687, 0.0259 }, 410 { -0.211, 0.5360, 0.0268, 0.214, -0.198, 0.6780 }, 411 { 0.566, -0.4440, 0.9040, -0.967, -0.740, 0.2250 }, 412 { 0.597, 0.1080, 0.8320, -0.514, -0.782, -0.4080 }, 413 { 0.823, -0.0452, 0.2710, -0.726, 0.998, 0.2750 }, 414 { -0.605, 0.2580, 0.4350, 0.608, -0.563, 0.0486 } 415 }; 416 checkUnsymmetricMatrix(MatrixUtils.createRealMatrix(randData2)); 417 } 418 419 @Test 420 @Ignore testRandomUnsymmetricMatrix()421 public void testRandomUnsymmetricMatrix() { 422 for (int run = 0; run < 100; run++) { 423 Random r = new Random(System.currentTimeMillis()); 424 425 // matrix size 426 int size = r.nextInt(20) + 4; 427 428 double[][] data = new double[size][size]; 429 for (int i = 0; i < size; i++) { 430 for (int j = 0; j < size; j++) { 431 data[i][j] = r.nextInt(100); 432 } 433 } 434 435 RealMatrix m = MatrixUtils.createRealMatrix(data); 436 checkUnsymmetricMatrix(m); 437 } 438 } 439 440 /** 441 * Tests the porting of a bugfix in Jama-1.0.3 (from changelog): 442 * 443 * Patched hqr2 method in Jama.EigenvalueDecomposition to avoid infinite loop; 444 * Thanks Frederic Devernay <frederic.devernay@m4x.org> 445 */ 446 @Test testMath1051()447 public void testMath1051() { 448 double[][] data = { 449 {0,0,0,0,0}, 450 {0,0,0,0,1}, 451 {0,0,0,1,0}, 452 {1,1,0,0,1}, 453 {1,0,1,0,1} 454 }; 455 456 RealMatrix m = MatrixUtils.createRealMatrix(data); 457 checkUnsymmetricMatrix(m); 458 } 459 460 @Test 461 @Ignore testNormalDistributionUnsymmetricMatrix()462 public void testNormalDistributionUnsymmetricMatrix() { 463 for (int run = 0; run < 100; run++) { 464 Random r = new Random(System.currentTimeMillis()); 465 NormalDistribution dist = new NormalDistribution(0.0, r.nextDouble() * 5); 466 467 // matrix size 468 int size = r.nextInt(20) + 4; 469 470 double[][] data = new double[size][size]; 471 for (int i = 0; i < size; i++) { 472 for (int j = 0; j < size; j++) { 473 data[i][j] = dist.sample(); 474 } 475 } 476 477 RealMatrix m = MatrixUtils.createRealMatrix(data); 478 checkUnsymmetricMatrix(m); 479 } 480 } 481 482 @Test testMath848()483 public void testMath848() { 484 double[][] data = { 485 { 0.1849449280, -0.0646971046, 0.0774755812, -0.0969651755, -0.0692648806, 0.3282344352, -0.0177423074, 0.2063136340}, 486 {-0.0742700134, -0.0289063030, -0.0017269460, -0.0375550146, -0.0487737922, -0.2616837868, -0.0821201295, -0.2530000167}, 487 { 0.2549910127, 0.0995733692, -0.0009718388, 0.0149282808, 0.1791878897, -0.0823182816, 0.0582629256, 0.3219545182}, 488 {-0.0694747557, -0.1880649148, -0.2740630911, 0.0720096468, -0.1800836914, -0.3518996425, 0.2486747833, 0.6257938167}, 489 { 0.0536360918, -0.1339297778, 0.2241579764, -0.0195327484, -0.0054103808, 0.0347564518, 0.5120802482, -0.0329902864}, 490 {-0.5933332356, -0.2488721082, 0.2357173629, 0.0177285473, 0.0856630593, -0.3567126300, -0.1600668126, -0.1010899621}, 491 {-0.0514349819, -0.0854319435, 0.1125050061, 0.0063453560, -0.2250000688, -0.2209343090, 0.1964623477, -0.1512329924}, 492 { 0.0197395947, -0.1997170581, -0.1425959019, -0.2749477910, -0.0969467073, 0.0603688520, -0.2826905192, 0.1794315473}}; 493 RealMatrix m = MatrixUtils.createRealMatrix(data); 494 checkUnsymmetricMatrix(m); 495 } 496 497 /** 498 * Checks that the eigen decomposition of a general (unsymmetric) matrix is valid by 499 * checking: A*V = V*D 500 */ checkUnsymmetricMatrix(final RealMatrix m)501 private void checkUnsymmetricMatrix(final RealMatrix m) { 502 try { 503 EigenDecomposition ed = new EigenDecomposition(m); 504 505 RealMatrix d = ed.getD(); 506 RealMatrix v = ed.getV(); 507 //RealMatrix vT = ed.getVT(); 508 509 RealMatrix x = m.multiply(v); 510 RealMatrix y = v.multiply(d); 511 512 double diffNorm = x.subtract(y).getNorm(); 513 Assert.assertTrue("The norm of (X-Y) is too large: " + diffNorm + ", matrix=" + m.toString(), 514 x.subtract(y).getNorm() < 1000 * Precision.EPSILON * FastMath.max(x.getNorm(), y.getNorm())); 515 516 RealMatrix invV = new LUDecomposition(v).getSolver().getInverse(); 517 double norm = v.multiply(d).multiply(invV).subtract(m).getNorm(); 518 Assert.assertEquals(0.0, norm, 1.0e-10); 519 } catch (Exception e) { 520 Assert.fail("Failed to create EigenDecomposition for matrix " + m.toString() + ", ex=" + e.toString()); 521 } 522 } 523 524 /** test eigenvectors */ 525 @Test testEigenvectors()526 public void testEigenvectors() { 527 EigenDecomposition ed; 528 ed = new EigenDecomposition(matrix); 529 for (int i = 0; i < matrix.getRowDimension(); ++i) { 530 double lambda = ed.getRealEigenvalue(i); 531 RealVector v = ed.getEigenvector(i); 532 RealVector mV = matrix.operate(v); 533 Assert.assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13); 534 } 535 } 536 537 /** test A = VDVt */ 538 @Test testAEqualVDVt()539 public void testAEqualVDVt() { 540 EigenDecomposition ed; 541 ed = new EigenDecomposition(matrix); 542 RealMatrix v = ed.getV(); 543 RealMatrix d = ed.getD(); 544 RealMatrix vT = ed.getVT(); 545 double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm(); 546 Assert.assertEquals(0, norm, 6.0e-13); 547 } 548 549 /** test that V is orthogonal */ 550 @Test testVOrthogonal()551 public void testVOrthogonal() { 552 RealMatrix v = new EigenDecomposition(matrix).getV(); 553 RealMatrix vTv = v.transpose().multiply(v); 554 RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension()); 555 Assert.assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13); 556 } 557 558 /** test diagonal matrix */ 559 @Test testDiagonal()560 public void testDiagonal() { 561 double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 }; 562 RealMatrix m = MatrixUtils.createRealDiagonalMatrix(diagonal); 563 EigenDecomposition ed; 564 ed = new EigenDecomposition(m); 565 Assert.assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15); 566 Assert.assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15); 567 Assert.assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15); 568 Assert.assertEquals(diagonal[3], ed.getRealEigenvalue(0), 2.0e-15); 569 } 570 571 /** 572 * Matrix with eigenvalues {8, -1, -1} 573 */ 574 @Test testRepeatedEigenvalue()575 public void testRepeatedEigenvalue() { 576 RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] { 577 {3, 2, 4}, 578 {2, 0, 2}, 579 {4, 2, 3} 580 }); 581 EigenDecomposition ed; 582 ed = new EigenDecomposition(repeated); 583 checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12); 584 checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12); 585 } 586 587 /** 588 * Matrix with eigenvalues {2, 0, 12} 589 */ 590 @Test testDistinctEigenvalues()591 public void testDistinctEigenvalues() { 592 RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] { 593 {3, 1, -4}, 594 {1, 3, -4}, 595 {-4, -4, 8} 596 }); 597 EigenDecomposition ed; 598 ed = new EigenDecomposition(distinct); 599 checkEigenValues((new double[] {2, 0, 12}), ed, 1E-12); 600 checkEigenVector((new double[] {1, -1, 0}), ed, 1E-12); 601 checkEigenVector((new double[] {1, 1, 1}), ed, 1E-12); 602 checkEigenVector((new double[] {-1, -1, 2}), ed, 1E-12); 603 } 604 605 /** 606 * Verifies operation on indefinite matrix 607 */ 608 @Test testZeroDivide()609 public void testZeroDivide() { 610 RealMatrix indefinite = MatrixUtils.createRealMatrix(new double [][] { 611 { 0.0, 1.0, -1.0 }, 612 { 1.0, 1.0, 0.0 }, 613 { -1.0,0.0, 1.0 } 614 }); 615 EigenDecomposition ed; 616 ed = new EigenDecomposition(indefinite); 617 checkEigenValues((new double[] {2, 1, -1}), ed, 1E-12); 618 double isqrt3 = 1/FastMath.sqrt(3.0); 619 checkEigenVector((new double[] {isqrt3,isqrt3,-isqrt3}), ed, 1E-12); 620 double isqrt2 = 1/FastMath.sqrt(2.0); 621 checkEigenVector((new double[] {0.0,-isqrt2,-isqrt2}), ed, 1E-12); 622 double isqrt6 = 1/FastMath.sqrt(6.0); 623 checkEigenVector((new double[] {2*isqrt6,-isqrt6,isqrt6}), ed, 1E-12); 624 } 625 626 /** 627 * Verifies operation on very small values. 628 * Matrix with eigenvalues {2e-100, 0, 12e-100} 629 */ 630 @Test testTinyValues()631 public void testTinyValues() { 632 final double tiny = 1e-100; 633 RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] { 634 {3, 1, -4}, 635 {1, 3, -4}, 636 {-4, -4, 8} 637 }); 638 distinct = distinct.scalarMultiply(tiny); 639 640 final EigenDecomposition ed = new EigenDecomposition(distinct); 641 checkEigenValues(MathArrays.scale(tiny, new double[] {2, 0, 12}), ed, 1e-12 * tiny); 642 checkEigenVector(new double[] {1, -1, 0}, ed, 1e-12); 643 checkEigenVector(new double[] {1, 1, 1}, ed, 1e-12); 644 checkEigenVector(new double[] {-1, -1, 2}, ed, 1e-12); 645 } 646 647 /** 648 * Verifies that the given EigenDecomposition has eigenvalues equivalent to 649 * the targetValues, ignoring the order of the values and allowing 650 * values to differ by tolerance. 651 */ checkEigenValues(double[] targetValues, EigenDecomposition ed, double tolerance)652 protected void checkEigenValues(double[] targetValues, 653 EigenDecomposition ed, double tolerance) { 654 double[] observed = ed.getRealEigenvalues(); 655 for (int i = 0; i < observed.length; i++) { 656 Assert.assertTrue(isIncludedValue(observed[i], targetValues, tolerance)); 657 Assert.assertTrue(isIncludedValue(targetValues[i], observed, tolerance)); 658 } 659 } 660 661 662 /** 663 * Returns true iff there is an entry within tolerance of value in 664 * searchArray. 665 */ isIncludedValue(double value, double[] searchArray, double tolerance)666 private boolean isIncludedValue(double value, double[] searchArray, 667 double tolerance) { 668 boolean found = false; 669 int i = 0; 670 while (!found && i < searchArray.length) { 671 if (FastMath.abs(value - searchArray[i]) < tolerance) { 672 found = true; 673 } 674 i++; 675 } 676 return found; 677 } 678 679 /** 680 * Returns true iff eigenVector is a scalar multiple of one of the columns 681 * of ed.getV(). Does not try linear combinations - i.e., should only be 682 * used to find vectors in one-dimensional eigenspaces. 683 */ checkEigenVector(double[] eigenVector, EigenDecomposition ed, double tolerance)684 protected void checkEigenVector(double[] eigenVector, 685 EigenDecomposition ed, double tolerance) { 686 Assert.assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance)); 687 } 688 689 /** 690 * Returns true iff there is a column that is a scalar multiple of column 691 * in searchMatrix (modulo tolerance) 692 */ isIncludedColumn(double[] column, RealMatrix searchMatrix, double tolerance)693 private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix, 694 double tolerance) { 695 boolean found = false; 696 int i = 0; 697 while (!found && i < searchMatrix.getColumnDimension()) { 698 double multiplier = 1.0; 699 boolean matching = true; 700 int j = 0; 701 while (matching && j < searchMatrix.getRowDimension()) { 702 double colEntry = searchMatrix.getEntry(j, i); 703 // Use the first entry where both are non-zero as scalar 704 if (FastMath.abs(multiplier - 1.0) <= FastMath.ulp(1.0) && FastMath.abs(colEntry) > 1E-14 705 && FastMath.abs(column[j]) > 1e-14) { 706 multiplier = colEntry / column[j]; 707 } 708 if (FastMath.abs(column[j] * multiplier - colEntry) > tolerance) { 709 matching = false; 710 } 711 j++; 712 } 713 found = matching; 714 i++; 715 } 716 return found; 717 } 718 719 @Before setUp()720 public void setUp() { 721 refValues = new double[] { 722 2.003, 2.002, 2.001, 1.001, 1.000, 0.001 723 }; 724 matrix = createTestMatrix(new Random(35992629946426l), refValues); 725 } 726 727 @After tearDown()728 public void tearDown() { 729 refValues = null; 730 matrix = null; 731 } 732 createTestMatrix(final Random r, final double[] eigenValues)733 static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) { 734 final int n = eigenValues.length; 735 final RealMatrix v = createOrthogonalMatrix(r, n); 736 final RealMatrix d = MatrixUtils.createRealDiagonalMatrix(eigenValues); 737 return v.multiply(d).multiply(v.transpose()); 738 } 739 createOrthogonalMatrix(final Random r, final int size)740 public static RealMatrix createOrthogonalMatrix(final Random r, final int size) { 741 742 final double[][] data = new double[size][size]; 743 744 for (int i = 0; i < size; ++i) { 745 final double[] dataI = data[i]; 746 double norm2 = 0; 747 do { 748 749 // generate randomly row I 750 for (int j = 0; j < size; ++j) { 751 dataI[j] = 2 * r.nextDouble() - 1; 752 } 753 754 // project the row in the subspace orthogonal to previous rows 755 for (int k = 0; k < i; ++k) { 756 final double[] dataK = data[k]; 757 double dotProduct = 0; 758 for (int j = 0; j < size; ++j) { 759 dotProduct += dataI[j] * dataK[j]; 760 } 761 for (int j = 0; j < size; ++j) { 762 dataI[j] -= dotProduct * dataK[j]; 763 } 764 } 765 766 // normalize the row 767 norm2 = 0; 768 for (final double dataIJ : dataI) { 769 norm2 += dataIJ * dataIJ; 770 } 771 final double inv = 1.0 / FastMath.sqrt(norm2); 772 for (int j = 0; j < size; ++j) { 773 dataI[j] *= inv; 774 } 775 776 } while (norm2 * size < 0.01); 777 } 778 779 return MatrixUtils.createRealMatrix(data); 780 781 } 782 } 783