1 package Jama;
2
3 /** Cholesky Decomposition.
4 <P>
5 For a symmetric, positive definite matrix A, the Cholesky decomposition
6 is an lower triangular matrix L so that A = L*L'.
7 <P>
8 If the matrix is not symmetric or positive definite, the constructor
9 returns a partial decomposition and sets an internal flag that may
10 be queried by the isSPD() method.
11 @version $Id: CholeskyDecomposition.java,v 1.1.1.1 2010/11/30 21:31:59 jeremy Exp $
12 */
13
14 public class CholeskyDecomposition implements java.io.Serializable {
15
16 /* ------------------------
17 Class variables
18 * ------------------------ */
19
20 /** Array for internal storage of decomposition.
21 @serial internal array storage.
22 */
23 private double[][] L;
24
25 /** Row and column dimension (square matrix).
26 @serial matrix dimension.
27 */
28 private int n;
29
30 /** Symmetric and positive definite flag.
31 @serial is symmetric and positive definite flag.
32 */
33 private boolean isspd;
34
35 /* ------------------------
36 Constructor
37 * ------------------------ */
38
39 /** Cholesky algorithm for symmetric and positive definite matrix.
40 @param Arg Square, symmetric matrix.
41 // Structure to access L and isspd flag.
42 */
43
44 public CholeskyDecomposition (Matrix Arg) {
45 // Initialize.
46 double[][] A = Arg.getArray();
47 n = Arg.getRowDimension();
48 L = new double[n][n];
49 isspd = (Arg.getColumnDimension() == n);
50 // Main loop.
51 for (int j = 0; j < n; j++) {
52 double[] Lrowj = L[j];
53 double d = 0.0;
54 for (int k = 0; k < j; k++) {
55 double[] Lrowk = L[k];
56 double s = 0.0;
57 for (int i = 0; i < k; i++) {
58 s += Lrowk[i]*Lrowj[i];
59 }
60 Lrowj[k] = s = (A[j][k] - s)/L[k][k];
61 d = d + s*s;
62 isspd = isspd & (A[k][j] == A[j][k]);
63 }
64 d = A[j][j] - d;
65 isspd = isspd & (d > 0.0);
66 L[j][j] = Math.sqrt(Math.max(d,0.0));
67 for (int k = j+1; k < n; k++) {
68 L[j][k] = 0.0;
69 }
70 }
71 }
72
73 /* ------------------------
74 Temporary, experimental code.
75 * ------------------------ *\
76
77 \** Right Triangular Cholesky Decomposition.
78 <P>
79 For a symmetric, positive definite matrix A, the Right Cholesky
80 decomposition is an upper triangular matrix R so that A = R'*R.
81 This constructor computes R with the Fortran inspired column oriented
82 algorithm used in LINPACK and MATLAB. In Java, we suspect a row oriented,
83 lower triangular decomposition is faster. We have temporarily included
84 this constructor here until timing experiments confirm this suspicion.
85 *\
86
87 \** Array for internal storage of right triangular decomposition. **\
88 private transient double[][] R;
89
90 \** Cholesky algorithm for symmetric and positive definite matrix.
91 @param A Square, symmetric matrix.
92 @param rightflag Actual value ignored.
93 @return Structure to access R and isspd flag.
94 *\
95
96 public CholeskyDecomposition (Matrix Arg, int rightflag) {
97 // Initialize.
98 double[][] A = Arg.getArray();
99 n = Arg.getColumnDimension();
100 R = new double[n][n];
101 isspd = (Arg.getColumnDimension() == n);
102 // Main loop.
103 for (int j = 0; j < n; j++) {
104 double d = 0.0;
105 for (int k = 0; k < j; k++) {
106 double s = A[k][j];
107 for (int i = 0; i < k; i++) {
108 s = s - R[i][k]*R[i][j];
109 }
110 R[k][j] = s = s/R[k][k];
111 d = d + s*s;
112 isspd = isspd & (A[k][j] == A[j][k]);
113 }
114 d = A[j][j] - d;
115 isspd = isspd & (d > 0.0);
116 R[j][j] = Math.sqrt(Math.max(d,0.0));
117 for (int k = j+1; k < n; k++) {
118 R[k][j] = 0.0;
119 }
120 }
121 }
122
123 \** Return upper triangular factor.
124 @return R
125 *\
126
127 public Matrix getR () {
128 return new Matrix(R,n,n);
129 }
130
131 \* ------------------------
132 End of temporary code.
133 * ------------------------ */
134
135 /* ------------------------
136 Public Methods
137 * ------------------------ */
138
139 /** Is the matrix symmetric and positive definite?
140 @return true if A is symmetric and positive definite.
141 */
142
143 public boolean isSPD () {
144 return isspd;
145 }
146
147 /** Return triangular factor.
148 @return L
149 */
150
151 public Matrix getL () {
152 return new Matrix(L,n,n);
153 }
154
155 /** Solve A*X = B
156 @param B A Matrix with as many rows as A and any number of columns.
157 @return X so that L*L'*X = B
158 @exception IllegalArgumentException Matrix row dimensions must agree.
159 @exception RuntimeException Matrix is not symmetric positive definite.
160 */
161
162 public Matrix solve (Matrix B) {
163 if (B.getRowDimension() != n) {
164 throw new IllegalArgumentException("Matrix row dimensions must agree.");
165 }
166 if (!isspd) {
167 throw new RuntimeException("Matrix is not symmetric positive definite.");
168 }
169
170 // Copy right hand side.
171 double[][] X = B.getArrayCopy();
172 int nx = B.getColumnDimension();
173
174 // Solve L*Y = B;
175 for (int k = 0; k < n; k++) {
176 for (int i = k+1; i < n; i++) {
177 for (int j = 0; j < nx; j++) {
178 X[i][j] -= X[k][j]*L[i][k];
179 }
180 }
181 for (int j = 0; j < nx; j++) {
182 X[k][j] /= L[k][k];
183 }
184 }
185
186 // Solve L'*X = Y;
187 for (int k = n-1; k >= 0; k--) {
188 for (int j = 0; j < nx; j++) {
189 X[k][j] /= L[k][k];
190 }
191 for (int i = 0; i < k; i++) {
192 for (int j = 0; j < nx; j++) {
193 X[i][j] -= X[k][j]*L[k][i];
194 }
195 }
196 }
197 return new Matrix(X,n,nx);
198 }
199 }