Author: luc
Date: Thu Sep 13 15:13:03 2012
New Revision: 1384363
URL: http://svn.apache.org/viewvc?rev=1384363&view=rev
Log:
Fixed an error in rectangular Cholesky decomposition.
JIRA: MATH-789
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java
(with props)
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1384363&r1=1384362&r2=1384363&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Thu Sep 13 15:13:03 2012
@@ -52,8 +52,11 @@ If the output is not quite correct, chec
<body>
<release version="3.1" date="TBD" description="
">
+ <action dev="lus" type="fix" issue="MATH-789">
+ Fixed an error in rectangular Cholesky decomposition.
+ </action>
<action dev="psteitz" type="update" issue="MATH-859">
- Clarified definition of isSupportXxxBoundInclusive in RealDistribution
+ Clarified definition of isSupportXxxBoundInclusive in
RealDistribution
interface, made code consistent with the definition, and deprecated
these methods, marking for removal in 4.0.
</action>
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java?rev=1384363&r1=1384362&r2=1384363&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java
Thu Sep 13 15:13:03 2012
@@ -62,11 +62,10 @@ public class RectangularCholeskyDecompos
public RectangularCholeskyDecomposition(RealMatrix matrix, double small)
throws NonPositiveDefiniteMatrixException {
- int order = matrix.getRowDimension();
- double[][] c = matrix.getData();
- double[][] b = new double[order][order];
+ final int order = matrix.getRowDimension();
+ final double[][] c = matrix.getData();
+ final 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;
@@ -76,21 +75,24 @@ public class RectangularCholeskyDecompos
for (boolean loop = true; loop;) {
// find maximal diagonal element
- swap[r] = r;
+ int swapR = 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;
+ int isr = index[swapR];
+ if (c[ii][ii] > c[isr][isr]) {
+ swapR = i;
}
}
// swap elements
- if (swap[r] != r) {
- int tmp = index[r];
- index[r] = index[swap[r]];
- index[swap[r]] = tmp;
+ if (swapR != r) {
+ final int tmpIndex = index[r];
+ index[r] = index[swapR];
+ index[swapR] = tmpIndex;
+ final double[] tmpRow = b[r];
+ b[r] = b[swapR];
+ b[swapR] = tmpRow;
}
// check diagonal element
@@ -118,17 +120,18 @@ public class RectangularCholeskyDecompos
} else {
// transform the matrix
- double sqrt = FastMath.sqrt(c[ir][ir]);
+ final double sqrt = FastMath.sqrt(c[ir][ir]);
b[r][r] = sqrt;
- double inverse = 1 / sqrt;
+ final double inverse = 1 / sqrt;
+ final double inverse2 = 1 / c[ir][ir];
for (int i = r + 1; i < order; ++i) {
- int ii = index[i];
- double e = inverse * c[ii][ir];
+ final int ii = index[i];
+ final double e = inverse * c[ii][ir];
b[i][r] = e;
- c[ii][ii] -= e * e;
+ c[ii][ii] -= c[ii][ir] * c[ii][ir] * inverse2;
for (int j = r + 1; j < i; ++j) {
- int ij = index[j];
- double f = c[ii][ij] - e * b[j][r];
+ final int ij = index[j];
+ final double f = c[ii][ij] - e * b[j][r];
c[ii][ij] = f;
c[ij][ii] = f;
}
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java?rev=1384363&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java
Thu Sep 13 15:13:03 2012
@@ -0,0 +1,112 @@
+/*
+ * 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.math3.linear;
+
+import org.junit.Test;
+import org.junit.Assert;
+
+public class RectangularCholeskyDecompositionTest {
+
+ @Test
+ public void testDecomposition3x3() {
+
+ RealMatrix m = MatrixUtils.createRealMatrix(new double[][] {
+ { 1, 9, 9 },
+ { 9, 225, 225 },
+ { 9, 225, 625 }
+ });
+
+ RectangularCholeskyDecomposition d =
+ new RectangularCholeskyDecomposition(m, 1.0e-6);
+
+ // as this decomposition permutes lines and columns, the root is NOT
triangular
+ // (in fact here it is the lower right part of the matrix which is
zero and
+ // the upper left non-zero)
+ Assert.assertEquals(0.8, d.getRootMatrix().getEntry(0, 2), 1.0e-15);
+ Assert.assertEquals(25.0, d.getRootMatrix().getEntry(2, 0), 1.0e-15);
+ Assert.assertEquals(0.0, d.getRootMatrix().getEntry(2, 2), 1.0e-15);
+
+ RealMatrix root = d.getRootMatrix();
+ RealMatrix rebuiltM = root.multiply(root.transpose());
+ Assert.assertEquals(0.0, m.subtract(rebuiltM).getNorm(), 1.0e-15);
+
+ }
+
+ @Test
+ public void testFullRank() {
+
+ RealMatrix base = MatrixUtils.createRealMatrix(new double[][] {
+ { 0.1159548705, 0., 0., 0. },
+ { 0.0896442724, 0.1223540781, 0., 0. },
+ { 0.0852155322, 4.558668e-3, 0.1083577299, 0. },
+ { 0.0905486674, 0.0213768077, 0.0128878333, 0.1014155693 }
+ });
+
+ RealMatrix m = base.multiply(base.transpose());
+
+ RectangularCholeskyDecomposition d =
+ new RectangularCholeskyDecomposition(m, 1.0e-10);
+
+ RealMatrix root = d.getRootMatrix();
+ RealMatrix rebuiltM = root.multiply(root.transpose());
+ Assert.assertEquals(0.0, m.subtract(rebuiltM).getNorm(), 1.0e-15);
+
+ // the pivoted Cholesky decomposition is *not* unique. Here, the root
is
+ // not equal to the original trianbular base matrix
+ Assert.assertTrue(root.subtract(base).getNorm() > 0.3);
+
+ }
+
+ @Test
+ public void testMath789() {
+
+ final RealMatrix m1 = MatrixUtils.createRealMatrix(new double[][]{
+ {0.013445532, 0.010394690, 0.009881156, 0.010499559},
+ {0.010394690, 0.023006616, 0.008196856, 0.010732709},
+ {0.009881156, 0.008196856, 0.019023866, 0.009210099},
+ {0.010499559, 0.010732709, 0.009210099, 0.019107243}
+ });
+ RealMatrix root1 = new RectangularCholeskyDecomposition(m1,
1.0e-10).getRootMatrix();
+ RealMatrix rebuiltM1 = root1.multiply(root1.transpose());
+ Assert.assertEquals(0.0, m1.subtract(rebuiltM1).getNorm(), 1.0e-16);
+
+ final RealMatrix m2 = MatrixUtils.createRealMatrix(new double[][]{
+ {0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.013445532, 0.010394690, 0.009881156, 0.010499559},
+ {0.0, 0.010394690, 0.023006616, 0.008196856, 0.010732709},
+ {0.0, 0.009881156, 0.008196856, 0.019023866, 0.009210099},
+ {0.0, 0.010499559, 0.010732709, 0.009210099, 0.019107243}
+ });
+ RealMatrix root2 = new RectangularCholeskyDecomposition(m2,
1.0e-10).getRootMatrix();
+ RealMatrix rebuiltM2 = root2.multiply(root2.transpose());
+ Assert.assertEquals(0.0, m2.subtract(rebuiltM2).getNorm(), 1.0e-16);
+
+ final RealMatrix m3 = MatrixUtils.createRealMatrix(new double[][]{
+ {0.013445532, 0.010394690, 0.0, 0.009881156, 0.010499559},
+ {0.010394690, 0.023006616, 0.0, 0.008196856, 0.010732709},
+ {0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.009881156, 0.008196856, 0.0, 0.019023866, 0.009210099},
+ {0.010499559, 0.010732709, 0.0, 0.009210099, 0.019107243}
+ });
+ RealMatrix root3 = new RectangularCholeskyDecomposition(m3,
1.0e-10).getRootMatrix();
+ RealMatrix rebuiltM3 = root3.multiply(root3.transpose());
+ Assert.assertEquals(0.0, m3.subtract(rebuiltM3).getNorm(), 1.0e-16);
+
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java?rev=1384363&r1=1384362&r2=1384363&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java
Thu Sep 13 15:13:03 2012
@@ -17,8 +17,13 @@
package org.apache.commons.math3.random;
+import java.util.Arrays;
+
+import org.apache.commons.math3.TestUtils;
+import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
+import org.apache.commons.math3.stat.correlation.StorelessCovariance;
import org.apache.commons.math3.stat.descriptive.moment.VectorialCovariance;
import org.apache.commons.math3.stat.descriptive.moment.VectorialMean;
import org.apache.commons.math3.util.FastMath;
@@ -82,9 +87,19 @@ public class CorrelatedRandomVectorGener
CorrelatedRandomVectorGenerator sg =
new CorrelatedRandomVectorGenerator(mean, covRM, 0.00001, rg);
+ double[] min = new double[mean.length];
+ Arrays.fill(min, Double.POSITIVE_INFINITY);
+ double[] max = new double[mean.length];
+ Arrays.fill(max, Double.NEGATIVE_INFINITY);
for (int i = 0; i < 10; i++) {
double[] generated = sg.nextVector();
- Assert.assertTrue(FastMath.abs(generated[0] - 1) > 0.1);
+ for (int j = 0; j < generated.length; ++j) {
+ min[j] = FastMath.min(min[j], generated[j]);
+ max[j] = FastMath.max(max[j], generated[j]);
+ }
+ }
+ for (int j = 0; j < min.length; ++j) {
+ Assert.assertTrue(max[j] - min[j] > 2.0);
}
}
@@ -123,4 +138,61 @@ public class CorrelatedRandomVectorGener
}
}
+
+ @Test
+ public void testSampleWithZeroCovariance() {
+ final double[][] covMatrix1 = new double[][]{
+ {0.013445532, 0.010394690, 0.009881156, 0.010499559},
+ {0.010394690, 0.023006616, 0.008196856, 0.010732709},
+ {0.009881156, 0.008196856, 0.019023866, 0.009210099},
+ {0.010499559, 0.010732709, 0.009210099, 0.019107243}
+ };
+
+ final double[][] covMatrix2 = new double[][]{
+ {0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.013445532, 0.010394690, 0.009881156, 0.010499559},
+ {0.0, 0.010394690, 0.023006616, 0.008196856, 0.010732709},
+ {0.0, 0.009881156, 0.008196856, 0.019023866, 0.009210099},
+ {0.0, 0.010499559, 0.010732709, 0.009210099, 0.019107243}
+ };
+
+ final double[][] covMatrix3 = new double[][]{
+ {0.013445532, 0.010394690, 0.0, 0.009881156, 0.010499559},
+ {0.010394690, 0.023006616, 0.0, 0.008196856, 0.010732709},
+ {0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.009881156, 0.008196856, 0.0, 0.019023866, 0.009210099},
+ {0.010499559, 0.010732709, 0.0, 0.009210099, 0.019107243}
+ };
+
+ testSampler(covMatrix1, 10000, 0.001);
+ testSampler(covMatrix2, 10000, 0.001);
+ testSampler(covMatrix3, 10000, 0.001);
+
+ }
+
+ private CorrelatedRandomVectorGenerator createSampler(double[][] cov) {
+ RealMatrix matrix = new Array2DRowRealMatrix(cov);
+ double small = 10e-12 * matrix.getNorm();
+ return new CorrelatedRandomVectorGenerator(
+ new double[cov.length],
+ matrix,
+ small,
+ new GaussianRandomGenerator(new JDKRandomGenerator()));
+ }
+
+ private void testSampler(final double[][] covMatrix, int samples, double
epsilon) {
+ CorrelatedRandomVectorGenerator sampler = createSampler(covMatrix);
+
+ StorelessCovariance cov = new StorelessCovariance(covMatrix.length);
+ for (int i = 0; i < samples; ++i) {
+ cov.increment(sampler.nextVector());
+ }
+
+ double[][] sampleCov = cov.getData();
+ for (int r = 0; r < covMatrix.length; ++r) {
+ TestUtils.assertEquals(covMatrix[r], sampleCov[r], epsilon);
+ }
+
+ }
+
}
