MatrixUtils { public static RealMatrix blockInverse(RealMatrix m, int splitIndex) { final int n = m.getRowDimension(); if (m.getColumnDimension() != n) { throw new NonSquareMatrixException(m.getRowDimension(), m.getColumnDimension()); } final int splitIndex1 = splitIndex + 1; final RealMatrix a = m.getSubMatrix(0, splitIndex, 0, splitIndex); final RealMatrix b = m.getSubMatrix(0, splitIndex, splitIndex1, n - 1); final RealMatrix c = m.getSubMatrix(splitIndex1, n - 1, 0, splitIndex); final RealMatrix d = m.getSubMatrix(splitIndex1, n - 1, splitIndex1, n - 1); final SingularValueDecomposition aDec = new SingularValueDecomposition(a); final DecompositionSolver aSolver = aDec.getSolver(); if (!aSolver.isNonSingular()) { throw new SingularMatrixException(); } final RealMatrix aInv = aSolver.getInverse(); final SingularValueDecomposition dDec = new SingularValueDecomposition(d); final DecompositionSolver dSolver = dDec.getSolver(); if (!dSolver.isNonSingular()) { throw new SingularMatrixException(); } final RealMatrix dInv = dSolver.getInverse(); final RealMatrix tmp1 = a.subtract(b.multiply(dInv).multiply(c)); final SingularValueDecomposition tmp1Dec = new SingularValueDecomposition(tmp1); final DecompositionSolver tmp1Solver = tmp1Dec.getSolver(); if (!tmp1Solver.isNonSingular()) { throw new SingularMatrixException(); } final RealMatrix result00 = tmp1Solver.getInverse(); final RealMatrix tmp2 = d.subtract(c.multiply(aInv).multiply(b)); final SingularValueDecomposition tmp2Dec = new SingularValueDecomposition(tmp2); final DecompositionSolver tmp2Solver = tmp2Dec.getSolver(); if (!tmp2Solver.isNonSingular()) { throw new SingularMatrixException(); } final RealMatrix result11 = tmp2Solver.getInverse(); final RealMatrix result01 = aInv.multiply(b).multiply(result11).scalarMultiply(-1); final RealMatrix result10 = dInv.multiply(c).multiply(result00).scalarMultiply(-1); final RealMatrix result = new Array2DRowRealMatrix(n, n); result.setSubMatrix(result00.getData(), 0, 0); result.setSubMatrix(result01.getData(), 0, splitIndex1); result.setSubMatrix(result10.getData(), splitIndex1, 0); result.setSubMatrix(result11.getData(), splitIndex1, splitIndex1); return result; } private  MatrixUtils(); static RealMatrix createRealMatrix(final int rows, final int columns); static FieldMatrix<T> createFieldMatrix(final Field<T> field,
                                                                               final int rows,
                                                                               final int columns); static RealMatrix createRealMatrix(double[][] data); static FieldMatrix<T> createFieldMatrix(T[][] data); static RealMatrix createRealIdentityMatrix(int dimension); static FieldMatrix<T> createFieldIdentityMatrix(final Field<T> field, final int dimension); static RealMatrix createRealDiagonalMatrix(final double[] diagonal); static FieldMatrix<T> createFieldDiagonalMatrix(final T[] diagonal); static RealVector createRealVector(double[] data); static FieldVector<T> createFieldVector(final T[] data); static RealMatrix createRowRealMatrix(double[] rowData); static FieldMatrix<T> createRowFieldMatrix(final T[] rowData); static RealMatrix createColumnRealMatrix(double[] columnData); static FieldMatrix<T> createColumnFieldMatrix(final T[] columnData); static void checkSymmetric(RealMatrix matrix,
                                      double eps); static boolean isSymmetric(RealMatrix matrix,
                                      double eps); static void checkMatrixIndex(final AnyMatrix m,
                                        final int row, final int column); static void checkRowIndex(final AnyMatrix m, final int row); static void checkColumnIndex(final AnyMatrix m, final int column); static void checkSubMatrixIndex(final AnyMatrix m,
                                           final int startRow, final int endRow,
                                           final int startColumn, final int endColumn); static void checkSubMatrixIndex(final AnyMatrix m,
                                           final int[] selectedRows,
                                           final int[] selectedColumns); static void checkAdditionCompatible(final AnyMatrix left, final AnyMatrix right); static void checkSubtractionCompatible(final AnyMatrix left, final AnyMatrix right); static void checkMultiplicationCompatible(final AnyMatrix left, final AnyMatrix right); static Array2DRowRealMatrix fractionMatrixToRealMatrix(final FieldMatrix<Fraction> m); static Array2DRowRealMatrix bigFractionMatrixToRealMatrix(final FieldMatrix<BigFraction> m); static void serializeRealVector(final RealVector vector,
                                           final ObjectOutputStream oos); static void deserializeRealVector(final Object instance,
                                             final String fieldName,
                                             final ObjectInputStream ois); static void serializeRealMatrix(final RealMatrix matrix,
                                           final ObjectOutputStream oos); static void deserializeRealMatrix(final Object instance,
                                             final String fieldName,
                                             final ObjectInputStream ois); static void solveLowerTriangularSystem(RealMatrix rm, RealVector b); static void solveUpperTriangularSystem(RealMatrix rm, RealVector b); static RealMatrix blockInverse(RealMatrix m,
                                          int splitIndex); static RealMatrix inverse(RealMatrix matrix); static RealMatrix inverse(RealMatrix matrix, double threshold); static final RealMatrixFormat DEFAULT_FORMAT; static final RealMatrixFormat OCTAVE_FORMAT; }
@Test public void testBlockInverse() { final double[][] data = { { -1, 0, 123, 4 }, { -56, 78.9, -0.1, -23.4 }, { 5.67, 8, -9, 1011 }, { 12, 345, -67.8, 9 }, }; final RealMatrix m = new Array2DRowRealMatrix(data); final int len = data.length; final double tol = 1e-14; for (int splitIndex = 0; splitIndex < 3; splitIndex++) { final RealMatrix mInv = MatrixUtils.blockInverse(m, splitIndex); final RealMatrix id = m.multiply(mInv); for (int i = 0; i < len; i++) { for (int j = 0; j < len; j++) { final double entry = id.getEntry(i, j); if (i == j) { Assert.assertEquals("[" + i + "][" + j + "]", 1, entry, tol); } else { Assert.assertEquals("[" + i + "][" + j + "]", 0, entry, tol); } } } } }
@Test(expected=SingularMatrixException.class) public void testBlockInverseNonInvertible() { final double[][] data = { { -1, 0, 123, 4 }, { -56, 78.9, -0.1, -23.4 }, { 5.67, 8, -9, 1011 }, { 5.67, 8, -9, 1011 }, }; MatrixUtils.blockInverse(new Array2DRowRealMatrix(data), 2); }