source: src/main/java/agents/org/apache/commons/math/linear/CholeskyDecompositionImpl.java

Last change on this file was 1, checked in by Wouter Pasman, 7 years ago

Initial import : Genius 9.0.0
File size: 12.5 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package agents.org.apache.commons.math.linear;
19
20import agents.org.apache.commons.math.MathRuntimeException;
21import agents.org.apache.commons.math.exception.util.LocalizedFormats;
22import agents.org.apache.commons.math.util.FastMath;
23
24
25/**
26 * Calculates the Cholesky decomposition of a matrix.
27 * <p>The Cholesky decomposition of a real symmetric positive-definite
28 * matrix A consists of a lower triangular matrix L with same size that
29 * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p>
30 *
31 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
32 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
33 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
34 * @since 2.0
35 */
36public class CholeskyDecompositionImpl implements CholeskyDecomposition {
37
38 /** Default threshold above which off-diagonal elements are considered too different
39 * and matrix not symmetric. */
40 public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
41
42 /** Default threshold below which diagonal elements are considered null
43 * and matrix not positive definite. */
44 public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
45
46 /** Row-oriented storage for L<sup>T</sup> matrix data. */
47 private double[][] lTData;
48
49 /** Cached value of L. */
50 private RealMatrix cachedL;
51
52 /** Cached value of LT. */
53 private RealMatrix cachedLT;
54
55 /**
56 * Calculates the Cholesky decomposition of the given matrix.
57 * <p>
58 * Calling this constructor is equivalent to call {@link
59 * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the
60 * thresholds set to the default values {@link
61 * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
62 * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
63 * </p>
64 * @param matrix the matrix to decompose
65 * @exception NonSquareMatrixException if matrix is not square
66 * @exception NotSymmetricMatrixException if matrix is not symmetric
67 * @exception NotPositiveDefiniteMatrixException if the matrix is not
68 * strictly positive definite
69 * @see #CholeskyDecompositionImpl(RealMatrix, double, double)
70 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
71 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
72 */
73 public CholeskyDecompositionImpl(final RealMatrix matrix)
74 throws NonSquareMatrixException,
75 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
76 this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
77 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
78 }
79
80 /**
81 * Calculates the Cholesky decomposition of the given matrix.
82 * @param matrix the matrix to decompose
83 * @param relativeSymmetryThreshold threshold above which off-diagonal
84 * elements are considered too different and matrix not symmetric
85 * @param absolutePositivityThreshold threshold below which diagonal
86 * elements are considered null and matrix not positive definite
87 * @exception NonSquareMatrixException if matrix is not square
88 * @exception NotSymmetricMatrixException if matrix is not symmetric
89 * @exception NotPositiveDefiniteMatrixException if the matrix is not
90 * strictly positive definite
91 * @see #CholeskyDecompositionImpl(RealMatrix)
92 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
93 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
94 */
95 public CholeskyDecompositionImpl(final RealMatrix matrix,
96 final double relativeSymmetryThreshold,
97 final double absolutePositivityThreshold)
98 throws NonSquareMatrixException,
99 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
100
101 if (!matrix.isSquare()) {
102 throw new NonSquareMatrixException(matrix.getRowDimension(),
103 matrix.getColumnDimension());
104 }
105
106 final int order = matrix.getRowDimension();
107 lTData = matrix.getData();
108 cachedL = null;
109 cachedLT = null;
110
111 // check the matrix before transformation
112 for (int i = 0; i < order; ++i) {
113
114 final double[] lI = lTData[i];
115
116 // check off-diagonal elements (and reset them to 0)
117 for (int j = i + 1; j < order; ++j) {
118 final double[] lJ = lTData[j];
119 final double lIJ = lI[j];
120 final double lJI = lJ[i];
121 final double maxDelta =
122 relativeSymmetryThreshold * FastMath.max(FastMath.abs(lIJ), FastMath.abs(lJI));
123 if (FastMath.abs(lIJ - lJI) > maxDelta) {
124 throw new NotSymmetricMatrixException();
125 }
126 lJ[i] = 0;
127 }
128 }
129
130 // transform the matrix
131 for (int i = 0; i < order; ++i) {
132
133 final double[] ltI = lTData[i];
134
135 // check diagonal element
136 if (ltI[i] < absolutePositivityThreshold) {
137 throw new NotPositiveDefiniteMatrixException();
138 }
139
140 ltI[i] = FastMath.sqrt(ltI[i]);
141 final double inverse = 1.0 / ltI[i];
142
143 for (int q = order - 1; q > i; --q) {
144 ltI[q] *= inverse;
145 final double[] ltQ = lTData[q];
146 for (int p = q; p < order; ++p) {
147 ltQ[p] -= ltI[q] * ltI[p];
148 }
149 }
150
151 }
152
153 }
154
155 /** {@inheritDoc} */
156 public RealMatrix getL() {
157 if (cachedL == null) {
158 cachedL = getLT().transpose();
159 }
160 return cachedL;
161 }
162
163 /** {@inheritDoc} */
164 public RealMatrix getLT() {
165
166 if (cachedLT == null) {
167 cachedLT = MatrixUtils.createRealMatrix(lTData);
168 }
169
170 // return the cached matrix
171 return cachedLT;
172
173 }
174
175 /** {@inheritDoc} */
176 public double getDeterminant() {
177 double determinant = 1.0;
178 for (int i = 0; i < lTData.length; ++i) {
179 double lTii = lTData[i][i];
180 determinant *= lTii * lTii;
181 }
182 return determinant;
183 }
184
185 /** {@inheritDoc} */
186 public DecompositionSolver getSolver() {
187 return new Solver(lTData);
188 }
189
190 /** Specialized solver. */
191 private static class Solver implements DecompositionSolver {
192
193 /** Row-oriented storage for L<sup>T</sup> matrix data. */
194 private final double[][] lTData;
195
196 /**
197 * Build a solver from decomposed matrix.
198 * @param lTData row-oriented storage for L<sup>T</sup> matrix data
199 */
200 private Solver(final double[][] lTData) {
201 this.lTData = lTData;
202 }
203
204 /** {@inheritDoc} */
205 public boolean isNonSingular() {
206 // if we get this far, the matrix was positive definite, hence non-singular
207 return true;
208 }
209
210 /** {@inheritDoc} */
211 public double[] solve(double[] b)
212 throws IllegalArgumentException, InvalidMatrixException {
213
214 final int m = lTData.length;
215 if (b.length != m) {
216 throw MathRuntimeException.createIllegalArgumentException(
217 LocalizedFormats.VECTOR_LENGTH_MISMATCH,
218 b.length, m);
219 }
220
221 final double[] x = b.clone();
222
223 // Solve LY = b
224 for (int j = 0; j < m; j++) {
225 final double[] lJ = lTData[j];
226 x[j] /= lJ[j];
227 final double xJ = x[j];
228 for (int i = j + 1; i < m; i++) {
229 x[i] -= xJ * lJ[i];
230 }
231 }
232
233 // Solve LTX = Y
234 for (int j = m - 1; j >= 0; j--) {
235 x[j] /= lTData[j][j];
236 final double xJ = x[j];
237 for (int i = 0; i < j; i++) {
238 x[i] -= xJ * lTData[i][j];
239 }
240 }
241
242 return x;
243
244 }
245
246 /** {@inheritDoc} */
247 public RealVector solve(RealVector b)
248 throws IllegalArgumentException, InvalidMatrixException {
249 try {
250 return solve((ArrayRealVector) b);
251 } catch (ClassCastException cce) {
252
253 final int m = lTData.length;
254 if (b.getDimension() != m) {
255 throw MathRuntimeException.createIllegalArgumentException(
256 LocalizedFormats.VECTOR_LENGTH_MISMATCH,
257 b.getDimension(), m);
258 }
259
260 final double[] x = b.getData();
261
262 // Solve LY = b
263 for (int j = 0; j < m; j++) {
264 final double[] lJ = lTData[j];
265 x[j] /= lJ[j];
266 final double xJ = x[j];
267 for (int i = j + 1; i < m; i++) {
268 x[i] -= xJ * lJ[i];
269 }
270 }
271
272 // Solve LTX = Y
273 for (int j = m - 1; j >= 0; j--) {
274 x[j] /= lTData[j][j];
275 final double xJ = x[j];
276 for (int i = 0; i < j; i++) {
277 x[i] -= xJ * lTData[i][j];
278 }
279 }
280
281 return new ArrayRealVector(x, false);
282
283 }
284 }
285
286 /** Solve the linear equation A &times; X = B.
287 * <p>The A matrix is implicit here. It is </p>
288 * @param b right-hand side of the equation A &times; X = B
289 * @return a vector X such that A &times; X = B
290 * @exception IllegalArgumentException if matrices dimensions don't match
291 * @exception InvalidMatrixException if decomposed matrix is singular
292 */
293 public ArrayRealVector solve(ArrayRealVector b)
294 throws IllegalArgumentException, InvalidMatrixException {
295 return new ArrayRealVector(solve(b.getDataRef()), false);
296 }
297
298 /** {@inheritDoc} */
299 public RealMatrix solve(RealMatrix b)
300 throws IllegalArgumentException, InvalidMatrixException {
301
302 final int m = lTData.length;
303 if (b.getRowDimension() != m) {
304 throw MathRuntimeException.createIllegalArgumentException(
305 LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
306 b.getRowDimension(), b.getColumnDimension(), m, "n");
307 }
308
309 final int nColB = b.getColumnDimension();
310 double[][] x = b.getData();
311
312 // Solve LY = b
313 for (int j = 0; j < m; j++) {
314 final double[] lJ = lTData[j];
315 final double lJJ = lJ[j];
316 final double[] xJ = x[j];
317 for (int k = 0; k < nColB; ++k) {
318 xJ[k] /= lJJ;
319 }
320 for (int i = j + 1; i < m; i++) {
321 final double[] xI = x[i];
322 final double lJI = lJ[i];
323 for (int k = 0; k < nColB; ++k) {
324 xI[k] -= xJ[k] * lJI;
325 }
326 }
327 }
328
329 // Solve LTX = Y
330 for (int j = m - 1; j >= 0; j--) {
331 final double lJJ = lTData[j][j];
332 final double[] xJ = x[j];
333 for (int k = 0; k < nColB; ++k) {
334 xJ[k] /= lJJ;
335 }
336 for (int i = 0; i < j; i++) {
337 final double[] xI = x[i];
338 final double lIJ = lTData[i][j];
339 for (int k = 0; k < nColB; ++k) {
340 xI[k] -= xJ[k] * lIJ;
341 }
342 }
343 }
344
345 return new Array2DRowRealMatrix(x, false);
346
347 }
348
349 /** {@inheritDoc} */
350 public RealMatrix getInverse() throws InvalidMatrixException {
351 return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
352 }
353
354 }
355
356}
Note: See TracBrowser for help on using the repository browser.