Singular Value Decomposition now computes either the compact SVD (using only positive singular values) or truncated SVD (using a user-specified maximal number of singular values).
Fixed Singular Value Decomposition solving of singular systems. JIRA: MATH-320, MATH-321 git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@894908 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
81ae49bafe
commit
b2f3f6db41
|
|
@ -159,27 +159,28 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio
|
|||
if (m >= n) {
|
||||
// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
|
||||
final RealMatrix e =
|
||||
eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1);
|
||||
eigenDecomposition.getV().getSubMatrix(0, n - 1, 0, p - 1);
|
||||
final double[][] eData = e.getData();
|
||||
final double[][] wData = new double[m][p];
|
||||
double[] ei1 = eData[0];
|
||||
for (int i = 0; i < p - 1; ++i) {
|
||||
for (int i = 0; i < p; ++i) {
|
||||
// compute W = B.E.S^(-1) where E is the eigenvectors matrix
|
||||
final double mi = mainBidiagonal[i];
|
||||
final double si = secondaryBidiagonal[i];
|
||||
final double[] ei0 = ei1;
|
||||
final double[] wi = wData[i];
|
||||
ei1 = eData[i + 1];
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
|
||||
if (i < n - 1) {
|
||||
ei1 = eData[i + 1];
|
||||
final double si = secondaryBidiagonal[i];
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
|
||||
}
|
||||
} else {
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wi[j] = mi * ei0[j] / singularValues[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
// last row
|
||||
final double lastMain = mainBidiagonal[p - 1];
|
||||
final double[] wr1 = wData[p - 1];
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wr1[j] = ei1[j] * lastMain / singularValues[j];
|
||||
}
|
||||
|
||||
for (int i = p; i < m; ++i) {
|
||||
wData[i] = new double[p];
|
||||
}
|
||||
|
|
@ -247,26 +248,26 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio
|
|||
// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
|
||||
// compute W = Bt.E.S^(-1) where E is the eigenvectors matrix
|
||||
final RealMatrix e =
|
||||
eigenDecomposition.getV().getSubMatrix(0, p - 1, 0, p - 1);
|
||||
eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1);
|
||||
final double[][] eData = e.getData();
|
||||
final double[][] wData = new double[n][p];
|
||||
double[] ei1 = eData[0];
|
||||
for (int i = 0; i < p - 1; ++i) {
|
||||
for (int i = 0; i < p; ++i) {
|
||||
final double mi = mainBidiagonal[i];
|
||||
final double si = secondaryBidiagonal[i];
|
||||
final double[] ei0 = ei1;
|
||||
final double[] wi = wData[i];
|
||||
ei1 = eData[i + 1];
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
|
||||
if (i < m - 1) {
|
||||
ei1 = eData[i + 1];
|
||||
final double si = secondaryBidiagonal[i];
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j];
|
||||
}
|
||||
} else {
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wi[j] = mi * ei0[j] / singularValues[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
// last row
|
||||
final double lastMain = mainBidiagonal[p - 1];
|
||||
final double[] wr1 = wData[p - 1];
|
||||
for (int j = 0; j < p; ++j) {
|
||||
wr1[j] = ei1[j] * lastMain / singularValues[j];
|
||||
}
|
||||
for (int i = p; i < n; ++i) {
|
||||
wData[i] = new double[p];
|
||||
}
|
||||
|
|
|
|||
|
|
@ -39,10 +39,18 @@ The <action> type attribute can be add,update,fix,remove.
|
|||
</properties>
|
||||
<body>
|
||||
<release version="2.1" date="TBD" description="TBD">
|
||||
<action dev="luc" type="add" issue="MATH-321" >
|
||||
Singular Value Decomposition now computes either the compact SVD (using only
|
||||
positive singular values) or truncated SVD (using a user-specified maximal
|
||||
number of singular values).
|
||||
</action>
|
||||
<action dev="luc" type="fix" issue="MATH-320" >
|
||||
Fixed Singular Value Decomposition solving of singular systems.
|
||||
</action>
|
||||
<action dev="psteitz" type="update" issue="MATH-239" due-to="Christian Semrau">
|
||||
Added MathUtils methods to compute gcd and lcm for long arguments.
|
||||
</actions>
|
||||
<action dev="psteitz" tyoe="update" issue="MATH-287" due-to="Matthew Rowles">
|
||||
</action>
|
||||
<action dev="psteitz" type="update" issue="MATH-287" due-to="Matthew Rowles">
|
||||
Added support for weighted univariate statistics.
|
||||
</action>
|
||||
<action dev="luc" type="fix" issue="MATH-326" due-to="Jake Mannix">
|
||||
|
|
|
|||
|
|
@ -138,6 +138,32 @@ public class SingularValueSolverTest {
|
|||
Assert.assertEquals(3.0, svd.getConditionNumber(), 1.0e-15);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testTruncated() {
|
||||
|
||||
RealMatrix rm = new Array2DRowRealMatrix(new double[][] {
|
||||
{ 1.0, 2.0, 3.0 }, { 2.0, 3.0, 4.0 }, { 3.0, 5.0, 7.0 }
|
||||
});
|
||||
double s439 = Math.sqrt(439.0);
|
||||
double[] reference = new double[] {
|
||||
Math.sqrt(3.0 * (21.0 + s439))
|
||||
};
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(rm, 1);
|
||||
|
||||
// check we get the expected theoretical singular values
|
||||
double[] singularValues = svd.getSingularValues();
|
||||
Assert.assertEquals(reference.length, singularValues.length);
|
||||
for (int i = 0; i < reference.length; ++i) {
|
||||
Assert.assertEquals(reference[i], singularValues[i], 4.0e-13);
|
||||
}
|
||||
|
||||
// check the truncated decomposition DON'T allows to recover the original matrix
|
||||
RealMatrix recomposed = svd.getU().multiply(svd.getS()).multiply(svd.getVT());
|
||||
Assert.assertTrue(recomposed.subtract(rm).getNorm() > 1.4);
|
||||
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testMath320A() {
|
||||
RealMatrix rm = new Array2DRowRealMatrix(new double[][] {
|
||||
|
|
@ -149,44 +175,38 @@ public class SingularValueSolverTest {
|
|||
};
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(rm);
|
||||
|
||||
// check we get the expected theoretical singular values
|
||||
double[] singularValues = svd.getSingularValues();
|
||||
Assert.assertEquals(reference.length, singularValues.length);
|
||||
for (int i = 0; i < reference.length; ++i) {
|
||||
Assert.assertEquals(reference[i], singularValues[i], 4.0e-13);
|
||||
}
|
||||
regularElements(svd.getU());
|
||||
regularElements(svd.getVT());
|
||||
// double[] b = new double[] { 5.0, 6.0, 7.0 };
|
||||
// double[] resSVD = svd.getSolver().solve(b);
|
||||
// Assert.assertEquals(rm.getColumnDimension(), resSVD.length);
|
||||
// System.out.println("resSVD = " + resSVD[0] + " " + resSVD[1] + " " + resSVD[2]);
|
||||
// double minResidual = Double.POSITIVE_INFINITY;
|
||||
// double d0Min = Double.NaN;
|
||||
// double d1Min = Double.NaN;
|
||||
// double d2Min = Double.NaN;
|
||||
// double h = 0.01;
|
||||
// int k = 100;
|
||||
// for (double d0 = -k * h; d0 <= k * h; d0 += h) {
|
||||
// for (double d1 = -k * h ; d1 <= k * h; d1 += h) {
|
||||
// for (double d2 = -k * h; d2 <= k * h; d2 += h) {
|
||||
// double[] f = rm.operate(new double[] { resSVD[0] + d0, resSVD[1] + d1, resSVD[2] + d2 });
|
||||
// double residual = Math.sqrt((f[0] - b[0]) * (f[0] - b[0]) +
|
||||
// (f[1] - b[1]) * (f[1] - b[1]) +
|
||||
// (f[2] - b[2]) * (f[2] - b[2]));
|
||||
// if (residual < minResidual) {
|
||||
// d0Min = d0;
|
||||
// d1Min = d1;
|
||||
// d2Min = d2;
|
||||
// minResidual = residual;
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
// System.out.println(d0Min + " " + d1Min + " " + d2Min + " -> " + minResidual);
|
||||
// Assert.assertEquals(0, d0Min, 1.0e-15);
|
||||
// Assert.assertEquals(0, d1Min, 1.0e-15);
|
||||
// Assert.assertEquals(0, d2Min, 1.0e-15);
|
||||
}
|
||||
|
||||
// check the decomposition allows to recover the original matrix
|
||||
RealMatrix recomposed = svd.getU().multiply(svd.getS()).multiply(svd.getVT());
|
||||
Assert.assertEquals(0.0, recomposed.subtract(rm).getNorm(), 5.0e-13);
|
||||
|
||||
// check we can solve a singular system
|
||||
double[] b = new double[] { 5.0, 6.0, 7.0 };
|
||||
double[] resSVD = svd.getSolver().solve(b);
|
||||
Assert.assertEquals(rm.getColumnDimension(), resSVD.length);
|
||||
|
||||
// check the solution really minimizes the residuals
|
||||
double svdMinResidual = residual(rm, resSVD, b);
|
||||
double epsilon = 2 * Math.ulp(svdMinResidual);
|
||||
double h = 0.1;
|
||||
int k = 3;
|
||||
for (double d0 = -k * h; d0 <= k * h; d0 += h) {
|
||||
for (double d1 = -k * h ; d1 <= k * h; d1 += h) {
|
||||
for (double d2 = -k * h; d2 <= k * h; d2 += h) {
|
||||
double[] x = new double[] { resSVD[0] + d0, resSVD[1] + d1, resSVD[2] + d2 };
|
||||
Assert.assertTrue((residual(rm, x, b) - svdMinResidual) > -epsilon);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testMath320B() {
|
||||
|
|
@ -195,18 +215,17 @@ public class SingularValueSolverTest {
|
|||
});
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(rm);
|
||||
regularElements(svd.getU());
|
||||
regularElements(svd.getVT());
|
||||
RealMatrix recomposed = svd.getU().multiply(svd.getS()).multiply(svd.getVT());
|
||||
Assert.assertEquals(0.0, recomposed.subtract(rm).getNorm(), 2.0e-15);
|
||||
}
|
||||
|
||||
private void regularElements(RealMatrix m) {
|
||||
for (int i = 0; i < m.getRowDimension(); ++i) {
|
||||
for (int j = 0; j < m.getColumnDimension(); ++j) {
|
||||
double mij = m.getEntry(i, j);
|
||||
Assert.assertFalse(Double.isInfinite(mij));
|
||||
Assert.assertFalse(Double.isNaN(mij));
|
||||
}
|
||||
private double residual(RealMatrix a, double[] x, double[] b) {
|
||||
double[] ax = a.operate(x);
|
||||
double sum = 0;
|
||||
for (int i = 0; i < ax.length; ++i) {
|
||||
sum += (ax[i] - b[i]) * (ax[i] - b[i]);
|
||||
}
|
||||
return Math.sqrt(sum);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue