From 621bbb8fc7ef562d86e18bbc28c5d4f25a68b442 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Thu, 13 Sep 2012 15:13:03 +0000 Subject: [PATCH] Fixed an error in rectangular Cholesky decomposition. JIRA: MATH-789 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1384363 13f79535-47bb-0310-9956-ffa450edef68 --- src/changes/changes.xml | 5 +- .../RectangularCholeskyDecomposition.java | 41 ++++--- .../RectangularCholeskyDecompositionTest.java | 112 ++++++++++++++++++ .../CorrelatedRandomVectorGeneratorTest.java | 74 +++++++++++- 4 files changed, 211 insertions(+), 21 deletions(-) create mode 100644 src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java diff --git a/src/changes/changes.xml b/src/changes/changes.xml index ec44dd50f..046450552 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -52,8 +52,11 @@ If the output is not quite correct, check for invisible trailing spaces! + + Fixed an error in rectangular Cholesky decomposition. + - 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. diff --git a/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java b/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java index 38584d420..aba7b9806 100644 --- a/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java +++ b/src/main/java/org/apache/commons/math3/linear/RectangularCholeskyDecomposition.java @@ -62,11 +62,10 @@ public class RectangularCholeskyDecomposition { 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 RectangularCholeskyDecomposition { 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 RectangularCholeskyDecomposition { } 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; } diff --git a/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java b/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java new file mode 100644 index 000000000..6fbfdfd3c --- /dev/null +++ b/src/test/java/org/apache/commons/math3/linear/RectangularCholeskyDecompositionTest.java @@ -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); + + } + +} diff --git a/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java b/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java index 167382a8f..eba52a387 100644 --- a/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java +++ b/src/test/java/org/apache/commons/math3/random/CorrelatedRandomVectorGeneratorTest.java @@ -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 CorrelatedRandomVectorGeneratorTest { 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 CorrelatedRandomVectorGeneratorTest { } } + + @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); + } + + } + }