Sparse solve vs BandedMatrix; time and allocation surprise

General Usage Performance
1

I’m trying to understand why the SparseSuite \ is so much worse than what I get with BandedMatrices.
My test case, which reflects the results of my application, uses two matrices with the same elements, but stored differently.

Running v 1.5 on a Mac.

n=20000;
D=rand(n,);
D1=rand(n-1,);
Dm1=rand(n-1,);
D2=rand(n-2,);
Dm2=rand(n-2,);
d0=Pair(0,D); d1=Pair(1,D1); d2=Pair(2,D2);
dm1=Pair(-1,Dm1); dm2=Pair(-2,Dm2);
FVS=spdiagm(dm2,dm1,d0,d1,d2);
FVB=BandedMatrix(dm2,dm1,d0,d1,d2);
b=rand(n,);

So I run the solvers FVS\b and FVB\b and I get this

julia> @btime $FVB\$b;
2.226 ms (11 allocations: 1.37 MiB)
julia> @btime $FVS\$b;
22.556 ms (73 allocations: 25.04 MiB)

I figured the generic spase solve would be worse, but not that much worse. Is there anything I shoud be doing that I’ve missed? Is there a direct call to a function inside SparseSuite that could help?

2

SparseMatrices are inherently slower because they don’t really work well with things like simd, and they also cause branch prediction issues. Both of those are really big deals.

1 Like
3

Nope, you didn’t miss anything. It’s that much of a difference.

4

This is a surprise to me. The C code for SuiteSparse was written by the best in the field

http://faculty.cse.tamu.edu/davis/suitesparse.html

and seems to work very will inside Matlab, but my Matlab experience with it was on pretty small problems and I was only comparing it to the previous sparse solvers in Matlab.

5

It doesn’t matter how good you are if you have no information. If you can specialize on properties of the system, you can do better. SuiteSparse is as good as a completely generic sparse algorithm can get, but completely sparse generic is too generic. Even good sparse algorithms work by trying to add structure.

4 Likes
6

Fair enough. SuiteSparse has always done what I’ve asked it to do and been fast enough. This is the first time I’ve done a comparison like this however.

7

It sounds like you think that SuiteSparse is disappointingly slow? Isn’t it rather the case that BandedMatrices is impressively fast?

(Is the glass 1% empty, or 99% full?)

5 Likes
8

I’ve never compared the LAPACK band solver vs SuiteSparse on the same problem. In hindsight, I should have expected this. Not disappointed in SparseSuite and I knew band solvers are fast. I did not know that the difference would be so large on what seems to be an easy problem for both.

Learn something every day…

1 Like
9

I tried this in Matlab (same computer and OS). Matlab does some analysis before a solve (looking for spd, bands, … ) and is smart enough to use the band solver and hence is as fast as Julia. The timings (assuming you believe tic and toc) are almost identical.

>> n=20000; em2=rand(n,1); em1=rand(n,1); em0=rand(n,1);
>> e1=rand(n,1); e2=rand(n,1);
>> A=spdiags([em1 em2 em0 e1 e2], -2:2, n,n);
>> b=rand(n,1);
>> tic; c=A\b; toc
Elapsed time is 0.002861 seconds.
>> tic; [l,u]=lu(A); toc
Elapsed time is 0.005762 seconds.

I do not know how to force Matlab to use SuiteSparse. Mathworks hired Tim Davis as a consultant to get it in there, so it would be interesting to see how/if they to better than Julia.

1 Like
10

Open an issue. The Julia one could probably detect this case and then run a fast path.

11

Where should I open that issue? Base? SuiteSparse?

There’s some subtlety happening. If you want to do qr! or lu!, you need to allocate a few extra bands. If all you’re doing is backslash, then that extra room needs some management.

12

Julia Base.

13

If the matrix is really simple and has a narrow band, of course the banded solver will win. Is there surprise there?

I have tested previously SuiteSparse with both Matlab and Julia on the same matrices. They were of equivalent speed.

14

Winner = no surprise
Margin of victory = some surprise

I’m moving a lot of my work from Matlab → Julia and finding new things all the time. Most of what I see is well explained if I RTFM. This one was bit different.

2 Likes
15

I wonder how much it might cost to detect the “bandedness” and reorganize the matrix. In other words, is the detection and reorganization done after a solve is called, or before. If the latter, this cost is not part of the timing.

16

Not completely sure. I call tic before anything happens and toc afterwards, so if tic/toc works the way I think it does, it times everything.

17

What I mean is: do you suppose the storage is selected when the matrix is created, or is it changed appropriately after a solve is called?

18

The bandedness cost would be O(n) (check each col’s greatest and least row element). The conversion is O(nd+nnz) since you need to allocate a zero banded matrix and fill in the non-zero entries.

19

Good news: to create a banded matrix from a sparse array can be reasonably amortized.
With Julia 1.6, BandedMatrices v0.15.17 :

@btime $FVB=BandedMatrix($FVS)
@btime $FVB\$b;
@btime $FVS\$b;
  1.739 ms (3 allocations: 781.38 KiB)                                            
  4.180 ms (11 allocations: 1.37 MiB)                                             
  48.322 ms (73 allocations: 25.04 MiB)
2 Likes
20

While the Julia BandedMatrices.jl package was written by the worst in the field :rofl:

EDIT: \ with BLAS types lowers to LAPACK’s banded LU, whose developers are also pretty good.

5 Likes