[MATH-789] Fixed rank calculation in case of dependant columns, added additional constructor that repaces small parameter.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1384945 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Thomas Neidhart 2012-09-14 21:55:46 +00:00
parent c8b8e61243
commit 51f446bf0c
3 changed files with 33 additions and 13 deletions

View File

@ -52,9 +52,28 @@ public class RectangularCholeskyDecomposition {
/**
* Decompose a symmetric positive semidefinite matrix.
* <p>
* <b>Note:</b> this constructor follows the linpack method to detect dependent
* columns by proceeding with the Cholesky algorithm until a nonpositive diagonal
* element is encountered.
*
* @see <a href="http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf">
* Analysis of the Cholesky Decomposition of a Semi-definite Matrix</a>
*
* @param matrix Symmetric positive semidefinite matrix.
* @param small Diagonal elements threshold under which column are
* @exception NonPositiveDefiniteMatrixException if the matrix is not
* positive semidefinite.
*/
public RectangularCholeskyDecomposition(RealMatrix matrix)
throws NonPositiveDefiniteMatrixException {
this(matrix, 0);
}
/**
* Decompose a symmetric positive semidefinite matrix.
*
* @param matrix Symmetric positive semidefinite matrix.
* @param small Diagonal elements threshold under which columns are
* considered to be dependent on previous ones and are discarded.
* @exception NonPositiveDefiniteMatrixException if the matrix is not
* positive semidefinite.
@ -97,7 +116,7 @@ public class RectangularCholeskyDecomposition {
// check diagonal element
int ir = index[r];
if (c[ir][ir] < small) {
if (c[ir][ir] <= small) {
if (r == 0) {
throw new NonPositiveDefiniteMatrixException(c[ir][ir], ir, small);
@ -114,7 +133,6 @@ public class RectangularCholeskyDecomposition {
// 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 {

View File

@ -81,9 +81,7 @@ public class RectangularCholeskyDecompositionTest {
{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);
composeAndTest(m1, 4);
final RealMatrix m2 = MatrixUtils.createRealMatrix(new double[][]{
{0.0, 0.0, 0.0, 0.0, 0.0},
@ -92,9 +90,7 @@ public class RectangularCholeskyDecompositionTest {
{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);
composeAndTest(m2, 4);
final RealMatrix m3 = MatrixUtils.createRealMatrix(new double[][]{
{0.013445532, 0.010394690, 0.0, 0.009881156, 0.010499559},
@ -103,10 +99,16 @@ public class RectangularCholeskyDecompositionTest {
{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);
composeAndTest(m3, 4);
}
private void composeAndTest(RealMatrix m, int expectedRank) {
RectangularCholeskyDecomposition r = new RectangularCholeskyDecomposition(m);
Assert.assertEquals(expectedRank, r.getRank());
RealMatrix root = r.getRootMatrix();
RealMatrix rebuiltMatrix = root.multiply(root.transpose());
Assert.assertEquals(0.0, m.subtract(rebuiltMatrix).getNorm(), 1.0e-16);
}
}

View File

@ -68,7 +68,7 @@ public class CorrelatedRandomVectorGeneratorTest {
@Test
public void testRank() {
Assert.assertEquals(3, generator.getRank());
Assert.assertEquals(2, generator.getRank());
}
@Test