From b4764661a3f6d7a74acfadd8fbe89b914b90811c Mon Sep 17 00:00:00 2001 From: Sebastien Brisard Date: Sat, 24 Sep 2011 04:54:28 +0000 Subject: [PATCH] Merged RectangularCholeskyDecomposition and RectangularCholeskyDecompositionImpl (see MATH-662). git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1175105 13f79535-47bb-0310-9956-ffa450edef68 --- .../RectangularCholeskyDecomposition.java | 127 +++++++++++++-- .../RectangularCholeskyDecompositionImpl.java | 152 ------------------ .../CorrelatedRandomVectorGenerator.java | 5 +- 3 files changed, 120 insertions(+), 164 deletions(-) delete mode 100644 src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java diff --git a/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java b/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java index bdb5bc0ac..24a91e174 100644 --- a/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java +++ b/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java @@ -17,11 +17,10 @@ package org.apache.commons.math.linear; +import org.apache.commons.math.util.FastMath; /** - * An interface to classes that implement an algorithm to calculate a - * rectangular variation of Cholesky decomposition of a real symmetric - * positive semidefinite matrix. + * Calculates the rectangular Cholesky decomposition of a matrix. *

The rectangular Cholesky decomposition of a real symmetric positive * semidefinite matrix A consists of a rectangular matrix B with the same * number of rows such that: A is almost equal to BBT, depending @@ -38,12 +37,118 @@ package org.apache.commons.math.linear; * linear systems, so it does not provide any {@link DecompositionSolver * decomposition solver}.

* - * @see CholeskyDecomposition - * @see org.apache.commons.math.random.CorrelatedRandomVectorGenerator + * @see MathWorld + * @see Wikipedia * @version $Id$ - * @since 3.0 + * @since 2.0 (changed to concrete class in 3.0) */ -public interface RectangularCholeskyDecomposition { +public class RectangularCholeskyDecomposition { + + /** Permutated Cholesky root of the symmetric positive semidefinite matrix. */ + private final RealMatrix root; + + /** Rank of the symmetric positive semidefinite matrix. */ + private int rank; + + /** + * Decompose a symmetric positive semidefinite matrix. + * + * @param matrix Symmetric positive semidefinite matrix. + * @param small Diagonal elements threshold under which column are + * considered to be dependent on previous ones and are discarded. + * @exception NonPositiveDefiniteMatrixException if the matrix is not + * positive semidefinite. + */ + public RectangularCholeskyDecomposition(RealMatrix matrix, double small) + throws NonPositiveDefiniteMatrixException { + + int order = matrix.getRowDimension(); + double[][] c = matrix.getData(); + double[][] b = new double[order][order]; + + int[] swap = new int[order]; + int[] index = new int[order]; + for (int i = 0; i < order; ++i) { + index[i] = i; + } + + int r = 0; + for (boolean loop = true; loop;) { + + // find maximal diagonal element + swap[r] = r; + for (int i = r + 1; i < order; ++i) { + int ii = index[i]; + int isi = index[swap[i]]; + if (c[ii][ii] > c[isi][isi]) { + swap[r] = i; + } + } + + + // swap elements + if (swap[r] != r) { + int tmp = index[r]; + index[r] = index[swap[r]]; + index[swap[r]] = tmp; + } + + // check diagonal element + int ir = index[r]; + if (c[ir][ir] < small) { + + if (r == 0) { + throw new NonPositiveDefiniteMatrixException(c[ir][ir], ir, small); + } + + // check remaining diagonal elements + for (int i = r; i < order; ++i) { + if (c[index[i]][index[i]] < -small) { + // there is at least one sufficiently negative diagonal element, + // the symmetric positive semidefinite matrix is wrong + throw new NonPositiveDefiniteMatrixException(c[index[i]][index[i]], i, small); + } + } + + // all remaining diagonal elements are close to zero, we consider we have + // found the rank of the symmetric positive semidefinite matrix + ++r; + loop = false; + + } else { + + // transform the matrix + double sqrt = FastMath.sqrt(c[ir][ir]); + b[r][r] = sqrt; + double inverse = 1 / sqrt; + for (int i = r + 1; i < order; ++i) { + int ii = index[i]; + double e = inverse * c[ii][ir]; + b[i][r] = e; + c[ii][ii] -= e * e; + for (int j = r + 1; j < i; ++j) { + int ij = index[j]; + double f = c[ii][ij] - e * b[j][r]; + c[ii][ij] = f; + c[ij][ii] = f; + } + } + + // prepare next iteration + loop = ++r < order; + } + } + + // build the root matrix + rank = r; + root = MatrixUtils.createRealMatrix(order, r); + for (int i = 0; i < order; ++i) { + for (int j = 0; j < r; ++j) { + root.setEntry(index[i], j, b[i][j]); + } + } + + } /** Get the root of the covariance matrix. * The root is the rectangular matrix B such that @@ -51,7 +156,9 @@ public interface RectangularCholeskyDecomposition { * @return root of the square matrix * @see #getRank() */ - RealMatrix getRootMatrix(); + public RealMatrix getRootMatrix() { + return root; + } /** Get the rank of the symmetric positive semidefinite matrix. * The r is the number of independent rows in the symmetric positive semidefinite @@ -60,6 +167,8 @@ public interface RectangularCholeskyDecomposition { * @return r of the square matrix. * @see #getRootMatrix() */ - int getRank(); + public int getRank() { + return rank; + } } diff --git a/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java deleted file mode 100644 index 989d5f4c0..000000000 --- a/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java +++ /dev/null @@ -1,152 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math.linear; - -import org.apache.commons.math.util.FastMath; - -/** - * Calculates the rectangular Cholesky decomposition of a matrix. - *

The rectangular Cholesky decomposition of a real symmetric positive - * semidefinite matrix A consists of a rectangular matrix B with the same - * number of rows such that: A is almost equal to BBT, depending - * on a user-defined tolerance. In a sense, this is the square root of A.

- * - * @see MathWorld - * @see Wikipedia - * @version $Id$ - * @since 2.0 - */ -public class RectangularCholeskyDecompositionImpl implements RectangularCholeskyDecomposition { - - /** Permutated Cholesky root of the symmetric positive semidefinite matrix. */ - private final RealMatrix root; - - /** Rank of the symmetric positive semidefinite matrix. */ - private int rank; - - /** - * Decompose a symmetric positive semidefinite matrix. - * - * @param matrix Symmetric positive semidefinite matrix. - * @param small Diagonal elements threshold under which column are - * considered to be dependent on previous ones and are discarded. - * @exception NonPositiveDefiniteMatrixException if the matrix is not - * positive semidefinite. - */ - public RectangularCholeskyDecompositionImpl(RealMatrix matrix, double small) - throws NonPositiveDefiniteMatrixException { - - int order = matrix.getRowDimension(); - double[][] c = matrix.getData(); - double[][] b = new double[order][order]; - - int[] swap = new int[order]; - int[] index = new int[order]; - for (int i = 0; i < order; ++i) { - index[i] = i; - } - - int r = 0; - for (boolean loop = true; loop;) { - - // find maximal diagonal element - swap[r] = r; - for (int i = r + 1; i < order; ++i) { - int ii = index[i]; - int isi = index[swap[i]]; - if (c[ii][ii] > c[isi][isi]) { - swap[r] = i; - } - } - - - // swap elements - if (swap[r] != r) { - int tmp = index[r]; - index[r] = index[swap[r]]; - index[swap[r]] = tmp; - } - - // check diagonal element - int ir = index[r]; - if (c[ir][ir] < small) { - - if (r == 0) { - throw new NonPositiveDefiniteMatrixException(c[ir][ir], ir, small); - } - - // check remaining diagonal elements - for (int i = r; i < order; ++i) { - if (c[index[i]][index[i]] < -small) { - // there is at least one sufficiently negative diagonal element, - // the symmetric positive semidefinite matrix is wrong - throw new NonPositiveDefiniteMatrixException(c[index[i]][index[i]], i, small); - } - } - - // all remaining diagonal elements are close to zero, we consider we have - // found the rank of the symmetric positive semidefinite matrix - ++r; - loop = false; - - } else { - - // transform the matrix - double sqrt = FastMath.sqrt(c[ir][ir]); - b[r][r] = sqrt; - double inverse = 1 / sqrt; - for (int i = r + 1; i < order; ++i) { - int ii = index[i]; - double e = inverse * c[ii][ir]; - b[i][r] = e; - c[ii][ii] -= e * e; - for (int j = r + 1; j < i; ++j) { - int ij = index[j]; - double f = c[ii][ij] - e * b[j][r]; - c[ii][ij] = f; - c[ij][ii] = f; - } - } - - // prepare next iteration - loop = ++r < order; - } - } - - // build the root matrix - rank = r; - root = MatrixUtils.createRealMatrix(order, r); - for (int i = 0; i < order; ++i) { - for (int j = 0; j < r; ++j) { - root.setEntry(index[i], j, b[i][j]); - } - } - - } - - /** {@inheritDoc} */ - public RealMatrix getRootMatrix() { - return root; - } - - /** {@inheritDoc} */ - public int getRank() { - return rank; - } - -} diff --git a/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java index ab9dca6b7..b3d82d9a0 100644 --- a/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java +++ b/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java @@ -20,7 +20,6 @@ package org.apache.commons.math.random; import org.apache.commons.math.exception.DimensionMismatchException; import org.apache.commons.math.linear.RealMatrix; import org.apache.commons.math.linear.RectangularCholeskyDecomposition; -import org.apache.commons.math.linear.RectangularCholeskyDecompositionImpl; /** * A {@link RandomVectorGenerator} that generates vectors with with @@ -95,7 +94,7 @@ public class CorrelatedRandomVectorGenerator this.mean = mean.clone(); final RectangularCholeskyDecomposition decomposition = - new RectangularCholeskyDecompositionImpl(covariance, small); + new RectangularCholeskyDecomposition(covariance, small); root = decomposition.getRootMatrix(); this.generator = generator; @@ -124,7 +123,7 @@ public class CorrelatedRandomVectorGenerator } final RectangularCholeskyDecomposition decomposition = - new RectangularCholeskyDecompositionImpl(covariance, small); + new RectangularCholeskyDecomposition(covariance, small); root = decomposition.getRootMatrix(); this.generator = generator;