Compute determinant(A), where A = L*L^T via Cholesky decomposition is a symmetric, positive definite matrix. - Java java.lang

Java examples for java.lang:Math Matrix

  1. HOME
  2. Java
  3. java.lang
  4. Math Matrix

Description

Compute determinant(A), where A = L*L^T via Cholesky decomposition is a symmetric, positive definite matrix.

Demo Code

/*//from ww  w.j av a  2  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{
    /**
     * <p>Compute determinant(A), where A = L*L^T via Cholesky decomposition
     * is a symmetric, positive definite matrix.
     * </p>
     * 
     * <p>This is an efficient computation, since det(A) = det(L*L^T) = det(L)*det(L^T)
     * and L is triangular so det(L) is just the product along the diagonal
     * </p>
     * 
     * @param A symmetric, positive definite matrix to take the determinant of
     * @return the determinant of A
     * @throws Exception when the matrix is not symmetric or positive definite
     */
    public static double determinantSymmPosDefMatrix(double[][] A)
            throws Exception {
        double[][] L = CholeskyDecomposition(A);

        return determinantViaCholeskyResult(L);
    }
    /**
     * <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>Compute determinant(A), where A = L*L^T via Cholesky decomposition
     * is a symmetric, positive definite matrix.
     * </p>
     * 
     * <p>This is an efficient computation, since det(A) = det(L*L^T) = det(L)*det(L^T)
     * and L is triangular so det(L) is just the product along the diagonal
     * </p>
     * 
     * @param L Cholesky decomposition L of the symmetric positive definite matrix A
     * @return det(A)
     * @see #CholeskyDecomposition(double[][])
     */
    public static double determinantViaCholeskyResult(double[][] L) {
        double detL = 1.0;
        int n = L.length;
        for (int i = 0; i < n; i++) {
            detL *= L[i][i];
        }
        return detL * detL;
    }
}

Related Tutorials