Author: luc Date: Mon Apr 25 15:28:12 2011 New Revision: 1096496 URL: http://svn.apache.org/viewvc?rev=1096496&view=rev Log: Added a "rectangular" Cholesky decomposition for positive semidefinite matrices.
JIRA: MATH-541 Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java (with props) commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java (with props) Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java commons/proper/math/trunk/src/site/xdoc/changes.xml Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java?rev=1096496&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java Mon Apr 25 15:28:12 2011 @@ -0,0 +1,66 @@ +/* + * 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.random.CorrelatedRandomVectorGenerator; + + +/** + * An interface to classes that implement an algorithm to calculate a + * rectangular variation of Cholesky decomposition of a real symmetric + * positive semidefinite matrix. + * <p>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 BB<sup>T</sup>, depending + * on a user-defined tolerance. In a sense, this is the square root of A.</p> + * <p>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.</p> + * <p>Rectangular Cholesky decomposition is <em>not</em> suited for solving + * linear systems, so it does not provide any {@link DecompositionSolver + * decomposition solver}.</p> + * @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 matrix <code>B</code> such that + * the covariance matrix is equal to <code>B.B<sup>T</sup></code> + * @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(); + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecomposition.java ------------------------------------------------------------------------------ svn:keywords = Author Date Id Revision Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java?rev=1096496&view=auto ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java (added) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java Mon Apr 25 15:28:12 2011 @@ -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. + * <p>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 BB<sup>T</sup>, depending + * on a user-defined tolerance. In a sense, this is the square root of A.</p> + * + * @see <a href="http://mathworld.wolfram.com/CholeskyDecomposition.html">MathWorld</a> + * @see <a href="http://en.wikipedia.org/wiki/Cholesky_decomposition">Wikipedia</a> + * @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; + } + +} Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java ------------------------------------------------------------------------------ svn:eol-style = native Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/RectangularCholeskyDecompositionImpl.java ------------------------------------------------------------------------------ svn:keywords = Author Date Id Revision Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java?rev=1096496&r1=1096495&r2=1096496&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java Mon Apr 25 15:28:12 2011 @@ -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 CorrelatedRandomVectorGener 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 CorrelatedRandomVectorGener } 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 CorrelatedRandomVectorGener 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,130 +139,24 @@ public class CorrelatedRandomVectorGener return generator; } - /** Get the root of the covariance matrix. - * The root is the rectangular matrix <code>B</code> such that - * the covariance matrix is equal to <code>B.B<sup>T</sup></code> - * @return root of the square matrix - * @see #getRank() - */ - public RealMatrix getRootMatrix() { - 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. + * matrix, it is also the number of columns of the root matrix. * @return rank of the square matrix. * @see #getRootMatrix() */ public int getRank() { - return rank; + return normalized.length; } - /** Decompose the original square matrix. - * <p>The decomposition is based on a Choleski decomposition - * where additional transforms are performed: - * <ul> - * <li>the rows of the decomposed matrix are permuted</li> - * <li>columns with the too small diagonal element are discarded</li> - * <li>the matrix is permuted</li> - * </ul> - * This means that rather than computing M = U<sup>T</sup>.U where U - * is an upper triangular matrix, this method computed M=B.B<sup>T</sup> - * where B is a rectangular matrix. - * @param covariance covariance matrix - * @param small diagonal elements threshold under which column are - * considered to be dependent on previous ones and are discarded - * @throws org.apache.commons.math.linear.NonPositiveDefiniteMatrixException - * if the covariance matrix is not strictly positive definite. + /** Get the root of the covariance matrix. + * The root is the rectangular matrix <code>B</code> such that + * the covariance matrix is equal to <code>B.B<sup>T</sup></code> + * @return root of the square matrix + * @see #getRank() */ - private void decompose(RealMatrix covariance, double small) { - int order = covariance.getRowDimension(); - double[][] c = covariance.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; - } - - rank = 0; - for (boolean loop = true; loop;) { - - // find maximal diagonal element - swap[rank] = rank; - for (int i = rank + 1; i < order; ++i) { - int ii = index[i]; - int isi = index[swap[i]]; - if (c[ii][ii] > c[isi][isi]) { - swap[rank] = i; - } - } - - - // swap elements - if (swap[rank] != rank) { - int tmp = index[rank]; - index[rank] = index[swap[rank]]; - index[swap[rank]] = tmp; - } - - // check diagonal element - int ir = index[rank]; - if (c[ir][ir] < small) { - - if (rank == 0) { - throw new NonPositiveDefiniteMatrixException(ir, small); - } - - // check remaining diagonal elements - for (int i = rank; i < order; ++i) { - if (c[index[i]][index[i]] < -small) { - // there is at least one sufficiently negative diagonal element, - // the covariance 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 covariance matrix - ++rank; - loop = false; - - } else { - - // transform the matrix - double sqrt = FastMath.sqrt(c[ir][ir]); - b[rank][rank] = sqrt; - double inverse = 1 / sqrt; - for (int i = rank + 1; i < order; ++i) { - int ii = index[i]; - double e = inverse * c[ii][ir]; - b[i][rank] = e; - c[ii][ii] -= e * e; - for (int j = rank + 1; j < i; ++j) { - int ij = index[j]; - double f = c[ii][ij] - e * b[j][rank]; - c[ii][ij] = f; - c[ij][ii] = f; - } - } - - // prepare next iteration - loop = ++rank < order; - } - } - - // build the root matrix - root = MatrixUtils.createRealMatrix(order, rank); - for (int i = 0; i < order; ++i) { - for (int j = 0; j < rank; ++j) { - root.setEntry(index[i], j, b[i][j]); - } - } - + public RealMatrix getRootMatrix() { + return root; } /** Generate a correlated random vector. @@ -269,7 +166,7 @@ public class CorrelatedRandomVectorGener public double[] nextVector() { // generate uncorrelated vector - for (int i = 0; i < rank; ++i) { + for (int i = 0; i < normalized.length; ++i) { normalized[i] = generator.nextNormalizedDouble(); } @@ -277,11 +174,13 @@ public class CorrelatedRandomVectorGener double[] correlated = new double[mean.length]; for (int i = 0; i < correlated.length; ++i) { correlated[i] = mean[i]; - for (int j = 0; j < rank; ++j) { + for (int j = 0; j < root.getColumnDimension(); ++j) { correlated[i] += root.getEntry(i, j) * normalized[j]; } } return correlated; + } + } Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1096496&r1=1096495&r2=1096496&view=diff ============================================================================== --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original) +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Mon Apr 25 15:28:12 2011 @@ -52,9 +52,12 @@ The <action> type attribute can be add,u If the output is not quite correct, check for invisible trailing spaces! --> <release version="3.0" date="TBD" description="TBD"> + <action dev="luc" type="add" issue="MATH-541" > + Added a "rectangular" Cholesky decomposition for positive semidefinite matrices. + </action> <action dev="luc" type="add" issue="MATH-563" > Added setters allowing to change the step size control parameters of adaptive - step size ODE integrators + step size ODE integrators. </action> <action dev="luc" type="add" issue="MATH-557" > Added a compareTo method to MathUtils that uses a number of ulps as a