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