Way to invert matrix?

[djsegal]

Is there a way to invert BandedBlockBandedMatrix? I'm getting the following error from:

A \ b

---

MethodError: no method matching lu!(::BlockBandedMatrices.BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}}, ::Val{true}; check=true)
Closest candidates are:
  lu!(!Matched::Union{DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, ReinterpretArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, ::Union{Val{false}, Val{true}}; check) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lu.jl:37
  lu!(!Matched::Union{Hermitian{T,S}, Symmetric{T,S}} where S where T, ::Union{Val{false}, Val{true}}; check) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lu.jl:45
  lu!(!Matched::Union{DenseArray{T,2}, ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T, ::Union{Val{false}, Val{true}}; check) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lu.jl:86
  ...

Stacktrace:
 [1] #lu#103(::Bool, ::Function, ::BlockBandedMatrices.BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}}, ::Val{true}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lu.jl:142
 [2] lu(::BlockBandedMatrices.BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}}, ::Val{true}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lu.jl:142 (repeats 2 times)
 [3] \(::BlockBandedMatrices.BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2}}}, ::Array{Float64,1}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/generic.jl:870
 [4] solve!(::LinearSystem) at /Users/dan/code/NINT/src/methods/solve.jl:2
 [5] solve!(::LaplaceProblem) at /Users/dan/code/NINT/src/methods/solve.jl:6
 [6] #LaplaceProblem#8(::Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:r_min, :r_max, :reference_value),Tuple{Float64,Float64,Int64}}}, ::Type, ::Int64, ::Int64, ::Int64) at /Users/dan/code/NINT/src/structs/laplace_problem.jl:58
 [7] Type at ./none:0 [inlined]
 [8] (::getfield(Main, Symbol("##971#974")))(::Int64, ::Int64, ::String) at ./In[255]:2
 [9] map(::Function, ::Observable{Int64}, ::Observable{Int64}, ::Observable{Any}) at /Users/dan/.julia/packages/Observables/qCJWB/src/Observables.jl:174
 [10] top-level scope at /Users/dan/.julia/packages/InteractBase/3SqBl/src/manipulate.jl:21
 [11] top-level scope at In[255]:1

[dlfivefifty]

This isn't implemented yet. Note that the LU factorisation of a BandedBlockBandedMatrix may have dense blocks, so when implemented it would consist of first converting to a BlockSkylineMatrix and using the block banded LU factorisation. This is possible to implement using LAPack's dense LU on the column blocks, but we haven't got around to it.

[djsegal]

Is there anything I could do to help?

[dlfivefifty]

If you want to try implementing it that would be great. The goal is to have lu!(::BlockSkylineMatrix) return a BlockSkylineMatrix. There are a few steps:

1. Copy over lu!(::StridedMatrix) to LazyArrays.jl to work with general strided matrices. That is, something like

_lu!(_, A, pivot, check) = LinearAlgebra.lu!(A, pivot; check=check) # Fall back to standard `lu!`
function _lu!(::AbstractColumnMajor, A::AbstractMatrix{<:BlasFloat}, pivot, check)
    if pivot === Val(false)
        return generic_lufact!(A, pivot; check = check)
    end
    lpt = LAPACK.getrf!(A)
    check && checknonsingular(lpt[3])
    return LU{T,typeof(A)}(lpt[1], lpt[2], lpt[3])
end
lu!(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true);
             check::Bool = true) = _lu!(MemoryLayout(A), A, pivot, check)

