Matrix Computing is too slow

Question

I am developing the game that named Lights Out. So for solving this, i have to compute the answer of AX = B in modules 2. So, for this reason i choose jscience library. In this game the size of A is 25x25 matrix, X and B are both 25x1 matrix. I wrote the code such below :

AllLightOut.java class :

public class AllLightOut {
    public static final int SIZE = 5;

    public static double[] Action(int i, int j) {
        double[] change = new double[SIZE * SIZE];
        int count = 0;

        for (double[] d : Switch(new double[SIZE][SIZE], i, j))
            for (double e : d)
                change[count++] = e;

        return change;
    }

    public static double[][] MatrixA() {
        double[][] mat = new double[SIZE * SIZE][SIZE * SIZE];

        for (int i = 0; i < SIZE; i++)
            for (int j = 0; j < SIZE; j++)
                mat[i * SIZE + j] = Action(i, j);

        return mat;
    }

    public static SparseVector<ModuloInteger> ArrayToDenseVectorModule2(
            double[] array) {
        List<ModuloInteger> list = new ArrayList<ModuloInteger>();

        for (int i = 0; i < array.length; i++) {
            if (array[i] == 0)
                list.add(ModuloInteger.ZERO);
            else
                list.add(ModuloInteger.ONE);
        }

        return SparseVector.valueOf(DenseVector.valueOf(list),
                ModuloInteger.ZERO);
    }

    public static SparseMatrix<ModuloInteger> MatrixAModule2() {
        double[][] mat = MatrixA();
        List<DenseVector<ModuloInteger>> list = new ArrayList<DenseVector<ModuloInteger>>();

        for (int i = 0; i < mat.length; i++) {
            List<ModuloInteger> l = new ArrayList<ModuloInteger>();
            for (int j = 0; j < mat[i].length; j++) {
                if (mat[i][j] == 0)
                    l.add(ModuloInteger.ZERO);
                else
                    l.add(ModuloInteger.ONE);
            }

            list.add(DenseVector.valueOf(l));
        }

        return SparseMatrix.valueOf(DenseMatrix.valueOf(list),
                ModuloInteger.ZERO);
    }

    public static double[][] Switch(double[][] action, int i, int j) {
        action[i][j] = action[i][j] == 1 ? 0 : 1;

        if (i > 0)
            action[i - 1][j] = action[i - 1][j] == 1 ? 0 : 1;

        if (i < action.length - 1)
            action[i + 1][j] = action[i + 1][j] == 1 ? 0 : 1;

        if (j > 0)
            action[i][j - 1] = action[i][j - 1] == 1 ? 0 : 1;

        if (j < action.length - 1)
            action[i][j + 1] = action[i][j + 1] == 1 ? 0 : 1;

        return action;
    }
}

And the main class is as follow :

public class Main {
    public static void main(String[] args) {
        double[] bVec = new double[] { 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
                1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0 };

        SparseMatrix<ModuloInteger> matA = AllLightOut.MatrixAModule2();
        SparseVector<ModuloInteger> matB = AllLightOut
                .ArrayToDenseVectorModule2(bVec);

        ModuloInteger.setModulus(LargeInteger.valueOf(2));
    Vector<ModuloInteger> matX = matA.solve(matB);

        System.out.println(matX);
    }
}

I ran this program for about 30 minutes, but it had not result. Does my code include a fatal error or wrong ? Why it takes too long ?

Thanks for your attention :)

EDIT The slowdown happening in this line Matrix<ModuloInteger> matX = matA.inverse();. Note that the JScience benchmark result, speed for this library is very high, but i don't know why my program ran too slow!

EDIT2 Please note that when i try to SIZE = 3, i get the answer truly. For example: MatA :

{{1, 1, 0, 1, 0, 0, 0, 0, 0},
{1, 1, 1, 0, 1, 0, 0, 0, 0},
{0, 1, 1, 0, 0, 1, 0, 0, 0},
{1, 0, 0, 1, 1, 0, 1, 0, 0},
{0, 1, 0, 1, 1, 1, 0, 1, 0},
{0, 0, 1, 0, 1, 1, 0, 0, 1},
{0, 0, 0, 1, 0, 0, 1, 1, 0},
{0, 0, 0, 0, 1, 0, 1, 1, 1},
{0, 0, 0, 0, 0, 1, 0, 1, 1}}

MatB :

{1, 1, 1, 1, 1, 1, 1, 0, 0}

MatC :

{0, 0, 1, 1, 0, 0, 0, 0, 0}

But when i try SIZE = 5, slowdown occurred.

Answer 1

The slowdown happening in this line Matrix<ModuloInteger> matX = matA.inverse();

That would be because the coefficient matrix matA is not invertible for SIZE == 5 (or 4, 9, 11, 14, 16, ...?).

I'm a bit surprised the library didn't detect that and throw an exception. If the library tries to invert the matrix in solve(), that would have the same consequences.

