Check bands for zeros in `lu!` and `qr!`

[ctkelley]

Hello again, it seems that qr! thinks it has a method for BandedMatrices, but the results are incorrect. lu! knows it has no such method and issues the usual complaint. Example below.

I am trying to figure out how to do a PR for the Float32 issue. This may take me awhile because I'm still not sure where the problem lives in your code and have never issued a PR for a project that I'm not seriously involved with.

-- Tim

julia> B=brand(Float64,10,10,2,2); julia> b=rand(10,); julia> c=B\b; julia> P=qr(B); julia> d=P\b; norm(d-c) 1.00688e-15 julia> C=copy(B); Q=qr!(C); e=Q\b; norm(e-c) 3.58602e+00

[ctkelley]

Sorry, lu! does think it has a method, but that is wrong too.

julia> Z=lu!(B) julia> f=Z\b; julia> norm(f-c) 8.02474e+01

[dlfivefifty]

qr! assumes that you've pre-allocated the matrices with extra bands, as otherwise in-place qr is not possible. Would you have a better design in mind? We could check that the extra bands are zero.

[ctkelley]

Thanks,

That makes sense. lu! for dense matrices allocates the extra storage for the pivots but I can see that as a real pain for the banded case. Was that documented somewhere? I missed it.

This is also a very good case for using qr instead of lu.

You need to double the bandwidths if I recall correctly. Is that right?

— Tim

On Tue, Aug 18, 2020 at 10:24 AM Sheehan Olver notifications@github.com wrote:

qr! assumes that you've pre-allocated the matrices with extra bands, as otherwise in-place qr is not possible. Would you have a better design in mind? We could check that the extra bands are zero.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-675509505, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX62UNJNCM2NE4I44GXTSBKFIPANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[dlfivefifty]

If A has bandwidths (l,u) then qr(A).factors has bandwidth (l,u+l)

[ctkelley]

So lu! needs 2l and 2u?

On Tue, Aug 18, 2020 at 10:50 AM Sheehan Olver notifications@github.com wrote:

If A has bandwidths (l,u) then qr(A).factors has bandwidth (l,u+l)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-675525896, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX67QBCSS2VYTXZFTOXLSBKIM5ANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[dlfivefifty]

I didn't develop the lu! code so I'm not sure I can answer. The LAPACK documentation should provide details: http://www.netlib.org/lapack/explore-html/d2/d2d/group__double_g_bcomputational_gad1efab86e6d869915e059286ecf1bcb1.html

[ctkelley]

No problem. The LAPACK manual is less than a meeter away from me.

On Tue, Aug 18, 2020 at 10:58 AM Sheehan Olver notifications@github.com wrote:

I didn't develop the lu! code so I'm not sure I can answer. The LAPACK documentation should provide details: http://www.netlib.org/lapack/explore-html/d2/d2d/group__double_g_bcomputational_gad1efab86e6d869915e059286ecf1bcb1.html

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-675530443, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX6462A2MH5VFYFAT2JTSBKJI3ANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[ctkelley]

You need to embed your matrix into one with bandwidths ll, and lu+2 with zeros on the unused bands. Then it works fine. So does lu!

You should look at the LAPACK manual, as @dlfivefifty suggests.

[ctkelley]

Take this example, please.

julia> using BandedMatrices julia> B=BandedMatrix{Float64}(Zeros(20,20),(2,6)); julia> D=rand(20,); D1=rand(19,); D2=rand(18,); D3=rand(17,); D4=rand(16,); julia> Dm1=rand(19,); Dm2=rand(18,); julia> B[band(0)].=D; B[band(1)].=D1; B[band(2)].=D2; B[band(3)].=D3; julia> B[band(4)].=D4; B[band(-1)].=Dm1; B[band(-2)].=Dm2; julia> b=rand(20,); C=copy(B); x1=B\b; julia> QRB=qr!(C); julia> x2=QRB\b; julia> norm(x1-x2) 8.86313e-14

[dlfivefifty]

A simple for loop checking that the extra bands are zero, and throwing an error if not, should fix this. The cost is negligible to the actual qr!, though we can bypass the check in qr(A) where the bands are initialised correctly.

[MikaelSlevinsky]

There's an argument here for not naming the current code qr!(A) since it doesn't work for all banded matrices, just those with the extra data assumption. Certain assumptions on A already result in an InexactError according to the base docs.

New in 1.5 is lu! for sparse matrices. According to the docs, due to type conversions, the signature takes the form

lu!(F::UmfpackLU, A::SparseMatrixCSC) uses the UMFPACK library that is part of SuiteSparse.
As this library only supports sparse matrices with Float64 or ComplexF64 elements,
lu! converts A into a copy that is of type SparseMatrixCSC{Float64} or SparseMatrixCSC{ComplexF64} as appropriate.

Actually calling this according to the example only changes F not A even though A is the "input" data.

