Author: tdunning
Date: Mon Feb 4 23:12:41 2013
New Revision: 1442428
URL: http://svn.apache.org/viewvc?rev=1442428&view=rev
Log:
MAHOUT-1148 - Re-implemented QR decomposition. On a 60x60 random matrix it is
about 8x faster now after the JIT warms up (takes 100 or so decompositions).
Added:
mahout/trunk/math/src/main/java/org/apache/mahout/math/OldQRDecomposition.java
- copied, changed from r1441391,
mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java
mahout/trunk/math/src/test/java/org/apache/mahout/math/OldQRDecompositionTest.java
- copied, changed from r1441391,
mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java
Modified:
mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java
mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java
Copied:
mahout/trunk/math/src/main/java/org/apache/mahout/math/OldQRDecomposition.java
(from r1441391,
mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java)
URL:
http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/OldQRDecomposition.java?p2=mahout/trunk/math/src/main/java/org/apache/mahout/math/OldQRDecomposition.java&p1=mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java&r1=1441391&r2=1442428&rev=1442428&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java
(original)
+++
mahout/trunk/math/src/main/java/org/apache/mahout/math/OldQRDecomposition.java
Mon Feb 4 23:12:41 2013
@@ -41,7 +41,7 @@ import java.util.Locale;
*/
/** partially deprecated until unit tests are in place. Until this time, this
class/interface is unsupported. */
-public class QRDecomposition {
+public class OldQRDecomposition {
/** Array for internal storage of decomposition. */
private final Matrix qr;
@@ -61,7 +61,7 @@ public class QRDecomposition {
* @throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
*/
- public QRDecomposition(Matrix a) {
+ public OldQRDecomposition(Matrix a) {
// Initialize.
qr = a.clone();
Modified:
mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java
URL:
http://svn.apache.org/viewvc/mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java?rev=1442428&r1=1442427&r2=1442428&view=diff
==============================================================================
--- mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java
(original)
+++ mahout/trunk/math/src/main/java/org/apache/mahout/math/QRDecomposition.java
Mon Feb 4 23:12:41 2013
@@ -23,8 +23,10 @@
*/
package org.apache.mahout.math;
+import com.google.common.collect.Lists;
import org.apache.mahout.math.function.Functions;
+import java.util.List;
import java.util.Locale;
@@ -40,91 +42,68 @@ import java.util.Locale;
returns <tt>false</tt>.
*/
-/** partially deprecated until unit tests are in place. Until this time, this
class/interface is unsupported. */
public class QRDecomposition {
-
- /** Array for internal storage of decomposition. */
- private final Matrix qr;
-
- /** Row and column dimensions. */
- private final int originalRows;
- private final int originalColumns;
-
- /** Array for internal storage of diagonal of R. */
- private final Vector rDiag;
+ private static final int N = 10;
+ private final Matrix q, r;
+ private final boolean fullRank;
+ private final int rows;
+ private final int columns;
/**
- * Constructs and returns a new QR decomposition object; computed by
Householder reflections; The decomposed matrices
- * can be retrieved via instance methods of the returned decomposition
object.
+ * Constructs and returns a new QR decomposition object; computed by
Householder reflections; The
+ * decomposed matrices can be retrieved via instance methods of the returned
decomposition
+ * object.
*
* @param a A rectangular matrix.
* @throws IllegalArgumentException if <tt>A.rows() < A.columns()</tt>.
*/
-
public QRDecomposition(Matrix a) {
- // Initialize.
- qr = a.clone();
- originalRows = a.numRows();
- originalColumns = a.numCols();
- rDiag = new DenseVector(originalColumns);
-
- // precompute and cache some views to avoid regenerating them time and
again
- Vector[] QRcolumnsPart = new Vector[originalColumns];
- for (int k = 0; k < originalColumns; k++) {
- QRcolumnsPart[k] = qr.viewColumn(k).viewPart(k, originalRows - k);
- }
-
- // Main loop.
- for (int k = 0; k < originalColumns; k++) {
- //DoubleMatrix1D QRcolk = QR.viewColumn(k).viewPart(k,m-k);
- // Compute 2-norm of k-th column without under/overflow.
- double nrm = 0;
- //if (k<m) nrm = QRcolumnsPart[k].aggregate(hypot,F.identity);
-
- for (int i = k; i < originalRows; i++) { // fixes bug reported by
[email protected]
- nrm = Algebra.hypot(nrm, qr.getQuick(i, k));
- }
-
-
- if (nrm != 0.0) {
- // Form k-th Householder vector.
- if (qr.getQuick(k, k) < 0) {
- nrm = -nrm;
- }
- QRcolumnsPart[k].assign(Functions.div(nrm));
- /*
- for (int i = k; i < m; i++) {
- QR[i][k] /= nrm;
+ rows = a.rowSize();
+ int min = Math.min(a.rowSize(), a.columnSize());
+ columns = a.columnSize();
+
+ Matrix qTmp = a.clone();
+
+ boolean fullRank = true;
+
+ r = new DenseMatrix(min, columns);
+
+ for (int i = 0; i < min; i++) {
+ Vector qi = qTmp.viewColumn(i);
+ double alpha = qi.norm(2);
+ if (Math.abs(alpha) > Double.MIN_VALUE) {
+ qi.assign(Functions.div(alpha));
+ } else {
+ if (Double.isInfinite(alpha) || Double.isNaN(alpha)) {
+ throw new ArithmeticException("Invalid intermediate result");
}
- */
-
- qr.setQuick(k, k, qr.getQuick(k, k) + 1);
+ fullRank = false;
+ }
+ r.set(i, i, alpha);
- // Apply transformation to remaining columns.
- for (int j = k + 1; j < originalColumns; j++) {
- Vector QRcolj = qr.viewColumn(j).viewPart(k, originalRows - k);
- double s = QRcolumnsPart[k].dot(QRcolj);
- /*
- // fixes bug reported by John Chambers
- DoubleMatrix1D QRcolj = QR.viewColumn(j).viewPart(k,m-k);
- double s = QRcolumnsPart[k].zDotProduct(QRcolumns[j]);
- double s = 0.0;
- for (int i = k; i < m; i++) {
- s += QR[i][k]*QR[i][j];
+ for (int j = i + 1; j < columns; j++) {
+ Vector qj = qTmp.viewColumn(j);
+ double norm = qj.norm(2);
+ if (Math.abs(norm) > Double.MIN_VALUE) {
+ double beta = qi.dot(qj);
+ r.set(i, j, beta);
+ if (j < min) {
+ qj.assign(qi, Functions.plusMult(-beta));
}
- */
- s = -s / qr.getQuick(k, k);
- //QRcolumnsPart[j].assign(QRcolumns[k], F.plusMult(s));
-
- for (int i = k; i < originalRows; i++) {
- qr.setQuick(i, j, qr.getQuick(i, j) + s * qr.getQuick(i, k));
+ } else {
+ if (Double.isInfinite(norm) || Double.isNaN(norm)) {
+ throw new ArithmeticException("Invalid intermediate result");
}
-
}
}
- rDiag.setQuick(k, -nrm);
}
+ if (columns > min) {
+ q = qTmp.viewPart(0, rows, 0, min).clone();
+ } else {
+ q = qTmp;
+ }
+ this.fullRank = fullRank;
}
/**
@@ -133,19 +112,6 @@ public class QRDecomposition {
* @return <tt>Q</tt>
*/
public Matrix getQ() {
- int columns = Math.min(originalColumns, originalRows);
- Matrix q = qr.like(originalRows, columns);
- for (int k = columns - 1; k >= 0; k--) {
- Vector QRcolk = qr.viewColumn(k).viewPart(k, originalRows - k);
- q.set(k, k, 1);
- for (int j = k; j < columns; j++) {
- if (qr.get(k, k) != 0) {
- Vector Qcolj = q.viewColumn(j).viewPart(k, originalRows - k);
- double s = -QRcolk.dot(Qcolj) / qr.get(k, k);
- Qcolj.assign(QRcolk, Functions.plusMult(s));
- }
- }
- }
return q;
}
@@ -155,19 +121,6 @@ public class QRDecomposition {
* @return <tt>R</tt>
*/
public Matrix getR() {
- int rows = Math.min(originalRows, originalColumns);
- Matrix r = qr.like(rows, originalColumns);
- for (int i = 0; i < rows; i++) {
- for (int j = 0; j < originalColumns; j++) {
- if (i < j) {
- r.setQuick(i, j, qr.getQuick(i, j));
- } else if (i == j) {
- r.setQuick(i, j, rDiag.getQuick(i));
- } else {
- r.setQuick(i, j, 0);
- }
- }
- }
return r;
}
@@ -177,12 +130,7 @@ public class QRDecomposition {
* @return true if <tt>R</tt>, and hence <tt>A</tt>, has full rank.
*/
public boolean hasFullRank() {
- for (int j = 0; j < originalColumns; j++) {
- if (rDiag.getQuick(j) == 0) {
- return false;
- }
- }
- return true;
+ return fullRank;
}
/**
@@ -193,12 +141,12 @@ public class QRDecomposition {
* @throws IllegalArgumentException if <tt>B.rows() != A.rows()</tt>.
*/
public Matrix solve(Matrix B) {
- if (B.numRows() != originalRows) {
+ if (B.numRows() != rows) {
throw new IllegalArgumentException("Matrix row dimensions must agree.");
}
- int columns = B.numCols();
- Matrix x = B.like(originalColumns, columns);
+ int cols = B.numCols();
+ Matrix x = B.like(columns, cols);
// this can all be done a bit more efficiently if we don't actually
// form explicit versions of Q^T and R but this code isn't soo bad
@@ -207,13 +155,13 @@ public class QRDecomposition {
Matrix y = qt.times(B);
Matrix r = getR();
- for (int k = Math.min(originalColumns, originalRows) - 1; k >= 0; k--) {
+ for (int k = Math.min(columns, rows) - 1; k >= 0; k--) {
// X[k,] = Y[k,] / R[k,k], note that X[k,] starts with 0 so += is same
as =
x.viewRow(k).assign(y.viewRow(k), Functions.plusMult(1 / r.get(k, k)));
// Y[0:(k-1),] -= R[0:(k-1),k] * X[k,]
Vector rColumn = r.viewColumn(k).viewPart(0, k);
- for (int c = 0; c < columns; c++) {
+ for (int c = 0; c < cols; c++) {
y.viewColumn(c).viewPart(0, k).assign(rColumn,
Functions.plusMult(-x.get(k, c)));
}
}
@@ -225,6 +173,39 @@ public class QRDecomposition {
*/
@Override
public String toString() {
- return String.format(Locale.ENGLISH, "QR(%d,%d,fullRank=%s)",
originalColumns, originalRows, hasFullRank());
+ return String.format(Locale.ENGLISH, "QR(%d x %d,fullRank=%s)", rows,
columns, hasFullRank());
+ }
+
+ public static void main(String[] args) {
+ Matrix a = new DenseMatrix(60, 60).assign(Functions.random());
+
+ int n = 0;
+ List<Integer> counts = Lists.newArrayList(10, 20, 50, 100, 200, 500, 1000,
2000, 5000);
+ for (int k : counts) {
+ double warmup = 0;
+ double other = 0;
+
+ n += k;
+ for (int i = 0; i < k; i++) {
+ QRDecomposition qr = new QRDecomposition(a);
+ warmup = Math.max(warmup,
qr.getQ().transpose().times(qr.getQ()).viewDiagonal().assign(Functions.plus(-1)).norm(1));
+ Matrix z = qr.getQ().times(qr.getR()).minus(a);
+ other = Math.max(other, z.aggregate(Functions.MIN, Functions.ABS));
+ }
+
+ double maxIdent = 0;
+ double maxError = 0;
+
+ long t0 = System.nanoTime();
+ for (int i = 0; i < N; i++) {
+ QRDecomposition qr = new QRDecomposition(a);
+
+ maxIdent = Math.max(maxIdent,
qr.getQ().transpose().times(qr.getQ()).viewDiagonal().assign(Functions.plus(-1)).norm(1));
+ Matrix z = qr.getQ().times(qr.getR()).minus(a);
+ maxError = Math.max(maxError, z.aggregate(Functions.MIN,
Functions.ABS));
+ }
+ System.out.printf("%d\t%.1f\t%g\t%g\t%g\n", n, (System.nanoTime() - t0)
/ 1e3 / N, maxIdent, maxError, warmup);
+// System.out.printf("%g, %g\n", maxIdent, maxError);
+ }
}
}
Copied:
mahout/trunk/math/src/test/java/org/apache/mahout/math/OldQRDecompositionTest.java
(from r1441391,
mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java)
URL:
http://svn.apache.org/viewvc/mahout/trunk/math/src/test/java/org/apache/mahout/math/OldQRDecompositionTest.java?p2=mahout/trunk/math/src/test/java/org/apache/mahout/math/OldQRDecompositionTest.java&p1=mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java&r1=1441391&r2=1442428&rev=1442428&view=diff
==============================================================================
---
mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java
(original)
+++
mahout/trunk/math/src/test/java/org/apache/mahout/math/OldQRDecompositionTest.java
Mon Feb 4 23:12:41 2013
@@ -21,7 +21,7 @@ import org.apache.mahout.math.function.D
import org.apache.mahout.math.function.Functions;
import org.junit.Test;
-public final class QRDecompositionTest extends MahoutTestCase {
+public final class OldQRDecompositionTest extends MahoutTestCase {
@Test
public void rank1() {
Matrix x = new DenseMatrix(3, 3);
@@ -29,7 +29,7 @@ public final class QRDecompositionTest e
x.viewRow(1).assign(new double[]{2, 4, 6});
x.viewRow(2).assign(new double[]{3, 6, 9});
- QRDecomposition qr = new QRDecomposition(x);
+ OldQRDecomposition qr = new OldQRDecomposition(x);
assertFalse(qr.hasFullRank());
assertEquals(0, new DenseVector(new double[]{3.741657, 7.483315,
11.22497}).aggregate(qr.getR().viewRow(0), Functions.PLUS, new
DoubleDoubleFunction() {
@Override
@@ -42,7 +42,7 @@ public final class QRDecompositionTest e
@Test
public void fullRankTall() {
Matrix x = matrix();
- QRDecomposition qr = new QRDecomposition(x);
+ OldQRDecomposition qr = new OldQRDecomposition(x);
assertTrue(qr.hasFullRank());
Matrix rRef = reshape(new double[]{
-2.99129686445138, 0, 0, 0, 0,
@@ -92,7 +92,7 @@ public final class QRDecompositionTest e
@Test
public void fullRankWide() {
Matrix x = matrix().transpose();
- QRDecomposition qr = new QRDecomposition(x);
+ OldQRDecomposition qr = new OldQRDecomposition(x);
assertFalse(qr.hasFullRank());
Matrix rActual = qr.getR();
Modified:
mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java
URL:
http://svn.apache.org/viewvc/mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java?rev=1442428&r1=1442427&r2=1442428&view=diff
==============================================================================
---
mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java
(original)
+++
mahout/trunk/math/src/test/java/org/apache/mahout/math/QRDecompositionTest.java
Mon Feb 4 23:12:41 2013
@@ -23,6 +23,21 @@ import org.junit.Test;
public final class QRDecompositionTest extends MahoutTestCase {
@Test
+ public void randomMatrix() {
+ Matrix a = new DenseMatrix(60, 60).assign(Functions.random());
+ QRDecomposition qr = new QRDecomposition(a);
+
+ // how close is Q to actually being orthornormal?
+ double maxIdent =
qr.getQ().transpose().times(qr.getQ()).viewDiagonal().assign(Functions.plus(-1)).norm(1);
+ assertEquals(0, maxIdent, 1e-13);
+
+ // how close is Q R to the original value of A?
+ Matrix z = qr.getQ().times(qr.getR()).minus(a);
+ double maxError = z.aggregate(Functions.MIN, Functions.ABS);
+ assertEquals(0, maxError, 1e-13);
+ }
+
+ @Test
public void rank1() {
Matrix x = new DenseMatrix(3, 3);
x.viewRow(0).assign(new double[]{1, 2, 3});
@@ -45,26 +60,28 @@ public final class QRDecompositionTest e
QRDecomposition qr = new QRDecomposition(x);
assertTrue(qr.hasFullRank());
Matrix rRef = reshape(new double[]{
- -2.99129686445138, 0, 0, 0, 0,
- -0.0282260628674372, -2.38850244769059, 0, 0, 0,
- 0.733739310355871, 1.48042000631646, 2.29051263117895, 0, 0,
- -0.0394082168269326, 0.282829484207801, -0.00438521041803086,
-2.90823198084203, 0,
- 0.923669647838536, 1.76679276072492, 0.637690104222683,
-0.225890909498753, -1.35732293800944},
- 5, 5);
+ -2.99129686445138, 0, 0, 0, 0,
+ -0.0282260628674372, -2.38850244769059, 0, 0, 0,
+ 0.733739310355871, 1.48042000631646, 2.29051263117895, 0, 0,
+ -0.0394082168269326, 0.282829484207801, -0.00438521041803086,
-2.90823198084203, 0,
+ 0.923669647838536, 1.76679276072492, 0.637690104222683,
-0.225890909498753, -1.35732293800944},
+ 5, 5);
Matrix r = qr.getR();
- assertEquals(rRef, r, 1.0e-8);
+
+ // check identity down to sign
+ assertEquals(0,
r.clone().assign(Functions.ABS).minus(rRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS,
Functions.IDENTITY), 1e-12);
Matrix qRef = reshape(new double[]{
- -0.165178287646573, 0.0510035857637869, 0.13985915987379,
-0.120173729496501,
- -0.453198314345324, 0.644400679630493, -0.503117990820608,
0.24968739845381,
- 0.323968339146224, -0.465266080134262, 0.276508948773268,
-0.687909700644343,
- 0.0544048888907195, -0.0166677718378263, 0.171309755790717,
0.310339001630029,
- 0.674790532821663, 0.0058166082200493, -0.381707516461884,
0.300504956413142,
- -0.105751091334003, 0.410450870871096, 0.31113446615821,
0.179338172684956,
- 0.361951807617901, 0.763921725548796, 0.380327892605634,
-0.287274944594054,
- 0.0311604042556675, 0.0386096858143961, 0.0387156960650472,
-0.232975755728917,
- 0.0358178276684149, 0.173105775703199, 0.327321867815603,
0.328671945345279,
- -0.36015879836344, -0.444261660176044, 0.09438499563253,
0.646216148583769
+ -0.165178287646573, 0.0510035857637869, 0.13985915987379,
-0.120173729496501,
+ -0.453198314345324, 0.644400679630493, -0.503117990820608,
0.24968739845381,
+ 0.323968339146224, -0.465266080134262, 0.276508948773268,
-0.687909700644343,
+ 0.0544048888907195, -0.0166677718378263, 0.171309755790717,
0.310339001630029,
+ 0.674790532821663, 0.0058166082200493, -0.381707516461884,
0.300504956413142,
+ -0.105751091334003, 0.410450870871096, 0.31113446615821,
0.179338172684956,
+ 0.361951807617901, 0.763921725548796, 0.380327892605634,
-0.287274944594054,
+ 0.0311604042556675, 0.0386096858143961, 0.0387156960650472,
-0.232975755728917,
+ 0.0358178276684149, 0.173105775703199, 0.327321867815603,
0.328671945345279,
+ -0.36015879836344, -0.444261660176044, 0.09438499563253,
0.646216148583769
}, 8, 5);
printMatrix("qRef", qRef);
@@ -72,15 +89,15 @@ public final class QRDecompositionTest e
Matrix q = qr.getQ();
printMatrix("q", q);
- assertEquals(qRef, q, 1.0e-8);
+ assertEquals(0,
q.clone().assign(Functions.ABS).minus(qRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS,
Functions.IDENTITY), 1e-12);
Matrix x1 = qr.solve(reshape(new double[]{
- -0.0178247686747641, 0.68631714634098, -0.335464858468858,
1.50249941751569,
- -0.669901640772149, -0.977025038942455, -1.18857546169856,
-1.24792900492054
+ -0.0178247686747641, 0.68631714634098, -0.335464858468858,
1.50249941751569,
+ -0.669901640772149, -0.977025038942455, -1.18857546169856,
-1.24792900492054
}, 8, 1));
Matrix xref = reshape(new double[]{
- -0.0127440093664874, 0.655825940180799, -0.100755415991702,
-0.0349559562697406,
- -0.190744297762028
+ -0.0127440093664874, 0.655825940180799, -0.100755415991702,
-0.0349559562697406,
+ -0.190744297762028
}, 5, 1);
printMatrix("x1", x1);
@@ -97,25 +114,26 @@ public final class QRDecompositionTest e
Matrix rActual = qr.getR();
Matrix rRef = reshape(new double[]{
- -2.42812464965842, 0, 0, 0, 0,
- 0.303587286111356, -2.91663643494775, 0, 0, 0,
- -0.201812474153156, -0.765485720168378, 1.09989373598954, 0, 0,
- 1.47980701097885, -0.637545820524326, -1.55519859337935,
0.844655127991726, 0,
- 0.0248883129453161, 0.00115010570270549, -0.236340588891252,
-0.092924118200147, 1.42910099545547,
- -1.1678472412429, 0.531245845248056, 0.351978196071514,
-1.03241474816555, -2.20223861735426,
- -0.887809959067632, 0.189731251982918, -0.504321849233586,
0.490484123999836, 1.21266692336743,
- -0.633888169775463, 1.04738559065986, 0.284041239547031,
0.578183510077156, -0.942314870832456
+ -2.42812464965842, 0, 0, 0, 0,
+ 0.303587286111356, -2.91663643494775, 0, 0, 0,
+ -0.201812474153156, -0.765485720168378, 1.09989373598954, 0, 0,
+ 1.47980701097885, -0.637545820524326, -1.55519859337935,
0.844655127991726, 0,
+ 0.0248883129453161, 0.00115010570270549, -0.236340588891252,
-0.092924118200147, 1.42910099545547,
+ -1.1678472412429, 0.531245845248056, 0.351978196071514,
-1.03241474816555, -2.20223861735426,
+ -0.887809959067632, 0.189731251982918, -0.504321849233586,
0.490484123999836, 1.21266692336743,
+ -0.633888169775463, 1.04738559065986, 0.284041239547031,
0.578183510077156, -0.942314870832456
}, 5, 8);
printMatrix("rRef", rRef);
printMatrix("rActual", rActual);
- assertEquals(rRef, rActual, 1.0e-8);
+ assertEquals(0,
rActual.clone().assign(Functions.ABS).minus(rRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS,
Functions.IDENTITY), 1e-12);
+// assertEquals(rRef, rActual, 1.0e-8);
Matrix qRef = reshape(new double[]{
- -0.203489262374627, 0.316761677948356, -0.784155643293468,
0.394321494579, -0.29641971170211,
- 0.0311283614803723, -0.34755265020736, 0.137138511478328,
0.848579887681972, 0.373287266507375,
- -0.39603700561249, -0.787812566647329, -0.377864833067864,
-0.275080943427399, 0.0636764674878229,
- 0.0763976893309043, -0.318551137554327, 0.286407036668598,
0.206004127289883, -0.876482672226889,
- 0.89159476695423, -0.238213616975551, -0.376141107880836,
-0.0794701657055114, 0.0227025098210165
+ -0.203489262374627, 0.316761677948356, -0.784155643293468,
0.394321494579, -0.29641971170211,
+ 0.0311283614803723, -0.34755265020736, 0.137138511478328,
0.848579887681972, 0.373287266507375,
+ -0.39603700561249, -0.787812566647329, -0.377864833067864,
-0.275080943427399, 0.0636764674878229,
+ 0.0763976893309043, -0.318551137554327, 0.286407036668598,
0.206004127289883, -0.876482672226889,
+ 0.89159476695423, -0.238213616975551, -0.376141107880836,
-0.0794701657055114, 0.0227025098210165
}, 5, 5);
Matrix q = qr.getQ();
@@ -123,16 +141,19 @@ public final class QRDecompositionTest e
printMatrix("qRef", qRef);
printMatrix("q", q);
- assertEquals(qRef, q, 1.0e-8);
+ assertEquals(0,
q.clone().assign(Functions.ABS).minus(qRef.clone().assign(Functions.ABS)).aggregate(Functions.PLUS,
Functions.IDENTITY), 1e-12);
+// assertEquals(qRef, q, 1.0e-8);
Matrix x1 = qr.solve(b());
Matrix xRef = reshape(new double[]{
- -0.182580239668147, -0.437233627652114, 0.138787653097464,
0.672934739896228, -0.131420217069083, 0, 0, 0
+ -0.182580239668147, -0.437233627652114, 0.138787653097464,
0.672934739896228, -0.131420217069083, 0, 0, 0
}, 8, 1);
printMatrix("xRef", xRef);
printMatrix("x", x1);
assertEquals(xRef, x1, 1.0e-8);
+
+ assertEquals(x, qr.getQ().times(qr.getR()), 1e-15);
}
private static void assertEquals(Matrix ref, Matrix actual, double epsilon) {
@@ -155,17 +176,17 @@ public final class QRDecompositionTest e
private static Matrix matrix() {
double[] values = {
- 0.494097293912641, -0.152566866170993, -0.418360266395271,
0.359475300232312,
- 1.35565069667582, -1.92759373242903, 1.50497526839076,
-0.746889132087904,
- -0.769136838293565, 1.10984954080986, -0.664389974392489,
1.6464660350229,
- -0.11715420616969, 0.0216221197371269, -0.394972730980765,
-0.748293157213142,
- 1.90402764664962, -0.638042862848559, -0.362336344669668,
-0.418261074380526,
- -0.494211543128429, 1.38828971158414, 0.597110366867923,
1.05341387608687,
- -0.957461740877418, -2.35528802598249, -1.03171458944128,
0.644319090271635,
- -0.0569108993041965, -0.14419465550881, -0.0456801828174936,
- 0.754694392571835, 0.719744008628535, -1.17873249802301,
-0.155887528905918,
- -1.5159868405466, 0.0918931582603128, 1.42179027361583,
-0.100495054250176,
- 0.0687986548485584
+ 0.494097293912641, -0.152566866170993, -0.418360266395271,
0.359475300232312,
+ 1.35565069667582, -1.92759373242903, 1.50497526839076,
-0.746889132087904,
+ -0.769136838293565, 1.10984954080986, -0.664389974392489,
1.6464660350229,
+ -0.11715420616969, 0.0216221197371269, -0.394972730980765,
-0.748293157213142,
+ 1.90402764664962, -0.638042862848559, -0.362336344669668,
-0.418261074380526,
+ -0.494211543128429, 1.38828971158414, 0.597110366867923,
1.05341387608687,
+ -0.957461740877418, -2.35528802598249, -1.03171458944128,
0.644319090271635,
+ -0.0569108993041965, -0.14419465550881, -0.0456801828174936,
+ 0.754694392571835, 0.719744008628535, -1.17873249802301,
-0.155887528905918,
+ -1.5159868405466, 0.0918931582603128, 1.42179027361583,
-0.100495054250176,
+ 0.0687986548485584
};
return reshape(values, 8, 5);
}
@@ -182,6 +203,7 @@ public final class QRDecompositionTest e
private static Matrix b() {
return reshape(new double[]
- {-0.0178247686747641, 0.68631714634098, -0.335464858468858,
1.50249941751569, -0.669901640772149}, 5, 1);
+ {-0.0178247686747641, 0.68631714634098, -0.335464858468858,
1.50249941751569, -0.669901640772149}, 5, 1);
}
}
+
