org.apache.commons.math3.linear.QRDecomposition Java Examples

/**
 * 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);
}
 

/**
 * Get the covariance matrix of the optimized parameters.
 * <br/>
 * Note that this operation involves the inversion of the
 * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
 * Jacobian matrix.
 * The {@code threshold} parameter is a way for the caller to specify
 * that the result of this computation should be considered meaningless,
 * and thus trigger an exception.
 *
 * @param threshold Singularity threshold.
 * @return the covariance matrix.
 * @throws org.apache.commons.math3.linear.SingularMatrixException
 * if the covariance matrix cannot be computed (singular problem).
 */
public double[][] getCovariances(double threshold) {
    // Set up the jacobian.
    updateJacobian();

    // Compute transpose(J)J, without building intermediate matrices.
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < rows; ++k) {
                sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    // Compute the covariances matrix.
    final DecompositionSolver solver
        = new QRDecomposition(MatrixUtils.createRealMatrix(jTj), threshold).getSolver();
    return solver.getInverse().getData();
}
 

/**
 * Get the covariance matrix of the optimized parameters.
 * <br/>
 * Note that this operation involves the inversion of the
 * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
 * Jacobian matrix.
 * The {@code threshold} parameter is a way for the caller to specify
 * that the result of this computation should be considered meaningless,
 * and thus trigger an exception.
 *
 * @param threshold Singularity threshold.
 * @return the covariance matrix.
 * @throws org.apache.commons.math3.linear.SingularMatrixException
 * if the covariance matrix cannot be computed (singular problem).
 */
public double[][] getCovariances(double threshold) {
    // Set up the jacobian.
    updateJacobian();

    // Compute transpose(J)J, without building intermediate matrices.
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < rows; ++k) {
                sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    // Compute the covariances matrix.
    final DecompositionSolver solver
        = new QRDecomposition(MatrixUtils.createRealMatrix(jTj), threshold).getSolver();
    return solver.getInverse().getData();
}
 

/**
 * Calculates the QR-decomposition of a matrix. The QR-decomposition of a
 * matrix A consists of two matrices Q and R that satisfy: A = QR, Q is
 * orthogonal (QTQ = I), and R is upper triangular. If A is m×n, Q is m×m
 * and R m×n.
 *
 * @param a Given matrix.
 * @return Result Q/R arrays.
 */
public static Array[] qr(Array a) {
    int m = a.getShape()[0];
    int n = a.getShape()[1];
    Array Qa = Array.factory(DataType.DOUBLE, new int[]{m, m});
    Array Ra = Array.factory(DataType.DOUBLE, a.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray_Double(a);
    RealMatrix matrix = new Array2DRowRealMatrix(aa, false);
    QRDecomposition decomposition = new QRDecomposition(matrix);
    RealMatrix Q = decomposition.getQ();
    RealMatrix R = decomposition.getR();
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            Qa.setDouble(i * m + j, Q.getEntry(i, j));
        }
    }
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            Ra.setDouble(i * n + j, R.getEntry(i, j));
        }
    }

    return new Array[]{Qa, Ra};
}
 

/**
 * Get the covariance matrix of the optimized parameters.
 * <br/>
 * Note that this operation involves the inversion of the
 * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the
 * Jacobian matrix.
 * The {@code threshold} parameter is a way for the caller to specify
 * that the result of this computation should be considered meaningless,
 * and thus trigger an exception.
 *
 * @param threshold Singularity threshold.
 * @return the covariance matrix.
 * @throws org.apache.commons.math3.linear.SingularMatrixException
 * if the covariance matrix cannot be computed (singular problem).
 */
public double[][] getCovariances(double threshold) {
    // Set up the jacobian.
    updateJacobian();

    // Compute transpose(J)J, without building intermediate matrices.
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < rows; ++k) {
                sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    // Compute the covariances matrix.
    final DecompositionSolver solver
        = new QRDecomposition(MatrixUtils.createRealMatrix(jTj), threshold).getSolver();
    return solver.getInverse().getData();
}
 

