/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import java.util.Arrays;
import java.util.Random;
import org.apache.commons.math.linear.EigenDecomposition;
import org.apache.commons.math.linear.EigenDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.TriDiagonalTransformer;
import org.apache.commons.math.util.MathUtils;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class EigenDecompositionImplTest extends TestCase {
private double[] refValues;
private RealMatrix matrix;
public EigenDecompositionImplTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(EigenDecompositionImplTest.class);
suite.setName("EigenDecompositionImpl Tests");
return suite;
}
public void testDimension1() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(new double[][] { { 1.5 } });
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(1.5, ed.getRealEigenvalue(0), 1.0e-15);
}
public void testDimension2() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(new double[][] {
{ 59.0, 12.0 },
{ 12.0, 66.0 }
});
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(75.0, ed.getRealEigenvalue(0), 1.0e-15);
assertEquals(50.0, ed.getRealEigenvalue(1), 1.0e-15);
}
public void testDimension3() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(new double[][] {
{ 39632.0, -4824.0, -16560.0 },
{ -4824.0, 8693.0, 7920.0 },
{ -16560.0, 7920.0, 17300.0 }
});
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(50000.0, ed.getRealEigenvalue(0), 3.0e-11);
assertEquals(12500.0, ed.getRealEigenvalue(1), 3.0e-11);
assertEquals( 3125.0, ed.getRealEigenvalue(2), 3.0e-11);
}
public void testDimension4WithSplit() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(new double[][] {
{ 0.784, -0.288, 0.000, 0.000 },
{ -0.288, 0.616, 0.000, 0.000 },
{ 0.000, 0.000, 0.164, -0.048 },
{ 0.000, 0.000, -0.048, 0.136 }
});
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
}
public void testDimension4WithoutSplit() {
RealMatrix matrix =
MatrixUtils.createRealMatrix(new double[][] {
{ 0.5608, -0.2016, 0.1152, -0.2976 },
{ -0.2016, 0.4432, -0.2304, 0.1152 },
{ 0.1152, -0.2304, 0.3088, -0.1344 },
{ -0.2976, 0.1152, -0.1344, 0.3872 }
});
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(1.0, ed.getRealEigenvalue(0), 1.0e-15);
assertEquals(0.4, ed.getRealEigenvalue(1), 1.0e-15);
assertEquals(0.2, ed.getRealEigenvalue(2), 1.0e-15);
assertEquals(0.1, ed.getRealEigenvalue(3), 1.0e-15);
}
/** test a matrix already in tridiagonal form. */
public void testTridiagonal() {
Random r = new Random(4366663527842l);
double[] ref = new double[30];
for (int i = 0; i < ref.length; ++i) {
if (i < 5) {
ref[i] = 2 * r.nextDouble() - 1;
} else {
ref[i] = 0.0001 * r.nextDouble() + 6;
}
}
Arrays.sort(ref);
TriDiagonalTransformer t =
new TriDiagonalTransformer(createTestMatrix(r, ref));
EigenDecomposition ed =
new EigenDecompositionImpl(t.getMainDiagonalRef(),
t.getSecondaryDiagonalRef(),
MathUtils.SAFE_MIN);
double[] eigenValues = ed.getRealEigenvalues();
assertEquals(ref.length, eigenValues.length);
for (int i = 0; i < ref.length; ++i) {
assertEquals(ref[ref.length - i - 1], eigenValues[i], 2.0e-14);
}
}
/** test dimensions */
public void testDimensions() {
final int m = matrix.getRowDimension();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(m, ed.getV().getRowDimension());
assertEquals(m, ed.getV().getColumnDimension());
assertEquals(m, ed.getD().getColumnDimension());
assertEquals(m, ed.getD().getColumnDimension());
assertEquals(m, ed.getVT().getRowDimension());
assertEquals(m, ed.getVT().getColumnDimension());
}
/** test eigenvalues */
public void testEigenvalues() {
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
double[] eigenValues = ed.getRealEigenvalues();
assertEquals(refValues.length, eigenValues.length);
for (int i = 0; i < refValues.length; ++i) {
assertEquals(refValues[i], eigenValues[i], 3.0e-15);
}
}
/** test eigenvalues for a big matrix. */
public void testBigMatrix() {
Random r = new Random(17748333525117l);
double[] bigValues = new double[200];
for (int i = 0; i < bigValues.length; ++i) {
bigValues[i] = 2 * r.nextDouble() - 1;
}
Arrays.sort(bigValues);
EigenDecomposition ed =
new EigenDecompositionImpl(createTestMatrix(r, bigValues), MathUtils.SAFE_MIN);
double[] eigenValues = ed.getRealEigenvalues();
assertEquals(bigValues.length, eigenValues.length);
for (int i = 0; i < bigValues.length; ++i) {
assertEquals(bigValues[bigValues.length - i - 1], eigenValues[i], 2.0e-14);
}
}
/** test eigenvectors */
public void testEigenvectors() {
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
for (int i = 0; i < matrix.getRowDimension(); ++i) {
double lambda = ed.getRealEigenvalue(i);
RealVector v = ed.getEigenvector(i);
RealVector mV = matrix.operate(v);
assertEquals(0, mV.subtract(v.mapMultiplyToSelf(lambda)).getNorm(), 1.0e-13);
}
}
/** test A = VDVt */
public void testAEqualVDVt() {
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
RealMatrix v = ed.getV();
RealMatrix d = ed.getD();
RealMatrix vT = ed.getVT();
double norm = v.multiply(d).multiply(vT).subtract(matrix).getNorm();
assertEquals(0, norm, 6.0e-13);
}
/** test that V is orthogonal */
public void testVOrthogonal() {
RealMatrix v = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN).getV();
RealMatrix vTv = v.transpose().multiply(v);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
}
/** test diagonal matrix */
public void testDiagonal() {
double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
EigenDecomposition ed = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN);
assertEquals(diagonal[0], ed.getRealEigenvalue(3), 2.0e-15);
assertEquals(diagonal[1], ed.getRealEigenvalue(2), 2.0e-15);
assertEquals(diagonal[2], ed.getRealEigenvalue(1), 2.0e-15);
assertEquals(diagonal[3], ed.getRealEigenvalue(0), 2.0e-15);
}
/**
* Matrix with eigenvalues {8, -1, -1}
*/
public void testRepeatedEigenvalue() {
RealMatrix repeated = MatrixUtils.createRealMatrix(new double[][] {
{3, 2, 4},
{2, 0, 2},
{4, 2, 3}
});
EigenDecomposition ed = new EigenDecompositionImpl(repeated, MathUtils.SAFE_MIN);
checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
}
/**
* Matrix with eigenvalues {2, 0, 12}
*/
public void testDistinctEigenvalues() {
RealMatrix distinct = MatrixUtils.createRealMatrix(new double[][] {
{3, 1, -4},
{1, 3, -4},
{-4, -4, 8}
});
EigenDecomposition ed = new EigenDecompositionImpl(distinct, MathUtils.SAFE_MIN);
checkEigenValues((new double[] {2, 0, 12}), ed, 1E-12);
checkEigenVector((new double[] {1, -1, 0}), ed, 1E-12);
checkEigenVector((new double[] {1, 1, 1}), ed, 1E-12);
checkEigenVector((new double[] {-1, -1, 2}), ed, 1E-12);
}
/**
* Verifies that the given EigenDecomposition has eigenvalues equivalent to
* the targetValues, ignoring the order of the values and allowing
* values to differ by tolerance.
*/
protected void checkEigenValues(double[] targetValues,
EigenDecomposition ed, double tolerance) {
double[] observed = ed.getRealEigenvalues();
for (int i = 0; i < observed.length; i++) {
assertTrue(isIncludedValue(observed[i], targetValues, tolerance));
assertTrue(isIncludedValue(targetValues[i], observed, tolerance));
}
}
/**
* Returns true iff there is an entry within tolerance of value in
* searchArray.
*/
private boolean isIncludedValue(double value, double[] searchArray,
double tolerance) {
boolean found = false;
int i = 0;
while (!found && i < searchArray.length) {
if (Math.abs(value - searchArray[i]) < tolerance) {
found = true;
}
i++;
}
return found;
}
/**
* Returns true iff eigenVector is a scalar multiple of one of the columns
* of ed.getV(). Does not try linear combinations - i.e., should only be
* used to find vectors in one-dimensional eigenspaces.
*/
protected void checkEigenVector(double[] eigenVector,
EigenDecomposition ed, double tolerance) {
assertTrue(isIncludedColumn(eigenVector, ed.getV(), tolerance));
}
/**
* Returns true iff there is a column that is a scalar multiple of column
* in searchMatrix (modulo tolerance)
*/
private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix,
double tolerance) {
boolean found = false;
int i = 0;
while (!found && i < searchMatrix.getColumnDimension()) {
double multiplier = 1.0;
boolean matching = true;
int j = 0;
while (matching && j < searchMatrix.getRowDimension()) {
double colEntry = searchMatrix.getEntry(j, i);
// Use the first entry where both are non-zero as scalar
if (Math.abs(multiplier - 1.0) <= Math.ulp(1.0) && Math.abs(colEntry) > 1E-14
&& Math.abs(column[j]) > 1e-14) {
multiplier = colEntry / column[j];
}
if (Math.abs(column[j] * multiplier - colEntry) > tolerance) {
matching = false;
}
j++;
}
found = matching;
i++;
}
return found;
}
@Override
public void setUp() {
refValues = new double[] {
2.003, 2.002, 2.001, 1.001, 1.000, 0.001
};
matrix = createTestMatrix(new Random(35992629946426l), refValues);
}
@Override
public void tearDown() {
refValues = null;
matrix = null;
}
static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
final int n = eigenValues.length;
final RealMatrix v = createOrthogonalMatrix(r, n);
final RealMatrix d = createDiagonalMatrix(eigenValues, n, n);
return v.multiply(d).multiply(v.transpose());
}
public static RealMatrix createOrthogonalMatrix(final Random r, final int size) {
final double[][] data = new double[size][size];
for (int i = 0; i < size; ++i) {
final double[] dataI = data[i];
double norm2 = 0;
do {
// generate randomly row I
for (int j = 0; j < size; ++j) {
dataI[j] = 2 * r.nextDouble() - 1;
}
// project the row in the subspace orthogonal to previous rows
for (int k = 0; k < i; ++k) {
final double[] dataK = data[k];
double dotProduct = 0;
for (int j = 0; j < size; ++j) {
dotProduct += dataI[j] * dataK[j];
}
for (int j = 0; j < size; ++j) {
dataI[j] -= dotProduct * dataK[j];
}
}
// normalize the row
norm2 = 0;
for (final double dataIJ : dataI) {
norm2 += dataIJ * dataIJ;
}
final double inv = 1.0 / Math.sqrt(norm2);
for (int j = 0; j < size; ++j) {
dataI[j] *= inv;
}
} while (norm2 * size < 0.01);
}
return MatrixUtils.createRealMatrix(data);
}
public static RealMatrix createDiagonalMatrix(final double[] diagonal,
final int rows, final int columns) {
final double[][] dData = new double[rows][columns];
for (int i = 0; i < Math.min(rows, columns); ++i) {
dData[i][i] = diagonal[i];
}
return MatrixUtils.createRealMatrix(dData);
}
}