From 20b0b4ecbc330659d6949e17b5ba0f4fa3e3fe2b Mon Sep 17 00:00:00 2001 From: StaticBeagle Date: Fri, 10 Jan 2025 20:15:36 -0600 Subject: [PATCH 1/2] Initial commit --- .../ComplexHouseholderTransformations.java | 2 +- .../HessembergDecompositionDense.java | 66 ++++ .../HouseholderTransformations.java | 341 ++++++++++++++++++ .../etk4j/util/DoubleArrays.java | 26 ++ 4 files changed, 434 insertions(+), 1 deletion(-) create mode 100644 src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java create mode 100644 src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HouseholderTransformations.java diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHouseholderTransformations.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHouseholderTransformations.java index 52dd8c1..9918354 100644 --- a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHouseholderTransformations.java +++ b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHouseholderTransformations.java @@ -58,7 +58,7 @@ public static Complex[] genc(ComplexMatrixDense A, int r1, int r2, int c) scale.multiplyEquals(t); } - t = Complex.fromReal(1).divide(scale).uminus();; + t = Complex.fromReal(1).divide(scale).uminus(); A.unsafeSet(r1, c, t); for (i = 0; i < ru; i++) { diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java new file mode 100644 index 0000000..63c78b6 --- /dev/null +++ b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java @@ -0,0 +1,66 @@ +package com.wildbitsfoundry.etk4j.math.linearalgebra; + +import com.wildbitsfoundry.etk4j.math.complex.Complex; +import com.wildbitsfoundry.etk4j.util.ComplexArrays; + +/** + * Zhess implements the unitary reduction to Hessenberg form + * by a unitary similarity transformation. Specifically, given + * a square matrix A, there is a unitary matrix U such that + * 
+ *      H = U^H AU
+ *
+ * is upper Hessenberg. + * Zhess represents U and H as Zmats. + * + * @version Pre-alpha + * @author G. W. Stewart + */ +public class HessembergDecompositionDense { + + /** The upper Hessenberg matrix */ + public MatrixDense H; + + /** The unitary matrix */ + public MatrixDense U; + + /** Creates a Zhess from a square Zmat. Throws a + * JampackException for nonsquare matrx. + * + * @param A A Zmat + * Thrown if A is not square. + */ + public HessembergDecompositionDense(MatrixDense A) { + + if (A.getRowCount() != A.getColumnCount()) { + //throw new JampackException("Matrix not square"); TODO + } + + H = new MatrixDense(A); + U = MatrixDense.Factory.identity(H.getRowCount()); + + double[] work = new double[H.getRowCount()]; + + for (int k = 0; k < H.getColumnCount() - 2; k++) { + double[] u = HouseholderTransformations.genc(H, k + 1, H.getRowCount() - 1, k); + HouseholderTransformations.ua(u, H, k + 1, H.getRowCount() - 1, k + 1, H.getColumnCount() - 1, work); + HouseholderTransformations.au(H, u, 0, H.getRowCount() - 1, k + 1, H.getColumnCount() - 1, work); + HouseholderTransformations.au(U, u, 0, U.getRowCount() - 1, k + 1, U.getColumnCount() - 1, work); + } + } + + public static void main(String[] args) { + double[][] matrix = { + {65, 35, 40, 69}, + {99, 64, 37, 2}, + {39, 48,35, 90}, + {30, 93, 87, 17} + }; + + MatrixDense A = MatrixDense.from2DArray(matrix); + HessembergDecompositionDense hess = new HessembergDecompositionDense(A); + System.out.println(hess.H); + System.out.println(); + System.out.println(hess.U); // TODO do we have getters for the complex case? + } +} diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HouseholderTransformations.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HouseholderTransformations.java new file mode 100644 index 0000000..4e37872 --- /dev/null +++ b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HouseholderTransformations.java @@ -0,0 +1,341 @@ +package com.wildbitsfoundry.etk4j.math.linearalgebra; + +import com.wildbitsfoundry.etk4j.math.complex.Complex; +import com.wildbitsfoundry.etk4j.util.ComplexArrays; +import com.wildbitsfoundry.etk4j.util.DoubleArrays; + +public class HouseholderTransformations { + + /** + * Generates a Householder transformation from within the part of + * column c of a ComplexMatrixDense (altered) extending from rows + * r1 to r2. The method overwrites the + * column with the result of applying the transformation. + * + * @param A The matrix from which the transformation + * is to be generated (altered) + * @param r1 The index of the row in which the generating column + * begins + * @param r2 The index of the row in which the generating column + * ends + * @param c The index of the generating column + * @return A Complex[] of length r2-r1+1 + * containing the Householder vector + * @throws RuntimeException Passed from below. + */ + public static double[] genc(MatrixDense A, int r1, int r2, int c) + throws RuntimeException { + + int i, ru; + double norm; + double s; + double scale; + double t; + double t1; + + ru = r2 - r1 + 1; + + double[] u = new double[r2 - r1 + 1]; + + for (i = r1; i <= r2; i++) { + u[i - r1] = A.unsafeGet(i, c); + A.unsafeSet(i, c, 0); + } + + norm = DoubleArrays.normFro(u); + + if (r1 == r2 || norm == 0) { + A.unsafeSet(r1, c, -u[0]); + u[0] = Math.sqrt(2); + return u; + } + + scale = 1 / norm; + + if (u[0] != 0) { + t = u[0]; + t1 = t; + t = t1 / Math.abs(t); + scale *= t; + } + + t = -1 / scale; + A.unsafeSet(r1, c, t); + + for (i = 0; i < ru; i++) { + u[i] *= scale; + } + + u[0] = u[0] + 1; + s = Math.sqrt(1 / u[0]); + + for (i = 0; i < ru; i++) { + u[i] = s * u[i]; + } + + return u; + } + + /** + * Generates a Householder transformation from within the part of row + * r of a ComplexMatrixDense (altered) extending from columns c1 to + * c2. The method overwrites the row with the result + * of applying the transformation. + * + * @param A The matrix from which the transformation + * is to be generated (altered) + * @param r The index of the generating row + * @param c1 The index of the column in which the generating row + * begins + * @param c2 The index of the column in which the generating row + * ends + * @return A Complex[] of length r2-r1+1 + * containing the Householder vector + * @throws RuntimeException Passed from below. + */ + public static Complex[] genr(ComplexMatrixDense A, int r, int c1, int c2) + throws RuntimeException { + + int j, cu; + double norm, s; + Complex scale; + Complex t; + Complex t1; + + cu = c2 - c1 + 1; + + Complex[] u = new Complex[cu]; + + for (j = c1; j <= c2; j++) { + u[j - c1] = A.unsafeGet(r, j).copy(); + A.unsafeSet(r, j, new Complex()); + } + + norm = ComplexArrays.normFro(u); + + if (c1 == c2 || norm == 0) { + A.unsafeSet(r, c1, u[0].uminus()); + u[0] = Complex.fromReal(Math.sqrt(2)); + return u; + } + + scale = new Complex(1 / norm, 0); + + if (u[0].real() != 0 || u[0].imag() != 0) { + t = u[0]; + t1 = t.conj(); + t = t1.divide(t.abs()); + scale.multiplyEquals(t); + } + + t = Complex.fromReal(1).divide(scale).uminus(); + A.unsafeSet(r, c1, t); + + for (j = 0; j < cu; j++) { + u[j].multiplyEquals(scale); + } + + u[0] = Complex.fromReal(u[0].real() + 1); + s = Math.sqrt(1 / u[0].real()); + + for (j = 0; j < cu; j++) { + u[j] = new Complex(s * u[j].real(), -s * u[j].imag()); + } + + return u; + } + + /** + * Premultiplies the Householder transformation contained in a + * Complex[] into a ComplexMatrixDense A[r1:r2,c1:c2] and overwrites + * ComplexMatrixDense A[r1:r2,c1:c2] with the results. If r1 > r2 + * or c1 > c2 the method does nothing. + * + * @param u The Householder vector + * @param A The ComplexMatrixDense to which the transformation + * is to be applied (altered) + * @param r1 The index of the first row to which the transformation + * is to be applied + * @param r2 The index of the last row to which the transformation + * is to be applied + * @param c1 The index of the first column to which the transformation + * is index of the to be applied + * @param c2 The index of the last column to which the transformation + * is to be applied + * @param v A work array of length at least c2-c1+1 + * @return The transformed ComplexMatrixDense A + * @throws RuntimeException Thrown if either u or v is too short. + */ + public static MatrixDense ua(double[] u, MatrixDense A, int r1, int r2, int c1, int c2, double[] v) + throws RuntimeException { + + int i, j; + + + if (r2 < r1 || c2 < c1) { + return A; + } + + if (r2 - r1 + 1 > u.length) { + throw new RuntimeException + ("Householder vector too short."); + } + + if (c2 - c1 + 1 > v.length) { + throw new RuntimeException + ("Work vector too short."); + } + + for (j = c1; j <= c2; j++) { + v[j - c1] = 0; + } + + for (i = r1; i <= r2; i++) { + for (j = c1; j <= c2; j++) { + double a = v[j - c1]; + double b = u[i - r1]; + double c = A.unsafeGet(i, j); + + v[j - c1] = a + b * c; + } + } + + for (i = r1; i <= r2; i++) { + for (j = c1; j <= c2; j++) { + double a = A.unsafeGet(i, j); + double b = u[i - r1]; + double c = v[j - c1]; + + + A.unsafeSet(i, j, a - b * c); + } + } + return A; + } + + /** + * Premultiplies the Householder transformation contained in a + * Complex[] into a ComplexMatrixDense A[r1:r2,c1:c2] and overwrites + * ComplexMatrixDense A[r1:r2,c1:c2] with the results. If r1 > r2 + * or c1 > c2 the method does nothing. + * + * @param u The Householder vector + * @param A The ComplexMatrixDense to which the transformation + * is to be applied (altered) + * @param r1 The index of the first row to which the transformation + * is to be applied + * @param r2 The index of the last row to which the transformation + * is to be applied + * @param c1 The index of the first column to which the transformation + * is index of the to be applied + * @param c2 The index of the last column to which the transformation + * is to be applied + * @return The transformed ComplexMatrixDense A + * @throws RuntimeException Passed from below. + */ + public static MatrixDense ua(double[] u, MatrixDense A, int r1, int r2, int c1, int c2) + throws RuntimeException { + + if (c1 > c2) { + return A; + } + + return ua(u, A, r1, r2, c1, c2, new double[c2 - c1 + 1]); + } + + + /** + * Postmultiplies the Householder transformation contained in a + * Complex[] into a ComplexMatrixDense A[r1:r2,c1:c2] and overwrites + * ComplexMatrixDense A[r1:r2,c1:c2] with the results. If r1 > r2 + * or c1 > c2 the method does nothing. + * + * @param u The Householder vector + * @param A The ComplexMatrixDense to which the transformation + * is to be applied (altered) + * @param r1 The index of the first row to which the transformation + * is to be applied + * @param r2 The index of the last row to which the transformation + * is to be applied + * @param c1 The index of the first column to which the transformation + * is index of the to be applied + * @param c2 The index of the last column to which the transformation + * is to be applied + * @param v A work array of length at least c2-c1+1 + * @return The transformed ComplexMatrixDense A + * @throws RuntimeException Thrown if either u or v is too short. + */ + public static MatrixDense au(MatrixDense A, double[] u, int r1, int r2, int c1, int c2, double[] v) + throws RuntimeException { + + int i, j, cu; + + if (r2 < r1 || c2 < c1) { + return A; + } + + if (c2 - c1 + 1 > u.length) { + throw new RuntimeException + ("Householder vector too short."); + } + + if (r2 - r1 + 1 > v.length) { + throw new RuntimeException + ("Work vector too short."); + } + + for (i = r1; i <= r2; i++) { + v[i - r1] = 0; + for (j = c1; j <= c2; j++) { + double a = v[i - r1]; + double b = A.unsafeGet(i, j); + double c = u[j - c1]; + + v[i - r1] = a + b * c; + } + } + for (i = r1; i <= r2; i++) { + for (j = c1; j <= c2; j++) { + double a = A.unsafeGet(i, j); + double b = v[i - r1]; + double c = u[j - c1]; + + A.unsafeSet(i, j, a - b * c); + } + } + return A; + } + + + /** + * Postmultiplies the Householder transformation contained in a + * Complex[] into a ComplexMatrixDense A[r1:r2,c1:c2] and overwrites + * ComplexMatrixDense A[r1:r2,c1:c2] with the results. If r1 > r2 + * or c1 > c2 the method does nothing. + * + * @param u The Householder vector + * @param A The ComplexMatrixDense to which the transformation + * is to be applied (altered) + * @param r1 The index of the first row to which the transformation + * is to be applied + * @param r2 The index of the last row to which the transformation + * is to be applied + * @param c1 The index of the first column to which the transformation + * is index of the to be applied + * @param c2 The index of the last column to which the transformation + * is to be applied + * @return The transformed ComplexMatrixDense A + * @throws RuntimeException Passed from below. + */ + + public static MatrixDense au(MatrixDense A, double[] u, int r1, int r2, int c1, int c2) + throws RuntimeException { + + if (r2 < r1) { + return A; + } + + return au(A, u, r1, r2, c1, c2, new double[r2 - r1 + 1]); + } + +} diff --git a/src/main/java/com/wildbitsfoundry/etk4j/util/DoubleArrays.java b/src/main/java/com/wildbitsfoundry/etk4j/util/DoubleArrays.java index f257f76..f467abe 100644 --- a/src/main/java/com/wildbitsfoundry/etk4j/util/DoubleArrays.java +++ b/src/main/java/com/wildbitsfoundry/etk4j/util/DoubleArrays.java @@ -910,4 +910,30 @@ public static double[] max(double[] a, double[] b) { } return result; } + + public static double normFro(double[] a) { + int i; + double fac, nrm, scale; + + int n = a.length; + + scale = 0.0; + for (i = 0; i < n; i++) { + scale = Math.max(scale, Math.abs(a[i])); + } + if (scale == 0) { + return 0.0; + } + if (scale < 1) { + scale = scale * 1.0e20; + } + scale = 1 / scale; + nrm = 0; + + for (i = 0; i < n; i++) { + fac = scale * a[i]; + nrm = nrm + fac * fac; + } + return Math.sqrt(nrm) / scale; + } } From d196afbf750464f177ab6cbdd122e2320d5e6fc1 Mon Sep 17 00:00:00 2001 From: StaticBeagle Date: Sat, 11 Jan 2025 01:19:05 -0600 Subject: [PATCH 2/2] Added unit tests --- .../ComplexHessembergDecompositionDense.java | 8 +++++++ .../HessembergDecompositionDense.java | 17 +++++--------- .../math/linearalgebra/MatrixDenseTest.java | 22 +++++++++++++++++++ 3 files changed, 35 insertions(+), 12 deletions(-) diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHessembergDecompositionDense.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHessembergDecompositionDense.java index 1069fef..b7091dd 100644 --- a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHessembergDecompositionDense.java +++ b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/ComplexHessembergDecompositionDense.java @@ -48,4 +48,12 @@ public ComplexHessembergDecompositionDense(ComplexMatrixDense A) { ComplexHouseholderTransformations.au(U, u, 0, U.getRowCount() - 1, k + 1, U.getColumnCount() - 1, work); } } + + public ComplexMatrixDense getH() { + return H; + } + + public ComplexMatrixDense getU() { + return U; + } } diff --git a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java index 63c78b6..3b135cf 100644 --- a/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java +++ b/src/main/java/com/wildbitsfoundry/etk4j/math/linearalgebra/HessembergDecompositionDense.java @@ -49,18 +49,11 @@ public HessembergDecompositionDense(MatrixDense A) { } } - public static void main(String[] args) { - double[][] matrix = { - {65, 35, 40, 69}, - {99, 64, 37, 2}, - {39, 48,35, 90}, - {30, 93, 87, 17} - }; + public MatrixDense getH() { + return H; + } - MatrixDense A = MatrixDense.from2DArray(matrix); - HessembergDecompositionDense hess = new HessembergDecompositionDense(A); - System.out.println(hess.H); - System.out.println(); - System.out.println(hess.U); // TODO do we have getters for the complex case? + public MatrixDense getU() { + return U; } } diff --git a/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java b/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java index a924b30..ecf3cf6 100644 --- a/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java +++ b/src/test/java/com/wildbitsfoundry/etk4j/math/linearalgebra/MatrixDenseTest.java @@ -80,6 +80,28 @@ public void testGetCol() { assertArrayEquals(row, m.getCol(0), 1e-12); } + @Test + public void testHessembergDecomposition() { + double[][] matrix = { + {65, 35, 40, 69}, + {99, 64, 37, 2}, + {39, 48,35, 90}, + {30, 93, 87, 17} + }; + + MatrixDense A = MatrixDense.from2DArray(matrix); + HessembergDecompositionDense hess = new HessembergDecompositionDense(A); + + double[] expectedH = {65.0, -64.17727312578465, -58.72784396012248, 4.279948356457808, -110.55315463612966, + 123.81148748159052, 17.486023979786594, 22.14752145860223, 0.0, 100.60648832932314, 36.36015423097008, + 26.39970581394131, 0.0, 0.0, 58.20341156591577, -44.17164171256066}; + assertArrayEquals(expectedH, hess.getH().getArray(), 1e-12); + double[] expectedU = {1.0, 0.0, 0.0, 0.0, 0.0, -0.8954968343132741, 0.39724801290719464, 0.2006973741138373, + 0.0, -0.3527714801840171, -0.3585884935028961, -0.8642722806477718, 0.0, -0.2713626770646286, + -0.8447534010399773, 0.46125263026644026}; + assertArrayEquals(expectedU, hess.getU().getArray(), 1e-12); + } + @Test public void allTests() { MatrixDense A, B, C, Z, O, I, R, S, X, SUB, M, T, SQ, DEF, SOL;