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