Add inv of Bidiagonal

[dlfivefifty]

This is the first example of the inverse of a banded matrix being a Semiseparable matrix. If we add multiplications of SemiseparableMatrix we can pretty easy get inverse of tridiagonal matrices too via the UL decomposition.

[andreasnoack]

I very excited about this. I thought this approach was difficult to get working in practice because the vectors ended up under- and overflowing. I'm working on a landscape ecology application where we can't avoid to compute an inverse of a block-tridiagonal matrix representing movements over a 2D grid. It should also have a semiseparable representation so I'm looking forward to the developments here since it might allow us to compute results for much larger landscapes.

[dlfivefifty]

@MikaelSlevinsky Do you know anything about avoiding under-/overflow in this context?

@andreasnoack I can see your point: if the first column becomes small (which is typical) and underflows we aren't going to be able to define it as is. An easy way to avoid this is by making it a more specialised type instead of just storing the low-rank vectors. In fact, I wonder if, with the right storage, multiplication of semi-separable matrices (as needed in doing U,L = ul(A); inv(L)*inv(U) for inverting a tridiagonal matrix) becomes cleaner...

I mostly added this just as notes for myself for now so will likely come back later and explore it further. To be honest I don't really know the literature that well; I know that semi-separable matrices are usually handled with a hierarchical discretisation but that is clearly overkill for the inverse of a bidiagonal matrix.

[dlfivefifty]

inverse of a block-tridiagonal matrix representing movements over a 2D grid

Sounds close to something else I'm interested in: using this for block-Cholesky/Gaussian elimination for BandedBlockBandedMatrix. Do you also need to algebra with the inverses: multiply them, add them, etc.?

[MikaelSlevinsky]

Yeah there's definitely possibility for overflow/underflow. This happens in practice for semiseparable discretizations of self-adjoint eigenvalue problems (e.g. spherical harmonic transforms). That's why I prefer the symmetric-definite and banded formulations.

[MikaelSlevinsky]

HSS format can alleviate this to an extent

[andreasnoack]

Thanks for the feedback. The problem I'm solving is (inv(I - A)*B*inv(I - A)) ./ inv(I - A) where A and B represent movements in a 2D grid so they are both block-tridiagonal with diagonal or tridiagonal off-diagonal blocks. A is furthermore substochastic which ensures that I - A is invertible. If inv(I - A) could be computed in a sparse representation then we'd be able to solve for much larger landscapes than we are currently able to handle.

HSS format can alleviate this to an extent

Do you have a Julia implementation of this format that I can try out?

[dlfivefifty]

It’s somewhat neglected but you can use https://github.com/JuliaMatrices/HierarchicalMatrices.jl

[dlfivefifty]

Here's an idea that should overcome under/overflow: store log.(ℓ) instead of just ℓ. I.e. replicate how loggamma is useful for products of gamma functions.

(This is getting rather close to Clenshaw–Olver Level-index Arithmetic...)

[dlfivefifty]

