Fix issue 105 for EigenvalueDecomposition

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;
 	}

src/test/java/org/encog/mathutil/matrices/decomposition/EigenvalueDecompositionTest.java

@@ -0,0 +1,177 @@
+/*
+ * 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.*;
+
+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_asymentric() {
+		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_asymentric() {
+		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_asymentric() {
+		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_asymentric() {
+		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);
+	}
+
+	@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 },
+				}));
+	}
+}