Lampisan commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r777988773
##########
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:
In the following calculations both w1 and w1_t are needed.
We tried to do it inplace (by double transposing), but that did not work.
--
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]