One option would be to use qr!(F::BT, A::BandedMatrix), where BT is either a banded QR factorization type or another banded matrix whose superdiagonals are extended.

[ctkelley]

No, please no. qr! does exactly what it does in LAPACK and LINPACK. You have to pad the matrix with a couple zero bands, but then it works correctly and is the same as the calls in C or Fortran. lu! does the same thing on BandedMatrices.

The general sparse case is very different. I am not clear on what the right thing to do is there.

[dlfivefifty]

Wait, there's a banded qr! in LAPACK??

[ctkelley]

https://www.netlib.org/lapack/improvement.html item 2.4

On Tue, Aug 25, 2020 at 3:05 AM Sheehan Olver notifications@github.com wrote:

Wait, there's a banded qr! in LAPACK??

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-679844916, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX6ZASTZLJYBXLQZ3QT3SCNPDHANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[dlfivefifty]

Must have been using old docs 🤦‍♂️

We aren't actually calling that function, instead using a Julia native banded QR. This should be changed

[dlfivefifty]

Is there also a banded * banded function I've missed in BLAS or LAPACK???

[ctkelley]

I don't think so. Not in BLAS3 anyhow. There's support for Banded * vector in BLAS2.

On Tue, Aug 25, 2020 at 5:45 AM Sheehan Olver notifications@github.com wrote:

Is there also a banded * banded function I've missed in BLAS or LAPACK???

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-679922738, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX6Y773AZBCJLIAXXRMLSCOB5HANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[dlfivefifty]

Yes we use BLAS.gbmv! https://github.com/JuliaMatrices/BandedMatrices.jl/blob/12fbc03d3798f4b1ed68f55e92c6151aa1e30b65/src/generic/matmul.jl#L23

Though it's probably not much faster than the Julia native code (BLAS2 rarely seems to make a difference).

Unfortunately we are using a hacked together "gbmm!", which I think I made in Julia v0.3 and never properly updated. At some point I was keen on the idea that the sub-components of a banded matrix with u = l have BLAS storage: that is, it can be thought of as a block-tridiagonal matrix whose diagonal blocks are dense column major arrays and whose off-diagonal blocks are triangular. Then I realised there does not appear to be a BLAS3 triangular * triangular routine either. (Probably not too hard to make a Julia one that chops up a triangular matrix in a hierarchical way... that is, as a block-triangular matrix with triangular diagonal blocks and dense off-diagonal block...)

[ctkelley]

I don't know if you have seen this thread on discourse or not

https://discourse.julialang.org/t/sparse-solve-vs-bandedmatrix-time-and-allocation-surprise/45119

@ChrisRackauckas suggests that I open an issue to see if spdiagm can produce a banded matrix if the bandwidths are small enough for that to make sense. I have been composing that message and waiting for the thread to play out. Would it make sense for BandedMatrices to get merged into base as well?

— Tim

On Tue, Aug 25, 2020 at 9:28 AM Sheehan Olver notifications@github.com wrote:

Yes we use BLAS.gbmv! https://github.com/JuliaMatrices/BandedMatrices.jl/blob/12fbc03d3798f4b1ed68f55e92c6151aa1e30b65/src/generic/matmul.jl#L23

Though it's probably not much faster than the Julia native code (BLAS2 rarely seems to make a difference).

Unfortunately we are using a hacked together "gbmm!", which I think I made in Julia v0.3 and never properly updated. At some point I was keen on the idea that the sub-components of a banded matrix with u = l have BLAS storage: that is, it can be thought of as a block-tridiagonal matrix whose diagonal blocks are dense column major arrays and whose off-diagonal blocks are triangular. Then I realised there does not appear to be a BLAS3 triangular * triangular routine either. (Probably not too hard to make a Julia one that chops up a triangular matrix in a hierarchical way... that is, as a block-triangular matrix with triangular diagonal blocks and dense off-diagonal block...)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-680024970, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX64KPN6S7VBN7SAUVOTSCO4AFANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[ChrisRackauckas]

I think that would go into Julialang/julia since IIRC you're looking for lu! and qr! on dense matrices to detect bandedness. The current code for the default factorize used by \ is https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/dense.jl#L1238-L1315 which as you can see specializes a lot but specifically doesn't include bandedness as something it can do, because BandedMatrices is a separate library from the stdlib. Usually there's no reason to add a package to the stdlib, but BandedMatrices.jl is one of the few things I think should get promoted so that the default linear algebra tools can exploit bandedness in algorithms and outputs more. @dlfivefifty what do you think of that?

[dlfivefifty]

