[GitHub] [systemds] Lampisan opened a new pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Lampisan opened a new pull request #1489:
URL: https://github.com/apache/systemds/pull/1489


   [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition
   
   This commit adds implementations of the following matrix algorithms in LibCommonsMath: Lanczos algorithm, QR decomposition, QR algorithm and Householder transformation (unused). Furthermore tests for the Lanczos and the QR algorithm are added.
   Lanczos algorithm is numerically unstable and therefore sometimes produces incorrect results.
   QR decomposition is working as intended.
   QR algorithm works for symmetric matrices. While it should theoretically also work for non-symmetric matrices, it currently does not produce all correct eigenvectors in some cases. Furthermore a runtime-exception is thrown if the matrix has complex eigenvalues or the algorithm hasn't converged.
   We are under the impression that our implementations of these algorithms are really slow compared to the existing implementation. Although we only tested this with a single thread this implementation is orders of magnitudes slower.
   
   
   Co-authored-by: Nikolaus Grogger <ni...@student.tugraz.at>


-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Baunsgaard commented on pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#issuecomment-1068909748


   
   Thanks for the contribution, i will merge shortly after minor cleanup!
   I close the PR to add the modifications in another ( #1561 ).


-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] j143 commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

j143 commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r775730328



##########
File path: src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -76,11 +79,17 @@ else if (opcode.equals("cholesky"))
 	public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode) {
 		if(opcode.equals("qr"))
 			return computeQR(in);
+		else if (opcode.equals("qr2"))
+			return computeQR2(in);
 		else if (opcode.equals("lu"))
 			return computeLU(in);
 		else if (opcode.equals("eigen"))
 			return computeEigen(in);
-		else if ( opcode.equals("svd"))
+		else if (opcode.equals("eigen_lanczos"))
+			return computeEigenLanczos(in);
+		else if (opcode.equals("eigen_qr"))
+			return computeEigenQR(in);
+		else if (opcode.equals("svd"))

Review comment:
       Hi @Baunsgaard 
   
   how about we document these functions and mention they are experimental (in the doc)?




-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Lampisan commented on pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Lampisan commented on pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#issuecomment-1004710225


   Hi @Baunsgaard,
   Thanks for your suggestions, we added a few questions to your replies.
   We will push the changes after resolving those.
   


-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Baunsgaard commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r775947098



##########
File path: src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java
##########
@@ -76,11 +79,17 @@ else if (opcode.equals("cholesky"))
 	public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode) {
 		if(opcode.equals("qr"))
 			return computeQR(in);
+		else if (opcode.equals("qr2"))
+			return computeQR2(in);
 		else if (opcode.equals("lu"))
 			return computeLU(in);
 		else if (opcode.equals("eigen"))
 			return computeEigen(in);
-		else if ( opcode.equals("svd"))
+		else if (opcode.equals("eigen_lanczos"))
+			return computeEigenLanczos(in);
+		else if (opcode.equals("eigen_qr"))
+			return computeEigenQR(in);
+		else if (opcode.equals("svd"))

Review comment:
       for now no, these operations are only available sideways calling directly into LibCommonsMath




-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Lampisan commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Lampisan commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r777987955



##########
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:
       We are not sure how to create an identity matrix using the existing MatrixBlock constructors.
   We could create a double[] and use that with MatrixBlock.init(), but this might not be the solution you have in mind.




-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Lampisan commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Lampisan commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r777990577



##########
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:
       If we use i=1 we ignore the first row/col completely.
   We changed it and it did not work, maybe we misunderstand




-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Lampisan edited a comment on pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Lampisan edited a comment on pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#issuecomment-1004710225


   Hi @Baunsgaard and @foebe,
   Thanks for your suggestions, we added a few questions to your replies.
   We will push the changes after resolving those.
   


-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

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: dev-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard closed pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Baunsgaard closed pull request #1489:
URL: https://github.com/apache/systemds/pull/1489


   


-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Lampisan commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

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: dev-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Baunsgaard commented on pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#issuecomment-1005110463


   > Hi @Baunsgaard and @foebe (Nikolaus), Thanks for your suggestions, we added a few questions to your replies. We will push the changes after resolving those.
   
   Thanks for the GitHub handle.
   
   The comments you have responded on are indeed insufficiently described,
   hopefully my new answers cover your questions.


-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

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: dev-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Baunsgaard commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r778323021



##########
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:
       Oh, i read it wrongly, i quickly saw it as a matrix you filled with a specific value.
   for now lets keep it as it is, unless you observe that there is some bottleneck related to Q_n in the first iteration of the algorithm.




-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Lampisan commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

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: dev-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Lampisan edited a comment on pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Lampisan edited a comment on pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#issuecomment-1004710225


   Hi @Baunsgaard and @foebe (Nikolaus),
   Thanks for your suggestions, we added a few questions to your replies.
   We will push the changes after resolving those.
   


-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org

[GitHub] [systemds] Baunsgaard commented on a change in pull request #1489: [SYSTEMDS-3167] Alternative Implementations of Eigen Decomposition

[GitBox <gi...@apache.org>]

Baunsgaard commented on a change in pull request #1489:
URL: https://github.com/apache/systemds/pull/1489#discussion_r778333985



##########
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:
       very unclear comment from my side, and completely wrong, it seems like i was getting worse over time while reviewing.
   
   My intention was to ignore the diagonal values in the matrix, since they don't have an effect.
   but the code already does this, and i was just a bit to quick in my comment.
   
   That said you could consider changing the implementation to not call rand, and instead generate the values yourself since the current implementation goes through the matrix twice.
   If you don't want to (this is only a test generator anyway) then leave a TODO comment.
   
   




-- 
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-unsubscribe@systemds.apache.org

For queries about this service, please contact Infrastructure at:
users@infra.apache.org