Just pushed some tests that demonstrate the formula for inv(::Tridiagonal) in terms of the UL decomposition: if inv(L) == tril(a*b') and inv(U) == triu(c*d') we get a beautifully simple formula (here Ai is the inverse of a tridiagonal matrix):

            b̃ = cumsum(b .* c) .* d
            c̃ = a .* cumsum(b .* c)
            @test tril(a*b̃') ≈ tril(Ai)
            @test triu(c̃*d') ≈ triu(Ai)

Not sure how this fits in to the log storage / level-index arithmetic proposal for avoiding under/overflow: one would need to compute cumsum(b .* c) .* d in that format and its not immediately clear how that would work.

[andreasnoack]

I tried out a simple version of the log idea for the symmetric tridiagonal case where the input is chosen such that all elements of the inverse are positive and it looks quite promising

struct Semiseparable2{T} <: AbstractMatrix{T}
    u::Vector{T}
    v::Vector{T}
end

Base.size(S::Semiseparable2) = (length(S.u), length(S.v))

function Base.getindex(S::Semiseparable2, i::Integer, j::Integer)
    if j >= i
        return exp(S.u[i] + S.v[j])
    else
        return exp(S.u[j] + S.v[i])
    end
end

function myinv2(S::SymTridiagonal)
    α = S.dv
    β = S.ev
    u = similar(α)
    v = similar(α)
    n = length(α)

    u[1] = 0
    u[2] = log(α[1]) - log(-β[1])
    for i in 2:(n - 1)
        u[i + 1] = u[i] + log(α[i] + β[i - 1]*exp(u[i - 1] - u[i])) - log(-β[i])
    end
    v[n]     = -u[n] - log(α[n] + β[n - 1]*exp(u[n - 1] - u[n]))
    v[n - 1] = v[n] + log(α[n]) - log(-β[n - 1])
    for i in (n - 2):-1:1
        v[i] = v[i + 1] + log(α[i + 1] + β[i + 1]*exp(v[i + 2] - v[i + 1])) - log(-β[i])
    end
    return Semiseparable2(u, v)
end

julia> n = 100;

julia> A = SymTridiagonal(ones(n), fill(-0.25, n-1));

julia> norm(inv(A) - myinv2(A))/norm(inv(A))
1.4938522770351065e-14

julia> n = 10_000_000;

julia> A = SymTridiagonal(ones(n), fill(-0.25, n-1));

julia> @time myinv2(A)
  1.039874 seconds (5 allocations: 152.588 MiB, 6.39% gc time)
10000000×10000000 Semiseparable2{Float64}:
 1.07181      0.28719      0.0769524    0.0206193    …  0.0          0.0          0.0
 0.28719      1.14876      0.30781      0.0824774       0.0          0.0          0.0
 0.0769524    0.30781      1.15429      0.30929         0.0          0.0          0.0
 0.0206193    0.0824774    0.30929      1.15468         0.0          0.0          0.0
 0.00552494   0.0220997    0.082874     0.309396        0.0          0.0          0.0
 0.0014804    0.00592161   0.022206     0.0829025    …  0.0          0.0          0.0
 0.000396673  0.00158669   0.00595009   0.0222137       0.0          0.0          0.0
 0.000106288  0.000425152  0.00159432   0.00595213      0.0          0.0          0.0
 2.84798e-5   0.000113919  0.000427197  0.00159487      0.0          0.0          0.0
 7.63114e-6   3.05246e-5   0.000114467  0.000427344     0.0          0.0          0.0
 2.04476e-6   8.17903e-6   3.06714e-5   0.000114506  …  0.0          0.0          0.0
 ⋮                                                   ⋱
 0.0          0.0          0.0          0.0          …  0.000114466  3.05242e-5   7.63105e-6
 0.0          0.0          0.0          0.0             0.000427192  0.000113918  2.84795e-5
 0.0          0.0          0.0          0.0             0.0015943    0.000425147  0.000106287
 0.0          0.0          0.0          0.0             0.00595002   0.00158667   0.000396668
 0.0          0.0          0.0          0.0             0.0222058    0.00592154   0.00148039
 0.0          0.0          0.0          0.0          …  0.0828731    0.0220995    0.00552487
 0.0          0.0          0.0          0.0             0.309287     0.0824764    0.0206191
 0.0          0.0          0.0          0.0             1.15427      0.307806     0.0769515
 0.0          0.0          0.0          0.0             0.307806     1.14875      0.287187
 0.0          0.0          0.0          0.0             0.0769515    0.287187     1.0718

[dlfivefifty]

Awesome! Will think about how best to incorporate that into the package. Maybe just use your code renamed as InvSymTridiagonal for now. (I'm hopeful these ideas translate to make an InvBandedMatrix but don't have time right now to work out the details.)

[dlfivefifty]

(Of course need to still deal with negative numbers.... signs can be stored separately )

[MikaelSlevinsky]

For the signs, can just promote to complex? (I mean, using exp/log is already way more expensive that calling ldiv! with a Bidiagonal.)

[dlfivefifty]

I don’t see what you mean by more expensive: once the inverse is computed as above its O(1) to evaluate an entry of the inverse, vs O(n) if you compute the entry using ldiv!

[MikaelSlevinsky]

I mean that if you're viewing a column you get n entries in O(n). (Sure, the O(1) formula is nice.)

[dlfivefifty]

Yes duh the LU or UL factorisations will be much faster for applying an inverse... But anyways the more serious case is the block-tridiagonal

[codecov[bot]]

Codecov Report

Merging #13 into master will increase coverage by 0.35%. The diff coverage is 93.75%.

@@            Coverage Diff             @@
##           master      #13      +/-   ##
==========================================
+ Coverage   93.22%   93.57%   +0.35%     
==========================================
  Files           3        4       +1     
  Lines         236      249      +13     
==========================================
+ Hits          220      233      +13     
  Misses         16       16              

Impacted Files	Coverage Δ
src/SemiseparableMatrices.jl	80.00% <ø> (ø)
src/SemiseparableMatrix.jl	89.47% <88.88%> (+4.85%)	:arrow_up:
src/invbanded.jl	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 8eac8f2...d4021e4. Read the comment docs.