This is an automated email from the ASF dual-hosted git repository.
baunsgaard pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/systemds.git
The following commit(s) were added to refs/heads/main by this push:
new 7b00190 [SYSTEMDS-3167] Internal Eigen Decomposition
7b00190 is described below
commit 7b001906d1cb39c480f0ab4f411d6d8b6a0e44be
Author: Palmisano Alexander <[email protected]>
AuthorDate: Sat Nov 20 11:23:34 2021 +0100
[SYSTEMDS-3167] Internal Eigen Decomposition
This commit adds two internal implementation of Eigen Decomposition.
One using the Lanczos algorithm, another using QR Decomposition.
Also added is an initial householder transformation that is intended
to be used in QR decomposition as an accelerator.
The tests verify the correctness of the decomposition through
reconstruction of the input symmetric matrix using:
A = V %*% D %*% t(V)
Future work is to improve the results in higher dimensions,
that currently return unstable decompositions and performance.
Co-authored-by: Palmisano Alexander <[email protected]>
Co-authored-by: Nikolaus Grogger <[email protected]>
Closes #1489
Closes #1561
---
scripts/builtin/pcaInverse.dml | 2 +-
scripts/builtin/pcaTransform.dml | 2 +-
.../sysds/runtime/matrix/data/LibCommonsMath.java | 369 ++++++++++++++++++---
src/test/java/org/apache/sysds/test/TestUtils.java | 35 ++
.../test/component/matrix/EigenDecompTest.java | 152 +++++++++
5 files changed, 521 insertions(+), 39 deletions(-)
diff --git a/scripts/builtin/pcaInverse.dml b/scripts/builtin/pcaInverse.dml
index 96983a3..ff8f662 100644
--- a/scripts/builtin/pcaInverse.dml
+++ b/scripts/builtin/pcaInverse.dml
@@ -20,7 +20,7 @@
#-------------------------------------------------------------
# Principal Component Analysis (PCA) for reconstruction of approximation of
the original data.
-# This methods allows to reconstruct an approximation of the original matrix,
and is usefull for
+# This methods allows to reconstruct an approximation of the original matrix,
and is useful for
# calculating how much information is lost in the PCA.
#
# INPUT PARAMETERS:
diff --git a/scripts/builtin/pcaTransform.dml b/scripts/builtin/pcaTransform.dml
index b7588c8..07bdbee 100644
--- a/scripts/builtin/pcaTransform.dml
+++ b/scripts/builtin/pcaTransform.dml
@@ -19,7 +19,7 @@
#
#-------------------------------------------------------------
-# Principal Component Analysis (PCA) for dimensionality reduction prediciton
+# Principal Component Analysis (PCA) for dimensionality reduction prediction
# This method is used to transpose data, which the PCA model was not trained
on. To validate how good
# The PCA is, and to apply in production.
#
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 a49d1ce..c71838c 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
@@ -33,6 +33,20 @@ import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.SingularValueDecomposition;
import org.apache.sysds.runtime.DMLRuntimeException;
import org.apache.sysds.runtime.data.DenseBlock;
+import org.apache.sysds.runtime.functionobjects.Multiply;
+import org.apache.sysds.runtime.functionobjects.Divide;
+import org.apache.sysds.runtime.functionobjects.SwapIndex;
+import org.apache.sysds.runtime.functionobjects.MinusMultiply;
+import org.apache.sysds.runtime.functionobjects.Builtin;
+import org.apache.sysds.runtime.instructions.InstructionUtils;
+import org.apache.sysds.runtime.matrix.operators.RightScalarOperator;
+import org.apache.sysds.runtime.matrix.operators.LeftScalarOperator;
+import org.apache.sysds.runtime.matrix.operators.ScalarOperator;
+import org.apache.sysds.runtime.matrix.operators.UnaryOperator;
+import org.apache.sysds.runtime.matrix.operators.TernaryOperator;
+import org.apache.sysds.runtime.matrix.operators.ReorgOperator;
+import org.apache.sysds.runtime.matrix.operators.AggregateBinaryOperator;
+import org.apache.sysds.runtime.matrix.operators.BinaryOperator;
import org.apache.sysds.runtime.util.DataConverter;
/**
@@ -72,15 +86,32 @@ public class LibCommonsMath
return computeCholesky(matrixInput);
return null;
}
-
+
public static MatrixBlock[] multiReturnOperations(MatrixBlock in,
String opcode) {
+ return multiReturnOperations(in, opcode, 1, 1);
+ }
+
+ public static MatrixBlock[] multiReturnOperations(MatrixBlock in,
String opcode, int threads, int num_iterations, double tol) {
+ if(opcode.equals("eigen_qr"))
+ return computeEigenQR(in, num_iterations, tol, threads);
+ else
+ return multiReturnOperations(in, opcode, threads, 1);
+ }
+
+ public static MatrixBlock[] multiReturnOperations(MatrixBlock in,
String opcode, int threads, long seed) {
if(opcode.equals("qr"))
return computeQR(in);
+ else if (opcode.equals("qr2"))
+ return computeQR2(in, threads);
else if (opcode.equals("lu"))
return computeLU(in);
else if (opcode.equals("eigen"))
return computeEigen(in);
- else if ( opcode.equals("svd"))
+ 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("svd"))
return computeSvd(in);
return null;
}
@@ -148,8 +179,10 @@ public class LibCommonsMath
* @return array of matrix blocks
*/
private static MatrixBlock[] computeLU(MatrixBlock in) {
- if ( in.getNumRows() != in.getNumColumns() ) {
- throw new DMLRuntimeException("LU Decomposition can
only be done on a square matrix. Input matrix is rectangular (rows=" +
in.getNumRows() + ", cols="+ in.getNumColumns() +")");
+ if(in.getNumRows() != in.getNumColumns()) {
+ throw new DMLRuntimeException(
+ "LU Decomposition can only be done on a square
matrix. Input matrix is rectangular (rows="
+ + in.getNumRows() + ", cols=" +
in.getNumColumns() + ")");
}
Array2DRowRealMatrix matrixInput =
DataConverter.convertToArray2DRowRealMatrix(in);
@@ -194,35 +227,10 @@ public class LibCommonsMath
RealMatrix eVectorsMatrix = eigendecompose.getV();
double[][] eVectors = eVectorsMatrix.getData();
double[] eValues = eigendecompose.getRealEigenvalues();
-
- //Sort the eigen values (and vectors) in increasing order (to
be compatible w/ LAPACK.DSYEVR())
- int n = eValues.length;
- for (int i = 0; i < n; i++) {
- int k = i;
- double p = eValues[i];
- for (int j = i + 1; j < n; j++) {
- if (eValues[j] < p) {
- k = j;
- p = eValues[j];
- }
- }
- if (k != i) {
- eValues[k] = eValues[i];
- eValues[i] = p;
- for (int j = 0; j < n; j++) {
- p = eVectors[j][i];
- eVectors[j][i] = eVectors[j][k];
- eVectors[j][k] = p;
- }
- }
- }
- MatrixBlock mbValues =
DataConverter.convertToMatrixBlock(eValues, true);
- MatrixBlock mbVectors =
DataConverter.convertToMatrixBlock(eVectors);
-
- return new MatrixBlock[] { mbValues, mbVectors };
+ return sortEVs(eValues, eVectors);
}
-
+
private static EigenDecomposition computeEigenRegularized(MatrixBlock
in) {
if( in == null || in.isEmptyBlock(false) )
throw new DMLRuntimeException("Invalid empty block");
@@ -244,7 +252,6 @@ public class LibCommonsMath
DataConverter.convertToArray2DRowRealMatrix(in2));
}
-
/**
* Performs Singular Value Decomposition. Calls Apache Commons Math SVD.
* X = U * Sigma * Vt, where X is the input matrix,
@@ -277,9 +284,10 @@ public class LibCommonsMath
* @return matrix block
*/
private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix
in) {
- if ( !in.isSquare() )
- throw new DMLRuntimeException("Input to inv() must be
square matrix -- given: a " + in.getRowDimension() + "x" +
in.getColumnDimension() + " matrix.");
-
+ if(!in.isSquare())
+ throw new DMLRuntimeException("Input to inv() must be
square matrix -- given: a " + in.getRowDimension()
+ + "x" + in.getColumnDimension() + " matrix.");
+
QRDecomposition qrdecompose = new QRDecomposition(in);
DecompositionSolver solver = qrdecompose.getSolver();
RealMatrix inverseMatrix = solver.getInverse();
@@ -295,11 +303,298 @@ public class LibCommonsMath
* @return matrix block
*/
private static MatrixBlock computeCholesky(Array2DRowRealMatrix in) {
- if ( !in.isSquare() )
- throw new DMLRuntimeException("Input to cholesky() must
be square matrix -- given: a " + in.getRowDimension() + "x" +
in.getColumnDimension() + " matrix.");
+ if(!in.isSquare())
+ throw new DMLRuntimeException("Input to cholesky() must
be square matrix -- given: a "
+ + in.getRowDimension() + "x" +
in.getColumnDimension() + " matrix.");
CholeskyDecomposition cholesky = new CholeskyDecomposition(in,
RELATIVE_SYMMETRY_THRESHOLD,
CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
RealMatrix rmL = cholesky.getL();
return DataConverter.convertToMatrixBlock(rmL.getData());
}
+
+ /**
+ * Creates a random normalized vector with dim elements.
+ *
+ * @param dim dimension of the created vector
+ * @param threads number of threads
+ * @param seed seed for the random MatrixBlock generation
+ * @return random normalized vector
+ */
+ private static MatrixBlock randNormalizedVect(int dim, int threads,
long seed) {
+ MatrixBlock v1 = MatrixBlock.randOperations(dim, 1, 1.0, 0.0,
1.0, "UNIFORM", seed);
+
+ double v1_sum = v1.sum();
+ ScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT), threads,
true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+
+ if(Math.abs(v1.sumSq() - 1.0) >= 1e-7)
+ throw new DMLRuntimeException("v1 not correctly
normalized (maybe try changing the seed)");
+
+ return v1;
+ }
+
+ /**
+ * Function to perform the Lanczos algorithm and then computes the
Eigendecomposition.
+ * Caution: Lanczos is not numerically stable (see
https://en.wikipedia.org/wiki/Lanczos_algorithm)
+ * Input must be a symmetric (and square) matrix.
+ *
+ * @param in matrix object
+ * @param threads number of threads
+ * @param seed seed for the random MatrixBlock generation
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in, int
threads, long seed) {
+ if(in.getNumRows() != in.getNumColumns()) {
+ throw new DMLRuntimeException(
+ "Lanczos algorithm and Eigen Decomposition can
only be done on a square matrix. "
+ + "Input matrix is rectangular (rows="
+ in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
+ }
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = randNormalizedVect(m, threads, seed);
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, m, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(threads);
+ ScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ v1.putInto(TV, 0, i, false);
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock alpha =
w1.aggregateBinaryOperations(v1.reorgOperations(op_t, new MatrixBlock(), 0, 0,
m), w1, op_mul_agg);
+ if(i < m - 1) {
+ w1 = w1.ternaryOperations(op_minus_mul, v1,
alpha, new MatrixBlock());
+ w1 = w1.ternaryOperations(op_minus_mul, v0,
beta, new MatrixBlock());
+ beta.setValue(0, 0, Math.sqrt(w1.sumSq()));
+ v0.copy(v1);
+ op_div_scalar =
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ w1.scalarOperations(op_div_scalar, v1);
+
+ T.setValue(i + 1, i, beta.getValue(0, 0));
+ T.setValue(i, i + 1, beta.getValue(0, 0));
+ }
+ T.setValue(i, i, alpha.getValue(0, 0));
+ }
+
+ MatrixBlock[] e = multiReturnOperations(T, "eigen");
+ TV.setNonZeros((long) m*m);
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ * TODO: use Householder transformation and implicit shifts to further
speed up QR decompositions
+ *
+ * @param in matrix object
+ * @param threads number of threads
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in, int threads) {
+ if(in.getNumRows() != in.getNumColumns()) {
+ throw new DMLRuntimeException("QR2 Decomposition can
only be done on a square matrix. "
+ + "Input matrix is rectangular (rows=" +
in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
+ }
+
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock();
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, true);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ ScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, threads);
+ ScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, threads);
+
+ for(int k = 0; k < m; k++) {
+ MatrixBlock z = A_n.slice(k, m - 1, k, k);
+ MatrixBlock uk = new MatrixBlock(m - k, 1, 0.0);
+ uk.copy(z);
+ uk.setValue(0, 0, uk.getValue(0, 0) +
Math.signum(z.getValue(0, 0)) * Math.sqrt(z.sumSq()));
+ op_div_scalar =
op_div_scalar.setConstant(Math.sqrt(uk.sumSq()));
+ uk = uk.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ MatrixBlock vk = new MatrixBlock(m, 1, 0.0);
+ vk.copy(k, m - 1, 0, 0, uk, true);
+
+ MatrixBlock vkt = vk.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock vkt2 = vkt.scalarOperations(op_mult_2, new
MatrixBlock());
+ MatrixBlock vkvkt2 = vk.aggregateBinaryOperations(vk,
vkt2, op_mul_agg);
+
+ A_n = A_n.binaryOperations(op_sub,
A_n.aggregateBinaryOperations(vkvkt2, A_n, op_mul_agg));
+ Q_n = Q_n.binaryOperations(op_sub,
Q_n.aggregateBinaryOperations(Q_n, vkvkt2, op_mul_agg));
+ }
+ // QR decomp: Q: Q_n; R: A_n
+ return new MatrixBlock[] {Q_n, A_n};
+ }
+
+ /**
+ * Function that computes the Eigen Decomposition using the QR
algorithm.
+ * Caution: check if the QR algorithm is converged, if not increase
iterations
+ * Caution: if the input matrix has complex eigenvalues results will be
incorrect
+ *
+ * @param in Input matrix
+ * @param threads number of threads
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
threads) {
+ return computeEigenQR(in, 100, 1e-10, threads);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations, double tol, int threads) {
+ if(in.getNumRows() != in.getNumColumns()) {
+ throw new DMLRuntimeException("Eigen Decomposition (QR)
can only be done on a square matrix. "
+ + "Input matrix is rectangular (rows=" +
in.getNumRows() + ", cols=" + in.getNumColumns() + ")");
+ }
+
+ int m = in.rlen;
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(threads);
+
+ MatrixBlock Q_prod = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_prod.setValue(i, i, 1.0);
+ }
+
+ for(int i = 0; i < num_iterations; i++) {
+ MatrixBlock[] QR = computeQR2(in, threads);
+ Q_prod = Q_prod.aggregateBinaryOperations(Q_prod,
QR[0], op_mul_agg);
+ in = QR[1].aggregateBinaryOperations(QR[1], QR[0],
op_mul_agg);
+ }
+
+ // Is converged if all values are below tol and the there only
is values on the diagonal.
+
+ double[] check = in.getDenseBlockValues();
+ double[] eval = new double[m];
+ for(int i = 0; i < m; i++)
+ eval[i] = check[i*m+i];
+
+ double[] evec = Q_prod.getDenseBlockValues();
+ return sortEVs(eval, evec);
+ }
+
+ /**
+ * Function to compute the Householder transformation of a Matrix.
+ *
+ * @param in Input Matrix
+ * @param threads number of threads
+ * @return transformed matrix
+ */
+ @SuppressWarnings("unused")
+ private static MatrixBlock computeHouseholder(MatrixBlock in, int
threads) {
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ for(int k = 0; k < m - 2; k++) {
+ MatrixBlock ajk = A_n.slice(0, m - 1, k, k);
+ for(int i = 0; i <= k; i++) {
+ ajk.setValue(i, 0, 0.0);
+ }
+ double alpha = Math.sqrt(ajk.sumSq());
+ double ak1k = A_n.getDouble(k + 1, k);
+ if(ak1k > 0.0)
+ alpha *= -1;
+ double r = Math.sqrt(0.5 * (alpha * alpha - ak1k *
alpha));
+ MatrixBlock v = new MatrixBlock(m, 1, 0.0);
+ v.copy(ajk);
+ v.setValue(k + 1, 0, ak1k - alpha);
+ ScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 2 * r, threads);
+ v = v.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ MatrixBlock P = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ P.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(threads);
+ BinaryOperator op_add =
InstructionUtils.parseExtendedBinaryOperator("+");
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+
+ MatrixBlock v_t = v.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ v_t = v_t.binaryOperations(op_add, v_t);
+ MatrixBlock v_v_t_2 = A_n.aggregateBinaryOperations(v,
v_t, op_mul_agg);
+ P = P.binaryOperations(op_sub, v_v_t_2);
+ A_n = A_n.aggregateBinaryOperations(P,
A_n.aggregateBinaryOperations(A_n, P, op_mul_agg), op_mul_agg);
+ }
+ return A_n;
+ }
+
+ /**
+ * Sort the eigen values (and vectors) in increasing order (to be
compatible w/ LAPACK.DSYEVR())
+ *
+ * @param eValues Eigenvalues
+ * @param eVectors Eigenvectors
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] sortEVs(double[] eValues, double[][]
eVectors) {
+ int n = eValues.length;
+ for(int i = 0; i < n; i++) {
+ int k = i;
+ double p = eValues[i];
+ for(int j = i + 1; j < n; j++) {
+ if(eValues[j] < p) {
+ k = j;
+ p = eValues[j];
+ }
+ }
+ if(k != i) {
+ eValues[k] = eValues[i];
+ eValues[i] = p;
+ for(int j = 0; j < n; j++) {
+ p = eVectors[j][i];
+ eVectors[j][i] = eVectors[j][k];
+ eVectors[j][k] = p;
+ }
+ }
+ }
+
+ MatrixBlock eval = DataConverter.convertToMatrixBlock(eValues,
true);
+ MatrixBlock evec = DataConverter.convertToMatrixBlock(eVectors);
+ return new MatrixBlock[] {eval, evec};
+ }
+
+ private static MatrixBlock[] sortEVs(double[] eValues, double[]
eVectors) {
+ int n = eValues.length;
+ for(int i = 0; i < n; i++) {
+ int k = i;
+ double p = eValues[i];
+ for(int j = i + 1; j < n; j++) {
+ if(eValues[j] < p) {
+ k = j;
+ p = eValues[j];
+ }
+ }
+ if(k != i) {
+ eValues[k] = eValues[i];
+ eValues[i] = p;
+ for(int j = 0; j < n; j++) {
+ p = eVectors[j*n+i];
+ eVectors[j*n+i] = eVectors[j*n+k];
+ eVectors[j*n+k] = p;
+ }
+ }
+ }
+
+ MatrixBlock eval = DataConverter.convertToMatrixBlock(eValues,
true);
+ MatrixBlock evec = new MatrixBlock(n, n, false);
+ evec.init(eVectors, n, n);
+ return new MatrixBlock[] {eval, evec};
+ }
}
diff --git a/src/test/java/org/apache/sysds/test/TestUtils.java
b/src/test/java/org/apache/sysds/test/TestUtils.java
index 85d3ee7..4729e78 100644
--- a/src/test/java/org/apache/sysds/test/TestUtils.java
+++ b/src/test/java/org/apache/sysds/test/TestUtils.java
@@ -1823,11 +1823,46 @@ public class TestUtils
return matrix;
}
+ /**
+ *
+ * Generates a test matrix with the specified parameters as a
MatrixBlock.
+ *
+ * @param rows number of rows
+ * @param cols number of columns
+ * @param min minimum value
+ * @param max maximum value
+ * @param sparsity sparsity
+ * @param seed seed
+ * @return random MatrixBlock
+ */
public static MatrixBlock generateTestMatrixBlock(int rows, int cols,
double min, double max, double sparsity, long seed){
return MatrixBlock.randOperations(rows, cols, sparsity, min,
max, "Uniform", seed);
}
/**
+ *
+ * Generates a symmetric test matrix with the specified parameters as a
MatrixBlock.
+ * TODO: generate the values without using the randOperation
+ *
+ * @param rows number of rows
+ * @param cols number of columns
+ * @param min minimum value
+ * @param max maximum value
+ * @param sparsity sparsity
+ * @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++) {
+ m.setValue(i,j, m.getValue(j,i));
+ }
+ }
+ return m;
+ }
+
+ /**
* Generates a test matrix with the specified parameters as a two
* dimensional array.
* Set seed to -1 to use the current time as seed.
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
new file mode 100644
index 0000000..f79b44f
--- /dev/null
+++ b/src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
@@ -0,0 +1,152 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you 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.
+ */
+
+package org.apache.sysds.test.component.matrix;
+
+import static org.junit.Assert.fail;
+
+import org.apache.commons.logging.Log;
+import org.apache.commons.logging.LogFactory;
+import org.apache.sysds.runtime.matrix.data.LibCommonsMath;
+import org.apache.sysds.runtime.matrix.data.LibMatrixMult;
+import org.apache.sysds.runtime.matrix.data.LibMatrixReorg;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.test.TestUtils;
+import org.junit.Ignore;
+import org.junit.Test;
+
+public class EigenDecompTest {
+
+ protected static final Log LOG =
LogFactory.getLog(EigenDecompTest.class.getName());
+
+ private enum type {
+ COMMONS, LANCZOS, QR,
+ }
+
+ @Test
+ 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});
+ testEigen(in, 1e-4, 1, type.LANCZOS);
+ }
+
+ @Test
+ public void testLanczosSmall() {
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10,
0.0, 1.0, 1.0, 1);
+ testEigen(in, 1e-4, 1, type.LANCZOS);
+ }
+
+ @Test
+ public void testLanczosMedium() {
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(12, 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);
+ testEigen(in, 1e-4, 1, type.LANCZOS);
+ }
+
+ @Test
+ @Ignore
+ public void testLanczosLargeMT() {
+ int threads = 10;
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100, 100,
0.0, 1.0, 1.0, 1);
+ testEigen(in, 1e-4, threads, type.LANCZOS);
+ }
+
+ @Test
+ 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});
+ testEigen(in, 1e-4, 1, type.QR);
+ }
+
+ @Test
+ public void testQREigenSymSmall() {
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10,
0.0, 1.0, 1.0, 1);
+ testEigen(in, 1e-3, 1, type.QR);
+ }
+
+ @Test
+ public void testQREigenSymSmallMT() {
+ int threads = 10;
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10,
0.0, 1.0, 1.0, 1);
+ testEigen(in, 1e-3, threads, type.QR);
+ }
+
+ @Test
+ @Ignore
+ public void testQREigenSymLarge() {
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 50,
0.0, 1.0, 1.0, 1);
+ testEigen(in, 1e-4, 1, type.QR);
+ }
+
+ @Test
+ @Ignore
+ public void testQREigenSymLargeMT() {
+ int threads = 10;
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50, 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)
{
+ try {
+ MatrixBlock[] m = null;
+ switch(t) {
+ case COMMONS:
+ m =
LibCommonsMath.multiReturnOperations(in, "eigen", threads, 1);
+ break;
+ case LANCZOS:
+ m =
LibCommonsMath.multiReturnOperations(in, "eigen_lanczos", threads, 1);
+ break;
+ case QR:
+ m =
LibCommonsMath.multiReturnOperations(in, "eigen_qr", threads, 1);
+ break;
+ }
+
+ isValidDecomposition(in, m[1], m[0], tol, t.toString());
+
+ }
+ catch(Exception e) {
+ e.printStackTrace();
+ fail(e.getMessage());
+ }
+ }
+
+ 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
+ // D is the eigen values in a diagonal matrix
+ // V is the eigen vectors
+
+ final int m = V.getNumColumns();
+ final MatrixBlock D = new MatrixBlock(m, m, false);
+ LibMatrixReorg.diag(vD, D);
+
+ MatrixBlock VD = new MatrixBlock();
+ LibMatrixMult.matrixMult(V, D, VD);
+
+ MatrixBlock VDtV = new MatrixBlock();
+ LibMatrixMult.matrixMult(VD, LibMatrixReorg.transpose(V), VDtV);
+
+ TestUtils.compareMatrices(A, VDtV, tol, message);
+ }
+}
