You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2011/04/25 17:28:12 UTC
svn commit: r1096496 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/linear/
main/java/org/apache/commons/math/random/ site/xdoc/
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