Compute matrix inversion of a symmetric, positive definite matrix by using the Cholesky Decomposition L of the matrix A. - Java java.lang

Description

Compute matrix inversion of a symmetric, positive definite matrix by using the Cholesky Decomposition L of the matrix A.

Demo Code

/*/*from ww w  .j a  v  a2  s. c  o  m*/
 *  Java Information Dynamics Toolkit (JIDT)
 *  Copyright (C) 2012, Joseph T. Lizier
 *  
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *  
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *  
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
import java.io.PrintStream;
import java.util.Arrays;
import java.util.Comparator;
import java.util.Vector;

public class Main{
    /**
     * Compute matrix inversion of a symmetric, positive definite matrix
     *  by using the Cholesky Decomposition L of the matrix A.
     * Since A = L L^T, then A^-1 = (L^T)^-1 L^-1, and the inverses
     *  of 
     * 
     * @param A matrix to be inverted
     * @return the inverse of A
     * @throws Exception when the matrix is not symmetric or positive definite
     * @see #CholeskyDecomposition(double[][])
     */
    public static double[][] invertSymmPosDefMatrix(double[][] A)
            throws Exception {
        // First do the Cholesky Decomposition:

        double[][] L = CholeskyDecomposition(A);

        return solveViaCholeskyResult(L, identityMatrix(A.length));
    }
    /**
     * <p>Make the Cholesky decomposition L of a given input matrix A,
     * where:
     *    <ol>
     *       <li>A is symmetric and positive definite (has full rank)</li>
     *       <li>A = L L^T (L^T is the transpose of L - here A has real
     *       entries only, though a Cholesky decomposition is possible
     *       with complex entries)</li>
     *       <li>L is a lower triangular matrix</li>
     *    </ol>
     * We perform the decomposition using the Cholesky?Banachiewicz
     *  algorithm, computing L from the top left, row by row (see wikipedia)
     * </p>
     * 
     * <p>This method has been adapted from the JAMA project (public domain software)
     * </p>
     * 
     * @param A input matrix
     * @return L 
     * @throws Exception when the matrix A is not square, is asymmetric, or
     *       not positive definite
     * @see {@link en.wikipedia.org/wiki/Cholesky_decomposition}
     * @see {@link http://mathworld.wolfram.com/CholeskyDecomposition.html}
     * @see {@link http://en.wikipedia.org/wiki/Positive-definite_matrix}
     * @see {@link http://math.nist.gov/javanumerics/jama/}
     * @see {@link http://www2.gsu.edu/~mkteer/npdmatri.html}
     */
    public static double[][] CholeskyDecomposition(double[][] A)
            throws Exception {
        int n = A.length;
        double[][] L = new double[n][n];
        // Loop over all rows:
        for (int j = 0; j < n; j++) {
            // Check length of row keeps this a square matrix:
            if (A[j].length != n) {
                throw new Exception(
                        "CholeskyDecomposition is only performed on square matrices");
            }
            double d = 0.0;
            for (int k = 0; k < j; k++) {
                double s = 0.0;
                for (int i = 0; i < k; i++) {
                    s += L[k][i] * L[j][i];
                }
                L[j][k] = s = (A[j][k] - s) / L[k][k];
                d = d + s * s;
                // Check that these matrix entries remain symmetric:
                if (A[k][j] != A[j][k]) {
                    throw new Exception(
                            "CholeskyDecomposition is only performed on symmetric matrices");
                }
            }
            d = A[j][j] - d;
            // Check the positive definite condition:
            if (d <= 0.0) {
                // Throw an error with some suggestions. The last suggestion is from my observations
                //  from a simple test with Matlab - I should find a reference for this ...
                throw new NonPositiveDefiniteMatrixException(
                        "CholeskyDecomposition is only performed on positive-definite matrices. "
                                + "Some reasons for non-positive-definite matrix are listed at http://www2.gsu.edu/~mkteer/npdmatri.html - "
                                + "note: a correlation matrix is non-positive-definite if you have more variables than observations");
            }
            L[j][j] = Math.sqrt(d);
            // Set the upper triangular part to all zeros:
            for (int k = j + 1; k < n; k++) {
                L[j][k] = 0.0;
            }
        }
        return L;
    }
    /**
     * <p>Solve A*X = B, where A = L*L^T via Cholesky decomposition.
     * </p>
     * 
     * <p>This method has been adapted from the JAMA project (public domain software)
     * </p>
     *
     * @param L Cholesky decomposition of the matrix A
     * @param B matrix with as many rows as A and any number of columns
     * @return X so that A*X = B
     * @see {@link http://math.nist.gov/javanumerics/jama/}
     * @see #CholeskyDecomposition(double[][])
     */
    public static double[][] solveViaCholeskyResult(double[][] L,
            double[][] B) {
        int aRows = L.length;
        if (aRows != B.length) {
            throw new IllegalArgumentException(
                    "Matrix row dimensions must agree.");
        }

        // Copy B matrix 
        double[][] X = MatrixUtils.arrayCopy(B);
        int bCols = B[0].length;

        // Solve L*Y = B;
        for (int k = 0; k < aRows; k++) {
            for (int j = 0; j < bCols; j++) {
                for (int i = 0; i < k; i++) {
                    X[k][j] -= X[i][j] * L[k][i];
                }
                X[k][j] /= L[k][k];
            }
        }

        // Solve L'*X = Y;
        for (int k = aRows - 1; k >= 0; k--) {
            for (int j = 0; j < bCols; j++) {
                for (int i = k + 1; i < aRows; i++) {
                    X[k][j] -= X[i][j] * L[i][k];
                }
                X[k][j] /= L[k][k];
            }
        }
        return X;
    }
    /**
     * Generate the identity matrix of the given size
     * 
     * @param size (size along one dimension)
     * @return two dimensional double array representing the identity matrix
     */
    public static double[][] identityMatrix(int size) {
        double[][] I = new double[size][size];
        for (int r = 0; r < size; r++) {
            I[r][r] = 1.0;
        }
        return I;
    }
    /**
     * Copies all rows and columns between two double arrays
     * 
     * @param src
     * @param dest
     */
    public static void arrayCopy(double[][] src, double[][] dest) {
        for (int r = 0; r < src.length; r++) {
            System.arraycopy(src[r], 0, dest[r], 0, src[r].length);
        }
    }
    /**
     * Copies all rows and columns between two double arrays
     * 
     * @param src
     * @param dest
     */
    public static double[][] arrayCopy(double[][] src) {
        double[][] dest = new double[src.length][];
        for (int r = 0; r < src.length; r++) {
            dest[r] = new double[src[r].length];
            System.arraycopy(src[r], 0, dest[r], 0, src[r].length);
        }
        return dest;
    }
    /**
     * Copies the required rows and columns between two 
     * double arrays
     * 
     * @param src
     * @param srcStartRow
     * @param srcStartCol
     * @param dest
     * @param destStartRow
     * @param destStartCol
     * @param rows
     * @param cols
     */
    public static void arrayCopy(double[][] src, int srcStartRow,
            int srcStartCol, double[][] dest, int destStartRow,
            int destStartCol, int rows, int cols) {

        for (int r = 0; r < rows; r++) {
            System.arraycopy(src[srcStartRow + r], srcStartCol,
                    dest[destStartRow + r], destStartCol, cols);
        }
    }
    /**
     * Copies all rows and columns between two int arrays
     * 
     * @param src
     * @param dest
     */
    public static void arrayCopy(int[][] src, int[][] dest) {
        for (int r = 0; r < src.length; r++) {
            System.arraycopy(src[r], 0, dest[r], 0, src[r].length);
        }
    }
    /**
     * Copies the required rows and columns between two 
     * double arrays
     * 
     * @param src
     * @param srcStartRow
     * @param srcStartCol
     * @param dest
     * @param destStartRow
     * @param destStartCol
     * @param rows
     * @param cols
     */
    public static void arrayCopy(int[][] src, int srcStartRow,
            int srcStartCol, int[][] dest, int destStartRow,
            int destStartCol, int rows, int cols) {

        for (int r = 0; r < rows; r++) {
            System.arraycopy(src[srcStartRow + r], srcStartCol,
                    dest[destStartRow + r], destStartCol, cols);
        }
    }
}

• Previous
• Next

Related Tutorials

• Compute the covariance matrix between all column pairs (variables) in the multivariate data set
• Returns the determinant of the input matrix.
• Compute determinant(A), where A = L*L^T via Cholesky decomposition is a symmetric, positive definite matrix.
• Works out the index of the k minimum values in the matrix in a given column
• Works out the maximum value in the matrix in a given column