Baunsgaard commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r778328559
##########
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:
My comment is hard to understand and not saying what i intended.
You only use this transposed form on the next line for a matrix
multiplication,
this can be seen as:
alpha = t(w1) %*% v1
you could consider using (this gives a equivalent output):
alpha = t(t(v1) %*% W1)
since the v1 is a vector it will not be allocated when transposing (the
underlying dense double[] is the same),
making it a metadata operation that cost nothing, unlike the transpose of W1
that is potentially costly.
The challenging part is to write the calls in such a way that they do not
unnecessarily allocate alpha, and v1 more than once.
--
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]
