Implementing Strassen's Algorithm

Question

I would like to multiply matrices very fast on a single core. I have looked around of the web and came across a few algorithms and found out Strassen's algorithm is the only one, that is actually implemented by people. I have looked on a few examples and came to the solution below. I made a simple benchmark which generates two randomly filled 500x500 matrices. Strassen's algorithm took 18 seconds, where the high school algorithm was done in 0.4 seconds. Other people where very promising after implementing the algorithm, so what is wrong with mine, how can I make it quicker?

// return C = A * B
private Matrix strassenTimes(Matrix B, int LEAFSIZE) {
    Matrix A = this;
    if (B.M != A.M || B.N != A.N) throw new RuntimeException("Illegal matrix dimensions.");

    if (N <= LEAFSIZE || M <= LEAFSIZE) {
        return A.times(B);
    }

    // make new sub-matrices
    int newAcols = (A.N + 1) / 2;
    int newArows = (A.M + 1) / 2;
    Matrix a11 = new Matrix(newArows, newAcols);
    Matrix a12 = new Matrix(newArows, newAcols);
    Matrix a21 = new Matrix(newArows, newAcols);
    Matrix a22 = new Matrix(newArows, newAcols);

    int newBcols = (B.N + 1) / 2;
    int newBrows = (B.M + 1) / 2;
    Matrix b11 = new Matrix(newBrows, newBcols);
    Matrix b12 = new Matrix(newBrows, newBcols);
    Matrix b21 = new Matrix(newBrows, newBcols);
    Matrix b22 = new Matrix(newBrows, newBcols);


    for (int i = 1; i <= newArows; i++) {
        for (int j = 1; j <= newAcols; j++) {
            a11.setElement(i, j, A.saveGet(i, j)); // top left
            a12.setElement(i, j, A.saveGet(i, j + newAcols)); // top right
            a21.setElement(i, j, A.saveGet(i + newArows, j)); // bottom left
            a22.setElement(i, j, A.saveGet(i + newArows, j + newAcols)); // bottom right
        }
    }

    for (int i = 1; i <= newBrows; i++) {
        for (int j = 1; j <= newBcols; j++) {
            b11.setElement(i, j, B.saveGet(i, j)); // top left
            b12.setElement(i, j, B.saveGet(i, j + newBcols)); // top right
            b21.setElement(i, j, B.saveGet(i + newBrows, j)); // bottom left
            b22.setElement(i, j, B.saveGet(i + newBrows, j + newBcols)); // bottom right
        }
    }

    Matrix aResult;
    Matrix bResult;

    aResult = a11.add(a22);
    bResult = b11.add(b22);
    Matrix p1 = aResult.strassenTimes(bResult, LEAFSIZE);

    aResult = a21.add(a22);
    Matrix p2 = aResult.strassenTimes(b11, LEAFSIZE);

    bResult = b12.minus(b22); // b12 - b22
    Matrix p3 = a11.strassenTimes(bResult, LEAFSIZE);

    bResult = b21.minus(b11); // b21 - b11
    Matrix p4 = a22.strassenTimes(bResult, LEAFSIZE);

    aResult = a11.add(a12); // a11 + a12
    Matrix p5 = aResult.strassenTimes(b22, LEAFSIZE);

    aResult = a21.minus(a11); // a21 - a11
    bResult = b11.add(b12); // b11 + b12
    Matrix p6 = aResult.strassenTimes(bResult, LEAFSIZE);

    aResult = a12.minus(a22); // a12 - a22
    bResult = b21.add(b22); // b21 + b22
    Matrix p7 = aResult.strassenTimes(bResult, LEAFSIZE);

    Matrix c12 = p3.add(p5); // c12 = p3 + p5
    Matrix c21 = p2.add(p4); // c21 = p2 + p4

    aResult = p1.add(p4); // p1 + p4
    bResult = aResult.add(p7); // p1 + p4 + p7
    Matrix c11 = bResult.minus(p5);

    aResult = p1.add(p3); // p1 + p3
    bResult = aResult.add(p6); // p1 + p3 + p6
    Matrix c22 = bResult.minus(p2);

    // Grouping the results obtained in a single matrix:
    int rows = c11.nrRows();
    int cols = c11.nrColumns();

    Matrix C = new Matrix(A.M, B.N);
    for (int i = 1; i <= A.M; i++) {
        for (int j = 1; j <= B.N; j++) {
            int el;
            if (i <= rows) {
                if (j <= cols) {
                    el = c11.get(i, j);
                } else {
                    el = c12.get(i, j - cols);
                }
            } else {
                if (j <= cols) {
                    el = c21.get(i - rows, j);
                } else {
                    el = c22.get(i - rows, j - rows);
                }
            }
            C.setElement(i, j, el);
        }
    }
    return C;
}

