Lampisan commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r777989544
##########
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:
We are not sure if we can parallelize that loop, since it is a iterative
algorithm.
Or is there another way of making it parallel?
--
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]
