New blog post on banded matrices and ordinary differential equations using BandedMatrices.jl



I just posted my first blog post, on using BandedMatrices.jl to solve ordinary differential equations via finite differences / spectral methods:

This is first of a series outlining the results of the NumFocus small development grant to add support for GPUs and other backends to BlockBandedMatrices.jl


I accidentally posted an old working draft. I’ve uploaded the correct version which has some minor changes.


This is a very nice article and introduction to BandedMatrices.jl, thanks! It is disheartening to hear about the OpenBlas issues with banded solvers since installing MKL with Julia still requires some extra work…

A few very minor comments/typos in case you are updating it at some point:

  1. It would be nice if you could make the width of your text and code boxes larger. Right now many of the code boxes require horizontal scrolling to read all the text. (At least on Chrome on my Mac.)
  2. Right at the start of section 2 you say j=1,2,\dots,n. I assume you meant, k=1,2,\dots,n.
  3. When plotting the coefficients your code says: .+ 1E-40 in a comment but actually uses 1E-100.
  4. Where you btime the lu(L) and lu(L̃) I assume the comment for the latter meant pivoted LU for sparse matrices?


Me too on the boxes and scrolling.


It would be nice if you could make the width of your text and code boxes larger. Right now many of the code boxes require horizontal scrolling to read all the text. (At least on Chrome on my Mac.)

Yes, the formatting still could use improvement. Part of the issue is the Weave.jl generates raw HTML, and only Blogger supports this from what I can tell. Blogger is a bit daunting and I haven’t figured out the “right” way to adjust the default layout. In particular, I don’t know why the column of text is so narrow.

Your corrections are correct, thanks.


I’ve increased the widths, does this improve things for you?


Very interesting! Do you have an idea of the trade-offs (in terms of bandwith size perhaps) between this approach and that of using iterative solvers with FFTs to do the matvecs?


New widths look much better to me!


Under practically all circumstances iterative methods are not competitive for ODEs, and this is certainly the case for singularly perturbed examples where one would need a very good preconditioner to overcome the ill-conditioning. I don’t know of such a preconditioner.

PDEs are a different story as the complexity of applying an operator is better than inverting an operator. I’ll talk about PDEs in the next post.


Great work. I like the package and your notes. I never thought about obtaining the fourth derivative matrix as a multiplication of two central derivative matrices — very nice!. Looking forward to the PDE notes.


I’m not sure. On my screen, there’s still more than twice as much horizontal empty space as there is scroll box. I guess if your tools don’t make a “responsive” page, there’s not much you can do… I wouldn’t want the text to fill all the space anyway.


Nice post and blog !
I wonder about the opportunity to add fixed size bandwidth (tri/penta/hepta/N-diag) specializations. We implement such N-diag matrices julia types for one of our’s PDE solver and I should compare the performance with this package. We usually use in-place LLT or LDT factorization to solve corresponding symmetric systems.


This sounds like StaticArray v regular array. I don’t think there’s much benefit though in knowing the bandwidths at compile time, as one is usually dealing with large dimensional problems where the actual computation swamps things.


Yes, In our cases which maybe not generally relevant, the inner loop operates through the bandwidth. When it is known statically, the loop is unrolled and this can lead to significant speed-up. This is the case for a cholesky factorization for instance. We define our own macros for this purpose but I guess that it maybe implemented efficiently with StaticArrays.


I forgot to mention that we deal with larges systems but small bandwidths (3-9).


If you are sure knowing that at compile time is useful then it’s a good argument. But it would need a new SBandedMatix type that knows the bandwidths in the template variables.


That is exactly what we have implemented.


Please make a PR if you think it’s of general use



I have tried your package on a 2d problem and it is quite slower than using sparse. My matrix line numbers is Nx * Ny, same for columns. It has 5 bands at 0, ±1, ±Nx. It seems that you store a big matrix to hold this banded structure and that may be the reason of the slowdown I am experiencing. Do you have any further advice please?

Thank you


I wouldn’t expect a banded matrix structure to do well there since you want something like “dense bandedness”. But @dlfivefifty have you considered this kind of “sparse bandedness”?