Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Check bands for zeros in lu! and qr! #197

Open
ctkelley opened this issue Aug 18, 2020 · 27 comments
Open

Check bands for zeros in lu! and qr! #197

ctkelley opened this issue Aug 18, 2020 · 27 comments

Comments

@ctkelley
Copy link

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
Copy link
Author

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
Copy link
Member

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
Copy link
Author

ctkelley commented Aug 18, 2020 via email

@dlfivefifty
Copy link
Member

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

@ctkelley
Copy link
Author

ctkelley commented Aug 18, 2020 via email

@dlfivefifty
Copy link
Member

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
Copy link
Author

ctkelley commented Aug 18, 2020 via email

@YingboMa YingboMa mentioned this issue Aug 23, 2020
Closed
@ctkelley
Copy link
Author

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
Copy link
Author

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 dlfivefifty changed the title qr! confused on BandedMatrices Check bands for zeros in lu! and qr! Aug 24, 2020
@dlfivefifty
Copy link
Member

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
Copy link
Member

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
Copy link
Author

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
Copy link
Member

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

@ctkelley
Copy link
Author

ctkelley commented Aug 25, 2020 via email

@dlfivefifty
Copy link
Member

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
Copy link
Member

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

@ctkelley
Copy link
Author

ctkelley commented Aug 25, 2020 via email

@dlfivefifty
Copy link
Member

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
Copy link
Author

ctkelley commented Aug 25, 2020 via email

@ChrisRackauckas
Copy link

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
Copy link
Member

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
Copy link
Member

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
Copy link
Member

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
Copy link
Member

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

@ctkelley
Copy link
Author

ctkelley commented Aug 25, 2020 via email

@ctkelley
Copy link
Author

ctkelley commented Aug 25, 2020 via email

@dlfivefifty
Copy link
Member

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:

JuliaArrays/OffsetArrays.jl#87

@dlfivefifty dlfivefifty mentioned this issue Nov 26, 2024
Closed
@lxvm lxvm mentioned this issue Aug 5, 2024
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@dlfivefifty @ChrisRackauckas @MikaelSlevinsky @ctkelley