2. If A is a BlockSkylineMatrix, then L = colblockbandwidth(A,1)[K]; V = view(A, Block.(K:K+L) is already a ColumnMajor, so LazyArrays.lu! should work. We therefore just need to go block-column-by-block-column, call lu!, (which modifies A in place), then use this to update the columns to the right of A.
3. The last part is lu(A), which needs to pad A by adding extra blocks above the diagonal to capture the fill-in.

[djsegal]

Would the solver described in the following discourse post be useful for any of this?

• https://discourse.julialang.org/t/way-to-speed-up-a-2x2-block-matrix-solver/20628

[dlfivefifty]

I believe doing a hierarchical solver with block matrix inversion can be used to reduce the complexity from O(n^2) to O(n^(3/2)), where n is the total number of unknowns. The right way to do this would be to leverage https://github.com/JuliaMatrices/HierarchicalMatrices.jl.

That said, the first step is a proper standard LU decomposition, which while it is O(n^2) can easily leverage LAPack, and hence get multithreading for free, and for moderate n would surely be much faster than hierarchical solvers. Performance testing would then guide how many levels to use in a hierarchical solver.

[djsegal]

Putting the hierarchical matrices aside, how exactly is the LU decomposition not leveraging LAPack?

Is this like a 2 line fix or something much deeper?

// also, this segue should have probably gone on BandedMatrices... sorry about that 😬

[dlfivefifty]

Not two lines, I outlined what needs to be done above.

[djsegal]

Sounds good. Just wanted to know if there was a difference now that we're focused on BandedMatrices and not BandedBlockBandedMatrices.

…

On Feb 13, 2019, at 5:57 AM, Sheehan Olver ***@***.***> wrote: Not two lines, I outlined what needs to be done above. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

[dlfivefifty]

Oh, BandedMatrices.jl already has a fast LU using LAPack, so I'm not sure I understand the point of block inversion except for BlockBandedMatrices.jl

[djsegal]

This was two separate issues.

At first, I wanted a black-box solution to quickly solving A\b using BandedBlockBanded matrices.

After realizing this wasn't possible, I setup my matrix so that it has 4 blocks:

submatrix	type
A	Banded
B	Sparse
C	Sparse
D	Dense

Then using blockwise inversion, I setup a package that should in theory quickly do this. However I'm not getting the speedup I thought I'd get.

That's why I posted to discourse and thought it might be related to this issue. Apologies again for the confusion.

[dlfivefifty]

This is the sort of thing SparseMatrixCSC excels at, as it automatically deduces the optimal pivoting strategy. Is there a reason sparse matrices won't work for you here?

[djsegal]

It seemed like you could get a win in the limit of A becoming the whole matrix? Shouldn't a diagonal solver beat a sparse array solver (assuming narrow bands)?

[dlfivefifty]

For a well-conditioned banded matrix, UMFPack is almost as fast as LAPack, especially with OpenBLAS which has very slow banded solvers. My blog post https://approximatelyfunctioning.blogspot.com/2018/12/banded-matrices-and-ordinary.html gives some timing examples but these are with MKL. But even if you get a 4x speedup using a banded solver, that gain will probably be lost in a blockwise inversion.

[rveltz]

Can't it be done for BlockTridiagonalMatrix first? There should be an analytical formula

[rveltz]

For example, for BlockTridiagonalMatrix, the LU decomposition is known analytically. The issue is that the block of the decomposition could be dense whereas, in principle, we only need to store their action. I was thinking about using LazyArrays for this. Would you recommand something else?

[dlfivefifty]

Are you imagining an increasingly complicated operator tree? I don't think that's viable, and if it is LazyArrays.jl won't help since it templates the arguments, meaning having complicated operator trees requires excessive compiling.

For BlockTridiagonalMatrix you wouldn't want this anyways since the blocks are dense, so its possible to do the linear algebra directly.

[rveltz]

My blocks are 1000 x 1000 so I am not sure I want to fill them in with a LU decomposition. Also the matrice A is from a 5 points finite differences approximations, a bit like your Laplacian example but not homogenous. I was thinking about gauss seidel, but it is tranpose(A) which is diagonally dominant, so the convergence is slow.

[rveltz]

I forgot to say that you were right, this is what I intended to do until I read your answer...

[dlfivefifty]

Do you need to calculate the factorisation fast, or just invert systems fast? That is, are you using the same matrix with multiple right hand side?

The LU factors should be practically banded-block-banded (they’ll have exponential decay in the bands) so one option is to calculate the factorisation using dense LU and then compress the results.

But really this is a case where iterative methods are the go to solution.

[rveltz]

I need to invert systems fast. More precisely, I need to solve (I-A) \ rhs1 and (I-A/2) \ rhs2. Right now, I am doing this with SparseArrays: one third of the time is spent to build A and 2/3 of the time to solve the two equations.

I tried iterative methods but I dont have a good preconditioner. The convergence is way too slow.

[dlfivefifty]

Right. What you are asking for may not currently be possible. I don't think hierarchical solvers will help much here either as the involve dense inverses, that is, for an n^2 x n^2 matrix we need to invert log(n) matrices of size n x n. After this inverse is calculated it can be compressed. So usually they are only used when the same inverse is reused multiple times.

I have some vague ideas on using the fact that the inverse of banded matrices is banded + semi-separable, hence can be represented using a small amount of memory and I believe form an algebra, but I don't have time to pursue this at the moment (I'm hoping for an interested student at some point).