/**
 * 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);
}
 

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an n x m - matrix given as double[n][m]</li>
 * <li>b is an m - vector given as double[m],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 *
 * @param matrixA The matrix A (left hand side of the linear equation).
 * @param b The vector (right hand of the linear equation).
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquation(final double[][] matrixA, final double[] b) {

	if(isSolverUseApacheCommonsMath) {
		final Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(matrixA);

		DecompositionSolver solver;
		if(matrix.getColumnDimension() == matrix.getRowDimension()) {
			solver = new LUDecomposition(matrix).getSolver();
		}
		else {
			solver = new QRDecomposition(new Array2DRowRealMatrix(matrixA)).getSolver();
		}

		// Using SVD - very slow
		//			solver = new SingularValueDecomposition(new Array2DRowRealMatrix(A)).getSolver();

		return solver.solve(new Array2DRowRealMatrix(b)).getColumn(0);
	}
	else {
		return org.jblas.Solve.solve(new org.jblas.DoubleMatrix(matrixA), new org.jblas.DoubleMatrix(b)).data;

		// For use of colt:
		// cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra();
		// return linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray();

		// For use of parallel colt:
		// return new cern.colt.matrix.tdouble.algo.decomposition.DenseDoubleLUDecomposition(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix2D(A)).solve(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix1D(b)).toArray();
	}
}
 

/**
 * Computes model parameters and parameter variance using a QR decomposition of the X matrix
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBetaQR(RealVector y, RealMatrix x) {
    final int n = x.getRowDimension();
    final int p = x.getColumnDimension();
    final int offset = hasIntercept() ? 1 : 0;
    final QRDecomposition decomposition = new QRDecomposition(x, threshold);
    final RealVector betaVector = decomposition.getSolver().solve(y);
    final RealVector residuals = y.subtract(x.operate(betaVector));
    this.rss = residuals.dotProduct(residuals);
    this.errorVariance = rss / (n - p);
    this.stdError = Math.sqrt(errorVariance);
    this.residuals = createResidualsFrame(residuals);
    final RealMatrix rAug = decomposition.getR().getSubMatrix(0, p - 1, 0, p - 1);
    final RealMatrix rInv = new LUDecomposition(rAug).getSolver().getInverse();
    final RealMatrix covMatrix = rInv.multiply(rInv.transpose()).scalarMultiply(errorVariance);
    final RealMatrix result = new Array2DRowRealMatrix(p, 2);
    if (hasIntercept()) {
        result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
        result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
    }
    for (int i = 0; i < getRegressors().size(); i++) {
        final int index = i + offset;
        final double variance = covMatrix.getEntry(index, index);
        result.setEntry(index, 1, variance);
        result.setEntry(index, 0, betaVector.getEntry(index));
    }
    return result;
}
 

/**
 * 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());
}
 

@Test(dataProvider = "styles")
public void testPseudoInverse(DataFrameAlgebra.Lib lib, boolean parallel) {
    DataFrameAlgebra.LIBRARY.set(lib);
    final DataFrame<Integer,String> source = DataFrame.read().csv("./src/test/resources/pca/svd/poppet-svd-eigenvectors.csv");
    Array.of(20, 77, 95, 135, 233, 245).forEach(count -> {
        final DataFrame<Integer,String> frame = source.cols().select(col -> col.ordinal() < count);
        final DataFrame<Integer,Integer> inverse = frame.inverse();
        final RealMatrix matrix = new QRDecomposition(toMatrix(frame)).getSolver().getInverse();
        assertEquals(inverse, matrix);
    });
}
 

/**
 * 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 };
}
 

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an n x m - matrix given as double[n][m]</li>
 * <li>b is an m - vector given as double[m],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 *
 * @param matrixA The matrix A (left hand side of the linear equation).
 * @param b The vector (right hand of the linear equation).
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquation(final double[][] matrixA, final double[] b) {

	if(isSolverUseApacheCommonsMath) {
		final Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(matrixA);

		DecompositionSolver solver;
		if(matrix.getColumnDimension() == matrix.getRowDimension()) {
			solver = new LUDecomposition(matrix).getSolver();
		}
		else {
			solver = new QRDecomposition(new Array2DRowRealMatrix(matrixA)).getSolver();
		}

		// Using SVD - very slow
		//			solver = new SingularValueDecomposition(new Array2DRowRealMatrix(A)).getSolver();

		return solver.solve(new Array2DRowRealMatrix(b)).getColumn(0);
	}
	else {
		return org.jblas.Solve.solve(new org.jblas.DoubleMatrix(matrixA), new org.jblas.DoubleMatrix(b)).data;

		// For use of colt:
		// cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra();
		// return linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray();

		// For use of parallel colt:
		// return new cern.colt.matrix.tdouble.algo.decomposition.DenseDoubleLUDecomposition(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix2D(A)).solve(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix1D(b)).toArray();
	}
}
 

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an n x m - matrix given as double[n][m]</li>
 * <li>b is an m - vector given as double[m],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 * using a standard Tikhonov regularization, i.e., we solve in the least square sense
 *   A* x = b*
 * where A* = (A^T, lambda I)^T and b* = (b^T , 0)^T.
 *
 * @param matrixA The matrix A (left hand side of the linear equation).
 * @param b The vector (right hand of the linear equation).
 * @param lambda The parameter lambda of the Tikhonov regularization. Lambda effectively measures which small numbers are considered zero.
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquationTikonov(final double[][] matrixA, final double[] b, final double lambda) {
	if(lambda == 0) {
		return solveLinearEquationLeastSquare(matrixA, b);
	}

	/*
	 * The copy of the array is inefficient, but the use cases for this method are currently limited.
	 * And SVD is an alternative to this method.
	 */
	final int rows = matrixA.length;
	final int cols = matrixA[0].length;
	final double[][] matrixRegularized = new double[rows+cols][cols];
	final double[] bRegularized = new double[rows+cols];					// Note the JVM initializes arrays to zero.
	for(int i=0; i<rows; i++) {
		System.arraycopy(matrixA[i], 0, matrixRegularized[i], 0, cols);
	}
	System.arraycopy(b, 0, bRegularized, 0, rows);

	for(int j=0; j<cols; j++) {
		final double[] matrixRow = matrixRegularized[rows+j];

		matrixRow[j] = lambda;
	}


	//		return solveLinearEquationLeastSquare(matrixRegularized, bRegularized);
	final DecompositionSolver solver = new QRDecomposition(new Array2DRowRealMatrix(matrixRegularized, false)).getSolver();
	return solver.solve(new ArrayRealVector(bRegularized, false)).toArray();
}
 

