commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject 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/
Date Mon, 25 Apr 2011 15:28:12 GMT
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



Mime
View raw message