From 79053bcebeb97b9d013273e06c296e82dd8fae6c Mon Sep 17 00:00:00 2001 From: Ian Flanigan Date: Mon, 28 May 2012 15:52:07 +0200 Subject: [PATCH 1/2] Added test for EgenvalueDecomposition --- .../EigenvalueDecompositionTest.java | 132 ++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java diff --git a/src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java b/src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java new file mode 100644 index 000000000..c106d603b --- /dev/null +++ b/src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java @@ -0,0 +1,132 @@ +/* + * Encog(tm) Core v3.1 - Java Version + * http://www.heatonresearch.com/encog/ + * http://code.google.com/p/encog-java/ + + * Copyright 2008-2012 Heaton Research, Inc. + * + * Licensed 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. + * + * For more information on Heaton Research copyrights, licenses + * and trademarks visit: + * http://www.heatonresearch.com/copyright + */ +package org.encog.mathutil.matrices.decomposition; + +import static org.junit.Assert.assertEquals; + +import org.encog.mathutil.matrices.Matrix; +import org.junit.Test; + +/** + * Simple tests for the EigenvalueDecomposition object. + */ +public class EigenvalueDecompositionTest { + + private Matrix matrix1 = new Matrix( + new double[][] { + { 2.0, 5.1 }, + { -0.5, 0.9 } + }); + + private Matrix matrix2 = new Matrix( + new double[][] { + { 2.0, 5.5 }, + { 5.5, -1.0 } + }); + + @Test + public void testEigenvalueDecomposition() { + new EigenvalueDecomposition(matrix1); + } + + @Test + public void testGetD_assymentric() { + EigenvalueDecomposition decomposition = + new EigenvalueDecomposition(matrix1); + Matrix d = decomposition.getD(); + assertEquals(2, d.getCols()); + assertEquals(2, d.getRows()); + assertEquals(1.45, d.get(0, 0), 1e-10); + assertEquals(1.4991664351, d.get(0, 1), 1e-10); + assertEquals(-1.4991664351, d.get(1, 0), 1e-10); + assertEquals(1.45, d.get(1, 1), 1e-10); + } + + @Test + public void testGetD_symentric() { + EigenvalueDecomposition decomposition = + new EigenvalueDecomposition(matrix2); + Matrix d = decomposition.getD(); + assertEquals(2, d.getCols()); + assertEquals(2, d.getRows()); + assertEquals(-5.2008771255, d.get(0, 0), 1e-10); + assertEquals(0.0, d.get(0, 1), 1e-10); + assertEquals(0.0, d.get(1, 0), 1e-10); + assertEquals(6.2008771255, d.get(1, 1), 1e-10); + } + + @Test + public void testGetImagEigenvalues_assymentric() { + EigenvalueDecomposition decomposition = + new EigenvalueDecomposition(matrix1); + double[] d = decomposition.getImagEigenvalues(); + assertEquals(2, d.length); + assertEquals(1.4991664351, d[0], 1e-10); + assertEquals(-1.4991664351, d[1], 1e-10); + } + + @Test + public void testGetImagEigenvalues_symentric() { + EigenvalueDecomposition decomposition = + new EigenvalueDecomposition(matrix2); + double[] d = decomposition.getImagEigenvalues(); + assertEquals(2, d.length); + assertEquals(0.0, d[0], 1e-10); + assertEquals(0.0, d[1], 1e-10); + } + + @Test + public void testGetRealEigenvalues_assymentric() { + EigenvalueDecomposition decomposition = + new EigenvalueDecomposition(matrix1); + double[] d = decomposition.getRealEigenvalues(); + assertEquals(2, d.length); + assertEquals(1.45, d[0], 1e-10); + assertEquals(1.45, d[1], 1e-10); + } + + @Test + public void testGetRealEigenvalues_symentric() { + EigenvalueDecomposition decomposition = + new EigenvalueDecomposition(matrix2); + double[] d = decomposition.getRealEigenvalues(); + assertEquals(2, d.length); + assertEquals(-5.2008771255, d[0], 1e-10); + assertEquals(6.2008771255, d[1], 1e-10); + } + + @Test + public void testGetV_assymentric() { + EigenvalueDecomposition decomposition = + new EigenvalueDecomposition(matrix1); + Matrix d = decomposition.getV(); + assertEquals(2, d.getCols()); + assertEquals(2, d.getRows()); + assertEquals(2.9983328701, d.get(0, 0), 1e-10); + assertEquals(-1.1, d.get(0, 1), 1e-10); + assertEquals(0.0, d.get(1, 0), 1e-10); + assertEquals(1.0, d.get(1, 1), 1e-10); + } + +} From 781178a732b71babc682090a8a489d3f2bb8693a Mon Sep 17 00:00:00 2001 From: Ian Flanigan Date: Mon, 28 May 2012 18:23:49 +0200 Subject: [PATCH 2/2] Optimize EigenvalueDecomposition --- .../EigenvalueDecomposition.java | 90 ++++++++++++------- .../EigenvalueDecompositionTest.java | 55 ++++++++++-- 2 files changed, 107 insertions(+), 38 deletions(-) diff --git a/src/main/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecomposition.java b/src/main/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecomposition.java index 99dee1993..99f2b8fd7 100644 --- a/src/main/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecomposition.java +++ b/src/main/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecomposition.java @@ -23,6 +23,8 @@ */ package org.encog.mathutil.matrices.decomposition; +import java.util.Arrays; + import org.encog.mathutil.EncogMath; import org.encog.mathutil.matrices.Matrix; @@ -55,7 +57,7 @@ public class EigenvalueDecomposition { /** * Symmetry flag. */ - private boolean issymmetric; + private final boolean issymmetric; /** * Arrays for internal storage of eigenvalues. @@ -105,18 +107,12 @@ public EigenvalueDecomposition(final Matrix matrix) { this.d = new double[this.n]; this.e = new double[this.n]; - this.issymmetric = true; - for (int j = 0; (j < this.n) & this.issymmetric; j++) { - for (int i = 0; (i < this.n) & this.issymmetric; i++) { - this.issymmetric = (a[i][j] == a[j][i]); - } - } + this.issymmetric = isSymmetric(a); if (this.issymmetric) { + // Copy matrix a to v. for (int i = 0; i < this.n; i++) { - for (int j = 0; j < this.n; j++) { - this.v[i][j] = a[i][j]; - } + System.arraycopy(a[i], 0, v[i], 0, this.n); } // Tridiagonalize. @@ -129,10 +125,9 @@ public EigenvalueDecomposition(final Matrix matrix) { this.h = new double[this.n][this.n]; this.ort = new double[this.n]; + // Copy matrix a to h. for (int j = 0; j < this.n; j++) { - for (int i = 0; i < this.n; i++) { - this.h[i][j] = a[i][j]; - } + System.arraycopy(a, 0, h, 0, n); } // Reduce to Hessenberg form. @@ -143,6 +138,45 @@ public EigenvalueDecomposition(final Matrix matrix) { } } + /** + * Returns whether the given arrays make a symmetric matrix. A symmetric + * matrix is defined as a square matrix that is identical when flipped + * around its diagonal. A matrix with no rows and no columns is defined to + * be symmetric. + * + * @param a the matrix to analyze. + * @return {@code true} iff the matrix is symmetric. Malformed matrices are + * considered asymetric. + */ + static boolean isSymmetric(final double[][] a) { + // TODO: Perhaps move this to the Matrix class. + // Note that we only need to analyze the positions on one side of the + // diagonal as the diagonal always stays the same. + final int len = a.length; + if (len == 0) { + return true; + } + + // Because we skip the first row, verify its length explicitly + if (a[0].length != len) { + return false; + } + + // Loop through all of the rows, skipping the first + for (int j = 1; j < len; j++) { + if (a[j].length != len) { + return false; + } + // Loop through all of the columns up to the diagonal + for (int i = 0; i < j; i++) { + if (a[i][j] != a[j][i]) { + return false; + } + } + } + return true; + } + // Symmetric tridiagonal QL algorithm. private void cdiv(final double xr, final double xi, final double yr, @@ -171,9 +205,7 @@ public Matrix getD() { final Matrix X = new Matrix(this.n, this.n); final double[][] D = X.getData(); for (int i = 0; i < this.n; i++) { - for (int j = 0; j < this.n; j++) { - D[i][j] = 0.0; - } + Arrays.fill(D[i], 0.0); D[i][i] = this.d[i]; if (this.e[i] > 0) { D[i][i + 1] = this.e[i]; @@ -647,9 +679,7 @@ private void hqr2() { for (int i = 0; i < nn; i++) { if ((i < low) | (i > high)) { - for (int j = i; j < nn; j++) { - this.v[i][j] = this.h[i][j]; - } + System.arraycopy(this.h, i, this.v, i, nn - i); } } @@ -731,10 +761,10 @@ private void orthes() { // Accumulate transformations (Algol's ortran). + // Fill v's diagonal with 1s. for (int i = 0; i < this.n; i++) { - for (int j = 0; j < this.n; j++) { - this.v[i][j] = (i == j ? 1.0 : 0.0); - } + Arrays.fill(this.v[i], 0.0); + this.v[i][i] = 1.0; } for (int m = high - 1; m >= low + 1; m--) { @@ -885,10 +915,8 @@ private void tred2() { // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding // Fortran subroutine in EISPACK. - for (int j = 0; j < this.n; j++) { - this.d[j] = this.v[this.n - 1][j]; - } - + System.arraycopy(this.v[this.n - 1], 0, this.d, 0, this.n); + // Householder reduction to tridiagonal form. for (int i = this.n - 1; i > 0; i--) { @@ -923,9 +951,7 @@ private void tred2() { this.e[i] = scale * g; h = h - f * g; this.d[i - 1] = f - g; - for (int j = 0; j < i; j++) { - this.e[j] = 0.0; - } + Arrays.fill(this.e, 0, i, 0.0); // Apply similarity transformation to remaining columns. @@ -985,10 +1011,8 @@ private void tred2() { this.v[k][i + 1] = 0.0; } } - for (int j = 0; j < this.n; j++) { - this.d[j] = this.v[this.n - 1][j]; - this.v[this.n - 1][j] = 0.0; - } + System.arraycopy(this.v[this.n - 1], 0, this.d, 0, this.n); + Arrays.fill(this.v[this.n - 1], 0.0); this.v[this.n - 1][this.n - 1] = 1.0; this.e[0] = 0.0; } diff --git a/src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java b/src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java index c106d603b..cde832216 100644 --- a/src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java +++ b/src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java @@ -23,7 +23,7 @@ */ package org.encog.mathutil.matrices.decomposition; -import static org.junit.Assert.assertEquals; +import static org.junit.Assert.*; import org.encog.mathutil.matrices.Matrix; import org.junit.Test; @@ -51,7 +51,7 @@ public void testEigenvalueDecomposition() { } @Test - public void testGetD_assymentric() { + public void testGetD_asymentric() { EigenvalueDecomposition decomposition = new EigenvalueDecomposition(matrix1); Matrix d = decomposition.getD(); @@ -77,7 +77,7 @@ public void testGetD_symentric() { } @Test - public void testGetImagEigenvalues_assymentric() { + public void testGetImagEigenvalues_asymentric() { EigenvalueDecomposition decomposition = new EigenvalueDecomposition(matrix1); double[] d = decomposition.getImagEigenvalues(); @@ -97,7 +97,7 @@ public void testGetImagEigenvalues_symentric() { } @Test - public void testGetRealEigenvalues_assymentric() { + public void testGetRealEigenvalues_asymentric() { EigenvalueDecomposition decomposition = new EigenvalueDecomposition(matrix1); double[] d = decomposition.getRealEigenvalues(); @@ -117,7 +117,7 @@ public void testGetRealEigenvalues_symentric() { } @Test - public void testGetV_assymentric() { + public void testGetV_asymentric() { EigenvalueDecomposition decomposition = new EigenvalueDecomposition(matrix1); Matrix d = decomposition.getV(); @@ -129,4 +129,49 @@ public void testGetV_assymentric() { assertEquals(1.0, d.get(1, 1), 1e-10); } + @Test + public void testIsSymmetric() { + assertTrue( + "An empty matrix should be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { })); + assertFalse( + "An matrix with a row but no values should not be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { { } })); + assertTrue( + "A singleton matrix should be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { { 0 } })); + assertFalse( + "A malformed matrix should not be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { + { 0, 1 } + })); + assertTrue( + "This simple matrix should be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { + { 0, 1 }, + { 1, 0 }, + })); + assertFalse( + "This simple matrix should not be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { + { 0, 1 }, + { 2, 0 }, + })); + assertTrue( + "This matrix should be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { + { 0, 1, 3, 8 }, + { 1, 5, 7, 1 }, + { 3, 7, 9, 0 }, + { 8, 1, 0, 6 }, + })); + assertFalse( + "This matrix should not be symmetric", + EigenvalueDecomposition.isSymmetric(new double[][] { + { 8, 1, 0, 6 }, + { 3, 7, 9, 0 }, + { 1, 5, 7, 1 }, + { 0, 1, 3, 8 }, + })); + } }