/**
 * Assert that the given matrix is unitary.
 * @param m
 */
public static void assertUnitaryMatrix(final RealMatrix m){
    final RealMatrix mInv = new QRDecomposition(m).getSolver().getInverse();
    final RealMatrix mT = m.transpose();

    for (int i = 0; i < mInv.getRowDimension(); i ++) {
        for (int j = 0; j < mInv.getColumnDimension(); j ++) {
            Assert.assertEquals(mInv.getEntry(i, j), mT.getEntry(i, j), 1e-7);
        }
    }
}
 

/** {@inheritDoc} */
public RealMatrix getCovariances(double threshold) {
    // Set up the Jacobian.
    final RealMatrix j = this.getJacobian();

    // Compute transpose(J)J.
    final RealMatrix jTj = j.transpose().multiply(j);

    // Compute the covariances matrix.
    final DecompositionSolver solver
            = new QRDecomposition(jTj, threshold).getSolver();
    return solver.getInverse();
}
 

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an n x m - matrix given as double[n][m]</li>
 * <li>b is an m - vector given as double[m],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 * using a standard Tikhonov regularization, i.e., we solve in the least square sense
 *   A* x = b*
 * where A* = (A^T, lambda I)^T and b* = (b^T , 0)^T.
 *
 * @param matrixA The matrix A (left hand side of the linear equation).
 * @param b The vector (right hand of the linear equation).
 * @param lambda The parameter lambda of the Tikhonov regularization. Lambda effectively measures which small numbers are considered zero.
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquationTikonov(final double[][] matrixA, final double[] b, final double lambda) {
	if(lambda == 0) {
		return solveLinearEquationLeastSquare(matrixA, b);
	}

	/*
	 * The copy of the array is inefficient, but the use cases for this method are currently limited.
	 * And SVD is an alternative to this method.
	 */
	final int rows = matrixA.length;
	final int cols = matrixA[0].length;
	final double[][] matrixRegularized = new double[rows+cols][cols];
	final double[] bRegularized = new double[rows+cols];					// Note the JVM initializes arrays to zero.
	for(int i=0; i<rows; i++) {
		System.arraycopy(matrixA[i], 0, matrixRegularized[i], 0, cols);
	}
	System.arraycopy(b, 0, bRegularized, 0, rows);

	for(int j=0; j<cols; j++) {
		final double[] matrixRow = matrixRegularized[rows+j];

		matrixRow[j] = lambda;
	}


	//		return solveLinearEquationLeastSquare(matrixRegularized, bRegularized);
	final DecompositionSolver solver = new QRDecomposition(new Array2DRowRealMatrix(matrixRegularized, false)).getSolver();
	return solver.solve(new ArrayRealVector(bRegularized, false)).toArray();
}
 

@Override
public QRDecompositionResult apply(DoubleMatrix x) {
  ArgChecker.notNull(x, "x");
  RealMatrix temp = CommonsMathWrapper.wrap(x);
  QRDecomposition qr = new QRDecomposition(temp);
  return new QRDecompositionCommonsResult(qr);
}
 