A consequence of the singularity of the coefficient matrix for some sizes is that not all puzzles for these sizes are solvable, and the others have multiple solutions.

Since we're calculating modulo 2, we can use bits or booleans to model our states/toggles, using XOR for addition and & for multiplication. I have cooked up a simple solver using Gaussian elimination, maybe it helps you (I haven't spent much time thinking about the design, so it's not pretty):

public class Lights{
    private static final int SIZE = 5;

    private static boolean[] toggle(int i, int j) {
        boolean[] action = new boolean[SIZE*SIZE];
        int idx = i*SIZE+j;
        action[idx] = true;
        if (j > 0)      action[idx-1] = true;
        if (j < SIZE-1) action[idx+1] = true;
        if (i > 0)      action[idx-SIZE] = true;
        if (i < SIZE-1) action[idx+SIZE] = true;
        return action;
    }
    private static boolean[][] matrixA() {
        boolean[][] mat = new boolean[SIZE*SIZE][];
        for(int i = 0; i < SIZE; ++i) {
            for(int j = 0; j < SIZE; ++j) {
                mat[i*SIZE+j] = toggle(i,j);
            }
        }
        return mat;
    }
    private static void rotateR(boolean[] a, int r) {
        r %= a.length;
        if (r < 0) r += a.length;
        if (r == 0) return;
        boolean[] tmp = new boolean[r];
        for(int i = 0; i < r; ++i) {
            tmp[i] = a[i];
        }
        for(int i = 0; i < a.length - r; ++i) {
            a[i] = a[i+r];
        }
        for(int i = 0; i < r; ++i) {
            a[i + a.length - r] = tmp[i];
        }
    }
    private static void rotateR(boolean[][] a, int r) {
        r %= a.length;
        if (r < 0) r += a.length;
        if (r == 0) return;
        boolean[][] tmp = new boolean[r][];
        for(int i = 0; i < r; ++i) {
            tmp[i] = a[i];
        }
        for(int i = 0; i < a.length - r; ++i) {
            a[i] = a[i+r];
        }
        for(int i = 0; i < r; ++i) {
            a[i + a.length - r] = tmp[i];
        }
    }
    private static int count(boolean[] a) {
        int c = 0;
        for(int i = 0; i < a.length; ++i) {
            if (a[i]) ++c;
        }
        return c;
    }
    private static void swapBits(boolean[] a, int i, int j) {
        boolean tmp = a[i];
        a[i] = a[j];
        a[j] = tmp;
    }
    private static void addBit(boolean[] a, int i, int j) {
        a[j] ^= a[i];
    }
    private static void swapRows(boolean[][] a, int i, int j) {
        boolean[] tmp = a[i];
        a[i] = a[j];
        a[j] = tmp;
    }
    private static void xorb(boolean[] a, boolean[] b) {
        for(int i = 0; i < a.length; ++i) {
            a[i] ^= b[i];
        }
    }
    private static boolean[] boolBits(int bits, long param) {
        boolean[] bitArr = new boolean[bits];
        for(int i = 0; i < bits; ++i) {
            if (((param >> i) & 1L) != 0) {
                bitArr[i] = true;
            }
        }
        return bitArr;
    }
    private static boolean[] solve(boolean[][] m, boolean[] b) {
        // Move first SIZE rows to bottom, so that on the diagonal
        // above the lowest SIZE rows, there are unit matrices
        rotateR(m, SIZE);
        // modify right hand side accordingly
        rotateR(b,SIZE);
        // clean first SIZE*(SIZE-1) columns
        for(int i = 0; i < SIZE*(SIZE-1); ++i) {
            for(int k = 0; k < SIZE*SIZE; ++k) {
                if (k == i) continue;
                if (m[k][i]) {
                    xorb(m[k], m[i]);
                    b[k] ^= b[i];
                }
            }
        }
        // Now we have a block matrix
        /*
         * E 0 0 ... 0 X
         * 0 E 0 ... 0 X
         * 0 0 E ... 0 X
         * ...
         * 0 0 ... E 0 X
         * 0 0 ... 0 E X
         * 0 0 ... 0 0 Y
         *
         */
        // Bring Y to row-echelon form
        int i = SIZE*(SIZE-1), j, k, mi = i;
        while(mi < SIZE*SIZE){
            // Try to find a row with mi-th bit set
            for(j = i; j < SIZE*SIZE; ++j) {
                if (m[j][mi]) break;
            }
            if (j < SIZE*SIZE) {
                // Found one
                if (j > i) {
                    swapRows(m,i,j);
                    swapBits(b,i,j);
                }
                for(k = 0; k < SIZE*SIZE; ++k) {
                    if (k == i) continue;
                    if (m[k][mi]) {
                        xorb(m[k], m[i]);
                        b[k] ^= b[i];
                    }
                }
                // cleaned up column, good row, next
                ++i;
            }
            // Look at next column
            ++mi;
        }
        printMat(m,b);
        boolean[] best = b;
        if (i < SIZE*SIZE) {
            // We have zero-rows in the matrix,
            // check whether the puzzle is solvable at all,
            // i.e. all corresponding bits in the rhs are 0
            for(j = i; j < SIZE*SIZE; ++j) {
                if (b[j]) {
                    System.out.println("Puzzle not solvable, some lights must remain lit.");
                    break;
//                     throw new IllegalArgumentException("Puzzle is not solvable!");
                }
            }
            // Pretending it were solvable if not
            if (j < SIZE*SIZE) {
                System.out.println("Pretending the puzzle were solvable...");
                for(; j < SIZE*SIZE; ++j) {
                    b[j] = false;
                }
            }
            // Okay, puzzle is solvable, but there are several solutions
            // Let's try to find the one with the least toggles.

            // We have the canonical solution with last bits all zero
            int toggles = count(b);
            System.out.println(toggles + " toggles in canonical solution");
            int freeBits = SIZE*SIZE - i;
            long max = 1L << freeBits;
            System.out.println(freeBits + " free bits");
            // Check all combinations of free bits whether they produce
            // something better
            for(long param = 1; param < max; ++param) {
                boolean[] base = boolBits(freeBits,param);
                boolean[] c = new boolean[SIZE*SIZE];
                for(k = 0; k < freeBits; ++k) {
                    c[i+k] = base[k];
                }
                for(k = 0; k < i; ++k) {
                    for(j = 0; j < freeBits; ++j) {
                        c[k] ^= base[j] && m[k][j+i];
                    }
                }
                xorb(c,b);
                int t = count(c);
                if (t < toggles) {
                    System.out.printf("Found new best for param %x, %d toggles\n",param,t);
                    printMat(m,c,b);
                    toggles = t;
                    best = c;
                } else {
                    System.out.printf("%d toggles for parameter %x\n", t, param);
                }
            }
        }
        return best;
    }
    private static boolean[] parseLights(int[] lights) {
        int lim = lights.length;
        if (SIZE*SIZE < lim) lim = SIZE*SIZE;
        boolean[] b = new boolean[SIZE*SIZE];
        for(int i = 0; i < lim; ++i) {
            b[i] = (lights[i] != 0);
        }
        return b;
    }
    private static void printToggles(boolean[] s) {
        for(int i = 0; i < s.length; ++i) {
            if (s[i]) {
                System.out.print("(" + (i/SIZE + 1) + ", " + (i%SIZE + 1) + "); ");
            }
        }
        System.out.println();
    }
    private static void printMat(boolean[][] a, boolean[] rhs) {
        for(int i = 0; i < SIZE*SIZE; ++i) {
            for(int j = 0; j < SIZE*SIZE; ++j) {
                System.out.print((a[i][j] ? "1 " : "0 "));
            }
            System.out.println("| " + (rhs[i] ? "1" : "0"));
        }
    }
    private static void printMat(boolean[][] a, boolean[] sol, boolean[] rhs) {
        for(int i = 0; i < SIZE*SIZE; ++i) {
            for(int j = 0; j < SIZE*SIZE; ++j) {
                System.out.print((a[i][j] ? "1 " : "0 "));
            }
            System.out.println("| " + (sol[i] ? "1" : "0") + " | " + (rhs[i] ? "1" : "0"));
        }
    }
    private static void printGrid(boolean[] g) {
        for(int i = 0; i < SIZE; ++i) {
            for(int j = 0; j < SIZE; ++j) {
                System.out.print(g[i*SIZE+j] ? "1" : "0");
            }
            System.out.println();
        }
    }
    public static void main(String[] args) {
        int[] initialLights = new int[] { 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0,
                1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0 };
        boolean[] b = parseLights(initialLights);
        boolean[] b2 = b.clone();
        boolean[][] coefficients = matrixA();
        boolean[] toggles = solve(coefficients, b);
        printGrid(b2);
        System.out.println("--------");
        boolean[][] check = matrixA();
        boolean[] verify = new boolean[SIZE*SIZE];
        for(int i = 0; i < SIZE*SIZE; ++i) {
            if (toggles[i]) {
                xorb(verify, check[i]);
            }
        }
        printGrid(verify);
        xorb(b2,verify);
        if (count(b2) > 0) {
            System.out.println("Aww, shuck, screwed up!");
            printGrid(b2);
        }
        printToggles(toggles);
    }
}

Answer 2

You almost never want to calculate the actual inverse of a matrix if it can be avoided. Such operations are problematic and highly time consuming. Looking at the docs for JScience have you considered using the solve method? Something along the lines of matX = matA.solve(matB) should give you what you're looking for and I doubt they're using an inverse to calculate that, although I haven't dug that far into JScience so it's not impossible.