private static Vector solve(Matrix Ai, Matrix Vi) { return new QRDecomposition(Ai).solve(Vi).viewColumn(0); }
/** * Returns a rough string rendition of a QR. */ @Override public String toString() { return String.format(Locale.ENGLISH, "QR(%d x %d,fullRank=%s)", rows, columns, hasFullRank()); } }
@Test public void fullRankWide() { Matrix x = matrix().transpose(); QRDecomposition qr = new QRDecomposition(x); assertTrue(qr.hasFullRank()); Matrix rActual = qr.getR(); }, 5, 5); Matrix q = qr.getQ(); Matrix x1 = qr.solve(b()); Matrix xRef = reshape(new double[]{ -0.182580239668147, -0.437233627652114, 0.138787653097464, 0.672934739896228, -0.131420217069083, 0, 0, 0 assertEquals(xRef, x1, 1.0e-8); assertEquals(x, qr.getQ().times(qr.getR()), 1.0e-15);
@Test public void rank1() { Matrix x = new DenseMatrix(3, 3); x.viewRow(0).assign(new double[]{1, 2, 3}); x.viewRow(1).assign(new double[]{2, 4, 6}); x.viewRow(2).assign(new double[]{3, 6, 9}); QRDecomposition qr = new QRDecomposition(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 public double apply(double arg1, double arg2) { return Math.abs(arg1) - Math.abs(arg2); } }), 1.0e-5); }
@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, 1.0e-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, 1.0e-13); }
@Test public void testProjection() { Vector v1 = new DenseVector(10).assign(Functions.random()); WeightedVector v2 = new WeightedVector(v1, v1, 31); assertEquals(v1.dot(v1), v2.getWeight(), 1.0e-13); assertEquals(31, v2.getIndex()); Matrix y = new DenseMatrix(10, 4).assign(Functions.random()); Matrix q = new QRDecomposition(y.viewPart(0, 10, 0, 3)).getQ(); Vector nullSpace = y.viewColumn(3).minus(q.times(q.transpose().times(y.viewColumn(3)))); WeightedVector v3 = new WeightedVector(q.viewColumn(0).plus(q.viewColumn(1)), nullSpace, 1); assertEquals(0, v3.getWeight(), 1.0e-13); Vector qx = q.viewColumn(0).plus(q.viewColumn(1)).normalize(); WeightedVector v4 = new WeightedVector(qx, q.viewColumn(0), 2); assertEquals(Math.sqrt(0.5), v4.getWeight(), 1.0e-13); WeightedVector v5 = WeightedVector.project(q.viewColumn(0), qx); assertEquals(Math.sqrt(0.5), v5.getWeight(), 1.0e-13); }
@Override public QR decompose(Matrix a) { return new QRDecomposition(a); } }, s1, a, "new");
@Test public void fullRankTall() { Matrix x = matrix(); QRDecomposition qr = new QRDecomposition(x); assertTrue(qr.hasFullRank()); Matrix rRef = reshape(new double[]{ -2.99129686445138, 0, 0, 0, 0, 0.923669647838536, 1.76679276072492, 0.637690104222683, -0.225890909498753, -1.35732293800944}, 5, 5); Matrix r = qr.getR(); Matrix q = qr.getQ(); printMatrix("q", q); Matrix x1 = qr.solve(reshape(new double[]{ -0.0178247686747641, 0.68631714634098, -0.335464858468858, 1.50249941751569, -0.669901640772149, -0.977025038942455, -1.18857546169856, -1.24792900492054
@Test public void rank1() { Matrix x = new DenseMatrix(3, 3); x.viewRow(0).assign(new double[]{1, 2, 3}); x.viewRow(1).assign(new double[]{2, 4, 6}); x.viewRow(2).assign(new double[]{3, 6, 9}); QRDecomposition qr = new QRDecomposition(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 public double apply(double arg1, double arg2) { return Math.abs(arg1) - Math.abs(arg2); } }), 1.0e-5); }
for (int i = 0; i < n; i++) y.viewRow(i).assign(s_o, Functions.MINUS); QRDecomposition qr = new QRDecomposition(y); Matrix qm = qr.getQ();
@Test public void fullRankTall() { Matrix x = matrix(); QRDecomposition qr = new QRDecomposition(x); assertTrue(qr.hasFullRank()); Matrix rRef = reshape(new double[]{ -2.99129686445138, 0, 0, 0, 0, 0.923669647838536, 1.76679276072492, 0.637690104222683, -0.225890909498753, -1.35732293800944}, 5, 5); Matrix r = qr.getR(); assertEquals(rRef, r, 1.0e-8); Matrix q = qr.getQ(); printMatrix("q", q); Matrix x1 = qr.solve(reshape(new double[]{ -0.0178247686747641, 0.68631714634098, -0.335464858468858, 1.50249941751569, -0.669901640772149, -0.977025038942455, -1.18857546169856, -1.24792900492054
private static Vector solve(Matrix A, Matrix y) { return new QRDecomposition(A).solve(y).viewColumn(0); }
/** * Least squares solution of <tt>A*X = B</tt>; <tt>returns X</tt>. * * @param B A matrix with as many rows as <tt>A</tt> and any number of columns. * @return <tt>X</tt> that minimizes the two norm of <tt>Q*R*X - B</tt>. * @throws IllegalArgumentException if <tt>B.rows() != A.rows()</tt>. */ public Matrix solve(Matrix B) { if (B.numRows() != originalRows) { throw new IllegalArgumentException("Matrix row dimensions must agree."); } int columns = B.numCols(); Matrix x = B.like(originalColumns, columns); // 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 // and it is much easier to understand Matrix qt = getQ().transpose(); Matrix y = qt.times(B); Matrix r = getR(); for (int k = Math.min(originalColumns, originalRows) - 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++) { y.viewColumn(c).viewPart(0, k).assign(rColumn, Functions.plusMult(-x.get(k, c))); } } return x; }
/** * Returns a rough string rendition of a QR. */ @Override public String toString() { return String.format(Locale.ENGLISH, "QR(%d x %d,fullRank=%s)", rows, columns, hasFullRank()); } }
@Test public void fullRankWide() { Matrix x = matrix().transpose(); QRDecomposition qr = new QRDecomposition(x); assertFalse(qr.hasFullRank()); Matrix rActual = qr.getR(); }, 5, 5); Matrix q = qr.getQ(); Matrix x1 = qr.solve(b()); Matrix xRef = reshape(new double[]{ -0.182580239668147, -0.437233627652114, 0.138787653097464, 0.672934739896228, -0.131420217069083, 0, 0, 0
/** * Predict the power law growth in number of unique samples from the first few data points. * Also check that the fitted growth coefficient is about right. * * @param m * @param currentIndex Total data points seen so far. Unique values should be log(currentIndex)*expectedCoefficient + offset. * @param expectedCoefficient What slope do we expect. * @return The predicted value for log(currentIndex) */ private static double predictSize(Matrix m, int currentIndex, double expectedCoefficient) { int rows = m.rowSize(); Matrix a = m.viewPart(0, rows, 1, 2); Matrix b = m.viewPart(0, rows, 0, 1); Matrix ata = a.transpose().times(a); Matrix atb = a.transpose().times(b); QRDecomposition s = new QRDecomposition(ata); Matrix r = s.solve(atb).transpose(); assertEquals(expectedCoefficient, r.get(0, 0), 0.2); return r.times(new DenseVector(new double[]{Math.log(currentIndex), 1})).get(0); }
/** * Returns a rough string rendition of a QR. */ @Override public String toString() { return String.format(Locale.ENGLISH, "QR(%d,%d,fullRank=%s)", originalColumns, originalRows, hasFullRank()); } }
private static Vector solve(Matrix Ai, Matrix Vi) { return new QRDecomposition(Ai).solve(Vi).viewColumn(0); }