public SequentialBigSvd(Matrix A, int p) { // Y = A * \Omega y = A.times(new RandomTrinaryMatrix(A.columnSize(), p)); // R'R = Y' Y cd1 = new CholeskyDecomposition(y.transpose().times(y)); // B = Q" A = (Y R^{-1} )' A b = cd1.solveRight(y).transpose().times(A); // L L' = B B' cd2 = new CholeskyDecomposition(b.times(b.transpose())); // U_0 D V_0' = L svd = new SingularValueDecomposition(cd2.getL()); }
public CholeskyDecomposition(Matrix a, boolean pivot) { int rows = a.rowSize(); L = new PivotedMatrix(new DenseMatrix(rows, rows)); // must be square Preconditions.checkArgument(rows == a.columnSize(), "Must be a Square Matrix"); if (pivot) { decomposeWithPivoting(a); } else { decompose(a); } }
@Test public void testRankDeficient() { Matrix A = rank4Matrix(); CholeskyDecomposition cd = new CholeskyDecomposition(A); PivotedMatrix Ax = new PivotedMatrix(A, cd.getPivot()); CholeskyDecomposition cd2 = new CholeskyDecomposition(Ax, false); assertEquals(0, cd2.getL().times(cd2.getL().transpose()).minus(Ax).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); assertEquals(0, cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); Assert.assertFalse(cd.isPositiveDefinite()); Matrix L = cd.getL(); Matrix Abar = L.times(L.transpose()); double error = A.minus(Abar).aggregate(Functions.MAX, Functions.ABS); Assert.assertEquals(0, error, 1.0e-10); }
@Test public void test2() { // Test matrix from Nicholas Higham's paper at http://eprints.ma.man.ac.uk/1199/01/covered/MIMS_ep2008_116.pdf double[][] values = new double[3][]; values[0] = new double[]{1, -1, 1}; values[1] = new double[]{-1, 1, -1}; values[2] = new double[]{1, -1, 2}; Matrix A = new DenseMatrix(values); // without pivoting CholeskyDecomposition cd = new CholeskyDecomposition(A, false); assertEquals(0, cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); // with pivoting cd = new CholeskyDecomposition(A); assertEquals(0, cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); }
CholeskyDecomposition cd = new CholeskyDecomposition(A, type); Matrix L = cd.getL(); Matrix q = cd.solveLeft(z); Matrix id = q.times(q.transpose()); for (int i = 0; i < id.columnSize(); i++) { q = cd.solveRight(z.transpose()); id = q.transpose().times(q); for (int i = 0; i < id.columnSize(); i++) {
public Matrix getU() { // U = (Y inv(R)) U_0 return cd1.solveRight(y).times(svd.getU()); }
CholeskyDecomposition cd = new CholeskyDecomposition(A, type); Matrix L = cd.getL(); Matrix q = cd.solveLeft(z); Matrix id = q.times(q.transpose()); for (int i = 0; i < id.columnSize(); i++) { q = cd.solveRight(z.transpose()); id = q.transpose().times(q); for (int i = 0; i < id.columnSize(); i++) {
@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}); CholeskyDecomposition rr = new CholeskyDecomposition(x.transpose().times(x), false); assertEquals(0, new DenseVector(new double[]{3.741657, 7.483315, 11.22497}).aggregate(rr.getL().transpose().viewRow(0), Functions.PLUS, new DoubleDoubleFunction() { @Override public double apply(double arg1, double arg2) { return Math.abs(arg1) - Math.abs(arg2); } }), 1.0e-5); assertEquals(0, rr.getL().viewPart(0, 3, 1, 2).aggregate(Functions.PLUS, Functions.ABS), 1.0e-9); }
public Matrix getV() { // V = (B' inv(L')) V_0 return cd2.solveRight(b.transpose()).times(svd.getV()); } }
public SequentialBigSvd(Matrix A, int p) { // Y = A * \Omega y = A.times(new RandomTrinaryMatrix(A.columnSize(), p)); // R'R = Y' Y cd1 = new CholeskyDecomposition(y.transpose().times(y)); // B = Q" A = (Y R^{-1} )' A b = cd1.solveRight(y).transpose().times(A); // L L' = B B' cd2 = new CholeskyDecomposition(b.times(b.transpose())); // U_0 D V_0' = L svd = new SingularValueDecomposition(cd2.getL()); }
@Test public void test2() { // Test matrix from Nicholas Higham's paper at http://eprints.ma.man.ac.uk/1199/01/covered/MIMS_ep2008_116.pdf double[][] values = new double[3][]; values[0] = new double[]{1, -1, 1}; values[1] = new double[]{-1, 1, -1}; values[2] = new double[]{1, -1, 2}; Matrix A = new DenseMatrix(values); // without pivoting CholeskyDecomposition cd = new CholeskyDecomposition(A, false); assertEquals(0, cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); // with pivoting cd = new CholeskyDecomposition(A); assertEquals(0, cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); }
@Test public void testRankDeficient() { Matrix A = rank4Matrix(); CholeskyDecomposition cd = new CholeskyDecomposition(A); PivotedMatrix Ax = new PivotedMatrix(A, cd.getPivot()); CholeskyDecomposition cd2 = new CholeskyDecomposition(Ax, false); assertEquals(0, cd2.getL().times(cd2.getL().transpose()).minus(Ax).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); assertEquals(0, cd.getL().times(cd.getL().transpose()).minus(A).aggregate(Functions.PLUS, Functions.ABS), 1.0e-10); Assert.assertFalse(cd.isPositiveDefinite()); Matrix L = cd.getL(); Matrix Abar = L.times(L.transpose()); double error = A.minus(Abar).aggregate(Functions.MAX, Functions.ABS); Assert.assertEquals(0, error, 1.0e-10); }
public Matrix getU() { // U = (Y inv(R)) U_0 return cd1.solveRight(y).times(svd.getU()); }
public CholeskyDecomposition(Matrix a, boolean pivot) { int rows = a.rowSize(); L = new PivotedMatrix(new DenseMatrix(rows, rows)); // must be square Preconditions.checkArgument(rows == a.columnSize(), "Must be a Square Matrix"); if (pivot) { decomposeWithPivoting(a); } else { decompose(a); } }
public SequentialBigSvd(Matrix A, int p) { // Y = A * \Omega y = A.times(new RandomTrinaryMatrix(A.columnSize(), p)); // R'R = Y' Y cd1 = new CholeskyDecomposition(y.transpose().times(y)); // B = Q" A = (Y R^{-1} )' A b = cd1.solveRight(y).transpose().times(A); // L L' = B B' cd2 = new CholeskyDecomposition(b.times(b.transpose())); // U_0 D V_0' = L svd = new SingularValueDecomposition(cd2.getL()); }
@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}); CholeskyDecomposition rr = new CholeskyDecomposition(x.transpose().times(x), false); assertEquals(0, new DenseVector(new double[]{3.741657, 7.483315, 11.22497}).aggregate(rr.getL().transpose().viewRow(0), Functions.PLUS, new DoubleDoubleFunction() { @Override public double apply(double arg1, double arg2) { return Math.abs(arg1) - Math.abs(arg2); } }), 1.0e-5); assertEquals(0, rr.getL().viewPart(0, 3, 1, 2).aggregate(Functions.PLUS, Functions.ABS), 1.0e-9); }
public Matrix getU() { // U = (Y inv(R)) U_0 return cd1.solveRight(y).times(svd.getU()); }
public CholeskyDecomposition(Matrix a, boolean pivot) { int rows = a.rowSize(); L = new PivotedMatrix(new DenseMatrix(rows, rows)); // must be square Preconditions.checkArgument(rows == a.columnSize()); if (pivot) { decomposeWithPivoting(a); } else { decompose(a); } }
r2 = new CholeskyDecomposition(y2); Matrix yI = aI.times(omega); Matrix aIJ = aI.viewPart(0, aI.rowSize(), j, Math.min(columnsPerSlice, aI.columnSize() - j)); Matrix bIJ = r2.solveRight(yI).transpose().times(aIJ); addToSavedCopy(bFile(tmpDir, j), bIJ); l2 = new CholeskyDecomposition(b2); svd = new SingularValueDecomposition(l2.getL());
public void computeV(File tmpDir, int ncols) throws IOException { // step 5, compute pieces of V for (int j = 0; j < ncols; j += columnsPerSlice) { File bPath = bFile(tmpDir, j); if (bPath.exists()) { MatrixWritable m = new MatrixWritable(); try (DataInputStream in = new DataInputStream(new FileInputStream(bPath))) { m.readFields(in); } m.set(l2.solveRight(m.get().transpose()).times(svd.getV())); try (DataOutputStream out = new DataOutputStream(new FileOutputStream( new File(tmpDir, String.format("V-%s", bPath.getName().replaceAll(".*-", "")))))) { m.write(out); } } } }