Hacker News Reader

Comment by vmarsy | Hacker News Reader

Writing the fastest matrix multiplication you can is pretty fun. It probably won't be as fast a BLAS or MKL, but if you get close to it, it's a rewarding experience.

> see what else can be done I recommend reading this paper by Kazushige Goto

After doing the easiest change of loop reordering for already significant perf improvements, the paper linked ("Anatomy on a high performance matrix multiplication") is really excellent.

It will walk you through all the optimizations like tiling for reducing L1, L2, L3, and TLB cache misses, and leverage vectorization.

Then for squeezing out even more perf I remember there's the things like loop prefetching, loop unrolling, which you would expect the compiler with the best optimizations flags even targeted for the native architecture to take care of automatically. But you realize that it's not necessarily true.

View on HN

If your particular matrix problem has special properties that aren't covered by blas, you can wring out an apparently impossible improvement. I once found a matrix problem (periodic cubic splines) which turned out to have a solution in O(N) - the matrix columns were related such that it could be replaced by a vector.

Right, there are many tricks for special cases, even when the matrix is not sparse. One example is Toeplitz-Matrix multiplication as described in https://math.mit.edu/icg/resources/teaching/18.085-spring201...

That paper is light on details but, for example, in R you can define the function for toeplitz multiplication as below. This has given the matrix multiplication part of my code a >100x speedup before.

     '%t*%' <- function(A,v){
        n <- nrow(A)
        x <- as.matrix(c(A[1,], 0, A[1,][n:2]))
        p <- c(v, rep(0, n))
        h <- as.vector(fft(p)*fft(x))
        out <- Re(pracma::ifft(h)[1:n])
       return(matrix(out, n))
     }

     all.equal(A %t*% v,  A %*% v) #TRUE

I had a problem of inverting an N x N upper bidiagonal matrix a while ago. I spent ten minutes fruitlessly googling for papers about efficient algorithms for it, then thought about it for thirty seconds, and realised you can do it straightforwardly in N divisions and N subtractions.

I can attest this too, I was recently working with some massive and sparse binary matrix, with only 0.1% elements have values, it was magnitudes faster to index and sum the elements than to do a dot product. Not even a Tesla GPU was close to the performance of a custom matrix operation we wrote.

A lot of matrix routines will have sparse matrix support. Though often the line between wanting a sparse matrix and wanting some kind of adjacency-list graph structure can be pretty fine...

Also, only tangentially related, but if you can change your problem formulation to make the matrix of interest sparse(r), that can be a huge win. Two things that have helped me in the past:

- Rows or columns that have mostly the same (nonzero) coefficient throughout. You can normally do some simple substitution to turn them into zeros.

- Rows or columns that are naturally very similar to each other. You can often replace `y` with `y-x` or replace `x[i]` with `x[i]-x[i+1]` to turn one (or `n-1`) columns "mostly empty."

Did you beat one of the sparse matrix packages? Or were you comparing with a library for dense matrices?

We tested it against Tensorflow and scipy sparse matrices, we beat those as well, mostly because those general purpose sparse matrices do not optimize for binary matrix (eliminates the need for multiplication)

Reply on news.ycombinator.com