001package algs91; // section 9.9
002import stdlib.*;
003/* ***********************************************************************
004 *  Compilation:  javac Cholesky.java
005 *  Execution:    java Cholesky
006 *
007 *  Compute XCholesky decomposition of symmetric positive definite
008 *  matrix A = LL^T.
009 *
010 *  % java Cholesky
011 *  2.00000  0.00000  0.00000
012 *  0.50000  2.17945  0.00000
013 *  0.50000  1.26179  3.62738
014 *
015 *************************************************************************/
016
017public class XCholesky {
018        //private static final double EPSILON = 1e-10;
019
020        // is symmetric
021        public static boolean isSymmetric(double[][] A) {
022                int N = A.length;
023                for (int i = 0; i < N; i++) {
024                        for (int j = 0; j < i; j++) {
025                                if (A[i][j] != A[j][i]) return false;
026                        }
027                }
028                return true;
029        }
030
031        // is symmetric
032        public static boolean isSquare(double[][] A) {
033                int N = A.length;
034                for (int i = 0; i < N; i++) {
035                        if (A[i].length != N) return false;
036                }
037                return true;
038        }
039
040
041        // return XCholesky factor L of psd matrix A = L L^T
042        public static double[][] cholesky(double[][] A) {
043                if (!isSquare(A)) {
044                        throw new Error("Matrix is not square");
045                }
046                if (!isSymmetric(A)) {
047                        throw new Error("Matrix is not symmetric");
048                }
049
050                int N  = A.length;
051                double[][] L = new double[N][N];
052
053                for (int i = 0; i < N; i++)  {
054                        for (int j = 0; j <= i; j++) {
055                                double sum = 0.0;
056                                for (int k = 0; k < j; k++) {
057                                        sum += L[i][k] * L[j][k];
058                                }
059                                if (i == j) L[i][i] = Math.sqrt(A[i][i] - sum);
060                                else        L[i][j] = 1.0 / L[j][j] * (A[i][j] - sum);
061                        }
062                        if (L[i][i] <= 0) {
063                                throw new Error("Matrix not positive definite");
064                        }
065                }
066                return L;
067        }
068
069
070        // sample client
071        public static void main(String[] args) {
072                int N = 3;
073                double[][] A = {
074                                { 4, 1,  1 },
075                                { 1, 5,  3 },
076                                { 1, 3, 15 }
077                };
078                double[][] L = cholesky(A);
079                for (int i = 0; i < N; i++) {
080                        for (int j = 0; j < N; j++) {
081                                StdOut.format("%8.5f ", L[i][j]);
082                        }
083                        StdOut.println();
084                }
085
086        }
087
088}