Author: luc
Date: Thu Dec 31 17:52:16 2009
New Revision: 894908
URL: http://svn.apache.org/viewvc?rev=894908&view=rev
Log:
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
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java?rev=894908&r1=894907&r2=894908&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java
Thu Dec 31 17:52:16 2009
@@ -159,27 +159,28 @@
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 @@
// 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];
}
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=894908&r1=894907&r2=894908&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Thu Dec 31 17:52:16 2009
@@ -39,10 +39,18 @@
</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">
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java?rev=894908&r1=894907&r2=894908&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/linear/SingularValueSolverTest.java
Thu Dec 31 17:52:16 2009
@@ -139,6 +139,32 @@
}
@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[][] {
{ 1.0, 2.0, 3.0 }, { 2.0, 3.0, 4.0 }, { 3.0, 5.0, 7.0 }
@@ -149,44 +175,38 @@
};
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 @@
});
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);
}
}
