Merge pull request #51 from StaticBeagle/adding-solve-to-matrix-sparse

Integrating decomposition and solve methods to matrix classes

Author: staticbeagle

README.md

@@ -36,12 +36,12 @@ Maven
 <dependency>
     <groupId>com.wildbitsfoundry</groupId>
     <artifactId>etk4j</artifactId>
-    <version>2.2.0</version>
+    <version>2.2.1</version>
 </dependency>
 ```
 Gradle
 ```bash
-implementation 'com.wildbitsfoundry:etk4j:2.2.0'
+implementation 'com.wildbitsfoundry:etk4j:2.2.1'
 ```
 
 Requirements

pom.xml

@@ -2,7 +2,7 @@
     <modelVersion>4.0.0</modelVersion>
     <groupId>com.wildbitsfoundry</groupId>
     <artifactId>etk4j</artifactId>
-    <version>2.2.0</version>
+    <version>2.2.1</version>
     <name>Engineering Toolkit for Java</name>
     <description>Tools and implementation of mathematical methods or backing math for engineering problems and
         applications.

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexLUDecompositionDense.java

@@ -4,16 +4,15 @@
 
 import java.util.Arrays;
 
-public class ComplexLUDecompositionDense {
+public class ComplexLUDecompositionDense extends ComplexLUDecomposition<ComplexMatrix> {
 	protected Complex[] _data;
-	protected final int _rows;
-	protected final int _cols;
 
 	protected int _pivotsign = 1;
 	protected int[] _pivot;
 
 	public ComplexLUDecompositionDense(ComplexMatrixDense matrix) {
-		final int rows = matrix.getRowCount();
+        super(matrix);
+        final int rows = matrix.getRowCount();
 		final int cols = matrix.getColumnCount();
 		Complex[] data = matrix.getArrayCopy();
 
@@ -67,28 +66,26 @@ public ComplexLUDecompositionDense(ComplexMatrixDense matrix) {
 			}
 		}
 		_data = data;
-		_rows = rows;
-		_cols = cols;
 	}
 
-	public boolean isNonSingular() {
-		for (int j = 0; j < _cols; ++j) {
-			if (_data[j * _cols + j].equals(new Complex())) {
-				return false;
+	@Override
+	public boolean isSingular() {
+		for (int j = 0; j < cols; ++j) {
+			if (_data[j * cols + j].equals(new Complex())) {
+				return true;
 			}
 		}
-		return true;
+		return false;
 	}
 
+
 	/**
 	 * Return lower triangular factor
 	 * 
 	 * @return L
 	 */
 
 	public ComplexMatrixDense getL() {
-		final int rows = _rows;
-		final int cols = _cols;
 		Complex[] L = new Complex[rows * cols];
 		for (int i = 0; i < rows; ++i) {
 			for (int j = 0; j < cols; ++j) {
@@ -111,8 +108,6 @@ public ComplexMatrixDense getL() {
 	 */
 
 	public ComplexMatrixDense getU() {
-		final int rows = _rows;
-		final int cols = _cols;
 		Complex[] U = new Complex[rows * cols];
 		for (int i = 0; i < cols; i++) {
 			for (int j = 0; j < cols; j++) {
@@ -126,6 +121,11 @@ public ComplexMatrixDense getU() {
 		return new ComplexMatrixDense(U, rows, cols);
 	}
 
+	@Override
+	public ComplexMatrix solve(double[] b) {
+		return solve(ComplexMatrixDense.fromRealMatrix(new MatrixDense(b, b.length)));
+	}
+
 	/**
 	 * Return pivot permutation vector
 	 * 
@@ -143,7 +143,6 @@ public int[] getPivot() {
 	 */
 
 	public double[] getPivotAsDouble() {
-		final int rows = _rows;
 		double[] vals = new double[rows];
 		for (int i = 0; i < rows; i++) {
 			vals[i] = (double) _pivot[i];
@@ -152,12 +151,12 @@ public double[] getPivotAsDouble() {
 	}
 
 	public Complex det() {
-		if (_rows != _cols) {
+		if (rows != cols) {
 			throw new IllegalArgumentException("Matrix must be square.");
 		}
 		Complex det = Complex.fromReal(_pivotsign);
-		for (int i = 0; i < _cols; i++) {
-			det.multiplyEquals(this._data[i * _cols + i]);
+		for (int i = 0; i < cols; i++) {
+			det.multiplyEquals(this._data[i * cols + i]);
 		}
 		return det;
 	}
@@ -190,10 +189,10 @@ public ComplexMatrixDense solve(MatrixDense B) {
 	 *                Matrix is singular.
 	 */
 	public ComplexMatrixDense solve(ComplexMatrixDense B) {
-		if (B.getRowCount() != _rows) {
+		if (B.getRowCount() != rows) {
 			throw new IllegalArgumentException("Matrix row dimensions must agree.");
 		}
-		if (!this.isNonSingular()) {
+		if (this.isSingular()) {
 			throw new RuntimeException("Matrix is singular.");
 		}
 
@@ -202,8 +201,6 @@ public ComplexMatrixDense solve(ComplexMatrixDense B) {
 		ComplexMatrixDense Xmat = B.subMatrix(_pivot, 0, nx - 1);
 		Complex[] X = Xmat.getArray();
 
-		final int cols = _cols;
-
 		// Solve L * Y = B(_pivot,:)
 		for (int k = 0; k < cols; ++k) {
 			for (int i = k + 1; i < cols; ++i) {
@@ -229,8 +226,8 @@ public ComplexMatrixDense solve(ComplexMatrixDense B) {
 	@Override
 	public String toString() {
 		StringBuilder sb = new StringBuilder();
-		for (int i = 0; i < _rows * _cols; ++i) {
-			if (i > 0 && i % _cols == 0) {
+		for (int i = 0; i < rows * cols; ++i) {
+			if (i > 0 && i % cols == 0) {
 				sb.append(System.lineSeparator());
 			}
 			sb.append(String.format("%.4f", _data[i])).append(" ");

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexMatrix.java

@@ -113,7 +113,7 @@ public boolean isSquare() {
      * @return The {@link CholeskyDecompositionDense} of the {@code Matrix}.
      * @see <a href="https://en.wikipedia.org/wiki/Cholesky_decomposition">Cholesky Decomposition</a>
      */
-    //public abstract ComplexCholeskyDecomposition<?> Chol();
+    public abstract ComplexCholeskyDecomposition<?> Chol();
     // endregion
 
     /***

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexMatrixDense.java

@@ -112,13 +112,41 @@ public double det() {
     }
 
     @Override
-    public ComplexLUDecomposition<?> LU() {
-        return null;
+    public ComplexLUDecompositionDense LU() {
+        return new ComplexLUDecompositionDense(this);
     }
 
     @Override
-    public ComplexQRDecomposition<?> QR() {
-        return null;
+    public ComplexQRDecompositionDense QR() {
+        return new ComplexQRDecompositionDense(this);
+    }
+
+    @Override
+    public ComplexCholeskyDecompositionDense Chol() {
+        return new ComplexCholeskyDecompositionDense(this);
+    }
+
+    public ComplexSingularValueDecompositionDense SVD() {
+        return new ComplexSingularValueDecompositionDense(this);
+    }
+
+    /**
+     * Eigenvalue decomposition. The {@link Matrix} is balanced ({@link MatrixDense#balance()}) prior to the decomposition.
+     *
+     * @return The {@link EigenvalueDecompositionDense} of the {@code Matrix}.
+     */
+    public ComplexEigenvalueDecompositionDense eig() {
+        return new ComplexEigenvalueDecompositionDense(this);
+    }
+
+    /**
+     * Schur Decomposition of the {@code Matrix}.
+     *
+     * @return The {@link SchurDecompositionDense} of the {@code Matrix}.
+     * @see <a href="https://en.wikipedia.org/wiki/Schur_decomposition">Schur Decomposition</a>
+     */
+    public ComplexSchurDecompositionDense Schur() {
+        return new ComplexSchurDecompositionDense(this);
     }
 
     public Complex[] getArray() {
@@ -285,7 +313,7 @@ public ComplexMatrixDense transpose() {
     }
 
     public ComplexMatrixDense inv() {
-        return this.solve(MatrixDense.Factory.identity(rows));
+        return this.solve(ComplexMatrixDense.Factory.identity(rows));
     }
 
     // region solve
@@ -509,10 +537,6 @@ public ComplexMatrixDense balance() {
         return new ComplexMatrixDense(data, rows, cols);
     }
 
-    public ComplexSingularValueDecompositionDense SVD() {
-        return new ComplexSingularValueDecompositionDense(this);
-    }
-
     public ComplexMatrixDense pinv() {
         int rows = this.rows;
         int cols = this.cols;

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDense.java

@@ -530,6 +530,16 @@ public CholeskyDecompositionDense Chol() {
     public SingularValueDecompositionDense SVD() {
         return new SingularValueDecompositionDense(this);
     }
+
+    /**
+     * Schur Decomposition of the {@code Matrix}.
+     *
+     * @return The {@link SchurDecompositionDense} of the {@code Matrix}.
+     * @see <a href="https://en.wikipedia.org/wiki/Schur_decomposition">Schur Decomposition</a>
+     */
+    public SchurDecompositionDense Schur() {
+        return new SchurDecompositionDense(this);
+    }
     // endregion
 
     /***

src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixSparse.java

@@ -63,7 +63,8 @@ public MatrixSparse(int rows, int cols, int arrayLength) {
 
     public MatrixSparse(MatrixSparse original) {
         this(original.rows, original.cols, original.nz_length);
-
+        nz_values = original.nz_values;
+        nz_rows = original.nz_rows;
         setTo(original);
     }
 
@@ -76,14 +77,13 @@ public MatrixSparse createLike() {
     }
 
     public void setTo(MatrixSparse original) {
-        MatrixSparse o = original;
-        reshape(o.rows, o.cols, o.nz_length);
-        this.nz_length = o.nz_length;
+        reshape(original.rows, original.cols, original.nz_length);
+        this.nz_length = original.nz_length;
 
-        System.arraycopy(o.nz_values, 0, nz_values, 0, nz_length);
-        System.arraycopy(o.nz_rows, 0, nz_rows, 0, nz_length);
-        System.arraycopy(o.col_idx, 0, col_idx, 0, cols + 1);
-        this.indicesSorted = o.indicesSorted;
+        System.arraycopy(original.nz_values, 0, nz_values, 0, nz_length);
+        System.arraycopy(original.nz_rows, 0, nz_rows, 0, nz_length);
+        System.arraycopy(original.col_idx, 0, col_idx, 0, cols + 1);
+        this.indicesSorted = original.indicesSorted;
     }
 
 //    public void print() {
@@ -182,6 +182,46 @@ public void set(int row, int col, double val) {
         unsafeSet(row, col, val);
     }
 
+    /**
+     * Solve system of linear equations. Three different algorithms are used depending on the shape of the matrix:
+     * <pre>
+     *     LU Decomposition if the matrix is squared.
+     *     QR if the matrix is thin in other words it has more rows than columns. (Overdetermined system)
+     *     Pseudo inverse * b if the matrix is short and wide in other words it has more columns than rows. (Under-determined system)
+     * </pre>
+     *
+     * @param b The solution {@link Matrix}.
+     * @return The solution to {@code Ax = b}
+     */
+    public MatrixSparse solve(MatrixSparse b) {
+        if (rows == cols) { // Matrix is Squared
+            return new LUDecompositionSparse(this).solve(b);
+        } else if (rows > cols) { // Matrix is tall and narrow (Overdetermined system)
+            return new QRDecompositionSparse(this).solve(b);
+        } else { // Matrix is short and wide (Under-determined system)
+            throw new UnsupportedOperationException("This operation is not supported for system with more columns than rows");
+        }
+    }
+
+    /**
+     * Solve system of linear equations. Three different algorithms are used depending on the shape of the matrix:
+     * <pre>
+     *     LU Decomposition if the matrix is squared.
+     *     QR if the matrix is thin in other words it has more rows than columns. (Overdetermined system)
+     *     Pseudo inverse * b if the matrix is short and wide in other words it has more columns than rows. (Under-determined system)
+     * </pre>
+     *
+     * @param b The solution {@link Matrix}.
+     * @return The solution to {@code Ax = b}
+     */
+    public MatrixSparse solve(double[] b) {
+        double[][] matrix = new double[b.length][1];
+        for(int i = 0; i < b.length; i++) {
+            matrix[i][0] = b[i];
+        }
+        return solve(from2DArray(matrix));
+    }
+
     @Override
     public double det() {
         throw new UnsupportedOperationException("Not implemented yet");
@@ -198,9 +238,7 @@ public QRDecompositionSparse QR() {
     }
 
     @Override
-    public CholeskyDecompositionSparse Chol() {
-        throw new UnsupportedOperationException("Not implemented yet");
-    }
+    public CholeskyDecompositionSparse Chol() { return new CholeskyDecompositionSparse(this); }
 
     @Override
     public boolean isEmpty() {

src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixSparseTest.java

@@ -127,6 +127,20 @@ public void testLUDecomposition() {
                 solution.get(1, 0), solution.get(2, 0)}, 1e-12);
     }
 
+    @Test
+    public void testSolveArrayRHS() {
+        double[][] matrix = {
+                {1, 4, 7},
+                {2, 5, 8},
+                {3, 6, 10},
+        };
+        MatrixSparse sparseCSC = MatrixSparse.from2DArray(matrix, ConstantsETK.DOUBLE_EPS);
+        double[] b = {1, 2, 0};
+        MatrixSparse solution = sparseCSC.solve(b);
+        assertArrayEquals(new double[]{-2, 6.000000000000005, -3.0000000000000027}, new double[]{solution.get(0, 0),
+                solution.get(1, 0), solution.get(2, 0)}, 1e-12);
+    }
+
     @Test
     public void testLUDecompositionSparseMatrixRHS() {
         double[][] matrix = {