From 61b0610048b2ac4cbde4f2868c9b167dcc2ae745 Mon Sep 17 00:00:00 2001 From: Alam Sanchez Date: Sun, 30 Jun 2024 23:45:52 +0200 Subject: [PATCH 1/6] Adapt EigenDecomposition implementation --- .../runtime/matrix/data/LibCommonsMath.java | 342 ++++++++++++++++++ .../java/org/apache/sysds/test/TestUtils.java | 8 +- .../component/matrix/EigenDecompTest.java | 83 ++++- 3 files changed, 413 insertions(+), 20 deletions(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java index 8757fda8228..967bc6914c4 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java @@ -27,10 +27,12 @@ import java.util.HashMap; import java.util.Map; import java.util.Map.Entry; +import java.util.Arrays; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.util.LocalizedFormats; import org.apache.commons.math3.linear.Array2DRowRealMatrix; import org.apache.commons.math3.linear.BlockRealMatrix; import org.apache.commons.math3.linear.CholeskyDecomposition; @@ -40,11 +42,13 @@ import org.apache.commons.math3.linear.QRDecomposition; import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.linear.SingularValueDecomposition; +import org.apache.commons.math3.util.FastMath; import org.apache.sysds.runtime.DMLRuntimeException; import org.apache.sysds.runtime.codegen.CodegenUtils; import org.apache.sysds.runtime.codegen.SpoofOperator.SideInput; import org.apache.sysds.runtime.compress.utils.IntArrayList; import org.apache.sysds.runtime.data.DenseBlock; +import org.apache.sysds.runtime.data.DenseBlockFactory; import org.apache.sysds.runtime.functionobjects.Builtin; import org.apache.sysds.runtime.functionobjects.Divide; import org.apache.sysds.runtime.functionobjects.MinusMultiply; @@ -61,6 +65,7 @@ import org.apache.sysds.runtime.matrix.operators.UnaryOperator; import org.apache.sysds.runtime.util.DataConverter; import org.apache.sysds.runtime.util.UtilFunctions; +import org.apache.commons.math3.util.Precision; /** * Library for matrix operations that need invocation of @@ -74,6 +79,8 @@ public class LibCommonsMath private static final Log LOG = LogFactory.getLog(LibCommonsMath.class.getName()); private static final double RELATIVE_SYMMETRY_THRESHOLD = 1e-6; private static final double EIGEN_LAMBDA = 1e-8; + private static final double EIGEN_EPS = Precision.EPSILON; + private static final int EIGEN_MAX_ITER = 30; private LibCommonsMath() { //prevent instantiation via private constructor @@ -146,6 +153,8 @@ else if (opcode.equals("eigen_lanczos")) return computeEigenLanczos(in, threads, seed); else if (opcode.equals("eigen_qr")) return computeEigenQR(in, threads); + else if (opcode.equals("eigen_symm")) + return computeEigenDecompositionSymm(in, threads); else if (opcode.equals("svd")) return computeSvd(in); else if (opcode.equals("fft")) @@ -209,6 +218,339 @@ private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) { return DataConverter.convertToMatrixBlock(solutionMatrix); } + + /** + * Computes the eigen decomposition of a symmetric matrix. + * + * @param in The input matrix to compute the eigen decomposition on. + * @param threads The number of threads to use for computation. + * @return An array of MatrixBlock objects containing the real eigen values and eigen vectors. + */ + public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock in, int threads) { + + // TODO: Verify matrix is symmetric + final double[] mainDiag = new double[in.rlen]; + final double[] secDiag = new double[in.rlen - 1]; + + final MatrixBlock householderVectors = transformToTridiagonal(in, mainDiag, secDiag, threads); + + // TODO: Consider using sparse blocks + final double[] hv = householderVectors.getDenseBlockValues(); + MatrixBlock houseHolderMatrix = getQ(hv, mainDiag, secDiag); + + MatrixBlock[] evResult = findEigenVectors(mainDiag, secDiag, houseHolderMatrix, EIGEN_MAX_ITER, threads); + + MatrixBlock realEigenValues = evResult[0]; + MatrixBlock eigenVectors = evResult[1]; + + realEigenValues.setNonZeros(realEigenValues.rlen * realEigenValues.rlen); + eigenVectors.setNonZeros(eigenVectors.rlen * eigenVectors.clen); + + return new MatrixBlock[] { realEigenValues, eigenVectors }; + } + + private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] main, double[] secondary, int threads) { + final int m = matrix.rlen; + + // MatrixBlock householderVectors = ; + MatrixBlock householderVectors = matrix.extractTriangular(new MatrixBlock(m, m, false), false, true, true); + if (householderVectors.isInSparseFormat()) { + householderVectors.sparseToDense(threads); + } + + final double[] z = new double[m]; + + // TODO: Consider sparse block case + final double[] hv = householderVectors.getDenseBlockValues(); + + for (int k = 0; k < m - 1; k++) { + final int k_kp1 = k * m + k + 1; + final int km = k * m; + + // zero-out a row and a column simultaneously + main[k] = hv[km + k]; + + // TODO: may or may not improve, seems it does not + // double xNormSqr = householderVectors.slice(k, k, k + 1, m - 1).sumSq(); + // double xNormSqr = LibMatrixMult.dotProduct(hv, hv, k_kp1, k_kp1, m - (k + 1)); + double xNormSqr = 0; + for (int j = k + 1; j < m; ++j) { + final double c = hv[k * m + j]; + xNormSqr += c * c; + } + + final double a = (hv[k_kp1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); + + secondary[k] = a; + + if (a != 0.0) { + // apply Householder transform from left and right simultaneously + hv[k_kp1] -= a; + + final double beta = -1 / (a * hv[k_kp1]); + double gamma = 0; + + // compute a = beta A v, where v is the Householder vector + // this loop is written in such a way + // 1) only the upper triangular part of the matrix is accessed + // 2) access is cache-friendly for a matrix stored in rows + Arrays.fill(z, k + 1, m, 0); + for (int i = k + 1; i < m; ++i) { + final double hKI = hv[km + i]; + double zI = hv[i * m + i] * hKI; + for (int j = i + 1; j < m; ++j) { + final double hIJ = hv[i * m + j]; + zI += hIJ * hv[k * m + j]; + z[j] += hIJ * hKI; + } + z[i] = beta * (z[i] + zI); + + gamma += z[i] * hv[km + i]; + } + + gamma *= beta / 2; + + // compute z = z - gamma v + // LibMatrixMult.vectMultiplyAdd(-gamma, hv, z, k_kp1, k + 1, m - (k + 1)); + for (int i = k + 1; i < m; ++i) { + z[i] -= gamma * hv[km + i]; + } + + // update matrix: A = A - v zT - z vT + // only the upper triangular part of the matrix is updated + for (int i = k + 1; i < m; ++i) { + final double hki = hv[km + i]; + for (int j = i; j < m; ++j) { + final double hkj = hv[km + j]; + + hv[i * m + j] -= (hki * z[j] + z[i] * hkj); + } + } + } + } + + main[m - 1] = hv[(m - 1) * m + m - 1]; + + return householderVectors; + } + + + /** + * Computes the orthogonal matrix Q using Householder transforms. + * The matrix Q is built by applying Householder transforms to the input vectors. + * + * @param hv The input vector containing the Householder vectors. + * @param main The main diagonal of the matrix. + * @param secondary The secondary diagonal of the matrix. + * @return The orthogonal matrix Q. + */ + public static MatrixBlock getQ(final double[] hv, double[] main, double[] secondary) { + final int m = main.length; + + // orthogonal matrix, most operations are column major (TODO) + DenseBlock qaB = DenseBlockFactory.createDenseBlock(m, m); + double[] qaV = qaB.valuesAt(0); + + // build up first part of the matrix by applying Householder transforms + for (int k = m - 1; k >= 1; --k) { + final int km = k * m; + final int km1m = (k - 1) * m; + + qaV[km + k] = 1.0; + if (hv[km1m + k] != 0.0) { + final double inv = 1.0 / (secondary[k - 1] * hv[km1m + k]); + + double beta = 1.0 / secondary[k - 1]; + + qaV[km + k] = 1 + beta * hv[km1m + k]; + + // TODO: may speedup vector operations + for (int i = k + 1; i < m; ++i) { + qaV[i * m + k] = beta * hv[km1m + i]; + } + + for (int j = k + 1; j < m; ++j) { + beta = 0; + for (int i = k + 1; i < m; ++i) { + beta += qaV[m * i + j] * hv[km1m + i]; + } + beta *= inv; + qaV[m * k + j] = hv[km1m + k] * beta; + + for (int i = k + 1; i < m; ++i) { + qaV[m * i + j] += beta * hv[km1m + i]; + } + } + } + } + + qaV[0] = 1.0; + MatrixBlock res = new MatrixBlock(m, m, qaB); + + return res; + } + + + /** + * Finds the eigen vectors corresponding to the given eigen values using the Householder transformation. + * (Dubrulle et al., 1971). + * + * @param main The main diagonal of the tridiagonal matrix. + * @param secondary The secondary diagonal of the tridiagonal matrix. + * @param hhMatrix The Householder matrix. + * @param maxIter The maximum number of iterations for convergence. + * @param threads The number of threads to use for computation. + * @return An array of two MatrixBlock objects: eigen values and eigen vectors. + * @throws MaxCountExceededException If the maximum number of iterations is exceeded and convergence fails. + */ + private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, final MatrixBlock hhMatrix, + final int maxIter, int threads) { + + DenseBlock hhDense = hhMatrix.getDenseBlock(); + + final int n = hhMatrix.rlen; + MatrixBlock eigenValues = new MatrixBlock(n, 1, main); + + double[] ev = eigenValues.denseBlock.valuesAt(0); + final double[] e = new double[n]; + + System.arraycopy(secondary, 0, e, 0, n - 1); + e[n - 1] = 0; + + // Determine the largest main and secondary value in absolute term. + double maxAbsoluteValue = 0; + for (int i = 0; i < n; i++) { + maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i])); + maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(e[i])); + } + // Make null any main and secondary value too small to be significant + if (maxAbsoluteValue != 0) { + for (int i = 0; i < n; i++) { + if (FastMath.abs(ev[i]) <= EIGEN_EPS * maxAbsoluteValue && ev[i] != 0.0) { + ev[i] = 0; + } + if (FastMath.abs(e[i]) <= EIGEN_EPS * maxAbsoluteValue && e[i] != 0.0) { + e[i] = 0; + } + } + } + + for (int j = 0; j < n; j++) { + int its = 0; + int m; + do { + for (m = j; m < n - 1; m++) { + final double delta = FastMath.abs(ev[m]) + + FastMath.abs(ev[m + 1]); + if (FastMath.abs(e[m]) + delta == delta) { + break; + } + } + if (m != j) { + if (its == maxIter) { + throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter); + } + + its++; + double q = (ev[j + 1] - ev[j]) / (2 * e[j]); + double t = FastMath.sqrt(1 + q * q); + if (q < 0.0) { + q = ev[m] - ev[j] + e[j] / (q - t); + } else { + q = ev[m] - ev[j] + e[j] / (q + t); + } + double u = 0.0; + double s = 1.0; + double c = 1.0; + int i; + for (i = m - 1; i >= j; i--) { + double p = s * e[i]; + double h = c * e[i]; + if (FastMath.abs(p) >= FastMath.abs(q)) { + c = q / p; + t = FastMath.sqrt(c * c + 1.0); + e[i + 1] = p * t; + s = 1.0 / t; + c *= s; + } else { + s = p / q; + t = FastMath.sqrt(s * s + 1.0); + e[i + 1] = q * t; + c = 1.0 / t; + s *= c; + } + if (e[i + 1] == 0.0) { + ev[i + 1] -= u; + e[m] = 0.0; + break; + } + q = ev[i + 1] - u; + t = (ev[i] - q) * s + 2.0 * c * h; + u = s * t; + ev[i + 1] = q + u; + q = c * t - h; + + for (int ia = 0; ia < n; ++ia) { + p = hhDense.get(ia, i + 1); + hhDense.set(ia, i + 1, s * hhDense.get(ia, i) + c * p); + hhDense.set(ia, i, c * hhDense.get(ia, i) - s * p); + } + } + + if (t == 0.0 && i >= j) { + continue; + } + ev[j] -= u; + e[j] = q; + e[m] = 0.0; + + } + } while (m != j); + } + + // Sort the eigen values (and vectors) in increase order + for (int i = 0; i < n; i++) { + int k = i; + double p = ev[i]; + for (int j = i + 1; j < n; j++) { + // reversed order from original implementation + if (ev[j] < p) { + k = j; + p = ev[j]; + } + } + if (k != i) { + ev[k] = ev[i]; + ev[i] = p; + for (int j = 0; j < n; j++) { + // TODO: an operation like SwapIndex be faster? + p = hhDense.get(j, i); + hhDense.set(j, i, hhDense.get(j, k)); + hhDense.set(j, k, p); + + } + } + } + + // Determine the largest eigen value in absolute term. + maxAbsoluteValue = 0; + for (int i = 0; i < n; i++) { + maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i])); + } + // Make null any eigen value too small to be significant + if (maxAbsoluteValue != 0.0) { + for (int i = 0; i < n; i++) { + if (FastMath.abs(ev[i]) < EIGEN_EPS * maxAbsoluteValue) { + ev[i] = 0; + } + } + } + + // MatrixBlock realEigenValues = new MatrixBlock(z.rlen, 1, ev); + + return new MatrixBlock[] { eigenValues, hhMatrix }; + } + /** * Function to perform QR decomposition on a given matrix. diff --git a/src/test/java/org/apache/sysds/test/TestUtils.java b/src/test/java/org/apache/sysds/test/TestUtils.java index 0c407cc71a4..c96138e50df 100644 --- a/src/test/java/org/apache/sysds/test/TestUtils.java +++ b/src/test/java/org/apache/sysds/test/TestUtils.java @@ -2084,10 +2084,10 @@ public static MatrixBlock generateTestMatrixBlock(int rows, int cols, double min * @param seed seed * @return random symmetric MatrixBlock */ - public static MatrixBlock generateTestMatrixBlockSym(int rows, int cols, double min, double max, double sparsity, long seed){ - MatrixBlock m = MatrixBlock.randOperations(rows, cols, sparsity, min, max, "Uniform", seed); - for(int i = 0; i < rows; i++) { - for(int j = i+1; j < cols; j++) { + public static MatrixBlock generateTestMatrixBlockSym(int size, double min, double max, double sparsity, long seed){ + MatrixBlock m = MatrixBlock.randOperations(size, size, sparsity, min, max, "Uniform", seed); + for(int i = 0; i < size; i++) { + for(int j = i+1; j < size; j++) { m.set(i,j, m.get(j,i)); } } diff --git a/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java b/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java index 6292a141389..8da74a4477a 100644 --- a/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java +++ b/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java @@ -37,31 +37,34 @@ public class EigenDecompTest { protected static final Log LOG = LogFactory.getLog(EigenDecompTest.class.getName()); private enum type { - COMMONS, LANCZOS, QR, + COMMONS, LANCZOS, QR, COMMONS_NATIVE } @Test + @Ignore public void testLanczosSimple() { - MatrixBlock in = new MatrixBlock(4, 4, new double[] {4, 1, -2, 2, 1, 2, 0, 1, -2, 0, 3, -2, 2, 1, -2, -1}); + MatrixBlock in = new MatrixBlock(4, 4, new double[] { 4, 1, -2, 2, 1, 2, 0, 1, -2, 0, 3, -2, 2, 1, -2, -1 }); testEigen(in, 1e-4, 1, type.LANCZOS); } @Test + @Ignore public void testLanczosSmall() { - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-4, 1, type.LANCZOS); } @Test + @Ignore public void testLanczosMedium() { - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(12, 12, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(12, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-4, 1, type.LANCZOS); } @Test @Ignore public void testLanczosLarge() { - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 100, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-4, 1, type.LANCZOS); } @@ -69,34 +72,75 @@ public void testLanczosLarge() { @Ignore public void testLanczosLargeMT() { int threads = 10; - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 100, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-4, threads, type.LANCZOS); } @Test + @Ignore public void testQREigenSimple() { MatrixBlock in = new MatrixBlock(4, 4, - new double[] {52, 30, 49, 28, 30, 50, 8, 44, 49, 8, 46, 16, 28, 44, 16, 22}); + new double[] { 52, 30, 49, 28, 30, 50, 8, 44, 49, 8, 46, 16, 28, 44, 16, 22 }); testEigen(in, 1e-4, 1, type.QR); } @Test + @Ignore public void testQREigenSymSmall() { - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-3, 1, type.QR); } @Test + @Ignore public void testQREigenSymSmallMT() { int threads = 10; - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-3, threads, type.QR); } + @Test + public void testEigenSymSmall() { + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1); + testEigen(in, 1e-4, 10, type.COMMONS); + } + + @Test + public void testEigenSymMedium() { + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(200, 0.0, 1.0, 1.0, 1); + testEigen(in, 1e-4, 10, type.COMMONS); + } + + @Test + public void testEigenSymLarge() { + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(1000, 0.0, 1.0, 1.0, 1); + testEigen(in, 1e-4, 10, type.COMMONS); + } + + @Test + public void testNativeEigenSymSmall() { + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1); + testEigen(in, 1e-4, 10, type.COMMONS_NATIVE); + } + + @Test + public void testNativeEigenSymMedium() { + int threads = 10; + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(200, 0.0, 1.0, 1.0, 1); + testEigen(in, 1e-4, threads, type.COMMONS_NATIVE); + } + + @Test + public void testNativeEigenSymLarge() { + int threads = 10; + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(1000, 0.0, 1.0, 1.0, 1); + testEigen(in, 1e-4, threads, type.COMMONS_NATIVE); + } + @Test @Ignore public void testQREigenSymLarge() { - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 50, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-4, 1, type.QR); } @@ -104,14 +148,14 @@ public void testQREigenSymLarge() { @Ignore public void testQREigenSymLargeMT() { int threads = 10; - MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 50, 0.0, 1.0, 1.0, 1); + MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 0.0, 1.0, 1.0, 1); testEigen(in, 1e-4, threads, type.QR); } - private void testEigen(MatrixBlock in, double tol, int threads, type t) { + private void testEigen(MatrixBlock in, double tol, int threads, type t, boolean validate) { try { MatrixBlock[] m = null; - switch(t) { + switch (t) { case COMMONS: m = LibCommonsMath.multiReturnOperations(in, "eigen", threads, 1); break; @@ -121,19 +165,26 @@ private void testEigen(MatrixBlock in, double tol, int threads, type t) { case QR: m = LibCommonsMath.multiReturnOperations(in, "eigen_qr", threads, 1); break; + case COMMONS_NATIVE: + m = LibCommonsMath.multiReturnOperations(in, "eigen_symm", threads, 1); + break; default: throw new NotImplementedException(); } - isValidDecomposition(in, m[1], m[0], tol, t.toString()); + if (validate) + isValidDecomposition(in, m[1], m[0], tol, t.toString()); - } - catch(Exception e) { + } catch (Exception e) { e.printStackTrace(); fail(e.getMessage()); } } + private void testEigen(MatrixBlock in, double tol, int threads, type t) { + testEigen(in, tol, threads, t, true); + } + private void isValidDecomposition(MatrixBlock A, MatrixBlock V, MatrixBlock vD, double tol, String message) { // Any eigen decomposition is valid if A = V %*% D %*% t(V) // A is the input of the eigen decomposition From ff1b77deda7d28e667668acb758843c20f474ac6 Mon Sep 17 00:00:00 2001 From: Alam Sanchez Date: Wed, 3 Jul 2024 17:58:02 +0200 Subject: [PATCH 2/6] Apply formatting --- .../runtime/matrix/data/LibCommonsMath.java | 354 +++++++++--------- 1 file changed, 178 insertions(+), 176 deletions(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java index 967bc6914c4..921c0024e14 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java @@ -221,13 +221,13 @@ private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) { /** * Computes the eigen decomposition of a symmetric matrix. - * - * @param in The input matrix to compute the eigen decomposition on. - * @param threads The number of threads to use for computation. - * @return An array of MatrixBlock objects containing the real eigen values and eigen vectors. - */ + * + * @param in The input matrix to compute the eigen decomposition on. + * @param threads The number of threads to use for computation. + * @return An array of MatrixBlock objects containing the real eigen values and eigen vectors. + */ public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock in, int threads) { - + // TODO: Verify matrix is symmetric final double[] mainDiag = new double[in.rlen]; final double[] secDiag = new double[in.rlen - 1]; @@ -243,19 +243,20 @@ public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock in, int th MatrixBlock realEigenValues = evResult[0]; MatrixBlock eigenVectors = evResult[1]; - realEigenValues.setNonZeros(realEigenValues.rlen * realEigenValues.rlen); - eigenVectors.setNonZeros(eigenVectors.rlen * eigenVectors.clen); + // TODO: Count number of zeros on construction + realEigenValues.setNonZeros(realEigenValues.denseBlock.countNonZeros()); - return new MatrixBlock[] { realEigenValues, eigenVectors }; + return new MatrixBlock[] {realEigenValues, eigenVectors}; } - private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] main, double[] secondary, int threads) { + private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] main, double[] secondary, + int threads) { final int m = matrix.rlen; // MatrixBlock householderVectors = ; MatrixBlock householderVectors = matrix.extractTriangular(new MatrixBlock(m, m, false), false, true, true); - if (householderVectors.isInSparseFormat()) { - householderVectors.sparseToDense(threads); + if(householderVectors.isInSparseFormat()) { + householderVectors.sparseToDense(threads); } final double[] z = new double[m]; @@ -263,7 +264,7 @@ private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] m // TODO: Consider sparse block case final double[] hv = householderVectors.getDenseBlockValues(); - for (int k = 0; k < m - 1; k++) { + for(int k = 0; k < m - 1; k++) { final int k_kp1 = k * m + k + 1; final int km = k * m; @@ -274,59 +275,59 @@ private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] m // double xNormSqr = householderVectors.slice(k, k, k + 1, m - 1).sumSq(); // double xNormSqr = LibMatrixMult.dotProduct(hv, hv, k_kp1, k_kp1, m - (k + 1)); double xNormSqr = 0; - for (int j = k + 1; j < m; ++j) { - final double c = hv[k * m + j]; - xNormSqr += c * c; + for(int j = k + 1; j < m; ++j) { + final double c = hv[k * m + j]; + xNormSqr += c * c; } final double a = (hv[k_kp1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); secondary[k] = a; - if (a != 0.0) { - // apply Householder transform from left and right simultaneously - hv[k_kp1] -= a; - - final double beta = -1 / (a * hv[k_kp1]); - double gamma = 0; - - // compute a = beta A v, where v is the Householder vector - // this loop is written in such a way - // 1) only the upper triangular part of the matrix is accessed - // 2) access is cache-friendly for a matrix stored in rows - Arrays.fill(z, k + 1, m, 0); - for (int i = k + 1; i < m; ++i) { - final double hKI = hv[km + i]; - double zI = hv[i * m + i] * hKI; - for (int j = i + 1; j < m; ++j) { - final double hIJ = hv[i * m + j]; - zI += hIJ * hv[k * m + j]; - z[j] += hIJ * hKI; - } - z[i] = beta * (z[i] + zI); + if(a != 0.0) { + // apply Householder transform from left and right simultaneously + hv[k_kp1] -= a; + + final double beta = -1 / (a * hv[k_kp1]); + double gamma = 0; + + // compute a = beta A v, where v is the Householder vector + // this loop is written in such a way + // 1) only the upper triangular part of the matrix is accessed + // 2) access is cache-friendly for a matrix stored in rows + Arrays.fill(z, k + 1, m, 0); + for(int i = k + 1; i < m; ++i) { + final double hKI = hv[km + i]; + double zI = hv[i * m + i] * hKI; + for(int j = i + 1; j < m; ++j) { + final double hIJ = hv[i * m + j]; + zI += hIJ * hv[k * m + j]; + z[j] += hIJ * hKI; + } + z[i] = beta * (z[i] + zI); - gamma += z[i] * hv[km + i]; - } + gamma += z[i] * hv[km + i]; + } - gamma *= beta / 2; + gamma *= beta / 2; - // compute z = z - gamma v - // LibMatrixMult.vectMultiplyAdd(-gamma, hv, z, k_kp1, k + 1, m - (k + 1)); - for (int i = k + 1; i < m; ++i) { - z[i] -= gamma * hv[km + i]; - } + // compute z = z - gamma v + // LibMatrixMult.vectMultiplyAdd(-gamma, hv, z, k_kp1, k + 1, m - (k + 1)); + for(int i = k + 1; i < m; ++i) { + z[i] -= gamma * hv[km + i]; + } - // update matrix: A = A - v zT - z vT - // only the upper triangular part of the matrix is updated - for (int i = k + 1; i < m; ++i) { - final double hki = hv[km + i]; - for (int j = i; j < m; ++j) { - final double hkj = hv[km + j]; + // update matrix: A = A - v zT - z vT + // only the upper triangular part of the matrix is updated + for(int i = k + 1; i < m; ++i) { + final double hki = hv[km + i]; + for(int j = i; j < m; ++j) { + final double hkj = hv[km + j]; - hv[i * m + j] -= (hki * z[j] + z[i] * hkj); + hv[i * m + j] -= (hki * z[j] + z[i] * hkj); + } } } - } } main[m - 1] = hv[(m - 1) * m + m - 1]; @@ -334,10 +335,9 @@ private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] m return householderVectors; } - /** - * Computes the orthogonal matrix Q using Householder transforms. - * The matrix Q is built by applying Householder transforms to the input vectors. + * Computes the orthogonal matrix Q using Householder transforms. The matrix Q is built by applying Householder + * transforms to the input vectors. * * @param hv The input vector containing the Householder vectors. * @param main The main diagonal of the matrix. @@ -352,36 +352,36 @@ public static MatrixBlock getQ(final double[] hv, double[] main, double[] second double[] qaV = qaB.valuesAt(0); // build up first part of the matrix by applying Householder transforms - for (int k = m - 1; k >= 1; --k) { + for(int k = m - 1; k >= 1; --k) { final int km = k * m; final int km1m = (k - 1) * m; qaV[km + k] = 1.0; - if (hv[km1m + k] != 0.0) { - final double inv = 1.0 / (secondary[k - 1] * hv[km1m + k]); - - double beta = 1.0 / secondary[k - 1]; + if(hv[km1m + k] != 0.0) { + final double inv = 1.0 / (secondary[k - 1] * hv[km1m + k]); - qaV[km + k] = 1 + beta * hv[km1m + k]; + double beta = 1.0 / secondary[k - 1]; - // TODO: may speedup vector operations - for (int i = k + 1; i < m; ++i) { - qaV[i * m + k] = beta * hv[km1m + i]; - } + qaV[km + k] = 1 + beta * hv[km1m + k]; - for (int j = k + 1; j < m; ++j) { - beta = 0; - for (int i = k + 1; i < m; ++i) { - beta += qaV[m * i + j] * hv[km1m + i]; + // TODO: may speedup vector operations + for(int i = k + 1; i < m; ++i) { + qaV[i * m + k] = beta * hv[km1m + i]; } - beta *= inv; - qaV[m * k + j] = hv[km1m + k] * beta; - for (int i = k + 1; i < m; ++i) { - qaV[m * i + j] += beta * hv[km1m + i]; + for(int j = k + 1; j < m; ++j) { + beta = 0; + for(int i = k + 1; i < m; ++i) { + beta += qaV[m * i + j] * hv[km1m + i]; + } + beta *= inv; + qaV[m * k + j] = hv[km1m + k] * beta; + + for(int i = k + 1; i < m; ++i) { + qaV[m * i + j] += beta * hv[km1m + i]; + } } } - } } qaV[0] = 1.0; @@ -390,21 +390,20 @@ public static MatrixBlock getQ(final double[] hv, double[] main, double[] second return res; } - /** - * Finds the eigen vectors corresponding to the given eigen values using the Householder transformation. - * (Dubrulle et al., 1971). + * Finds the eigen vectors corresponding to the given eigen values using the Householder transformation. (Dubrulle + * et al., 1971). * * @param main The main diagonal of the tridiagonal matrix. * @param secondary The secondary diagonal of the tridiagonal matrix. * @param hhMatrix The Householder matrix. * @param maxIter The maximum number of iterations for convergence. * @param threads The number of threads to use for computation. - * @return An array of two MatrixBlock objects: eigen values and eigen vectors. + * @return An array of two MatrixBlock objects: eigen values and eigen vectors. * @throws MaxCountExceededException If the maximum number of iterations is exceeded and convergence fails. */ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, final MatrixBlock hhMatrix, - final int maxIter, int threads) { + final int maxIter, int threads) { DenseBlock hhDense = hhMatrix.getDenseBlock(); @@ -419,138 +418,141 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, // Determine the largest main and secondary value in absolute term. double maxAbsoluteValue = 0; - for (int i = 0; i < n; i++) { + for(int i = 0; i < n; i++) { maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i])); maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(e[i])); } // Make null any main and secondary value too small to be significant - if (maxAbsoluteValue != 0) { - for (int i = 0; i < n; i++) { - if (FastMath.abs(ev[i]) <= EIGEN_EPS * maxAbsoluteValue && ev[i] != 0.0) { - ev[i] = 0; - } - if (FastMath.abs(e[i]) <= EIGEN_EPS * maxAbsoluteValue && e[i] != 0.0) { - e[i] = 0; - } + if(maxAbsoluteValue != 0) { + for(int i = 0; i < n; i++) { + if(FastMath.abs(ev[i]) <= EIGEN_EPS * maxAbsoluteValue && ev[i] != 0.0) { + ev[i] = 0; + } + if(FastMath.abs(e[i]) <= EIGEN_EPS * maxAbsoluteValue && e[i] != 0.0) { + e[i] = 0; + } } } - for (int j = 0; j < n; j++) { + for(int j = 0; j < n; j++) { int its = 0; int m; do { - for (m = j; m < n - 1; m++) { - final double delta = FastMath.abs(ev[m]) + - FastMath.abs(ev[m + 1]); - if (FastMath.abs(e[m]) + delta == delta) { - break; - } - } - if (m != j) { - if (its == maxIter) { - throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter); + for(m = j; m < n - 1; m++) { + final double delta = FastMath.abs(ev[m]) + FastMath.abs(ev[m + 1]); + if(FastMath.abs(e[m]) + delta == delta) { + break; + } } + if(m != j) { + if(its == maxIter) { + throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter); + } - its++; - double q = (ev[j + 1] - ev[j]) / (2 * e[j]); - double t = FastMath.sqrt(1 + q * q); - if (q < 0.0) { - q = ev[m] - ev[j] + e[j] / (q - t); - } else { - q = ev[m] - ev[j] + e[j] / (q + t); - } - double u = 0.0; - double s = 1.0; - double c = 1.0; - int i; - for (i = m - 1; i >= j; i--) { - double p = s * e[i]; - double h = c * e[i]; - if (FastMath.abs(p) >= FastMath.abs(q)) { - c = q / p; - t = FastMath.sqrt(c * c + 1.0); - e[i + 1] = p * t; - s = 1.0 / t; - c *= s; - } else { - s = p / q; - t = FastMath.sqrt(s * s + 1.0); - e[i + 1] = q * t; - c = 1.0 / t; - s *= c; - } - if (e[i + 1] == 0.0) { - ev[i + 1] -= u; + its++; + double q = (ev[j + 1] - ev[j]) / (2 * e[j]); + double t = FastMath.sqrt(1 + q * q); + if(q < 0.0) { + q = ev[m] - ev[j] + e[j] / (q - t); + } + else { + q = ev[m] - ev[j] + e[j] / (q + t); + } + double u = 0.0; + double s = 1.0; + double c = 1.0; + int i; + for(i = m - 1; i >= j; i--) { + double p = s * e[i]; + double h = c * e[i]; + if(FastMath.abs(p) >= FastMath.abs(q)) { + c = q / p; + t = FastMath.sqrt(c * c + 1.0); + e[i + 1] = p * t; + s = 1.0 / t; + c *= s; + } + else { + s = p / q; + t = FastMath.sqrt(s * s + 1.0); + e[i + 1] = q * t; + c = 1.0 / t; + s *= c; + } + if(e[i + 1] == 0.0) { + ev[i + 1] -= u; + e[m] = 0.0; + break; + } + q = ev[i + 1] - u; + t = (ev[i] - q) * s + 2.0 * c * h; + u = s * t; + ev[i + 1] = q + u; + q = c * t - h; + + for(int ia = 0; ia < n; ++ia) { + p = hhDense.get(ia, i + 1); + hhDense.set(ia, i + 1, s * hhDense.get(ia, i) + c * p); + hhDense.set(ia, i, c * hhDense.get(ia, i) - s * p); + } + } + + if(t == 0.0 && i >= j) { + continue; + } + ev[j] -= u; + e[j] = q; e[m] = 0.0; - break; - } - q = ev[i + 1] - u; - t = (ev[i] - q) * s + 2.0 * c * h; - u = s * t; - ev[i + 1] = q + u; - q = c * t - h; - - for (int ia = 0; ia < n; ++ia) { - p = hhDense.get(ia, i + 1); - hhDense.set(ia, i + 1, s * hhDense.get(ia, i) + c * p); - hhDense.set(ia, i, c * hhDense.get(ia, i) - s * p); - } - } - if (t == 0.0 && i >= j) { - continue; } - ev[j] -= u; - e[j] = q; - e[m] = 0.0; - } - } while (m != j); + while(m != j); } // Sort the eigen values (and vectors) in increase order - for (int i = 0; i < n; i++) { + for(int i = 0; i < n; i++) { int k = i; double p = ev[i]; - for (int j = i + 1; j < n; j++) { - // reversed order from original implementation - if (ev[j] < p) { - k = j; - p = ev[j]; - } + for(int j = i + 1; j < n; j++) { + // reversed order from original implementation + if(ev[j] < p) { + k = j; + p = ev[j]; + } } - if (k != i) { - ev[k] = ev[i]; - ev[i] = p; - for (int j = 0; j < n; j++) { - // TODO: an operation like SwapIndex be faster? - p = hhDense.get(j, i); - hhDense.set(j, i, hhDense.get(j, k)); - hhDense.set(j, k, p); + if(k != i) { + ev[k] = ev[i]; + ev[i] = p; + for(int j = 0; j < n; j++) { + // TODO: an operation like SwapIndex be faster? + p = hhDense.get(j, i); + hhDense.set(j, i, hhDense.get(j, k)); + hhDense.set(j, k, p); - } + } } } // Determine the largest eigen value in absolute term. maxAbsoluteValue = 0; - for (int i = 0; i < n; i++) { + for(int i = 0; i < n; i++) { maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i])); } // Make null any eigen value too small to be significant - if (maxAbsoluteValue != 0.0) { - for (int i = 0; i < n; i++) { - if (FastMath.abs(ev[i]) < EIGEN_EPS * maxAbsoluteValue) { - ev[i] = 0; - } + int zeros = 0; + if(maxAbsoluteValue != 0.0) { + for(int i = 0; i < n; i++) { + if(FastMath.abs(ev[i]) < EIGEN_EPS * maxAbsoluteValue) { + ev[i] = 0; + zeros++; + } } } - // MatrixBlock realEigenValues = new MatrixBlock(z.rlen, 1, ev); + eigenValues.setNonZeros(n - zeros); - return new MatrixBlock[] { eigenValues, hhMatrix }; + return new MatrixBlock[] {eigenValues, hhMatrix}; } - /** * Function to perform QR decomposition on a given matrix. From 617b7c90053e238a24ed49bad0383d5cd38eb6ec Mon Sep 17 00:00:00 2001 From: Alam Sanchez Date: Wed, 3 Jul 2024 19:10:13 +0200 Subject: [PATCH 3/6] Fix tests --- .../org/apache/sysds/runtime/matrix/data/LibCommonsMath.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java index 921c0024e14..0a5e0d721be 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java @@ -244,7 +244,7 @@ public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock in, int th MatrixBlock eigenVectors = evResult[1]; // TODO: Count number of zeros on construction - realEigenValues.setNonZeros(realEigenValues.denseBlock.countNonZeros()); + eigenVectors.setNonZeros(eigenVectors.denseBlock.countNonZeros()); return new MatrixBlock[] {realEigenValues, eigenVectors}; } From 36525ba9ae9d6be43fd8e97218adc1fdecbd14db Mon Sep 17 00:00:00 2001 From: Alam Sanchez Date: Mon, 8 Jul 2024 11:34:28 +0200 Subject: [PATCH 4/6] Apply transposes for row major operations --- .../runtime/matrix/data/LibCommonsMath.java | 205 +++++++++--------- 1 file changed, 98 insertions(+), 107 deletions(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java index 0a5e0d721be..48402aa3663 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java @@ -27,7 +27,6 @@ import java.util.HashMap; import java.util.Map; import java.util.Map.Entry; -import java.util.Arrays; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; @@ -65,7 +64,6 @@ import org.apache.sysds.runtime.matrix.operators.UnaryOperator; import org.apache.sysds.runtime.util.DataConverter; import org.apache.sysds.runtime.util.UtilFunctions; -import org.apache.commons.math3.util.Precision; /** * Library for matrix operations that need invocation of @@ -79,7 +77,8 @@ public class LibCommonsMath private static final Log LOG = LogFactory.getLog(LibCommonsMath.class.getName()); private static final double RELATIVE_SYMMETRY_THRESHOLD = 1e-6; private static final double EIGEN_LAMBDA = 1e-8; - private static final double EIGEN_EPS = Precision.EPSILON; + // Machine epsilon 2^-53 + private static final double EIGEN_MACHINE_EPS = 0x1.0p-53; private static final int EIGEN_MAX_ITER = 30; private LibCommonsMath() { @@ -148,7 +147,8 @@ else if (opcode.equals("qr2")) else if (opcode.equals("lu")) return computeLU(in); else if (opcode.equals("eigen")) - return computeEigen(in); + return computeEigenDecompositionSymm(in, threads); + // return computeEigen(in); else if (opcode.equals("eigen_lanczos")) return computeEigenLanczos(in, threads, seed); else if (opcode.equals("eigen_qr")) @@ -243,8 +243,7 @@ public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock in, int th MatrixBlock realEigenValues = evResult[0]; MatrixBlock eigenVectors = evResult[1]; - // TODO: Count number of zeros on construction - eigenVectors.setNonZeros(eigenVectors.denseBlock.countNonZeros()); + eigenVectors = LibMatrixReorg.transposeInPlace(eigenVectors, threads); return new MatrixBlock[] {realEigenValues, eigenVectors}; } @@ -253,51 +252,44 @@ private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] m int threads) { final int m = matrix.rlen; - // MatrixBlock householderVectors = ; MatrixBlock householderVectors = matrix.extractTriangular(new MatrixBlock(m, m, false), false, true, true); if(householderVectors.isInSparseFormat()) { householderVectors.sparseToDense(threads); } - final double[] z = new double[m]; - - // TODO: Consider sparse block case final double[] hv = householderVectors.getDenseBlockValues(); for(int k = 0; k < m - 1; k++) { - final int k_kp1 = k * m + k + 1; - final int km = k * m; + final int rowKp1 = k * m + k + 1; + final int rowK = k * m; // zero-out a row and a column simultaneously - main[k] = hv[km + k]; + main[k] = hv[rowK + k]; - // TODO: may or may not improve, seems it does not - // double xNormSqr = householderVectors.slice(k, k, k + 1, m - 1).sumSq(); - // double xNormSqr = LibMatrixMult.dotProduct(hv, hv, k_kp1, k_kp1, m - (k + 1)); double xNormSqr = 0; for(int j = k + 1; j < m; ++j) { final double c = hv[k * m + j]; xNormSqr += c * c; } - final double a = (hv[k_kp1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); + final double a = (hv[rowKp1] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); secondary[k] = a; if(a != 0.0) { // apply Householder transform from left and right simultaneously - hv[k_kp1] -= a; + hv[rowKp1] -= a; - final double beta = -1 / (a * hv[k_kp1]); + final double beta = -1 / (a * hv[rowKp1]); double gamma = 0; // compute a = beta A v, where v is the Householder vector // this loop is written in such a way // 1) only the upper triangular part of the matrix is accessed // 2) access is cache-friendly for a matrix stored in rows - Arrays.fill(z, k + 1, m, 0); + final double[] z = new double[m]; for(int i = k + 1; i < m; ++i) { - final double hKI = hv[km + i]; + final double hKI = hv[rowK + i]; double zI = hv[i * m + i] * hKI; for(int j = i + 1; j < m; ++j) { final double hIJ = hv[i * m + j]; @@ -306,24 +298,22 @@ private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] m } z[i] = beta * (z[i] + zI); - gamma += z[i] * hv[km + i]; + gamma += z[i] * hv[rowK + i]; } gamma *= beta / 2; // compute z = z - gamma v - // LibMatrixMult.vectMultiplyAdd(-gamma, hv, z, k_kp1, k + 1, m - (k + 1)); for(int i = k + 1; i < m; ++i) { - z[i] -= gamma * hv[km + i]; + z[i] -= gamma * hv[rowK + i]; } // update matrix: A = A - v zT - z vT // only the upper triangular part of the matrix is updated for(int i = k + 1; i < m; ++i) { - final double hki = hv[km + i]; + final double hki = hv[rowK + i]; for(int j = i; j < m; ++j) { - final double hkj = hv[km + j]; - + final double hkj = hv[rowK + j]; hv[i * m + j] -= (hki * z[j] + z[i] * hkj); } } @@ -347,38 +337,39 @@ private static MatrixBlock transformToTridiagonal(MatrixBlock matrix, double[] m public static MatrixBlock getQ(final double[] hv, double[] main, double[] secondary) { final int m = main.length; - // orthogonal matrix, most operations are column major (TODO) DenseBlock qaB = DenseBlockFactory.createDenseBlock(m, m); + double[] qaV = qaB.valuesAt(0); // build up first part of the matrix by applying Householder transforms for(int k = m - 1; k >= 1; --k) { - final int km = k * m; - final int km1m = (k - 1) * m; + final int rowK = k * m; + final int rowKm1 = (k - 1) * m; - qaV[km + k] = 1.0; - if(hv[km1m + k] != 0.0) { - final double inv = 1.0 / (secondary[k - 1] * hv[km1m + k]); + qaV[rowK + k] = 1.0; + if(hv[rowKm1 + k] != 0.0) { + final double inv = 1.0 / (secondary[k - 1] * hv[rowKm1 + k]); double beta = 1.0 / secondary[k - 1]; - qaV[km + k] = 1 + beta * hv[km1m + k]; + qaV[rowK + k] = 1 + beta * hv[rowKm1 + k]; - // TODO: may speedup vector operations for(int i = k + 1; i < m; ++i) { - qaV[i * m + k] = beta * hv[km1m + i]; + qaV[rowK + i] = beta * hv[rowKm1 + i]; } for(int j = k + 1; j < m; ++j) { + beta = 0; for(int i = k + 1; i < m; ++i) { - beta += qaV[m * i + j] * hv[km1m + i]; + beta += qaV[m * j + i] * hv[rowKm1 + i]; } beta *= inv; - qaV[m * k + j] = hv[km1m + k] * beta; + + qaV[m * j + k] = hv[rowKm1 + k] * beta; for(int i = k + 1; i < m; ++i) { - qaV[m * i + j] += beta * hv[km1m + i]; + qaV[m * j + i] += beta * hv[rowKm1 + i]; } } } @@ -386,6 +377,8 @@ public static MatrixBlock getQ(final double[] hv, double[] main, double[] second qaV[0] = 1.0; MatrixBlock res = new MatrixBlock(m, m, qaB); + // Arbitrarily set non zero count, will set actual count at the end of the algorithm + res.setNonZeros(m * m); return res; } @@ -402,10 +395,10 @@ public static MatrixBlock getQ(final double[] hv, double[] main, double[] second * @return An array of two MatrixBlock objects: eigen values and eigen vectors. * @throws MaxCountExceededException If the maximum number of iterations is exceeded and convergence fails. */ - private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, final MatrixBlock hhMatrix, + private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, MatrixBlock hhMatrix, final int maxIter, int threads) { - DenseBlock hhDense = hhMatrix.getDenseBlock(); + double[] hhvalues = hhMatrix.getDenseBlock().valuesAt(0); final int n = hhMatrix.rlen; MatrixBlock eigenValues = new MatrixBlock(n, 1, main); @@ -422,13 +415,14 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i])); maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(e[i])); } + // Make null any main and secondary value too small to be significant if(maxAbsoluteValue != 0) { for(int i = 0; i < n; i++) { - if(FastMath.abs(ev[i]) <= EIGEN_EPS * maxAbsoluteValue && ev[i] != 0.0) { + if(FastMath.abs(ev[i]) <= EIGEN_MACHINE_EPS * maxAbsoluteValue && ev[i] != 0.0) { ev[i] = 0; } - if(FastMath.abs(e[i]) <= EIGEN_EPS * maxAbsoluteValue && e[i] != 0.0) { + if(FastMath.abs(e[i]) <= EIGEN_MACHINE_EPS * maxAbsoluteValue && e[i] != 0.0) { e[i] = 0; } } @@ -436,85 +430,83 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, for(int j = 0; j < n; j++) { int its = 0; - int m; - do { + int m = -1; + while(m != j) { for(m = j; m < n - 1; m++) { final double delta = FastMath.abs(ev[m]) + FastMath.abs(ev[m + 1]); if(FastMath.abs(e[m]) + delta == delta) { break; } } - if(m != j) { - if(its == maxIter) { - throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter); - } + if(m == j) + break; + if(its == maxIter) { + throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter); + } - its++; - double q = (ev[j + 1] - ev[j]) / (2 * e[j]); - double t = FastMath.sqrt(1 + q * q); - if(q < 0.0) { - q = ev[m] - ev[j] + e[j] / (q - t); + its++; + double q = (ev[j + 1] - ev[j]) / (2 * e[j]); + double t = FastMath.sqrt(1 + q * q); + if(q < 0.0) { + q = ev[m] - ev[j] + e[j] / (q - t); + } + else { + q = ev[m] - ev[j] + e[j] / (q + t); + } + double u = 0.0; + double s = 1.0; + double c = 1.0; + int i; + for(i = m - 1; i >= j; i--) { + double p = s * e[i]; + double h = c * e[i]; + if(FastMath.abs(p) >= FastMath.abs(q)) { + c = q / p; + t = FastMath.sqrt(c * c + 1.0); + e[i + 1] = p * t; + s = 1.0 / t; + c *= s; } else { - q = ev[m] - ev[j] + e[j] / (q + t); + s = p / q; + t = FastMath.sqrt(s * s + 1.0); + e[i + 1] = q * t; + c = 1.0 / t; + s *= c; } - double u = 0.0; - double s = 1.0; - double c = 1.0; - int i; - for(i = m - 1; i >= j; i--) { - double p = s * e[i]; - double h = c * e[i]; - if(FastMath.abs(p) >= FastMath.abs(q)) { - c = q / p; - t = FastMath.sqrt(c * c + 1.0); - e[i + 1] = p * t; - s = 1.0 / t; - c *= s; - } - else { - s = p / q; - t = FastMath.sqrt(s * s + 1.0); - e[i + 1] = q * t; - c = 1.0 / t; - s *= c; - } - if(e[i + 1] == 0.0) { - ev[i + 1] -= u; - e[m] = 0.0; - break; - } - q = ev[i + 1] - u; - t = (ev[i] - q) * s + 2.0 * c * h; - u = s * t; - ev[i + 1] = q + u; - q = c * t - h; - - for(int ia = 0; ia < n; ++ia) { - p = hhDense.get(ia, i + 1); - hhDense.set(ia, i + 1, s * hhDense.get(ia, i) + c * p); - hhDense.set(ia, i, c * hhDense.get(ia, i) - s * p); - } + if(e[i + 1] == 0.0) { + ev[i + 1] -= u; + e[m] = 0.0; + break; } - - if(t == 0.0 && i >= j) { - continue; + q = ev[i + 1] - u; + t = (ev[i] - q) * s + 2.0 * c * h; + u = s * t; + ev[i + 1] = q + u; + q = c * t - h; + + for(int ia = 0; ia < n; ++ia) { + p = hhvalues[(i + 1) * n + ia]; + hhvalues[(i + 1) * n + ia] = s * hhvalues[i * n + ia] + c * p; + hhvalues[i * n + ia] = c * hhvalues[i * n + ia] - s * p; } - ev[j] -= u; - e[j] = q; - e[m] = 0.0; + } + if(t == 0.0 && i >= j) { + continue; } + ev[j] -= u; + e[j] = q; + e[m] = 0.0; } - while(m != j); + } - // Sort the eigen values (and vectors) in increase order + // Sort eigen values and vectors in decreasing order for(int i = 0; i < n; i++) { int k = i; double p = ev[i]; for(int j = i + 1; j < n; j++) { - // reversed order from original implementation if(ev[j] < p) { k = j; p = ev[j]; @@ -524,11 +516,9 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, ev[k] = ev[i]; ev[i] = p; for(int j = 0; j < n; j++) { - // TODO: an operation like SwapIndex be faster? - p = hhDense.get(j, i); - hhDense.set(j, i, hhDense.get(j, k)); - hhDense.set(j, k, p); - + p = hhvalues[i * n + j]; + hhvalues[i * n + j] = hhvalues[k * n + j]; + hhvalues[k * n + j] = p; } } } @@ -542,7 +532,7 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, int zeros = 0; if(maxAbsoluteValue != 0.0) { for(int i = 0; i < n; i++) { - if(FastMath.abs(ev[i]) < EIGEN_EPS * maxAbsoluteValue) { + if(FastMath.abs(ev[i]) < EIGEN_MACHINE_EPS * maxAbsoluteValue) { ev[i] = 0; zeros++; } @@ -550,6 +540,7 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, } eigenValues.setNonZeros(n - zeros); + hhMatrix.setNonZeros(hhMatrix.denseBlock.countNonZeros()); return new MatrixBlock[] {eigenValues, hhMatrix}; } From e4b6278ff061dfd24839fa6e9088754afcdfb1fb Mon Sep 17 00:00:00 2001 From: Alam Sanchez Date: Mon, 8 Jul 2024 12:31:36 +0200 Subject: [PATCH 5/6] Refactor matrix shift step --- .../runtime/matrix/data/LibCommonsMath.java | 134 ++++++++++-------- 1 file changed, 76 insertions(+), 58 deletions(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java index 48402aa3663..0aff7b731b3 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java @@ -220,15 +220,18 @@ private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) { } /** - * Computes the eigen decomposition of a symmetric matrix. + * Computes the eigen decomposition of a symmetric matrix using the Implicit QL Algorithm (Dubrulle et al., 1971). * * @param in The input matrix to compute the eigen decomposition on. * @param threads The number of threads to use for computation. * @return An array of MatrixBlock objects containing the real eigen values and eigen vectors. */ public static MatrixBlock[] computeEigenDecompositionSymm(MatrixBlock in, int threads) { + if ( in.getNumRows() != in.getNumColumns() ) { + throw new DMLRuntimeException("Eigen Decomposition can only be done on a square and symmetric matrix. " + + "Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")"); + } - // TODO: Verify matrix is symmetric final double[] mainDiag = new double[in.rlen]; final double[] secDiag = new double[in.rlen - 1]; @@ -384,12 +387,11 @@ public static MatrixBlock getQ(final double[] hv, double[] main, double[] second } /** - * Finds the eigen vectors corresponding to the given eigen values using the Householder transformation. (Dubrulle - * et al., 1971). + * Finds the eigen vectors corresponding to the given eigen values using the Householder transformation. * * @param main The main diagonal of the tridiagonal matrix. * @param secondary The secondary diagonal of the tridiagonal matrix. - * @param hhMatrix The Householder matrix. + * @param hhMatrix The Householder matrix (Q). * @param maxIter The maximum number of iterations for convergence. * @param threads The number of threads to use for computation. * @return An array of two MatrixBlock objects: eigen values and eigen vectors. @@ -401,6 +403,7 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, double[] hhvalues = hhMatrix.getDenseBlock().valuesAt(0); final int n = hhMatrix.rlen; + MatrixBlock eigenValues = new MatrixBlock(n, 1, main); double[] ev = eigenValues.denseBlock.valuesAt(0); @@ -440,11 +443,70 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, } if(m == j) break; - if(its == maxIter) { + if(its == maxIter) throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, maxIter); - } its++; + formMatrixShift(e, ev, hhvalues, j, m); + } + + } + + // Sort eigen values and vectors in decreasing order + for(int i = 0; i < n; i++) { + int k = i; + double p = ev[i]; + for(int j = i + 1; j < n; j++) { + if(ev[j] < p) { + k = j; + p = ev[j]; + } + } + if(k != i) { + ev[k] = ev[i]; + ev[i] = p; + for(int j = 0; j < n; j++) { + p = hhvalues[i * n + j]; + hhvalues[i * n + j] = hhvalues[k * n + j]; + hhvalues[k * n + j] = p; + } + } + } + + // Determine the largest eigen value in absolute term. + maxAbsoluteValue = 0; + for(int i = 0; i < n; i++) { + maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i])); + } + // Make null any eigen value too small to be significant + int zeros = 0; + if(maxAbsoluteValue != 0.0) { + for(int i = 0; i < n; i++) { + if(FastMath.abs(ev[i]) < EIGEN_MACHINE_EPS * maxAbsoluteValue) { + ev[i] = 0; + zeros++; + } + } + } + + eigenValues.setNonZeros(n - zeros); + hhMatrix.setNonZeros(hhMatrix.denseBlock.countNonZeros()); + + return new MatrixBlock[] {eigenValues, hhMatrix}; + } + + /** + * Performs a matrix shift operation on the given arrays. Implements 'imtqll' procedure (Dubrulle et al., 1971). + * + * @param e The array to store eigenvalues. + * @param ev The array (dense block) to store eigenvectors. + * @param hhValues The array of Householder vectors. + * @param j The starting index for the matrix shift operation. + * @param m The ending index for the matrix shift operation. + */ + private static void formMatrixShift(double[] e, double[] ev, double[] hhValues, int j, int m) { + final int n = e.length; + double q = (ev[j + 1] - ev[j]) / (2 * e[j]); double t = FastMath.sqrt(1 + q * q); if(q < 0.0) { @@ -453,6 +515,7 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, else { q = ev[m] - ev[j] + e[j] / (q + t); } + double u = 0.0; double s = 1.0; double c = 1.0; @@ -486,63 +549,18 @@ private static MatrixBlock[] findEigenVectors(double[] main, double[] secondary, q = c * t - h; for(int ia = 0; ia < n; ++ia) { - p = hhvalues[(i + 1) * n + ia]; - hhvalues[(i + 1) * n + ia] = s * hhvalues[i * n + ia] + c * p; - hhvalues[i * n + ia] = c * hhvalues[i * n + ia] - s * p; + p = hhValues[(i + 1) * n + ia]; + hhValues[(i + 1) * n + ia] = s * hhValues[i * n + ia] + c * p; + hhValues[i * n + ia] = c * hhValues[i * n + ia] - s * p; } } - if(t == 0.0 && i >= j) { - continue; - } + if(t == 0.0 && i >= j) + return; + ev[j] -= u; e[j] = q; e[m] = 0.0; - } - - } - - // Sort eigen values and vectors in decreasing order - for(int i = 0; i < n; i++) { - int k = i; - double p = ev[i]; - for(int j = i + 1; j < n; j++) { - if(ev[j] < p) { - k = j; - p = ev[j]; - } - } - if(k != i) { - ev[k] = ev[i]; - ev[i] = p; - for(int j = 0; j < n; j++) { - p = hhvalues[i * n + j]; - hhvalues[i * n + j] = hhvalues[k * n + j]; - hhvalues[k * n + j] = p; - } - } - } - - // Determine the largest eigen value in absolute term. - maxAbsoluteValue = 0; - for(int i = 0; i < n; i++) { - maxAbsoluteValue = FastMath.max(maxAbsoluteValue, FastMath.abs(ev[i])); - } - // Make null any eigen value too small to be significant - int zeros = 0; - if(maxAbsoluteValue != 0.0) { - for(int i = 0; i < n; i++) { - if(FastMath.abs(ev[i]) < EIGEN_MACHINE_EPS * maxAbsoluteValue) { - ev[i] = 0; - zeros++; - } - } - } - - eigenValues.setNonZeros(n - zeros); - hhMatrix.setNonZeros(hhMatrix.denseBlock.countNonZeros()); - - return new MatrixBlock[] {eigenValues, hhMatrix}; } /** From 2978f8c53ff50e4e50f7984de827617d18aedd68 Mon Sep 17 00:00:00 2001 From: Alam Sanchez Date: Mon, 8 Jul 2024 23:47:53 +0200 Subject: [PATCH 6/6] cleanup --- .../org/apache/sysds/runtime/matrix/data/LibCommonsMath.java | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java index 0aff7b731b3..c082acddd9b 100644 --- a/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java +++ b/src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java @@ -147,8 +147,7 @@ else if (opcode.equals("qr2")) else if (opcode.equals("lu")) return computeLU(in); else if (opcode.equals("eigen")) - return computeEigenDecompositionSymm(in, threads); - // return computeEigen(in); + return computeEigen(in); else if (opcode.equals("eigen_lanczos")) return computeEigenLanczos(in, threads, seed); else if (opcode.equals("eigen_qr"))