Baunsgaard commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r775948322
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -33,6 +33,9 @@
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.*;
Review comment:
don't do wildcard imports, only import things that are in use.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
Review comment:
this isSym check is taking much of your time,
In the computeEigen we have skipped this check, i think it is fair to do the
same here.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
Review comment:
use higher abstraction as the field type (ScalarOperator)
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
Review comment:
i think we have a in place transpose of vectors. make sure you use it
here to avoid allocation
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
+ int m = in.rlen;
+
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock Q_prod = new MatrixBlock(m, m, 0.0);
Review comment:
use other constructor
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
+ int m = in.rlen;
+
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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);
+ Q_prod = Q_prod.aggregateBinaryOperations(Q_prod,
QR[0], op_mul_agg);
+ in = in.aggregateBinaryOperations(QR[1], QR[0],
op_mul_agg);
+ }
+
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < m; j++) {
+ if(i != j) {
+ if(Math.abs(in.getValue(i, j)) > tol)
+ throw new
DMLRuntimeException("QR Eigen Decomposition not converged or contains complex
EVs"
+ + " (position i = " + i
+ ", j = " + j + ")");
+ }
+ }
+ }
+
+ // If complex evals then A is in real Schur-form
+ // 2x2 blocks on diagonal of A correspond to a matrix that has
the same complex evals
+ double[] eval = new double[m];
+ for(int i = 0; i < m; i++) {
+ eval[i] = in.getValue(i, i);
+ }
+
+ double[][] evec =
DataConverter.convertToArray2DRowRealMatrix(Q_prod).getData();
+ return sortEVs(eval, evec);
+ }
+
+ /**
+ * Function to compute the Householder transformation of a Matrix.
+ *
+ * @param in Input Matrix
+ * @return transformed matrix
+ */
+ private static MatrixBlock computeHouseholder(MatrixBlock in) {
+ int num_Threads = 1;
+ 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);
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 2 * r, num_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(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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) {
Review comment:
consider if you can modify this method to not use the double [][] of
eVectors, but instead use a the already materialized double[] underlying the
MatrixBlock.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
Review comment:
give the number of threads as argument to the method.
and add tests with 1 and more threads.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
Review comment:
this can be done in place on v1
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
Review comment:
you can allocate the matrix with a specific value in one of the other
constructors.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
Review comment:
i like the Coffee but allow the method to have a seed argument (you can
always add coffee to the seed if you want to have a little easteregg)!
i would like it if you isolate the entire allocation of v1 into another
method, this method should return the normalized v1.
This allows us to consider if we can find better ways of generating this,
and profilers will make it obvious if this is a bottleneck.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
Review comment:
tolerance in argument
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
Review comment:
add threads argument
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
Review comment:
remove commented code
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
+ int m = in.rlen;
+
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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++) {
Review comment:
add todo or make this loop parallel.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
+ int m = in.rlen;
+
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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);
+ Q_prod = Q_prod.aggregateBinaryOperations(Q_prod,
QR[0], op_mul_agg);
+ in = in.aggregateBinaryOperations(QR[1], QR[0],
op_mul_agg);
+ }
+
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < m; j++) {
+ if(i != j) {
+ if(Math.abs(in.getValue(i, j)) > tol)
Review comment:
call getDenseBlockValues to allow you to iterate through a double[]
array directly instead.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
Review comment:
Use higher abstraction for the field type
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
+ int m = in.rlen;
+
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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);
+ Q_prod = Q_prod.aggregateBinaryOperations(Q_prod,
QR[0], op_mul_agg);
+ in = in.aggregateBinaryOperations(QR[1], QR[0],
op_mul_agg);
+ }
+
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < m; j++) {
+ if(i != j) {
+ if(Math.abs(in.getValue(i, j)) > tol)
+ throw new
DMLRuntimeException("QR Eigen Decomposition not converged or contains complex
EVs"
+ + " (position i = " + i
+ ", j = " + j + ")");
+ }
+ }
+ }
+
+ // If complex evals then A is in real Schur-form
+ // 2x2 blocks on diagonal of A correspond to a matrix that has
the same complex evals
+ double[] eval = new double[m];
+ for(int i = 0; i < m; i++) {
+ eval[i] = in.getValue(i, i);
+ }
+
+ double[][] evec =
DataConverter.convertToArray2DRowRealMatrix(Q_prod).getData();
Review comment:
if we can avoid getting the values out as 2dRowRealMatrix, then it would
be nice.
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
+ int m = in.rlen;
+
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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);
+ Q_prod = Q_prod.aggregateBinaryOperations(Q_prod,
QR[0], op_mul_agg);
+ in = in.aggregateBinaryOperations(QR[1], QR[0],
op_mul_agg);
+ }
+
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < m; j++) {
+ if(i != j) {
+ if(Math.abs(in.getValue(i, j)) > tol)
+ throw new
DMLRuntimeException("QR Eigen Decomposition not converged or contains complex
EVs"
+ + " (position i = " + i
+ ", j = " + j + ")");
+ }
+ }
+ }
+
+ // If complex evals then A is in real Schur-form
+ // 2x2 blocks on diagonal of A correspond to a matrix that has
the same complex evals
+ double[] eval = new double[m];
+ for(int i = 0; i < m; i++) {
+ eval[i] = in.getValue(i, i);
+ }
+
+ double[][] evec =
DataConverter.convertToArray2DRowRealMatrix(Q_prod).getData();
+ return sortEVs(eval, evec);
+ }
+
+ /**
+ * Function to compute the Householder transformation of a Matrix.
+ *
+ * @param in Input Matrix
+ * @return transformed matrix
+ */
+ private static MatrixBlock computeHouseholder(MatrixBlock in) {
+ int num_Threads = 1;
+ 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);
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 2 * r, num_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(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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 boolean isSym(MatrixBlock in) {
Review comment:
1. consider if you want to use it the one place it is in use i commented
that you probably should not. if you decide to throw it away ignore 2.
2. The convertToArray2DRowReadlMatrix call makes this exceedingly slow,
instead use the getDenseDoubleValues() to get a row major double[] that you can
compare in.
##########
File path:
src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
##########
@@ -0,0 +1,112 @@
+/*
+ * 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 org.apache.sysds.runtime.matrix.data.LibCommonsMath;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.runtime.util.DataConverter;
+import org.apache.sysds.test.TestUtils;
+import org.junit.Test;
+import org.junit.Assert;
+
+public class EigenDecompTest {
+
+ @Test public void testLanczosSimple() {
+ double tol = 1e-4;
+
+ MatrixBlock in = new MatrixBlock(4, 4, false);
+ double[] a = { 4, 1, -2, 2,
+ 1, 2, 0, 1,
+ -2, 0, 3, -2,
+ 2, 1, -2, -1};
Review comment:
format
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
+ double tol = 1e-10;
+ int m = in.rlen;
+
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_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);
+ Q_prod = Q_prod.aggregateBinaryOperations(Q_prod,
QR[0], op_mul_agg);
+ in = in.aggregateBinaryOperations(QR[1], QR[0],
op_mul_agg);
+ }
+
+ for(int i = 0; i < m; i++) {
+ for(int j = 0; j < m; j++) {
+ if(i != j) {
+ if(Math.abs(in.getValue(i, j)) > tol)
+ throw new
DMLRuntimeException("QR Eigen Decomposition not converged or contains complex
EVs"
+ + " (position i = " + i
+ ", j = " + j + ")");
+ }
+ }
+ }
+
+ // If complex evals then A is in real Schur-form
+ // 2x2 blocks on diagonal of A correspond to a matrix that has
the same complex evals
+ double[] eval = new double[m];
+ for(int i = 0; i < m; i++) {
+ eval[i] = in.getValue(i, i);
Review comment:
here also leverage the double array from the comment above.
##########
File path: src/test/java/org/apache/sysds/test/TestUtils.java
##########
@@ -1802,6 +1802,16 @@ public static MatrixBlock generateTestMatrixBlock(int
rows, int cols, double min
return MatrixBlock.randOperations(rows, cols, sparsity, min,
max, "Uniform", seed);
}
+ 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++) {
Review comment:
consider starting from i = 1
##########
File path:
src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -295,11 +281,283 @@ private static MatrixBlock
computeMatrixInverse(Array2DRowRealMatrix in) {
* @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());
}
+
+ /**
+ * 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenLanczos(MatrixBlock in) {
+ 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() + ")");
+ }
+ if(!isSym(in)) {
+ throw new DMLRuntimeException("Lanczos algorithm can
only be done on a symmetric matrix.");
+ }
+
+ int num_Threads = 1;
+
+ int m = in.getNumRows();
+ MatrixBlock v0 = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock v1 = MatrixBlock.randOperations(m, 1, 1.0, 0.0,
1.0, "UNIFORM", 0xC0FFEE);
+
+ // normalize v1
+ double v1_sum = v1.sum();
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), v1_sum, num_Threads);
+ v1 = v1.scalarOperations(op_div_scalar, new MatrixBlock());
+ UnaryOperator op_sqrt = new
UnaryOperator(Builtin.getBuiltinFnObject(Builtin.BuiltinCode.SQRT),
num_Threads, true);
+ v1 = v1.unaryOperations(op_sqrt, new MatrixBlock());
+ if(v1.sumSq() != 1.0)
+ throw new DMLRuntimeException("Lanczos algorithm: v1
not correctly normalized");
+
+ MatrixBlock T = new MatrixBlock(m, m, 0.0);
+ MatrixBlock TV = new MatrixBlock(m, 1, 0.0);
+ MatrixBlock w1;
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ TernaryOperator op_minus_mul = new
TernaryOperator(MinusMultiply.getFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+
+ MatrixBlock beta = new MatrixBlock(1, 1, 0.0);
+ for(int i = 0; i < m; i++) {
+ if(i == 0)
+ TV.copy(v1);
+ else
+ TV = TV.append(v1, new MatrixBlock(), true);
+
+ w1 = in.aggregateBinaryOperations(in, v1, op_mul_agg);
+ MatrixBlock w1_t = w1.reorgOperations(op_t, new
MatrixBlock(), 0, 0, m);
+ MatrixBlock alpha =
w1_t.aggregateBinaryOperations(w1_t, v1, 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 = (RightScalarOperator)
op_div_scalar.setConstant(beta.getDouble(0, 0));
+ v1 = w1.scalarOperations(op_div_scalar, new
MatrixBlock());
+
+ 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");
+ e[1] = TV.aggregateBinaryOperations(TV, e[1], op_mul_agg);
+ return e;
+ }
+
+ /**
+ * Function to perform the QR decomposition.
+ * Input must be a square matrix.
+ *
+ * @param in matrix object
+ * @return array of matrix blocks [Q, R]
+ */
+ private static MatrixBlock[] computeQR2(MatrixBlock in) {
+ 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 num_Threads = 1;
+ int m = in.rlen;
+
+ MatrixBlock A_n = new MatrixBlock(m, m, 0.0);
+ A_n.copy(in);
+
+ MatrixBlock Q_n = new MatrixBlock(m, m, 0.0);
+ for(int i = 0; i < m; i++) {
+ Q_n.setValue(i, i, 1.0);
+ }
+
+ ReorgOperator op_t = new
ReorgOperator(SwapIndex.getSwapIndexFnObject(), num_Threads);
+ AggregateBinaryOperator op_mul_agg =
InstructionUtils.getMatMultOperator(num_Threads);
+ BinaryOperator op_sub =
InstructionUtils.parseExtendedBinaryOperator("-");
+ RightScalarOperator op_div_scalar = new
RightScalarOperator(Divide.getDivideFnObject(), 1, num_Threads);
+ LeftScalarOperator op_mult_2 = new
LeftScalarOperator(Multiply.getMultiplyFnObject(), 2, num_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 = (RightScalarOperator)
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 vkvkt = vk.aggregateBinaryOperations(vk,
vkt, op_mul_agg);
+ MatrixBlock vkvkt2 = vkvkt.scalarOperations(op_mult_2,
new MatrixBlock());
+
+ 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));
+ }
+ // Q_n= Q
+ // A_n = R
+ 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
+ * @return array of matrix blocks
+ */
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in) {
+ return computeEigenQR(in, 300);
+ }
+
+ private static MatrixBlock[] computeEigenQR(MatrixBlock in, int
num_iterations) {
+ 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 num_Threads = 1;
Review comment:
same, threads in argument
##########
File path: src/test/java/org/apache/sysds/test/TestUtils.java
##########
@@ -1802,6 +1802,16 @@ public static MatrixBlock generateTestMatrixBlock(int
rows, int cols, double min
return MatrixBlock.randOperations(rows, cols, sparsity, min,
max, "Uniform", seed);
}
+ public static MatrixBlock generateTestMatrixBlockSym(int rows, int
cols, double min, double max, double sparsity, long seed){
Review comment:
add java doc comment here, to make it easy to find this method,(add
keywords or something) we have a bunch of tests using symmetric matrices and
all implement their own, and it would be nice if everyone use the same
generator (this could be it!).
##########
File path:
src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
##########
@@ -0,0 +1,112 @@
+/*
+ * 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 org.apache.sysds.runtime.matrix.data.LibCommonsMath;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.runtime.util.DataConverter;
+import org.apache.sysds.test.TestUtils;
+import org.junit.Test;
+import org.junit.Assert;
+
+public class EigenDecompTest {
+
+ @Test public void testLanczosSimple() {
Review comment:
keyword line above
##########
File path:
src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
##########
@@ -0,0 +1,112 @@
+/*
+ * 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 org.apache.sysds.runtime.matrix.data.LibCommonsMath;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.runtime.util.DataConverter;
+import org.apache.sysds.test.TestUtils;
+import org.junit.Test;
+import org.junit.Assert;
+
+public class EigenDecompTest {
+
+ @Test public void testLanczosSimple() {
+ double tol = 1e-4;
+
+ MatrixBlock in = new MatrixBlock(4, 4, false);
+ double[] a = { 4, 1, -2, 2,
+ 1, 2, 0, 1,
+ -2, 0, 3, -2,
+ 2, 1, -2, -1};
+ in.init(a, 4, 4);
+ MatrixBlock[] m1 = LibCommonsMath.multiReturnOperations(in,
"eigen");
+ MatrixBlock[] m2 = LibCommonsMath.multiReturnOperations(in,
"eigen_lanczos");
+
+ TestUtils.compareMatrices(m1[0], m2[0], tol, "Result of
eigenvalues of new eigen_lanczos function wrong");
+ testEvecValues(m1[1], m2[1], tol);
+ }
+
+ @Test public void testLanczosRandom() {
Review comment:
put Test keyword on the line above.
##########
File path:
src/test/java/org/apache/sysds/test/component/matrix/EigenDecompTest.java
##########
@@ -0,0 +1,112 @@
+/*
+ * 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 org.apache.sysds.runtime.matrix.data.LibCommonsMath;
+import org.apache.sysds.runtime.matrix.data.MatrixBlock;
+import org.apache.sysds.runtime.util.DataConverter;
+import org.apache.sysds.test.TestUtils;
+import org.junit.Test;
+import org.junit.Assert;
+
+public class EigenDecompTest {
+
+ @Test public void testLanczosSimple() {
+ double tol = 1e-4;
+
+ MatrixBlock in = new MatrixBlock(4, 4, false);
+ double[] a = { 4, 1, -2, 2,
+ 1, 2, 0, 1,
+ -2, 0, 3, -2,
+ 2, 1, -2, -1};
+ in.init(a, 4, 4);
+ MatrixBlock[] m1 = LibCommonsMath.multiReturnOperations(in,
"eigen");
+ MatrixBlock[] m2 = LibCommonsMath.multiReturnOperations(in,
"eigen_lanczos");
+
+ TestUtils.compareMatrices(m1[0], m2[0], tol, "Result of
eigenvalues of new eigen_lanczos function wrong");
+ testEvecValues(m1[1], m2[1], tol);
+ }
+
+ @Test public void testLanczosRandom() {
+ double tol = 1e-4;
+
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10,
0.0, 1.0, 1.0, 1);
+ // MatrixBlock in = TestUtils.generateTestMatrixBlockSym(100,
100, 0.0, 1.0, 1.0, 1); // fails
+ long t1 = System.nanoTime();
+ MatrixBlock[] m1 = LibCommonsMath.multiReturnOperations(in,
"eigen");
+ long t2 = System.nanoTime();
+ MatrixBlock[] m2 = LibCommonsMath.multiReturnOperations(in,
"eigen_lanczos");
+ long t3 = System.nanoTime();
+
+ System.out.println(
+ "time eigen: " + (t2 - t1) + " time Lanczos: " + (t3 -
t2) + " Lanczos speedup: " + ((double) (t2 - t1) / (t3 - t2)));
+ TestUtils.compareMatrices(m1[0], m2[0], tol, "Result of
eigenvalues of new eigen_lanczos function wrong");
+ testEvecValues(m1[1], m2[1], tol);
+ }
+
+ @Test public void testQREigenSimple() {
+ double tol = 1e-4;
+
+ MatrixBlock in = new MatrixBlock(4, 4, false);
+ double[] a = { 52, 30, 49, 28,
+ 30, 50, 8, 44,
+ 49, 8, 46, 16,
+ 28, 44, 16, 22};
+ in.init(a, 4, 4);
+
+ MatrixBlock[] m1 = LibCommonsMath.multiReturnOperations(in,
"eigen");
+ MatrixBlock[] m2 = LibCommonsMath.multiReturnOperations(in,
"eigen_qr");
+
+ TestUtils.compareMatrices(m1[0], m2[0], tol, "Result of
eigenvalues of new eigen_qr function wrong");
+ testEvecValues(m1[1], m2[1], tol);
+ }
+
+ @Test public void testQREigenRandom() {
+ double tol = 1e-4;
+
+ MatrixBlock in = TestUtils.generateTestMatrixBlockSym(10, 10,
0.0, 1.0, 1.0, 1);
+ // MatrixBlock in = TestUtils.generateTestMatrixBlock(5, 5,
0.0, 1.0, 1.0, 5); // fails evals correct evecs wrong (evec corresponding to
largest eval correct)
+ // MatrixBlock in = TestUtils.generateTestMatrixBlock(10, 10,
0.0, 1.0, 1.0, 2); // fails complex EVs
+ // MatrixBlock in = TestUtils.generateTestMatrixBlockSym(50,
50, 0.0, 1.0, 1.0, 1); // fails not converged
Review comment:
instead of leaving this outcommented, add more tests, each calling
another method that execute the test with the given argument MatrixBlock "in"
--
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]