I think adding it completely to StdLib is a nonstandard for the time being: I need ArrayLayouts.jl to support views of banded matrices behaving like banded matrices, so one would need to first talk about incorporating (a replacement for) ArrayLayouts.jl into StdLib to do trait-based dispatch. (This almost made it into Julia v1.0 but we couldn't quite get agreement on what traits should be supported.)

That said, what can easily be done is incorporate the notion of bandedness into StdLib. For example, ArrayLayouts.jl has no notion of a banded matrix but already supports optimal complexity algorithms via iterating only over colsupport and rowsupport (the non-zero entries of the corresponding column and row...could be renamed colrange and rowrange since in practice the "support" is always a range). E.g., the following triangular lmul! works as-is for triangular banded matrices via (i + 1:m) ∩ rowsupport(Adata,i):

https://github.com/JuliaMatrices/ArrayLayouts.jl/blob/27736556b9dbdc1a8725b1fab1dc208c96b873a2/src/triangular.jl#L25

Then BandedMatrices.jl could play the corresponding roll as the "canonical banded matrix package", like OffsetArrays.jl does for non-1-based arrays.

[MikaelSlevinsky]

https://www.netlib.org/lapack/improvement.html item 2.4

Looking at the associated working note (http://www.netlib.org/lapack/lawnspdf/lawn203.pdf), the new-to-that-version computational routines, e.g. xLARFB.f, scan Householder reflectors to be applied to a dense n x n matrix for trailing zeros, saving some cost for low profile matrices. Scanning for zeros costs O(n^2).

This package, however, uses the banded storage scheme which is not the same as storing an n x n banded matrix densely.

As annoying as it may seem, you can't perform an in-place QR on a banded matrix without the extra data assumption. Calling that code qr!(A::BandedMatrix) is only going to bring more folks to the already filed issues.

A simple for loop checking that the extra bands are zero, and throwing an error if not, should fix this. The cost is negligible to the actual qr!, though we can bypass the check in qr(A) where the bands are initialised correctly.

Checking for zeros in pre-allocated upper bands is not the only plausible thing to do: if the data were allocated as junk (undef), or even if they're normwise small.

[dlfivefifty]

It's impossible to tell if it's junk, and normwise small is probably too clever; we can just add a flag to ignore bands that users can use for that.

[MikaelSlevinsky]

Oops, yeah, in isbits type floating-point it's impossible to tell if undef is junk

[ctkelley]

I don't think it's a burden on the user (aka me) to be aware of the need for the extra zero bands and populate them correctly. An automatic check, if it does not cost much, would be fine.

On Tue, Aug 25, 2020 at 10:20 AM Mikael Slevinsky notifications@github.com wrote:

Oops, yeah, in isbits type floating-point it's impossible to tell if undef is junk

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-680055120, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX67VHE4GHPGHXDXSN23SCPCEDANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[ctkelley]

I'm a bit confused now. Can/should something be added to base that supports banded matrices and solvers like lu, qr, lu!, and qr!?

— Tim

On Tue, Aug 25, 2020 at 10:06 AM Sheehan Olver notifications@github.com wrote:

I think adding it completely to StdLib is a nonstandard for the time being: I need ArrayLayouts.jl to support views of banded matrices behaving like banded matrices, so one would need to first talk about incorporating (a replacement for) ArrayLayouts.jl into StdLib to do trait-based dispatch. (This almost made it into Julia v1.0 but we couldn't quite get agreement on what traits should be supported.)

That said, what can easily be done is incorporate the notion of bandedness into StdLib. For example, ArrayLayouts.jl has no notion of a banded matrix but already supports optimal complexity algorithms via iterating only over colsupport and rowsupport (the non-zero entries of the corresponding column and row...could be renamed colrange and rowrange since in practice the "support" is always a range). E.g., the following triangular lmul! works as-is for triangular banded matrices via (i + 1:m) ∩ rowsupport(Adata,i):

https://github.com/JuliaMatrices/ArrayLayouts.jl/blob/27736556b9dbdc1a8725b1fab1dc208c96b873a2/src/triangular.jl#L25

Then BandedMatrices.jl could play the corresponding roll as the "canonical banded matrix package", like OffsetArrays.jl does for non-1-based arrays.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/JuliaMatrices/BandedMatrices.jl/issues/197#issuecomment-680046851, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOEX67U6RLXWFNIAQK5FV3SCPAPBANCNFSM4QDN27QQ .

-- C. T. Kelley Department of Mathematics, Box 8205 SAS Hall 2311 Stinson Drive North Carolina State University Raleigh, NC 27695-8205 (919) 515-7163, (919) 513-7336 (FAX) Tim_Kelley@ncsu.edu https://ctk.math.ncsu.edu

[dlfivefifty]

First a point of terminology: note it should be "added to StdLib": Base doesn't know anything about linear algebra, but the StdLib package LinearAlgebra.jl does.

Support could be added without adding all of BandedMatrices.jl: there's a notion of "sanctioned type piracy" so Base could in theory detect a matrix is banded and call a bandedlu that defaults to standard lu, unless a user has first done using BandedMatrices which will commit type piracy and overload bandedlu to use its functionality. As a precedent, see OffsetArrays.jl:

https://github.com/JuliaArrays/OffsetArrays.jl/issues/87