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
This commit is contained in:
Luc Maisonobe
parent
8469bc496c
commit
621bbb8fc7
|
|
@ -52,8 +52,11 @@ If the output is not quite correct, check for invisible trailing spaces!
|
|||
<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>
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -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);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue