Wouter-vw commented on code in PR #2038:
URL: https://github.com/apache/systemds/pull/2038#discussion_r1668652719
##########
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java:
##########
@@ -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) {
Review Comment:
Yes there is a reason. From page 263 of "Matrix Computations" by Gene H.
Golub and Charles F. Van Loan: "Recall that the leading (j − 1)-by-(j − 1)
portion of Qj is the identity. Thus, at the beginning of backward accumulation,
Q is “mostly the identity” and it gradually becomes full as the iteration
progresses. This pattern can be exploited to reduce the number of required
flops. In contrast, Q is full in forward accumulation after the first step. For
this reason, backward accumulation is cheaper and the strategy of choice."
--
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
To unsubscribe, e-mail: [email protected]
For queries about this service, please contact Infrastructure at:
[email protected]
