commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Peter <ppl...@gmail.com>
Subject [PATCH] Cholesky decomposition
Date Mon, 13 Dec 2004 21:31:15 GMT
This method implements the Cholesky decomposition of a RealMatrix. It
can be easily modified to be incorporated under the RealMatrix class and
to decompose BigMatrixes as well.
I hope you place it in the next release.

Peter P. Lupo
Computer Science barchelor student. UFRJ-Brazil
http://www.dcc.ufrj.br
http://www.ufrj.br

    /**
     * Computes the cholesky decomposition of matrix. The return is a lower triangular matrix
that may be empty and also coincide with the original matrix. The method only uses the upper
right part of the matrix, which must be square, symmetric and positive-definite by definition.
     * @param matrix
     *            The matrix to decompose.
     * @return
     *            The matrix with the result, which is lower triangular and has the same dimensions
of the original matrix.
     * @exception InvalidMatrixException
     *                if matrix is not square, symmetric or positive-definite.
     */
    public static RealMatrix choleskyDecomposition(RealMatrix matrix) {
        if (!matrix.isSquare()) throw new InvalidMatrixException("Not a
square/symmetric matrix."); //"matrix" is not square
        if (!matrix.equals(matrix.transpose())) throw new
InvalidMatrixException("Not a symmetric matrix."); //"matrix" is not square
        double[][] matrixData = matrix.copy().getData();
        double[][] cholesky = new
double[matrix.getRowDimension()][matrix.getColumnDimension()];
        for (int i = 0; i < matrix.getColumnDimension(); i++) {
            double[] ithRowOfL = cholesky[i];
            double[] ithRowOfA = matrixData[i];
            for (int j = i; j < matrix.getRowDimension(); j++) {
                double[] jthRowOfL = cholesky[j];
                double sum = ithRowOfA[j];
                for (int k = i - 1; k >= 0; k--){
                    sum -= ithRowOfL[k] * jthRowOfL[k];
                }
                if (i == j) {
                    if (sum <= 0)throw new InvalidMatrixException("Not a
positive definite matrix."); //"matrix", with rounding errors is not
positive definite
                    if (cholesky != null){
                        cholesky[j][i] = Math.sqrt(sum);
                    }
                } else {
                    if (cholesky != null){
                        cholesky[j][i] = sum / cholesky[i][i];
                    }
                    if (i < j){
                        cholesky[i][j] = 0;
                    }
                }
            }
        }
        return new RealMatrixImpl(cholesky);
    }


---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message