1 #ifndef JAMA_CHOLESKY_H
2 #define JAMA_CHOLESKY_H
3
4 #include "math.h"
5 /* needed for sqrt() below. */
6
7
8 namespace JAMA
9 {
10
11 using namespace TNT;
12
13 /**
14 <P>
15 For a symmetric, positive definite matrix A, this function
16 computes the Cholesky factorization, i.e. it computes a lower
17 triangular matrix L such that A = L*L'.
18 If the matrix is not symmetric or positive definite, the function
19 computes only a partial decomposition. This can be tested with
20 the is_spd() flag.
21
22 <p>Typical usage looks like:
23 <pre>
24 Array2D<double> A(n,n);
25 Array2D<double> L;
26
27 ...
28
29 Cholesky<double> chol(A);
30
31 if (chol.is_spd())
32 L = chol.getL();
33
34 else
35 cout << "factorization was not complete.\n";
36
37 </pre>
38
39
40 <p>
41 (Adapted from JAMA, a Java Matrix Library, developed by jointly
42 by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
43
44 */
45
46 template <class Real>
47 class Cholesky
48 {
49 Array2D<Real> L_; // lower triangular factor
50 int isspd; // 1 if matrix to be factored was SPD
51
52 public:
53
54 Cholesky();
55 Cholesky(const Array2D<Real> &A);
56 Array2D<Real> getL() const;
57 Array1D<Real> solve(const Array1D<Real> &B);
58 Array2D<Real> solve(const Array2D<Real> &B);
59 int is_spd() const;
60
61 };
62
63 template <class Real>
Cholesky()64 Cholesky<Real>::Cholesky() : L_(0,0), isspd(0) {}
65
66 /**
67 @return 1, if original matrix to be factored was symmetric
68 positive-definite (SPD).
69 */
70 template <class Real>
is_spd()71 int Cholesky<Real>::is_spd() const
72 {
73 return isspd;
74 }
75
76 /**
77 @return the lower triangular factor, L, such that L*L'=A.
78 */
79 template <class Real>
getL()80 Array2D<Real> Cholesky<Real>::getL() const
81 {
82 return L_;
83 }
84
85 /**
86 Constructs a lower triangular matrix L, such that L*L'= A.
87 If A is not symmetric positive-definite (SPD), only a
88 partial factorization is performed. If is_spd()
89 evalutate true (1) then the factorizaiton was successful.
90 */
91 template <class Real>
Cholesky(const Array2D<Real> & A)92 Cholesky<Real>::Cholesky(const Array2D<Real> &A)
93 {
94
95
96 int m = A.dim1();
97 int n = A.dim2();
98
99 isspd = (m == n);
100
101 if (m != n)
102 {
103 L_ = Array2D<Real>(0,0);
104 return;
105 }
106
107 L_ = Array2D<Real>(n,n);
108
109
110 // Main loop.
111 for (int j = 0; j < n; j++)
112 {
113 double d = 0.0;
114 for (int k = 0; k < j; k++)
115 {
116 Real s = 0.0;
117 for (int i = 0; i < k; i++)
118 {
119 s += L_[k][i]*L_[j][i];
120 }
121 L_[j][k] = s = (A[j][k] - s)/L_[k][k];
122 d = d + s*s;
123 isspd = isspd && (A[k][j] == A[j][k]);
124 }
125 d = A[j][j] - d;
126 isspd = isspd && (d > 0.0);
127 L_[j][j] = sqrt(d > 0.0 ? d : 0.0);
128 for (int k = j+1; k < n; k++)
129 {
130 L_[j][k] = 0.0;
131 }
132 }
133 }
134
135 /**
136
137 Solve a linear system A*x = b, using the previously computed
138 cholesky factorization of A: L*L'.
139
140 @param B A Matrix with as many rows as A and any number of columns.
141 @return x so that L*L'*x = b. If b is nonconformat, or if A
142 was not symmetric posidtive definite, a null (0x0)
143 array is returned.
144 */
145 template <class Real>
solve(const Array1D<Real> & b)146 Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
147 {
148 int n = L_.dim1();
149 if (b.dim1() != n)
150 return Array1D<Real>();
151
152
153 Array1D<Real> x = b.copy();
154
155
156 // Solve L*y = b;
157 for (int k = 0; k < n; k++)
158 {
159 for (int i = 0; i < k; i++)
160 x[k] -= x[i]*L_[k][i];
161 x[k] /= L_[k][k];
162
163 }
164
165 // Solve L'*X = Y;
166 for (int k = n-1; k >= 0; k--)
167 {
168 for (int i = k+1; i < n; i++)
169 x[k] -= x[i]*L_[i][k];
170 x[k] /= L_[k][k];
171 }
172
173 return x;
174 }
175
176
177 /**
178
179 Solve a linear system A*X = B, using the previously computed
180 cholesky factorization of A: L*L'.
181
182 @param B A Matrix with as many rows as A and any number of columns.
183 @return X so that L*L'*X = B. If B is nonconformat, or if A
184 was not symmetric posidtive definite, a null (0x0)
185 array is returned.
186 */
187 template <class Real>
solve(const Array2D<Real> & B)188 Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
189 {
190 int n = L_.dim1();
191 if (B.dim1() != n)
192 return Array2D<Real>();
193
194
195 Array2D<Real> X = B.copy();
196 int nx = B.dim2();
197
198 // Cleve's original code
199 #if 0
200 // Solve L*Y = B;
201 for (int k = 0; k < n; k++) {
202 for (int i = k+1; i < n; i++) {
203 for (int j = 0; j < nx; j++) {
204 X[i][j] -= X[k][j]*L_[k][i];
205 }
206 }
207 for (int j = 0; j < nx; j++) {
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 for (int j = 0; j < nx; j++) {
215 X[k][j] /= L_[k][k];
216 }
217 for (int i = 0; i < k; i++) {
218 for (int j = 0; j < nx; j++) {
219 X[i][j] -= X[k][j]*L_[k][i];
220 }
221 }
222 }
223 #endif
224
225
226 // Solve L*y = b;
227 for (int j=0; j< nx; j++)
228 {
229 for (int k = 0; k < n; k++)
230 {
231 for (int i = 0; i < k; i++)
232 X[k][j] -= X[i][j]*L_[k][i];
233 X[k][j] /= L_[k][k];
234 }
235 }
236
237 // Solve L'*X = Y;
238 for (int j=0; j<nx; j++)
239 {
240 for (int k = n-1; k >= 0; k--)
241 {
242 for (int i = k+1; i < n; i++)
243 X[k][j] -= X[i][j]*L_[i][k];
244 X[k][j] /= L_[k][k];
245 }
246 }
247
248
249
250 return X;
251 }
252
253
254 }
255 // namespace JAMA
256
257 #endif
258 // JAMA_CHOLESKY_H
259