Baunsgaard commented on a change in pull request #1489: URL: https://github.com/apache/systemds/pull/1489#discussion_r778330493
########## 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: unclear comment from my side again, i was referring to using the num_threads argument and make the instructions parallel. this would entail adding the num_threads to the computeQR2 call. -- 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: dev-unsubscr...@systemds.apache.org For queries about this service, please contact Infrastructure at: us...@infra.apache.org