← Back to branch summary

~ubuntu-branches/ubuntu/maverick/commons-math/maverick

« back to all changes in this revision

Viewing changes to src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java

  • Committer: Bazaar Package Importer
  • Author(s): Damien Raude-Morvan
  • Date: 2009-08-22 01:13:25 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20090822011325-hi4peq1ua5weguwn
Tags: 2.0-1
* New upstream release.
* Set Maintainer field to Debian Java Team
* Add myself as Uploaders
* Switch to Quilt patch system:
  - Refresh all patchs
  - Remove B-D on dpatch, Add B-D on quilt
  - Include patchsys-quilt.mk in debian/rules
* Bump Standards-Version to 3.8.3:
  - Add a README.source to describe patch system
* Maven POMs:
  - Add a Build-Depends-Indep dependency on maven-repo-helper
  - Use mh_installpom and mh_installjar to install the POM and the jar to the
    Maven repository
* Use default-jdk/jre:
  - Depends on java5-runtime-headless
  - Build-Depends on default-jdk
  - Use /usr/lib/jvm/default-java as JAVA_HOME
* Move api documentation to /usr/share/doc/libcommons-math-java/api
* Build-Depends on junit4 instead of junit

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java
 
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
 
 
18
package org.apache.commons.math.linear;
 
19
 
 
20
import org.apache.commons.math.MathRuntimeException;
 
21
 
 
22
 
 
23
/**
 
24
 * Calculates the Cholesky decomposition of a matrix.
 
25
 * <p>The Cholesky decomposition of a real symmetric positive-definite
 
26
 * matrix A consists of a lower triangular matrix L with same size that
 
27
 * satisfy: A = LL<sup>T</sup>Q = I). In a sense, this is the square root of A.</p>
 
28
 *
 
29
 * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a>
 
30
 * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a>
 
31
 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
 
32
 * @since 2.0
 
33
 */
 
34
public class CholeskyDecompositionImpl implements CholeskyDecomposition {
 
35
 
 
36
    /** Default threshold above which off-diagonal elements are considered too different
 
37
     * and matrix not symmetric. */
 
38
    public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
 
39
 
 
40
    /** Default threshold below which diagonal elements are considered null
 
41
     * and matrix not positive definite. */
 
42
    public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
 
43
 
 
44
    /** Row-oriented storage for L<sup>T</sup> matrix data. */
 
45
    private double[][] lTData;
 
46
 
 
47
    /** Cached value of L. */
 
48
    private RealMatrix cachedL;
 
49
 
 
50
    /** Cached value of LT. */
 
51
    private RealMatrix cachedLT;
 
52
 
 
53
    /**
 
54
     * Calculates the Cholesky decomposition of the given matrix.
 
55
     * <p>
 
56
     * Calling this constructor is equivalent to call {@link
 
57
     * #CholeskyDecompositionImpl(RealMatrix, double, double)} with the
 
58
     * thresholds set to the default values {@link
 
59
     * #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} and {@link
 
60
     * #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
 
61
     * </p>
 
62
     * @param matrix the matrix to decompose
 
63
     * @exception NonSquareMatrixException if matrix is not square
 
64
     * @exception NotSymmetricMatrixException if matrix is not symmetric
 
65
     * @exception NotPositiveDefiniteMatrixException if the matrix is not
 
66
     * strictly positive definite
 
67
     * @see #CholeskyDecompositionImpl(RealMatrix, double, double)
 
68
     * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
 
69
     * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
 
70
     */
 
71
    public CholeskyDecompositionImpl(final RealMatrix matrix)
 
72
        throws NonSquareMatrixException,
 
73
               NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
 
74
        this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
 
75
             DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
 
76
    }
 
77
 
 
78
    /**
 
79
     * Calculates the Cholesky decomposition of the given matrix.
 
80
     * @param matrix the matrix to decompose
 
81
     * @param relativeSymmetryThreshold threshold above which off-diagonal
 
82
     * elements are considered too different and matrix not symmetric
 
83
     * @param absolutePositivityThreshold threshold below which diagonal
 
84
     * elements are considered null and matrix not positive definite
 
85
     * @exception NonSquareMatrixException if matrix is not square
 
86
     * @exception NotSymmetricMatrixException if matrix is not symmetric
 
87
     * @exception NotPositiveDefiniteMatrixException if the matrix is not
 
88
     * strictly positive definite
 
89
     * @see #CholeskyDecompositionImpl(RealMatrix)
 
90
     * @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
 
91
     * @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
 
92
     */
 
93
    public CholeskyDecompositionImpl(final RealMatrix matrix,
 
94
                                     final double relativeSymmetryThreshold,
 
95
                                     final double absolutePositivityThreshold)
 
96
        throws NonSquareMatrixException,
 
97
               NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
 
98
 
 
99
        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 &times; X = B.
 
285
         * <p>The A matrix is implicit here. It is </p>
 
286
         * @param b right-hand side of the equation A &times; X = B
 
287
         * @return a vector X such that A &times; 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
}