src/main/java/org/apache/sysds/runtime/matrix/data/LibCommonsMath.java - systemds - Git at Google

 /*
  * 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.runtime.matrix.data;

 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
 import org.apache.commons.math3.linear.BlockRealMatrix;
 import org.apache.commons.math3.linear.CholeskyDecomposition;
 import org.apache.commons.math3.linear.DecompositionSolver;
 import org.apache.commons.math3.linear.EigenDecomposition;
 import org.apache.commons.math3.linear.LUDecomposition;
 import org.apache.commons.math3.linear.QRDecomposition;
 import org.apache.commons.math3.linear.RealMatrix;
 import org.apache.commons.math3.linear.SingularValueDecomposition;
 import org.apache.sysds.runtime.DMLRuntimeException;
 import org.apache.sysds.runtime.util.DataConverter;

 /**
  * Library for matrix operations that need invocation of
  * Apache Commons Math library.
  *
  * This library currently supports following operations:
  * matrix inverse, matrix decompositions (QR, LU, Eigen), solve
  */
 public class LibCommonsMath
 {
 	static final double RELATIVE_SYMMETRY_THRESHOLD = 1e-10;

 	private LibCommonsMath() {
 		//prevent instantiation via private constructor
 	}

 	public static boolean isSupportedUnaryOperation( String opcode ) {
 		return ( opcode.equals("inverse") || opcode.equals("cholesky") );
 	}

 	public static boolean isSupportedMultiReturnOperation( String opcode ) {
 		return ( opcode.equals("qr") || opcode.equals("lu") || opcode.equals("eigen") || opcode.equals("svd") );
 	}

 	public static boolean isSupportedMatrixMatrixOperation( String opcode ) {
 		return ( opcode.equals("solve") );
 	}

 	public static MatrixBlock unaryOperations(MatrixBlock inj, String opcode) {
 		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(inj);
 		if(opcode.equals("inverse"))
 			return computeMatrixInverse(matrixInput);
 		else if (opcode.equals("cholesky"))
 			return computeCholesky(matrixInput);
 		return null;
 	}

 	public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode) {
 		if(opcode.equals("qr"))
 			return computeQR(in);
 		else if (opcode.equals("lu"))
 			return computeLU(in);
 		else if (opcode.equals("eigen"))
 			return computeEigen(in);
 		else if ( opcode.equals("svd"))
 			return computeSvd(in);
 		return null;
 	}

 	public static MatrixBlock matrixMatrixOperations(MatrixBlock in1, MatrixBlock in2, String opcode) {
 		if(opcode.equals("solve")) {
 			if (in1.getNumRows() != in1.getNumColumns())
 				throw new DMLRuntimeException("The A matrix, in solve(A,b) should have squared dimensions.");
 			return computeSolve(in1, in2);
 		}
 		return null;
 	}

 	/**
 	 * Function to solve a given system of equations.
 	 *
 	 * @param in1 matrix object 1
 	 * @param in2 matrix object 2
 	 * @return matrix block
 	 */
 	private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) {
 		//convert to commons math BlockRealMatrix instead of Array2DRowRealMatrix
 		//to avoid unnecessary conversion as QR internally creates a BlockRealMatrix
 		BlockRealMatrix matrixInput = DataConverter.convertToBlockRealMatrix(in1);
 		BlockRealMatrix vectorInput = DataConverter.convertToBlockRealMatrix(in2);

 		/*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
 		DecompositionSolver lusolver = ludecompose.getSolver();
 		RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/

 		// Setup a solver based on QR Decomposition
 		QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
 		DecompositionSolver solver = qrdecompose.getSolver();
 		// Invoke solve
 		RealMatrix solutionMatrix = solver.solve(vectorInput);

 		return DataConverter.convertToMatrixBlock(solutionMatrix);
 	}

 	/**
 	 * Function to perform QR decomposition on a given matrix.
 	 *
 	 * @param in matrix object
 	 * @return array of matrix blocks
 	 */
 	private static MatrixBlock[] computeQR(MatrixBlock in) {
 		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

 		// Perform QR decomposition
 		QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
 		RealMatrix H = qrdecompose.getH();
 		RealMatrix R = qrdecompose.getR();

 		// Read the results into native format
 		MatrixBlock mbH = DataConverter.convertToMatrixBlock(H.getData());
 		MatrixBlock mbR = DataConverter.convertToMatrixBlock(R.getData());

 		return new MatrixBlock[] { mbH, mbR };
 	}

 	/**
 	 * Function to perform LU decomposition on a given matrix.
 	 *
 	 * @param in matrix object
 	 * @return array of matrix blocks
 	 */
 	private static MatrixBlock[] computeLU(MatrixBlock in) {
 		if ( in.getNumRows() != in.getNumColumns() ) {
 			throw new DMLRuntimeException("LU Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
 		}

 		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

 		// Perform LUP decomposition
 		LUDecomposition ludecompose = new LUDecomposition(matrixInput);
 		RealMatrix P = ludecompose.getP();
 		RealMatrix L = ludecompose.getL();
 		RealMatrix U = ludecompose.getU();

 		// Read the results into native format
 		MatrixBlock mbP = DataConverter.convertToMatrixBlock(P.getData());
 		MatrixBlock mbL = DataConverter.convertToMatrixBlock(L.getData());
 		MatrixBlock mbU = DataConverter.convertToMatrixBlock(U.getData());

 		return new MatrixBlock[] { mbP, mbL, mbU };
 	}

 	/**
 	 * Function to perform Eigen decomposition on a given matrix.
 	 * Input must be a symmetric matrix.
 	 *
 	 * @param in matrix object
 	 * @return array of matrix blocks
 	 */
 	private static MatrixBlock[] computeEigen(MatrixBlock in) {
 		if ( in.getNumRows() != in.getNumColumns() ) {
 			throw new DMLRuntimeException("Eigen Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
 		}

 		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

 		EigenDecomposition eigendecompose = new EigenDecomposition(matrixInput);
 		RealMatrix eVectorsMatrix = eigendecompose.getV();
 		double[][] eVectors = eVectorsMatrix.getData();
 		double[] eValues = eigendecompose.getRealEigenvalues();

 		//Sort the eigen values (and vectors) in increasing order (to be compatible w/ LAPACK.DSYEVR())
 		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 mbValues = DataConverter.convertToMatrixBlock(eValues, true);
 		MatrixBlock mbVectors = DataConverter.convertToMatrixBlock(eVectors);

 		return new MatrixBlock[] { mbValues, mbVectors };
 	}


 	/**
 	 * Performs Singular Value Decomposition. Calls Apache Commons Math SVD.
 	 * X = U * Sigma * Vt, where X is the input matrix,
 	 * U is the left singular matrix, Sigma is the singular values matrix returned as a
 	 * column matrix and Vt is the transpose of the right singular matrix V.
 	 * However, the returned array has  { U, Sigma, V}
 	 *
 	 * @param in Input matrix
 	 * @return An array containing U, Sigma & V
 	 */
 	private static MatrixBlock[] computeSvd(MatrixBlock in) {
 		Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

 		SingularValueDecomposition svd = new SingularValueDecomposition(matrixInput);
 		double[] sigma = svd.getSingularValues();
 		RealMatrix u = svd.getU();
 		RealMatrix v = svd.getV();
 		MatrixBlock U = DataConverter.convertToMatrixBlock(u.getData());
 		MatrixBlock Sigma = DataConverter.convertToMatrixBlock(sigma, true);
 		Sigma = LibMatrixReorg.diag(Sigma, new MatrixBlock(Sigma.rlen, Sigma.rlen, true));
 		MatrixBlock V = DataConverter.convertToMatrixBlock(v.getData());

 		return new MatrixBlock[] { U, Sigma, V };
 	}

 	/**
 	 * Function to compute matrix inverse via matrix decomposition.
 	 *
 	 * @param in commons-math3 Array2DRowRealMatrix
 	 * @return matrix block
 	 */
 	private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) {
 		if ( !in.isSquare() )
 			throw new DMLRuntimeException("Input to inv() must be square matrix -- given: a " + in.getRowDimension() + "x" + in.getColumnDimension() + " matrix.");

 		QRDecomposition qrdecompose = new QRDecomposition(in);
 		DecompositionSolver solver = qrdecompose.getSolver();
 		RealMatrix inverseMatrix = solver.getInverse();

 		return DataConverter.convertToMatrixBlock(inverseMatrix.getData());
 	}

 	/**
 	 * Function to compute Cholesky decomposition of the given input matrix.
 	 * The input must be a real symmetric positive-definite matrix.
 	 *
 	 * @param in commons-math3 Array2DRowRealMatrix
 	 * @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.");
 		CholeskyDecomposition cholesky = new CholeskyDecomposition(in, RELATIVE_SYMMETRY_THRESHOLD,
 			CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
 		RealMatrix rmL = cholesky.getL();
 		return DataConverter.convertToMatrixBlock(rmL.getData());
 	}
 }
	/*
	* 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.runtime.matrix.data;

	import org.apache.commons.math3.linear.Array2DRowRealMatrix;
	import org.apache.commons.math3.linear.BlockRealMatrix;
	import org.apache.commons.math3.linear.CholeskyDecomposition;
	import org.apache.commons.math3.linear.DecompositionSolver;
	import org.apache.commons.math3.linear.EigenDecomposition;
	import org.apache.commons.math3.linear.LUDecomposition;
	import org.apache.commons.math3.linear.QRDecomposition;
	import org.apache.commons.math3.linear.RealMatrix;
	import org.apache.commons.math3.linear.SingularValueDecomposition;
	import org.apache.sysds.runtime.DMLRuntimeException;
	import org.apache.sysds.runtime.util.DataConverter;

	/**
	* Library for matrix operations that need invocation of
	* Apache Commons Math library.
	*
	* This library currently supports following operations:
	* matrix inverse, matrix decompositions (QR, LU, Eigen), solve
	*/
	public class LibCommonsMath
	{
	static final double RELATIVE_SYMMETRY_THRESHOLD = 1e-10;

	private LibCommonsMath() {
	//prevent instantiation via private constructor
	}

	public static boolean isSupportedUnaryOperation( String opcode ) {
	return ( opcode.equals("inverse") \|\| opcode.equals("cholesky") );
	}

	public static boolean isSupportedMultiReturnOperation( String opcode ) {
	return ( opcode.equals("qr") \|\| opcode.equals("lu") \|\| opcode.equals("eigen") \|\| opcode.equals("svd") );
	}

	public static boolean isSupportedMatrixMatrixOperation( String opcode ) {
	return ( opcode.equals("solve") );
	}

	public static MatrixBlock unaryOperations(MatrixBlock inj, String opcode) {
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(inj);
	if(opcode.equals("inverse"))
	return computeMatrixInverse(matrixInput);
	else if (opcode.equals("cholesky"))
	return computeCholesky(matrixInput);
	return null;
	}

	public static MatrixBlock[] multiReturnOperations(MatrixBlock in, String opcode) {
	if(opcode.equals("qr"))
	return computeQR(in);
	else if (opcode.equals("lu"))
	return computeLU(in);
	else if (opcode.equals("eigen"))
	return computeEigen(in);
	else if ( opcode.equals("svd"))
	return computeSvd(in);
	return null;
	}

	public static MatrixBlock matrixMatrixOperations(MatrixBlock in1, MatrixBlock in2, String opcode) {
	if(opcode.equals("solve")) {
	if (in1.getNumRows() != in1.getNumColumns())
	throw new DMLRuntimeException("The A matrix, in solve(A,b) should have squared dimensions.");
	return computeSolve(in1, in2);
	}
	return null;
	}

	/**
	* Function to solve a given system of equations.
	*
	* @param in1 matrix object 1
	* @param in2 matrix object 2
	* @return matrix block
	*/
	private static MatrixBlock computeSolve(MatrixBlock in1, MatrixBlock in2) {
	//convert to commons math BlockRealMatrix instead of Array2DRowRealMatrix
	//to avoid unnecessary conversion as QR internally creates a BlockRealMatrix
	BlockRealMatrix matrixInput = DataConverter.convertToBlockRealMatrix(in1);
	BlockRealMatrix vectorInput = DataConverter.convertToBlockRealMatrix(in2);

	/*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
	DecompositionSolver lusolver = ludecompose.getSolver();
	RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/

	// Setup a solver based on QR Decomposition
	QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
	DecompositionSolver solver = qrdecompose.getSolver();
	// Invoke solve
	RealMatrix solutionMatrix = solver.solve(vectorInput);

	return DataConverter.convertToMatrixBlock(solutionMatrix);
	}

	/**
	* Function to perform QR decomposition on a given matrix.
	*
	* @param in matrix object
	* @return array of matrix blocks
	*/
	private static MatrixBlock[] computeQR(MatrixBlock in) {
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

	// Perform QR decomposition
	QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
	RealMatrix H = qrdecompose.getH();
	RealMatrix R = qrdecompose.getR();

	// Read the results into native format
	MatrixBlock mbH = DataConverter.convertToMatrixBlock(H.getData());
	MatrixBlock mbR = DataConverter.convertToMatrixBlock(R.getData());

	return new MatrixBlock[] { mbH, mbR };
	}

	/**
	* Function to perform LU decomposition on a given matrix.
	*
	* @param in matrix object
	* @return array of matrix blocks
	*/
	private static MatrixBlock[] computeLU(MatrixBlock in) {
	if ( in.getNumRows() != in.getNumColumns() ) {
	throw new DMLRuntimeException("LU Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
	}

	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

	// Perform LUP decomposition
	LUDecomposition ludecompose = new LUDecomposition(matrixInput);
	RealMatrix P = ludecompose.getP();
	RealMatrix L = ludecompose.getL();
	RealMatrix U = ludecompose.getU();

	// Read the results into native format
	MatrixBlock mbP = DataConverter.convertToMatrixBlock(P.getData());
	MatrixBlock mbL = DataConverter.convertToMatrixBlock(L.getData());
	MatrixBlock mbU = DataConverter.convertToMatrixBlock(U.getData());

	return new MatrixBlock[] { mbP, mbL, mbU };
	}

	/**
	* Function to perform Eigen decomposition on a given matrix.
	* Input must be a symmetric matrix.
	*
	* @param in matrix object
	* @return array of matrix blocks
	*/
	private static MatrixBlock[] computeEigen(MatrixBlock in) {
	if ( in.getNumRows() != in.getNumColumns() ) {
	throw new DMLRuntimeException("Eigen Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
	}

	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

	EigenDecomposition eigendecompose = new EigenDecomposition(matrixInput);
	RealMatrix eVectorsMatrix = eigendecompose.getV();
	double[][] eVectors = eVectorsMatrix.getData();
	double[] eValues = eigendecompose.getRealEigenvalues();

	//Sort the eigen values (and vectors) in increasing order (to be compatible w/ LAPACK.DSYEVR())
	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 mbValues = DataConverter.convertToMatrixBlock(eValues, true);
	MatrixBlock mbVectors = DataConverter.convertToMatrixBlock(eVectors);

	return new MatrixBlock[] { mbValues, mbVectors };
	}


	/**
	* Performs Singular Value Decomposition. Calls Apache Commons Math SVD.
	* X = U * Sigma * Vt, where X is the input matrix,
	* U is the left singular matrix, Sigma is the singular values matrix returned as a
	* column matrix and Vt is the transpose of the right singular matrix V.
	* However, the returned array has { U, Sigma, V}
	*
	* @param in Input matrix
	* @return An array containing U, Sigma & V
	*/
	private static MatrixBlock[] computeSvd(MatrixBlock in) {
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

	SingularValueDecomposition svd = new SingularValueDecomposition(matrixInput);
	double[] sigma = svd.getSingularValues();
	RealMatrix u = svd.getU();
	RealMatrix v = svd.getV();
	MatrixBlock U = DataConverter.convertToMatrixBlock(u.getData());
	MatrixBlock Sigma = DataConverter.convertToMatrixBlock(sigma, true);
	Sigma = LibMatrixReorg.diag(Sigma, new MatrixBlock(Sigma.rlen, Sigma.rlen, true));
	MatrixBlock V = DataConverter.convertToMatrixBlock(v.getData());

	return new MatrixBlock[] { U, Sigma, V };
	}

	/**
	* Function to compute matrix inverse via matrix decomposition.
	*
	* @param in commons-math3 Array2DRowRealMatrix
	* @return matrix block
	*/
	private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) {
	if ( !in.isSquare() )
	throw new DMLRuntimeException("Input to inv() must be square matrix -- given: a " + in.getRowDimension() + "x" + in.getColumnDimension() + " matrix.");

	QRDecomposition qrdecompose = new QRDecomposition(in);
	DecompositionSolver solver = qrdecompose.getSolver();
	RealMatrix inverseMatrix = solver.getInverse();

	return DataConverter.convertToMatrixBlock(inverseMatrix.getData());
	}

	/**
	* Function to compute Cholesky decomposition of the given input matrix.
	* The input must be a real symmetric positive-definite matrix.
	*
	* @param in commons-math3 Array2DRowRealMatrix
	* @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.");
	CholeskyDecomposition cholesky = new CholeskyDecomposition(in, RELATIVE_SYMMETRY_THRESHOLD,
	CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
	RealMatrix rmL = cholesky.getL();
	return DataConverter.convertToMatrixBlock(rmL.getData());
	}
	}