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 import org.apache.commons.math.exception.util.LocalizedFormats;
022 import org.apache.commons.math.util.FastMath;
023
024
025 /**
026 * Calculates the Cholesky decomposition of a matrix.
027 * <p>The Cholesky decomposition of a real symmetric positive-definite
028 * matrix A consists of a lower triangular matrix L with same size that
029 * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p>
030 *
031 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
032 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
033 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 ao??t 2010) $
034 * @since 2.0
035 */
036 public class CholeskyDecompositionImpl implements CholeskyDecomposition {
037
038 /** Default threshold above which off-diagonal elements are considered too different
039 * and matrix not symmetric. */
040 public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
041
042 /** Default threshold below which diagonal elements are considered null
043 * and matrix not positive definite. */
044 public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
045
046 /** Row-oriented storage for L<sup>T</sup> matrix data. */
047 private double[][] lTData;
048
049 /** Cached value of L. */
050 private RealMatrix cachedL;
051
052 /** Cached value of LT. */
053 private RealMatrix cachedLT;
054
055 /**
056 * Calculates the Cholesky decomposition of the given matrix.
057 * <p>
058 * Calling this constructor is equivalent to call {@link
059 * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the
060 * thresholds set to the default values {@link
061 * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
062 * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
063 * </p>
064 * @param matrix the matrix to decompose
065 * @exception NonSquareMatrixException if matrix is not square
066 * @exception NotSymmetricMatrixException if matrix is not symmetric
067 * @exception NotPositiveDefiniteMatrixException if the matrix is not
068 * strictly positive definite
069 * @see #CholeskyDecompositionImpl(RealMatrix, double, double)
070 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
071 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
072 */
073 public CholeskyDecompositionImpl(final RealMatrix matrix)
074 throws NonSquareMatrixException,
075 NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
076 this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
077 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
078 }
079
080 /**
081 * Calculates the Cholesky decomposition of the given matrix.
082 * @param matrix the matrix to decompose
083 * @param relativeSymmetryThreshold threshold above which off-diagonal
084 * elements are considered too different and matrix not symmetric
085 * @param absolutePositivityThreshold threshold below which diagonal
086 * elements are considered null and matrix not positive definite
087 * @exception NonSquareMatrixException if matrix is not square
088 * @exception NotSymmetricMatrixException if matrix is not symmetric
089 * @exception NotPositiveDefiniteMatrixException if the matrix is not
090 * strictly positive definite
091 * @see #CholeskyDecompositionImpl(RealMatrix)
092 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
093 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
094 */
095 public CholeskyDecompositionImpl(final RealMatrix matrix,
096 final double relativeSymmetryThreshold,
097 final double absolutePositivityThreshold)
098 throws NonSquareMatrixException,
099 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 × X = B.
287 * <p>The A matrix is implicit here. It is </p>
288 * @param b right-hand side of the equation A × X = B
289 * @return a vector X such that A × 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 }