6 ms·

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 mat

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.

yonkshi 8y ago
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.
- repsilat 8y ago
  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."
- greglindahl 8y ago
  Did you beat one of the sparse matrix packages? Or were you comparing with a library for dense matrices?
  - yonkshi 8y ago
    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)
twic 8y ago
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.
ryanmonroe 8y ago
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-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

yonkshi 8y ago

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.

repsilat 8y ago
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."
greglindahl 8y ago
Did you beat one of the sparse matrix packages? Or were you comparing with a library for dense matrices?
- yonkshi 8y ago
  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)

twic 8y ago

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.

ryanmonroe 8y ago

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-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