/**
 * 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 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());
}
 

/**
 * Creates an instance.
 * 
 * @param qr The result of the QR decomposition, not null
 */
public QRDecompositionCommonsResult(QRDecomposition qr) {
  ArgChecker.notNull(qr, "qr");
  _q = CommonsMathWrapper.unwrap(qr.getQ());
  _r = CommonsMathWrapper.unwrap(qr.getR());
  _qTranspose = _q.transpose();
  _solver = qr.getSolver();
}
 

/**
 * Check that the given matrix is unitary.
 */
public static boolean isUnitaryMatrix(final RealMatrix m){
    //Note can't use MatrixUtils.inverse because m may not be square
    final RealMatrix mInv = new QRDecomposition(m).getSolver().getInverse();
    final RealMatrix mT = m.transpose();

    for (int i = 0; i < mInv.getRowDimension(); i ++) {
        for (int j = 0; j < mInv.getColumnDimension(); j ++) {
            if (Math.abs(mInv.getEntry(i, j) - mT.getEntry(i, j)) > 1.0e-7){
                return false;
            }
        }
    }
    return true;
}
 

/** {@inheritDoc} */
public RealMatrix getCovariances(double threshold) {
    // Set up the Jacobian.
    final RealMatrix j = this.getJacobian();

    // Compute transpose(J)J.
    final RealMatrix jTj = j.transpose().multiply(j);

    // Compute the covariances matrix.
    final DecompositionSolver solver
            = new QRDecomposition(jTj, threshold).getSolver();
    return solver.getInverse();
}
 

/**
 * {@inheritDoc}
 * <p>This implementation computes and caches the QR decomposition of the X matrix.</p>
 */
@Override
public void newSampleData(double[] data, int nobs, int nvars) {
    super.newSampleData(data, nobs, nvars);
    qr = new QRDecomposition(getX());
}
 

/**
 * Calculates pseudo inverse of a matrix using QR decomposition
 * @param arr the array to invert
 * @return the pseudo inverted matrix
 */
public static INDArray pinvert(INDArray arr, boolean inPlace) {

    // TODO : do it natively instead of relying on commons-maths

    RealMatrix realMatrix = CheckUtil.convertToApacheMatrix(arr);
    QRDecomposition decomposition = new QRDecomposition(realMatrix, 0);
    DecompositionSolver solver = decomposition.getSolver();

    if (!solver.isNonSingular()) {
        throw new IllegalArgumentException("invalid array: must be singular matrix");
    }

    RealMatrix pinvRM = solver.getInverse();

    INDArray pseudoInverse = CheckUtil.convertFromApacheMatrix(pinvRM, arr.dataType());

    if (inPlace)
        arr.assign(pseudoInverse);
    return pseudoInverse;

}
 

/**
 * {@inheritDoc}
 * <p>This implementation computes and caches the QR decomposition of the X matrix
 * once it is successfully loaded.</p>
 */
@Override
protected void newXSampleData(double[][] x) {
    super.newXSampleData(x);
    qr = new QRDecomposition(getX());
}
 

/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {

    final ConvergenceChecker<PointVectorValuePair> checker
        = getConvergenceChecker();

    // iterate until convergence is reached
    PointVectorValuePair current = null;
    int iter = 0;
    for (boolean converged = false; !converged;) {
        ++iter;

        // evaluate the objective function and its jacobian
        PointVectorValuePair previous = current;
        updateResidualsAndCost();
        updateJacobian();
        current = new PointVectorValuePair(point, objective);

        final double[] targetValues = getTargetRef();
        final double[] residualsWeights = getWeightRef();

        // build the linear problem
        final double[]   b = new double[cols];
        final double[][] a = new double[cols][cols];
        for (int i = 0; i < rows; ++i) {

            final double[] grad   = weightedResidualJacobian[i];
            final double weight   = residualsWeights[i];
            final double residual = objective[i] - targetValues[i];

            // compute the normal equation
            final double wr = weight * residual;
            for (int j = 0; j < cols; ++j) {
                b[j] += wr * grad[j];
            }

            // build the contribution matrix for measurement i
            for (int k = 0; k < cols; ++k) {
                double[] ak = a[k];
                double wgk = weight * grad[k];
                for (int l = 0; l < cols; ++l) {
                    ak[l] += wgk * grad[l];
                }
            }
        }

        try {
            // solve the linearized least squares problem
            RealMatrix mA = new BlockRealMatrix(a);
            DecompositionSolver solver = useLU ?
                    new LUDecomposition(mA).getSolver() :
                    new QRDecomposition(mA).getSolver();
            final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
            // update the estimated parameters
            for (int i = 0; i < cols; ++i) {
                point[i] += dX[i];
            }
        } catch (SingularMatrixException e) {
            throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
        }

        // check convergence
        if (checker != null) {
            if (previous != null) {
                converged = checker.converged(iter, previous, current);
            }
        }
    }
    // we have converged
    return current;
}
 

