# Sparse solve vs BandedMatrix; time and allocation surprise

**URL:** <https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119>\
**Category:** Performance\
**Created:** [August 17, 2020, 4:48pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119 "2020-08-17T16:48:24Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 17, 2020, 4:48pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/1 "2020-08-17T16:48:24Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 17, 2020, 5:38pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/2 "2020-08-17T17:38:37Z")

</div>

`SparseMatrice`s 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 17, 2020, 5:40pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/3 "2020-08-17T17:40:54Z")

</div>

> [@ctkelley](#):
>
> 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?

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

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 17, 2020, 7:36pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/4 "2020-08-17T19:36:04Z")

</div>

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](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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 17, 2020, 7:39pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/5 "2020-08-17T19:39:38Z")

</div>

> [@ctkelley](#):
>
> The C code for SuiteSparse was written by the best in the field

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.

> **[Algorithm efficiency comes from problem information - Stochastic Lifestyle](https://www.stochasticlifestyle.com/algorithm-efficiency-comes-problem-information/)**
>
> This is a high level post about algorithms (especially mathematical, scientific, and data analysis algorithms) which I hope can help people who are not researchers or numerical software developers better understand how to choose and evaluate...

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 17, 2020, 7:55pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/6 "2020-08-17T19:55:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 17, 2020, 8:09pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/7 "2020-08-17T20:09:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 17, 2020, 8:18pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/8 "2020-08-17T20:18:27Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 18, 2020, 6:38pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/9 "2020-08-18T18:38:32Z")

</div>

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.

```julia
>> 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 18, 2020, 9:19pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/10 "2020-08-18T21:19:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 18, 2020, 9:32pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/11 "2020-08-18T21:32:02Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 18, 2020, 9:55pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/12 "2020-08-18T21:55:24Z")

</div>

> [@ctkelley](#):
>
> Where should I open that issue? Base? SuiteSparse?

Julia Base.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [August 19, 2020, 12:22am UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/13 "2020-08-19T00:22:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 19, 2020, 11:01am UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/14 "2020-08-19T11:01:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [August 19, 2020, 2:42pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/15 "2020-08-19T14:42:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [August 19, 2020, 2:51pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/16 "2020-08-19T14:51:42Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [August 19, 2020, 3:03pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/17 "2020-08-19T15:03:11Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 19, 2020, 3:06pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/18 "2020-08-19T15:06:24Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [August 19, 2020, 3:32pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/19 "2020-08-19T15:32:29Z")

</div>

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

```julia
@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)     

```

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [August 25, 2020, 2:06pm UTC](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119/20 "2020-08-25T14:06:14Z")

</div>

> [@ctkelley](#):
>
> This is a surprise to me. The C code for SuiteSparse was written by the best in the field

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

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

[Next page](https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119.md?page=2)