You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2012/09/23 21:34:03 UTC
svn commit: r1389129 - in /commons/proper/math/trunk/src: changes/
main/java/org/apache/commons/math3/linear/
test/java/org/apache/commons/math3/linear/
Author: tn
Date: Sun Sep 23 19:34:02 2012
New Revision: 1389129
URL: http://svn.apache.org/viewvc?rev=1389129&view=rev
Log:
[MATH-848] Fixed Schur transformation for certain input matrices, changed index parameter names to indicate their purpose.
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/SchurTransformerTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1389129&r1=1389128&r2=1389129&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Sun Sep 23 19:34:02 2012
@@ -52,6 +52,9 @@ If the output is not quite correct, chec
<body>
<release version="3.1" date="TBD" description="
">
+ <action dev="tn" type="fix" issue="MATH-848">
+ Fixed transformation to a Schur matrix for certain input matrices.
+ </action>
<action dev="erans" type="add" issue="MATH-860">
Added matrix "block inversion" (in "o.a.c.m.linear.MatrixUtils").
</action>
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java?rev=1389129&r1=1389128&r2=1389129&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/linear/SchurTransformer.java Sun Sep 23 19:34:02 2012
@@ -140,69 +140,66 @@ class SchurTransformer {
// Outer loop over eigenvalue index
int iteration = 0;
- int idx = n - 1;
- while (idx >= 0) {
+ int iu = n - 1;
+ while (iu >= 0) {
// Look for single small sub-diagonal element
- final int l = findSmallSubDiagonalElement(idx, norm);
+ final int il = findSmallSubDiagonalElement(iu, norm);
// Check for convergence
- if (l == idx) {
+ if (il == iu) {
// One root found
- matrixT[idx][idx] = matrixT[idx][idx] + shift.exShift;
- idx--;
+ matrixT[iu][iu] = matrixT[iu][iu] + shift.exShift;
+ iu--;
iteration = 0;
- } else if (l == idx - 1) {
+ } else if (il == iu - 1) {
// Two roots found
- shift.w = matrixT[idx][idx - 1] * matrixT[idx - 1][idx];
- double p = (matrixT[idx - 1][idx - 1] - matrixT[idx][idx]) / 2.0;
- double q = p * p + shift.w;
- double z = FastMath.sqrt(FastMath.abs(q));
- matrixT[idx][idx] = matrixT[idx][idx] + shift.exShift;
- matrixT[idx - 1][idx - 1] = matrixT[idx - 1][idx - 1] + shift.exShift;
- shift.x = matrixT[idx][idx];
+ double p = (matrixT[iu - 1][iu - 1] - matrixT[iu][iu]) / 2.0;
+ double q = p * p + matrixT[iu][iu - 1] * matrixT[iu - 1][iu];
+ matrixT[iu][iu] += shift.exShift;
+ matrixT[iu - 1][iu - 1] += shift.exShift;
if (q >= 0) {
+ double z = FastMath.sqrt(FastMath.abs(q));
if (p >= 0) {
z = p + z;
} else {
z = p - z;
}
- shift.x = matrixT[idx][idx - 1];
- double s = FastMath.abs(shift.x) + FastMath.abs(z);
- p = shift.x / s;
+ final double x = matrixT[iu][iu - 1];
+ final double s = FastMath.abs(x) + FastMath.abs(z);
+ p = x / s;
q = z / s;
- double r = FastMath.sqrt(p * p + q * q);
+ final double r = FastMath.sqrt(p * p + q * q);
p = p / r;
q = q / r;
// Row modification
- for (int j = idx - 1; j < n; j++) {
- z = matrixT[idx - 1][j];
- matrixT[idx - 1][j] = q * z + p * matrixT[idx][j];
- matrixT[idx][j] = q * matrixT[idx][j] - p * z;
+ for (int j = iu - 1; j < n; j++) {
+ z = matrixT[iu - 1][j];
+ matrixT[iu - 1][j] = q * z + p * matrixT[iu][j];
+ matrixT[iu][j] = q * matrixT[iu][j] - p * z;
}
// Column modification
- for (int i = 0; i <= idx; i++) {
- z = matrixT[i][idx - 1];
- matrixT[i][idx - 1] = q * z + p * matrixT[i][idx];
- matrixT[i][idx] = q * matrixT[i][idx] - p * z;
+ for (int i = 0; i <= iu; i++) {
+ z = matrixT[i][iu - 1];
+ matrixT[i][iu - 1] = q * z + p * matrixT[i][iu];
+ matrixT[i][iu] = q * matrixT[i][iu] - p * z;
}
// Accumulate transformations
for (int i = 0; i <= n - 1; i++) {
- z = matrixP[i][idx - 1];
- matrixP[i][idx - 1] = q * z + p * matrixP[i][idx];
- matrixP[i][idx] = q * matrixP[i][idx] - p * z;
+ z = matrixP[i][iu - 1];
+ matrixP[i][iu - 1] = q * z + p * matrixP[i][iu];
+ matrixP[i][iu] = q * matrixP[i][iu] - p * z;
}
}
- idx -= 2;
+ iu -= 2;
iteration = 0;
} else {
// No convergence yet
-
- computeShift(l, idx, iteration, shift);
+ computeShift(il, iu, iteration, shift);
// stop transformation after too many iterations
if (++iteration > MAX_ITERATIONS) {
@@ -210,43 +207,11 @@ class SchurTransformer {
MAX_ITERATIONS);
}
- // Look for two consecutive small sub-diagonal elements
- int m = idx - 2;
-
// the initial houseHolder vector for the QR step
final double[] hVec = new double[3];
- while (m >= l) {
- double z = matrixT[m][m];
- hVec[2] = shift.x - z;
- double s = shift.y - z;
- hVec[0] = (hVec[2] * s - shift.w) / matrixT[m + 1][m] + matrixT[m][m + 1];
- hVec[1] = matrixT[m + 1][m + 1] - z - hVec[2] - s;
- hVec[2] = matrixT[m + 2][m + 1];
- s = FastMath.abs(hVec[0]) + FastMath.abs(hVec[1]) + FastMath.abs(hVec[2]);
-
- if (m == l) {
- break;
- }
-
- for (int i = 0; i < hVec.length; i++) {
- hVec[i] /= s;
- }
-
- final double lhs = FastMath.abs(matrixT[m][m - 1]) *
- (FastMath.abs(hVec[1]) + FastMath.abs(hVec[2]));
-
- final double rhs = FastMath.abs(hVec[0]) *
- (FastMath.abs(matrixT[m - 1][m - 1]) + FastMath.abs(z) +
- FastMath.abs(matrixT[m + 1][m + 1]));
-
- if (lhs < epsilon * rhs) {
- break;
- }
- m--;
- }
-
- performDoubleQRStep(l, m, idx, shift, hVec);
+ final int im = initQRStep(il, iu, shift, hVec);
+ performDoubleQRStep(il, im, iu, shift, hVec);
}
}
}
@@ -278,7 +243,7 @@ class SchurTransformer {
int l = startIdx;
while (l > 0) {
double s = FastMath.abs(matrixT[l - 1][l - 1]) + FastMath.abs(matrixT[l][l]);
- if (Precision.equals(s, 0.0, epsilon)) {
+ if (s == 0.0) {
s = norm;
}
if (FastMath.abs(matrixT[l][l - 1]) < epsilon * s) {
@@ -312,8 +277,9 @@ class SchurTransformer {
for (int i = 0; i <= idx; i++) {
matrixT[i][i] -= shift.x;
}
- double s = FastMath.abs(matrixT[idx][idx - 1]) + FastMath.abs(matrixT[idx - 1][idx - 2]);
- shift.x = shift.y = 0.75 * s;
+ final double s = FastMath.abs(matrixT[idx][idx - 1]) + FastMath.abs(matrixT[idx - 1][idx - 2]);
+ shift.x = 0.75 * s;
+ shift.y = 0.75 * s;
shift.w = -0.4375 * s * s;
}
@@ -321,7 +287,7 @@ class SchurTransformer {
if (iteration == 30) {
double s = (shift.y - shift.x) / 2.0;
s = s * s + shift.w;
- if (Precision.compareTo(s, 0.0d, epsilon) > 0) {
+ if (s > 0.0) {
s = FastMath.sqrt(s);
if (shift.y < shift.x) {
s = -s;
@@ -337,15 +303,53 @@ class SchurTransformer {
}
/**
+ * Initialize the householder vectors for the QR step.
+ *
+ * @param il the index of the small sub-diagonal element
+ * @param iu the current eigenvalue index
+ * @param shift shift information holder
+ * @param hVec the initial houseHolder vector
+ * @return the start index for the QR step
+ */
+ private int initQRStep(int il, final int iu, final ShiftInfo shift, double[] hVec) {
+ // Look for two consecutive small sub-diagonal elements
+ int im = iu - 2;
+ while (im >= il) {
+ final double z = matrixT[im][im];
+ final double r = shift.x - z;
+ double s = shift.y - z;
+ hVec[0] = (r * s - shift.w) / matrixT[im + 1][im] + matrixT[im][im + 1];
+ hVec[1] = matrixT[im + 1][im + 1] - z - r - s;
+ hVec[2] = matrixT[im + 2][im + 1];
+
+ if (im == il) {
+ break;
+ }
+
+ final double lhs = FastMath.abs(matrixT[im][im - 1]) * (FastMath.abs(hVec[1]) + FastMath.abs(hVec[2]));
+ final double rhs = FastMath.abs(hVec[0]) * (FastMath.abs(matrixT[im - 1][im - 1]) +
+ FastMath.abs(z) +
+ FastMath.abs(matrixT[im + 1][im + 1]));
+
+ if (lhs < epsilon * rhs) {
+ break;
+ }
+ im--;
+ }
+
+ return im;
+ }
+
+ /**
* Perform a double QR step involving rows l:idx and columns m:n
*
- * @param l the index of the small sub-diagonal element
- * @param m the start index for the QR step
- * @param idx the current eigenvalue index
+ * @param il the index of the small sub-diagonal element
+ * @param im the start index for the QR step
+ * @param iu the current eigenvalue index
* @param shift shift information holder
* @param hVec the initial houseHolder vector
*/
- private void performDoubleQRStep(final int l, final int m, final int idx,
+ private void performDoubleQRStep(final int il, final int im, final int iu,
final ShiftInfo shift, final double[] hVec) {
final int n = matrixT.length;
@@ -353,9 +357,9 @@ class SchurTransformer {
double q = hVec[1];
double r = hVec[2];
- for (int k = m; k <= idx - 1; k++) {
- boolean notlast = k != idx - 1;
- if (k != m) {
+ for (int k = im; k <= iu - 1; k++) {
+ boolean notlast = k != (iu - 1);
+ if (k != im) {
p = matrixT[k][k - 1];
q = matrixT[k + 1][k - 1];
r = notlast ? matrixT[k + 2][k - 1] : 0.0;
@@ -366,17 +370,17 @@ class SchurTransformer {
r = r / shift.x;
}
}
- if (Precision.equals(shift.x, 0.0, epsilon)) {
+ if (shift.x == 0.0) {
break;
}
double s = FastMath.sqrt(p * p + q * q + r * r);
if (p < 0.0) {
s = -s;
}
- if (!Precision.equals(s, 0.0, epsilon)) {
- if (k != m) {
+ if (s != 0.0) {
+ if (k != im) {
matrixT[k][k - 1] = -s * shift.x;
- } else if (l != m) {
+ } else if (il != im) {
matrixT[k][k - 1] = -matrixT[k][k - 1];
}
p = p + s;
@@ -398,7 +402,7 @@ class SchurTransformer {
}
// Column modification
- for (int i = 0; i <= FastMath.min(idx, k + 3); i++) {
+ for (int i = 0; i <= FastMath.min(iu, k + 3); i++) {
p = shift.x * matrixT[i][k] + shift.y * matrixT[i][k + 1];
if (notlast) {
p = p + z * matrixT[i][k + 2];
@@ -423,9 +427,9 @@ class SchurTransformer {
} // k loop
// clean up pollution due to round-off errors
- for (int i = m+2; i <= idx; i++) {
+ for (int i = im + 2; i <= iu; i++) {
matrixT[i][i-2] = 0.0;
- if (i > m+2) {
+ if (i > im + 2) {
matrixT[i][i-3] = 0.0;
}
}
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java?rev=1389129&r1=1389128&r2=1389129&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/EigenDecompositionTest.java Sun Sep 23 19:34:02 2012
@@ -408,6 +408,21 @@ public class EigenDecompositionTest {
}
}
+ @Test
+ public void testMath848() {
+ double[][] data = {
+ { 0.1849449280, -0.0646971046, 0.0774755812, -0.0969651755, -0.0692648806, 0.3282344352, -0.0177423074, 0.2063136340},
+ {-0.0742700134, -0.0289063030, -0.0017269460, -0.0375550146, -0.0487737922, -0.2616837868, -0.0821201295, -0.2530000167},
+ { 0.2549910127, 0.0995733692, -0.0009718388, 0.0149282808, 0.1791878897, -0.0823182816, 0.0582629256, 0.3219545182},
+ {-0.0694747557, -0.1880649148, -0.2740630911, 0.0720096468, -0.1800836914, -0.3518996425, 0.2486747833, 0.6257938167},
+ { 0.0536360918, -0.1339297778, 0.2241579764, -0.0195327484, -0.0054103808, 0.0347564518, 0.5120802482, -0.0329902864},
+ {-0.5933332356, -0.2488721082, 0.2357173629, 0.0177285473, 0.0856630593, -0.3567126300, -0.1600668126, -0.1010899621},
+ {-0.0514349819, -0.0854319435, 0.1125050061, 0.0063453560, -0.2250000688, -0.2209343090, 0.1964623477, -0.1512329924},
+ { 0.0197395947, -0.1997170581, -0.1425959019, -0.2749477910, -0.0969467073, 0.0603688520, -0.2826905192, 0.1794315473}};
+ RealMatrix m = MatrixUtils.createRealMatrix(data);
+ checkUnsymmetricMatrix(m);
+ }
+
/**
* Checks that the eigen decomposition of a general (unsymmetric) matrix is valid by
* checking: A*V = V*D
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/SchurTransformerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/SchurTransformerTest.java?rev=1389129&r1=1389128&r2=1389129&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/SchurTransformerTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/linear/SchurTransformerTest.java Sun Sep 23 19:34:02 2012
@@ -130,6 +130,22 @@ public class SchurTransformerTest {
}
}
+ @Test
+ public void testMath848() {
+ double[][] data = {
+ { 0.1849449280, -0.0646971046, 0.0774755812, -0.0969651755, -0.0692648806, 0.3282344352, -0.0177423074, 0.2063136340},
+ {-0.0742700134, -0.0289063030, -0.0017269460, -0.0375550146, -0.0487737922, -0.2616837868, -0.0821201295, -0.2530000167},
+ { 0.2549910127, 0.0995733692, -0.0009718388, 0.0149282808, 0.1791878897, -0.0823182816, 0.0582629256, 0.3219545182},
+ {-0.0694747557, -0.1880649148, -0.2740630911, 0.0720096468, -0.1800836914, -0.3518996425, 0.2486747833, 0.6257938167},
+ { 0.0536360918, -0.1339297778, 0.2241579764, -0.0195327484, -0.0054103808, 0.0347564518, 0.5120802482, -0.0329902864},
+ {-0.5933332356, -0.2488721082, 0.2357173629, 0.0177285473, 0.0856630593, -0.3567126300, -0.1600668126, -0.1010899621},
+ {-0.0514349819, -0.0854319435, 0.1125050061, 0.0063453560, -0.2250000688, -0.2209343090, 0.1964623477, -0.1512329924},
+ { 0.0197395947, -0.1997170581, -0.1425959019, -0.2749477910, -0.0969467073, 0.0603688520, -0.2826905192, 0.1794315473}};
+ RealMatrix m = MatrixUtils.createRealMatrix(data);
+ RealMatrix s = checkAEqualPTPt(m);
+ checkSchurForm(s);
+ }
+
///////////////////////////////////////////////////////////////////////////
// Test helpers
///////////////////////////////////////////////////////////////////////////
@@ -143,7 +159,7 @@ public class SchurTransformerTest {
RealMatrix result = p.multiply(t).multiply(pT);
double norm = result.subtract(matrix).getNorm();
- Assert.assertEquals(0, norm, 1.0e-10);
+ Assert.assertEquals(0, norm, 1.0e-9);
return t;
}