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}