@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); }
Matrix qt = getQ().transpose(); Matrix y = qt.times(B);
@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); }
Matrix q = qr.getQ(); printMatrix("q", q);
Matrix qt = getQ().transpose(); Matrix y = qt.times(B);
/** * 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; }
Matrix qm = qr.getQ();
Matrix q = qr.getQ(); printMatrix("q", q);
}, 5, 5); Matrix q = qr.getQ();