001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.linear;
019
020 import org.apache.commons.math.MathRuntimeException;
021
022
023 /**
024 * Calculates the Cholesky decomposition of a matrix.
025 * <p>The Cholesky decomposition of a real symmetric positive-definite
026 * matrix A consists of a lower triangular matrix L with same size that
027 * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p>
028 *
029 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
030 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
031 * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
032 * @since 2.0
033 */
034 public class CholeskyDecompositionImpl implements CholeskyDecomposition {
035
036 /** Default threshold above which off-diagonal elements are considered too different
037 * and matrix not symmetric. */
038 public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
039
040 /** Default threshold below which diagonal elements are considered null
041 * and matrix not positive definite. */
042 public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
043
044 /** Row-oriented storage for L<sup>T</sup> matrix data. */
045 private double[][] lTData;
046
047 /** Cached value of L. */
048 private RealMatrix cachedL;
049
050 /** Cached value of LT. */
051 private RealMatrix cachedLT;
052
053 /**
054 * Calculates the Cholesky decomposition of the given matrix.
055 * <p>
056 * Calling this constructor is equivalent to call {@link
057 * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the
058 * thresholds set to the default values {@link
059 * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
060 * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
061 * </p>
062 * @param matrix the matrix to decompose
063 * @exception NonSquareMatrixException if matrix is not square
064 * @exception NotSymmetricMatrixException if matrix is not symmetric
065 * @exception NotPositiveDefiniteMatrixException if the matrix is not
066 * strictly positive definite
067 * @see #CholeskyDecompositionImpl(RealMatrix, double, double)
068 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
069 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
070 */
071 public CholeskyDecompositionImpl(final RealMatrix matrix)
072 throws NonSquareMatrixException,
073 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
074 this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
075 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
076 }
077
078 /**
079 * Calculates the Cholesky decomposition of the given matrix.
080 * @param matrix the matrix to decompose
081 * @param relativeSymmetryThreshold threshold above which off-diagonal
082 * elements are considered too different and matrix not symmetric
083 * @param absolutePositivityThreshold threshold below which diagonal
084 * elements are considered null and matrix not positive definite
085 * @exception NonSquareMatrixException if matrix is not square
086 * @exception NotSymmetricMatrixException if matrix is not symmetric
087 * @exception NotPositiveDefiniteMatrixException if the matrix is not
088 * strictly positive definite
089 * @see #CholeskyDecompositionImpl(RealMatrix)
090 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
091 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
092 */
093 public CholeskyDecompositionImpl(final RealMatrix matrix,
094 final double relativeSymmetryThreshold,
095 final double absolutePositivityThreshold)
096 throws NonSquareMatrixException,
097 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
098
099 if (!matrix.isSquare()) {
100 throw new NonSquareMatrixException(matrix.getRowDimension(),
101 matrix.getColumnDimension());
102 }
103
104 final int order = matrix.getRowDimension();
105 lTData = matrix.getData();
106 cachedL = null;
107 cachedLT = null;
108
109 // check the matrix before transformation
110 for (int i = 0; i < order; ++i) {
111
112 final double[] lI = lTData[i];
113
114 // check off-diagonal elements (and reset them to 0)
115 for (int j = i + 1; j < order; ++j) {
116 final double[] lJ = lTData[j];
117 final double lIJ = lI[j];
118 final double lJI = lJ[i];
119 final double maxDelta =
120 relativeSymmetryThreshold * Math.max(Math.abs(lIJ), Math.abs(lJI));
121 if (Math.abs(lIJ - lJI) > maxDelta) {
122 throw new NotSymmetricMatrixException();
123 }
124 lJ[i] = 0;
125 }
126 }
127
128 // transform the matrix
129 for (int i = 0; i < order; ++i) {
130
131 final double[] ltI = lTData[i];
132
133 // check diagonal element
134 if (ltI[i] < absolutePositivityThreshold) {
135 throw new NotPositiveDefiniteMatrixException();
136 }
137
138 ltI[i] = Math.sqrt(ltI[i]);
139 final double inverse = 1.0 / ltI[i];
140
141 for (int q = order - 1; q > i; --q) {
142 ltI[q] *= inverse;
143 final double[] ltQ = lTData[q];
144 for (int p = q; p < order; ++p) {
145 ltQ[p] -= ltI[q] * ltI[p];
146 }
147 }
148
149 }
150
151 }
152
153 /** {@inheritDoc} */
154 public RealMatrix getL() {
155 if (cachedL == null) {
156 cachedL = getLT().transpose();
157 }
158 return cachedL;
159 }
160
161 /** {@inheritDoc} */
162 public RealMatrix getLT() {
163
164 if (cachedLT == null) {
165 cachedLT = MatrixUtils.createRealMatrix(lTData);
166 }
167
168 // return the cached matrix
169 return cachedLT;
170
171 }
172
173 /** {@inheritDoc} */
174 public double getDeterminant() {
175 double determinant = 1.0;
176 for (int i = 0; i < lTData.length; ++i) {
177 double lTii = lTData[i][i];
178 determinant *= lTii * lTii;
179 }
180 return determinant;
181 }
182
183 /** {@inheritDoc} */
184 public DecompositionSolver getSolver() {
185 return new Solver(lTData);
186 }
187
188 /** Specialized solver. */
189 private static class Solver implements DecompositionSolver {
190
191 /** Row-oriented storage for L<sup>T</sup> matrix data. */
192 private final double[][] lTData;
193
194 /**
195 * Build a solver from decomposed matrix.
196 * @param lTData row-oriented storage for L<sup>T</sup> matrix data
197 */
198 private Solver(final double[][] lTData) {
199 this.lTData = lTData;
200 }
201
202 /** {@inheritDoc} */
203 public boolean isNonSingular() {
204 // if we get this far, the matrix was positive definite, hence non-singular
205 return true;
206 }
207
208 /** {@inheritDoc} */
209 public double[] solve(double[] b)
210 throws IllegalArgumentException, InvalidMatrixException {
211
212 final int m = lTData.length;
213 if (b.length != m) {
214 throw MathRuntimeException.createIllegalArgumentException(
215 "vector length mismatch: got {0} but expected {1}",
216 b.length, m);
217 }
218
219 final double[] x = b.clone();
220
221 // Solve LY = b
222 for (int j = 0; j < m; j++) {
223 final double[] lJ = lTData[j];
224 x[j] /= lJ[j];
225 final double xJ = x[j];
226 for (int i = j + 1; i < m; i++) {
227 x[i] -= xJ * lJ[i];
228 }
229 }
230
231 // Solve LTX = Y
232 for (int j = m - 1; j >= 0; j--) {
233 x[j] /= lTData[j][j];
234 final double xJ = x[j];
235 for (int i = 0; i < j; i++) {
236 x[i] -= xJ * lTData[i][j];
237 }
238 }
239
240 return x;
241
242 }
243
244 /** {@inheritDoc} */
245 public RealVector solve(RealVector b)
246 throws IllegalArgumentException, InvalidMatrixException {
247 try {
248 return solve((ArrayRealVector) b);
249 } catch (ClassCastException cce) {
250
251 final int m = lTData.length;
252 if (b.getDimension() != m) {
253 throw MathRuntimeException.createIllegalArgumentException(
254 "vector length mismatch: got {0} but expected {1}",
255 b.getDimension(), m);
256 }
257
258 final double[] x = b.getData();
259
260 // Solve LY = b
261 for (int j = 0; j < m; j++) {
262 final double[] lJ = lTData[j];
263 x[j] /= lJ[j];
264 final double xJ = x[j];
265 for (int i = j + 1; i < m; i++) {
266 x[i] -= xJ * lJ[i];
267 }
268 }
269
270 // Solve LTX = Y
271 for (int j = m - 1; j >= 0; j--) {
272 x[j] /= lTData[j][j];
273 final double xJ = x[j];
274 for (int i = 0; i < j; i++) {
275 x[i] -= xJ * lTData[i][j];
276 }
277 }
278
279 return new ArrayRealVector(x, false);
280
281 }
282 }
283
284 /** Solve the linear equation A × X = B.
285 * <p>The A matrix is implicit here. It is </p>
286 * @param b right-hand side of the equation A × X = B
287 * @return a vector X such that A × X = B
288 * @exception IllegalArgumentException if matrices dimensions don't match
289 * @exception InvalidMatrixException if decomposed matrix is singular
290 */
291 public ArrayRealVector solve(ArrayRealVector b)
292 throws IllegalArgumentException, InvalidMatrixException {
293 return new ArrayRealVector(solve(b.getDataRef()), false);
294 }
295
296 /** {@inheritDoc} */
297 public RealMatrix solve(RealMatrix b)
298 throws IllegalArgumentException, InvalidMatrixException {
299
300 final int m = lTData.length;
301 if (b.getRowDimension() != m) {
302 throw MathRuntimeException.createIllegalArgumentException(
303 "dimensions mismatch: got {0}x{1} but expected {2}x{3}",
304 b.getRowDimension(), b.getColumnDimension(), m, "n");
305 }
306
307 final int nColB = b.getColumnDimension();
308 double[][] x = b.getData();
309
310 // Solve LY = b
311 for (int j = 0; j < m; j++) {
312 final double[] lJ = lTData[j];
313 final double lJJ = lJ[j];
314 final double[] xJ = x[j];
315 for (int k = 0; k < nColB; ++k) {
316 xJ[k] /= lJJ;
317 }
318 for (int i = j + 1; i < m; i++) {
319 final double[] xI = x[i];
320 final double lJI = lJ[i];
321 for (int k = 0; k < nColB; ++k) {
322 xI[k] -= xJ[k] * lJI;
323 }
324 }
325 }
326
327 // Solve LTX = Y
328 for (int j = m - 1; j >= 0; j--) {
329 final double lJJ = lTData[j][j];
330 final double[] xJ = x[j];
331 for (int k = 0; k < nColB; ++k) {
332 xJ[k] /= lJJ;
333 }
334 for (int i = 0; i < j; i++) {
335 final double[] xI = x[i];
336 final double lIJ = lTData[i][j];
337 for (int k = 0; k < nColB; ++k) {
338 xI[k] -= xJ[k] * lIJ;
339 }
340 }
341 }
342
343 return new Array2DRowRealMatrix(x, false);
344
345 }
346
347 /** {@inheritDoc} */
348 public RealMatrix getInverse() throws InvalidMatrixException {
349 return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
350 }
351
352 }
353
354 }