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.
+ *The difference with respect to the regular {@link CholeskyDecomposition} + * is that rows/columns may be permuted (hence the rectangular shape instead + * of the traditional triangular shape) and there is a threshold to ignore + * small diagonal elements. This is used for example to generate {@link + * CorrelatedRandomVectorGenerator correlated random n-dimensions vectors} + * in a p-dimension subspace (p < n). In other words, it allows generating + * random vectors from a covariance matrix that is only positive semidefinite, + * and not positive definite.
+ *Rectangular Cholesky decomposition is not suited for solving + * linear systems, so it does not provide any {@link DecompositionSolver + * decomposition solver}.
+ * @see CholeskyDecomposition + * @see CorrelatedRandomVectorGenerator + * @version $Revision$ $Date$ + * @since 3.0 + */ +public interface RectangularCholeskyDecomposition { + + /** Get the root of the covariance matrix. + * The root is the rectangular matrixB such that
+ * the covariance matrix is equal to B.BT
+ * @return root of the square matrix
+ * @see #getRank()
+ */
+ RealMatrix getRootMatrix();
+
+ /** Get the rank of the symmetric positive semidefinite matrix.
+ * The r is the number of independent rows in the symmetric positive semidefinite
+ * matrix, it is also the number of columns of the rectangular
+ * matrix of the decomposition.
+ * @return r of the square matrix.
+ * @see #getRootMatrix()
+ */
+ int getRank();
+
+}
diff --git a/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java b/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java
new file mode 100644
index 000000000..94ab5d924
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java
@@ -0,0 +1,152 @@
+/*
+ * 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 $Revision$ $Date$ + * @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(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(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 264bfd296..96929dbf4 100644 --- a/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java +++ b/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java @@ -18,10 +18,9 @@ package org.apache.commons.math.random; import org.apache.commons.math.exception.DimensionMismatchException; -import org.apache.commons.math.linear.NonPositiveDefiniteMatrixException; -import org.apache.commons.math.linear.MatrixUtils; import org.apache.commons.math.linear.RealMatrix; -import org.apache.commons.math.util.FastMath; +import org.apache.commons.math.linear.RectangularCholeskyDecomposition; +import org.apache.commons.math.linear.RectangularCholeskyDecompositionImpl; /** * A {@link RandomVectorGenerator} that generates vectors with with @@ -68,10 +67,8 @@ public class CorrelatedRandomVectorGenerator private final NormalizedRandomGenerator generator; /** Storage for the normalized vector. */ private final double[] normalized; - /** Permutated Cholesky root of the covariance matrix. */ - private RealMatrix root; - /** Rank of the covariance matrix. */ - private int rank; + /** Root of the covariance matrix. */ + private final RealMatrix root; /** * Builds a correlated random vector generator from its mean @@ -97,10 +94,13 @@ public class CorrelatedRandomVectorGenerator } this.mean = mean.clone(); - decompose(covariance, small); + final RectangularCholeskyDecomposition decomposition = + new RectangularCholeskyDecompositionImpl(covariance, small); + root = decomposition.getRootMatrix(); this.generator = generator; - normalized = new double[rank]; + normalized = new double[decomposition.getRank()]; + } /** @@ -123,10 +123,13 @@ public class CorrelatedRandomVectorGenerator mean[i] = 0; } - decompose(covariance, small); + final RectangularCholeskyDecomposition decomposition = + new RectangularCholeskyDecompositionImpl(covariance, small); + root = decomposition.getRootMatrix(); this.generator = generator; - normalized = new double[rank]; + normalized = new double[decomposition.getRank()]; + } /** Get the underlying normalized components generator. @@ -136,6 +139,16 @@ public class CorrelatedRandomVectorGenerator return generator; } + /** Get the rank of the covariance matrix. + * The rank is the number of independent rows in the covariance + * matrix, it is also the number of columns of the root matrix. + * @return rank of the square matrix. + * @see #getRootMatrix() + */ + public int getRank() { + return normalized.length; + } + /** Get the root of the covariance matrix. * The root is the rectangular matrixB such that
* the covariance matrix is equal to B.BT
@@ -146,122 +159,6 @@ public class CorrelatedRandomVectorGenerator
return root;
}
- /** Get the rank of the covariance matrix.
- * The rank is the number of independent rows in the covariance
- * matrix, it is also the number of columns of the rectangular
- * matrix of the decomposition.
- * @return rank of the square matrix.
- * @see #getRootMatrix()
- */
- public int getRank() {
- return rank;
- }
-
- /** Decompose the original square matrix.
- * The decomposition is based on a Choleski decomposition - * where additional transforms are performed: - *