/**
 * {@inheritDoc}
 * <p>This implementation computes and caches the QR decomposition of the X matrix.</p>
 */
@Override
public void newSampleData(double[] data, int nobs, int nvars) {
    super.newSampleData(data, nobs, nvars);
    qr = new QRDecomposition(getX());
}
 

/**
 * {@inheritDoc}
 * <p>This implementation computes and caches the QR decomposition of the X matrix
 * once it is successfully loaded.</p>
 */
@Override
protected void newXSampleData(double[][] x) {
    super.newXSampleData(x);
    qr = new QRDecomposition(getX());
}
 

/** {@inheritDoc} */
@Override
public PointVectorValuePair doOptimize() {
    final ConvergenceChecker<PointVectorValuePair> checker
        = getConvergenceChecker();

    // Computation will be useless without a checker (see "for-loop").
    if (checker == null) {
        throw new NullArgumentException();
    }

    final double[] targetValues = getTarget();
    final int nR = targetValues.length; // Number of observed data.

    final RealMatrix weightMatrix = getWeight();
    // Diagonal of the weight matrix.
    final double[] residualsWeights = new double[nR];
    for (int i = 0; i < nR; i++) {
        residualsWeights[i] = weightMatrix.getEntry(i, i);
    }

    final double[] currentPoint = getStartPoint();
    final int nC = currentPoint.length;

    // iterate until convergence is reached
    PointVectorValuePair current = null;
    int iter = 0;
    for (boolean converged = false; !converged;) {
        ++iter;

        // evaluate the objective function and its jacobian
        PointVectorValuePair previous = current;
        // Value of the objective function at "currentPoint".
        final double[] currentObjective = computeObjectiveValue(currentPoint);
        final double[] currentResiduals = computeResiduals(currentObjective);
        final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint);
        current = new PointVectorValuePair(currentPoint, currentObjective);

        // build the linear problem
        final double[]   b = new double[nC];
        final double[][] a = new double[nC][nC];
        for (int i = 0; i < nR; ++i) {

            final double[] grad   = weightedJacobian.getRow(i);
            final double weight   = residualsWeights[i];
            final double residual = currentResiduals[i];

            // compute the normal equation
            final double wr = weight * residual;
            for (int j = 0; j < nC; ++j) {
                b[j] += wr * grad[j];
            }

            // build the contribution matrix for measurement i
            for (int k = 0; k < nC; ++k) {
                double[] ak = a[k];
                double wgk = weight * grad[k];
                for (int l = 0; l < nC; ++l) {
                    ak[l] += wgk * grad[l];
                }
            }
        }

        try {
            // solve the linearized least squares problem
            RealMatrix mA = new BlockRealMatrix(a);
            DecompositionSolver solver = useLU ?
                    new LUDecomposition(mA).getSolver() :
                    new QRDecomposition(mA).getSolver();
            final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
            // update the estimated parameters
            for (int i = 0; i < nC; ++i) {
                currentPoint[i] += dX[i];
            }
        } catch (SingularMatrixException e) {
            throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
        }

        // Check convergence.
        if (previous != null) {
            converged = checker.converged(iter, previous, current);
            if (converged) {
                cost = computeCost(currentResiduals);
                // Update (deprecated) "point" field.
                point = current.getPoint();
                return current;
            }
        }
    }
    // Must never happen.
    throw new MathInternalError();
}
 

/**
 * Calculates pseudo inverse of a matrix using QR decomposition
 * @param arr the array to invert
 * @return the pseudo inverted matrix
 */
public static INDArray pinvert(INDArray arr, boolean inPlace) {

    // TODO : do it natively instead of relying on commons-maths

    RealMatrix realMatrix = CheckUtil.convertToApacheMatrix(arr);
    QRDecomposition decomposition = new QRDecomposition(realMatrix, 0);
    DecompositionSolver solver = decomposition.getSolver();

    if (!solver.isNonSingular()) {
        throw new IllegalArgumentException("invalid array: must be singular matrix");
    }

    RealMatrix pinvRM = solver.getInverse();

    INDArray pseudoInverse = CheckUtil.convertFromApacheMatrix(pinvRM);

    if (inPlace)
        arr.assign(pseudoInverse);
    return pseudoInverse;

}