``` 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 ``` ```package algs9; // section 9.9 import stdlib.*; /* *********************************************************************** * Compilation: javac Cholesky.java * Execution: java Cholesky * * Compute XCholesky decomposition of symmetric positive definite * matrix A = LL^T. * * % java Cholesky * 2.00000 0.00000 0.00000 * 0.50000 2.17945 0.00000 * 0.50000 1.26179 3.62738 * *************************************************************************/ public class XCholesky { //private static final double EPSILON = 1e-10; // is symmetric public static boolean isSymmetric(double[][] A) { int N = A.length; for (int i = 0; i < N; i++) { for (int j = 0; j < i; j++) { if (A[i][j] != A[j][i]) return false; } } return true; } // is symmetric public static boolean isSquare(double[][] A) { int N = A.length; for (int i = 0; i < N; i++) { if (A[i].length != N) return false; } return true; } // return XCholesky factor L of psd matrix A = L L^T public static double[][] cholesky(double[][] A) { if (!isSquare(A)) { throw new Error("Matrix is not square"); } if (!isSymmetric(A)) { throw new Error("Matrix is not symmetric"); } int N = A.length; double[][] L = new double[N][N]; for (int i = 0; i < N; i++) { for (int j = 0; j <= i; j++) { double sum = 0.0; for (int k = 0; k < j; k++) { sum += L[i][k] * L[j][k]; } if (i == j) L[i][i] = Math.sqrt(A[i][i] - sum); else L[i][j] = 1.0 / L[j][j] * (A[i][j] - sum); } if (L[i][i] <= 0) { throw new Error("Matrix not positive definite"); } } return L; } // sample client public static void main(String[] args) { int N = 3; double[][] A = { { 4, 1, 1 }, { 1, 5, 3 }, { 1, 3, 15 } }; double[][] L = cholesky(A); for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) { StdOut.format("%8.5f ", L[i][j]); } StdOut.println(); } } } ```