The little benchmark has the following code:

int AM, AN, BM, BN;
AM = 500;
AN = BM = 500;
BN = 500;
Matrix a = new Matrix(AM, AN);
Matrix b = new Matrix(BM, BN);

Random random = new Random();

for (int i = 1; i <= AM; i++) {
    for (int j = 1; j <= AN; j++) {
        a.setElement(i, j, random.nextInt(20));
    }
}
for (int i = 1; i <= BM; i++) {
    for (int j = 1; j <= BN; j++) {
        b.setElement(i, j, random.nextInt(20));
    }
}

System.out.println("strassen: A x B");
long tijd = System.currentTimeMillis();
Matrix c = a.strassenTimes(b);
System.out.println("time = " + (System.currentTimeMillis() - tijd));

System.out.println("normal: A x B");
tijd = System.currentTimeMillis();
Matrix d = a.times(b);
System.out.println("time = " + (System.currentTimeMillis() - tijd));

System.out.println("nr of different elements = " + c.compare(d));

With the following results:

strassen: A x B
time = 18372
normal: A x B
time = 308
nr of different elements = 0

I know it's a low of code, but I would be very happy if you guys help me out ;)

EDIT 1: For the sake of completeness I add some methods that is used by the above code.

public int get(int r, int c) {
    if (c > nrColumns() || r > nrRows() || c <= 0 || r <= 0) {
        throw new ArrayIndexOutOfBoundsException("matrix is of size (" +
                nrRows() + ", " + nrColumns() + "), but tries to set element(" + r + ", " + c + ")");
    }

    return content[r - 1][c - 1];
}

private int saveGet(int r, int c) {
    if (c > nrColumns() || r > nrRows() || c <= 0 || r <= 0) {
        return 0;
    }

    return content[r - 1][c - 1];
}

public void setElement(int r, int c, int n) {
    if (c > nrColumns() || r > nrRows() || c <= 0 || r <= 0) {
        throw new ArrayIndexOutOfBoundsException("matrix is of size (" +
                nrRows() + ", " + nrColumns() + "), but tries to set element(" + r + ", " + c + ")");
    }
    content[r - 1][c - 1] = n;
}

// return C = A + B
public Matrix add(Matrix B) {
    Matrix A = this;
    if (B.M != A.M || B.N != A.N) throw new RuntimeException("Illegal matrix dimensions.");
    Matrix C = new Matrix(M, N);
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            C.content[i][j] = A.content[i][j] + B.content[i][j];
        }
    }
    return C;
}

Answer 1

I should have choosen another leaf size for Strassen´s algorithm. Therefore I did a little experiment. It seems that leaf size 256 works best with the code included in the question. Below a plot with different leaf sizes with each time a random matrix of size 1025 x 1025.

I have compared Strassen´s algorithm with leaf size 256 with the trivial algorithm for matrix multiplication, to see if it´s actually an improvement. It turned out to be an improvement, see below the results on random matrices of different sizes (in steps of 10 and repeated 50 times for each size).

Below the code for the trivial algorithm for matrix multiplication:

// return C = A * B
public Matrix times(Matrix B) {
    Matrix A = this;
    if (A.N != B.M) throw new RuntimeException("Illegal matrix dimensions.");
    Matrix C = new Matrix(A.M, B.N);
    for (int i = 0; i < C.M; i++) {
        for (int j = 0; j < C.N; j++) {
            for (int k = 0; k < A.N; k++) {
                C.content[i][j] += (A.content[i][k] * B.content[k][j]);
            }
        }
    }
    return C;
}

It still think there can be done other improvements on the implementation, but it turned out that the leaf size is a very important factor. All experiments are done with a machine running on Ubuntu 14.04 with the following specifications:

CPU:    Intel(R) Core(TM) i7-2600K CPU @ 3.40GHz
Memory: 2 x 4GB DDR3 1333 MHz