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.exception.DimensionMismatchException;
021 import org.apache.commons.math.util.FastMath;
022
023
024 /**
025 * Calculates the Cholesky decomposition of a matrix.
026 * <p>The Cholesky decomposition of a real symmetric positive-definite
027 * matrix A consists of a lower triangular matrix L with same size such
028 * that: A = LL<sup>T</sup>. In a sense, this is the square root of A.</p>
029 *
030 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
031 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
032 * @version $Id: CholeskyDecompositionImpl.java 1131229 2011-06-03 20:49:25Z luc $
033 * @since 2.0
034 */
035 public class CholeskyDecompositionImpl implements CholeskyDecomposition {
036 /**
037 * Default threshold above which off-diagonal elements are considered too different
038 * and matrix not symmetric.
039 */
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 */
045 public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
046 /** Row-oriented storage for L<sup>T</sup> matrix data. */
047 private double[][] lTData;
048 /** Cached value of L. */
049 private RealMatrix cachedL;
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 * @throws NonSquareMatrixException if the matrix is not square.
064 * @throws NonSymmetricMatrixException if the matrix is not symmetric.
065 * @throws NonPositiveDefiniteMatrixException 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 this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
073 DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
074 }
075
076 /**
077 * Calculates the Cholesky decomposition of the given matrix.
078 * @param matrix the matrix to decompose
079 * @param relativeSymmetryThreshold threshold above which off-diagonal
080 * elements are considered too different and matrix not symmetric
081 * @param absolutePositivityThreshold threshold below which diagonal
082 * elements are considered null and matrix not positive definite
083 * @throws NonSquareMatrixException if the matrix is not square.
084 * @throws NonSymmetricMatrixException if the matrix is not symmetric.
085 * @throws NonPositiveDefiniteMatrixException if the matrix is not
086 * strictly positive definite.
087 * @see #CholeskyDecompositionImpl(RealMatrix)
088 * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
089 * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
090 */
091 public CholeskyDecompositionImpl(final RealMatrix matrix,
092 final double relativeSymmetryThreshold,
093 final double absolutePositivityThreshold) {
094 if (!matrix.isSquare()) {
095 throw new NonSquareMatrixException(matrix.getRowDimension(),
096 matrix.getColumnDimension());
097 }
098
099 final int order = matrix.getRowDimension();
100 lTData = matrix.getData();
101 cachedL = null;
102 cachedLT = null;
103
104 // check the matrix before transformation
105 for (int i = 0; i < order; ++i) {
106 final double[] lI = lTData[i];
107
108 // check off-diagonal elements (and reset them to 0)
109 for (int j = i + 1; j < order; ++j) {
110 final double[] lJ = lTData[j];
111 final double lIJ = lI[j];
112 final double lJI = lJ[i];
113 final double maxDelta =
114 relativeSymmetryThreshold * FastMath.max(FastMath.abs(lIJ), FastMath.abs(lJI));
115 if (FastMath.abs(lIJ - lJI) > maxDelta) {
116 throw new NonSymmetricMatrixException(i, j, relativeSymmetryThreshold);
117 }
118 lJ[i] = 0;
119 }
120 }
121
122 // transform the matrix
123 for (int i = 0; i < order; ++i) {
124
125 final double[] ltI = lTData[i];
126
127 // check diagonal element
128 if (ltI[i] < absolutePositivityThreshold) {
129 throw new NonPositiveDefiniteMatrixException(i, absolutePositivityThreshold);
130 }
131
132 ltI[i] = FastMath.sqrt(ltI[i]);
133 final double inverse = 1.0 / ltI[i];
134
135 for (int q = order - 1; q > i; --q) {
136 ltI[q] *= inverse;
137 final double[] ltQ = lTData[q];
138 for (int p = q; p < order; ++p) {
139 ltQ[p] -= ltI[q] * ltI[p];
140 }
141 }
142 }
143 }
144
145 /** {@inheritDoc} */
146 public RealMatrix getL() {
147 if (cachedL == null) {
148 cachedL = getLT().transpose();
149 }
150 return cachedL;
151 }
152
153 /** {@inheritDoc} */
154 public RealMatrix getLT() {
155
156 if (cachedLT == null) {
157 cachedLT = MatrixUtils.createRealMatrix(lTData);
158 }
159
160 // return the cached matrix
161 return cachedLT;
162 }
163
164 /** {@inheritDoc} */
165 public double getDeterminant() {
166 double determinant = 1.0;
167 for (int i = 0; i < lTData.length; ++i) {
168 double lTii = lTData[i][i];
169 determinant *= lTii * lTii;
170 }
171 return determinant;
172 }
173
174 /** {@inheritDoc} */
175 public DecompositionSolver getSolver() {
176 return new Solver(lTData);
177 }
178
179 /** Specialized solver. */
180 private static class Solver implements DecompositionSolver {
181 /** Row-oriented storage for L<sup>T</sup> matrix data. */
182 private final double[][] lTData;
183
184 /**
185 * Build a solver from decomposed matrix.
186 * @param lTData row-oriented storage for L<sup>T</sup> matrix data
187 */
188 private Solver(final double[][] lTData) {
189 this.lTData = lTData;
190 }
191
192 /** {@inheritDoc} */
193 public boolean isNonSingular() {
194 // if we get this far, the matrix was positive definite, hence non-singular
195 return true;
196 }
197
198 /** {@inheritDoc} */
199 public double[] solve(double[] b) {
200 final int m = lTData.length;
201 if (b.length != m) {
202 throw new DimensionMismatchException(b.length, m);
203 }
204
205 final double[] x = b.clone();
206
207 // Solve LY = b
208 for (int j = 0; j < m; j++) {
209 final double[] lJ = lTData[j];
210 x[j] /= lJ[j];
211 final double xJ = x[j];
212 for (int i = j + 1; i < m; i++) {
213 x[i] -= xJ * lJ[i];
214 }
215 }
216
217 // Solve LTX = Y
218 for (int j = m - 1; j >= 0; j--) {
219 x[j] /= lTData[j][j];
220 final double xJ = x[j];
221 for (int i = 0; i < j; i++) {
222 x[i] -= xJ * lTData[i][j];
223 }
224 }
225
226 return x;
227 }
228
229 /** {@inheritDoc} */
230 public RealVector solve(RealVector b) {
231 try {
232 return solve((ArrayRealVector) b);
233 } catch (ClassCastException cce) {
234
235 final int m = lTData.length;
236 if (b.getDimension() != m) {
237 throw new DimensionMismatchException(b.getDimension(), m);
238 }
239
240 final double[] x = b.getData();
241
242 // Solve LY = b
243 for (int j = 0; j < m; j++) {
244 final double[] lJ = lTData[j];
245 x[j] /= lJ[j];
246 final double xJ = x[j];
247 for (int i = j + 1; i < m; i++) {
248 x[i] -= xJ * lJ[i];
249 }
250 }
251
252 // Solve LTX = Y
253 for (int j = m - 1; j >= 0; j--) {
254 x[j] /= lTData[j][j];
255 final double xJ = x[j];
256 for (int i = 0; i < j; i++) {
257 x[i] -= xJ * lTData[i][j];
258 }
259 }
260
261 return new ArrayRealVector(x, false);
262 }
263 }
264
265 /** Solve the linear equation A × X = B.
266 * <p>The A matrix is implicit here. It is </p>
267 * @param b right-hand side of the equation A × X = B
268 * @return a vector X such that A × X = B
269 * @throws DimensionMismatchException if the matrices dimensions do not match.
270 * @throws SingularMatrixException if the decomposed matrix is singular.
271 */
272 public ArrayRealVector solve(ArrayRealVector b) {
273 return new ArrayRealVector(solve(b.getDataRef()), false);
274 }
275
276 /** Solve the linear equation A × X = B for matrices A.
277 * <p>The A matrix is implicit, it is provided by the underlying
278 * decomposition algorithm.</p>
279 * @param b right-hand side of the equation A × X = B
280 * @param reuseB if true, the b array will be reused and returned,
281 * instead of being copied
282 * @return a matrix X that minimizes the two norm of A × X - B
283 * @throws org.apache.commons.math.exception.DimensionMismatchException
284 * if the matrices dimensions do not match.
285 * @throws SingularMatrixException
286 * if the decomposed matrix is singular.
287 */
288 private double[][] solve(double[][] b, boolean reuseB) {
289 final int m = lTData.length;
290 if (b.length != m) {
291 throw new DimensionMismatchException(b.length, m);
292 }
293
294 final int nColB = b[0].length;
295 final double[][] x;
296 if (reuseB) {
297 x = b;
298 } else {
299 x = new double[b.length][nColB];
300 for (int i = 0; i < b.length; ++i) {
301 System.arraycopy(b[i], 0, x[i], 0, nColB);
302 }
303 }
304
305 // Solve LY = b
306 for (int j = 0; j < m; j++) {
307 final double[] lJ = lTData[j];
308 final double lJJ = lJ[j];
309 final double[] xJ = x[j];
310 for (int k = 0; k < nColB; ++k) {
311 xJ[k] /= lJJ;
312 }
313 for (int i = j + 1; i < m; i++) {
314 final double[] xI = x[i];
315 final double lJI = lJ[i];
316 for (int k = 0; k < nColB; ++k) {
317 xI[k] -= xJ[k] * lJI;
318 }
319 }
320 }
321
322 // Solve LTX = Y
323 for (int j = m - 1; j >= 0; j--) {
324 final double lJJ = lTData[j][j];
325 final double[] xJ = x[j];
326 for (int k = 0; k < nColB; ++k) {
327 xJ[k] /= lJJ;
328 }
329 for (int i = 0; i < j; i++) {
330 final double[] xI = x[i];
331 final double lIJ = lTData[i][j];
332 for (int k = 0; k < nColB; ++k) {
333 xI[k] -= xJ[k] * lIJ;
334 }
335 }
336 }
337
338 return x;
339
340 }
341
342 /** {@inheritDoc} */
343 public double[][] solve(double[][] b) {
344 return solve(b, false);
345 }
346
347 /** {@inheritDoc} */
348 public RealMatrix solve(RealMatrix b) {
349 return new Array2DRowRealMatrix(solve(b.getData(), true), false);
350 }
351
352 /** {@inheritDoc} */
353 public RealMatrix getInverse() {
354 return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
355 }
356 }
357 }