4 ms·
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 experi
by vmarsy 8y ago
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.
- comicjk 8y agoIf 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.
- yonkshi 8y agoI 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.
- repsilat 8y agoA 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."
- greglindahl 8y agoDid you beat one of the sparse matrix packages? Or were you comparing with a library for dense matrices?
- yonkshi 8y agoWe 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)
- twic 8y agoI 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.
- ryanmonroe 8y agoRight, 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-spring2015/toeplitz.pdf 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
- deleted 8y ago[deleted]