Sparse solve vs BandedMatrix; time and allocation surprise

General Usage Performance
21

Is this a reasonable way to open the issue?

–
@ChrisRackauckas suggests that I raise this issue. It came up on

`Sparse solve vs BandedMatrix; time and allocation surprise

Is it reasonable for spdiagm to return a banded matrix and use the LAPACK
band solvers if the bandwidth is sufficiently narrow?

Right now you get a general sparse matrix and SuiteSparse does the solve.

Show the example from discourse

Matlab seems to do the right thing

Show them the other example

22

Type parameter for whether bidiagonal matrix is low(:L) bidiagonal or up(:U) bidiagonal. · Issue #33402 · JuliaLang/julia · GitHub is another case where banded matrices could/should be used in Base as the output type of an operation to allow for further optimizations, but currently the output is just a pure sparse matrix. @dlfivefifty mentioned in Check bands for zeros in `lu!` and `qr!` · Issue #197 · JuliaLinearAlgebra/BandedMatrices.jl · GitHub that it can’t be added to the StdLib yet because of ArrayLayouts.jl, but I think there’s more than a few cases where Base should be using BandedMatrix types in some places and so the only way to get the efficient algorithm would be to have some form of a BandedMatrix in Base, or to just copy some of the methods over and have the best possible algorithm in Base.

Banded matrices just show up as such a common form that I think we’re losing out by not having them in our pool of go-to things that Base is allowed to give you and specialize on.

8 Likes
23

So yes, let’s get an issue in Base to discuss this more. @dlfivefifty would you want to lead this to discuss what would need to be done in order to make that possible?

24

:+1:

25

@dlfivefifty, are you taking ownership of this issue? That would be fine with me and you have more gravitas on this topic than I do.

— Tim

26

No. That would be type-unstable.

A factor-of-10 slowdown between a generic sparse-direct solve, which assumes no matrix structure, and a specialized banded-matrix solve seems totally unsurprising to me. If you have a structured matrix, you should use a structured-matrix type to take advantage of it.

2 Likes
27

https://github.com/JuliaLang/julia/issues/37258

2 Likes
28

Nicely done. It’ll be interesting to see where this goes.

29

In a problem, till now, I was using sparse matrices (of size ~= (1_000_000, 1_000_000)) for
matrix-vector \ operation.

Now, I have realized that my matrix has a descent banded structure so I want to gain benefit from this.
this blog and the current thread of discussion, says that I would probably get a performance improvement.

In the flow of a program, I have a SparseMatrixCSC type matrix which needs to be converted to BandedMatrix type before using with \ operator.
I wrote the following piece of code for converting SparseMatrixCSC to BandedMatrix.

function banded_matrix(A::SparseMatrixCSC)
    I, J, V = findnz(A)
    vals = Dict{Int, Vector{Float64}}()
    for (ai, aj, av) in zip(I, J, V)
        d_ij = ai-aj
        if d_ij in keys(vals)
            push!(vals[d_ij], av)
        else
            vals[d_ij] = [av,]
        end
    end
    A = BandedMatrix(
        Tuple((k => v for (k, v) in vals)),
        size(A),
    )
    return A
end

Running this is giving StackOverflowError.

@dlfivefifty, Please let me know if there is any method already available in BandedMatrices.jl for this sparse to banded matrix conversion.

Thanks,

30

No there’s not yet a method. Probably the best approach would be to implement bandwidths(::SparseMatrixCSC), create a zero banded matrix with those bandwidths, then just iterate through the nonzero entries calling banded_setindex!

2 Likes