A Julia repository for linear algebra with infinite banded and block-banded matrices

This currently supports the infinite-dimensional QR factorization for banded matrices, also known as the adaptive QR decomposition as the entries of the QR decomposition are determined lazily.

As a simple example, consider the Bessel recurrence relationship:

```
julia> using InfiniteLinearAlgebra, InfiniteArrays, BandedMatrices, FillArrays, SpecialFunctions
julia> z = 10_000; # the bigger z the longer before we see convergence
julia> A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
∞×∞ BandedMatrix{Float64,ApplyArray{Float64,2,typeof(*),Tuple{Array{Float64,2},ApplyArray{Float64,2,typeof(vcat),Tuple{Transpose{Float64,InfiniteArrays.InfStepRange{Float64,Float64}},Ones{Float64,2,Tuple{Base.OneTo{Int64},InfiniteArrays.OneToInf{Int64}}},Ones{Float64,2,Tuple{Base.OneTo{Int64},InfiniteArrays.OneToInf{Int64}}}}}}},InfiniteArrays.OneToInf{Int64}} with indices OneToInf()×OneToInf():
0.0 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
1.0 -0.0002 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 1.0 -0.0004 1.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 1.0 -0.0006 1.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 1.0 -0.0008 1.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 1.0 -0.001 1.0 ⋅ ⋅ …
⋅ ⋅ ⋅ ⋅ ⋅ 1.0 -0.0012 1.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 -0.0014 1.0
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 -0.0016
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋮ ⋮ ⋱
```

The first row corresponds to specifying an initial condition. Thus we can determine the Bessel functions via solving the recurrence:

```
julia> A \ [besselj(1,z); Zeros(∞)]
∞-element LazyArrays.CachedArray{Float64,1,Array{Float64,1},Zeros{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}} with indices OneToInf():
-0.007096160353406478
0.0036474507555295833
0.007096889843557584
-0.0036446119995921654
-0.007099076610757337
0.0036389327383035616
0.00710271554349564
-0.003630409479651369
-0.007107798116767152
0.0036190370026645442
0.007114312383371949
-0.0036048083778978026
-0.0071222429618033245
0.0035877149947894766
0.007131571020789777
⋮
julia> J[1000] - besselj(999,z) # matches besselj to high (relative) accuracy
-6.8252695162307475e-15
julia> J[11_000] - besselj(11_000-1, z)
3.3730094946097293e-143
```

We're even faster than SpecialFunctions.jl for constructing a range of Bessel functions:

```
julia> @time [besselj(k-1, z) for k=0:11_000-1];
0.188690 seconds (77.20 k allocations: 3.295 MiB)
julia> @time J = A \ Vcat([besselj(1,z)], Zeros(∞));
0.006354 seconds (90.93 k allocations: 6.791 MiB)
```

This currently supports the infinite-dimensional QL factorization for perturbations of Toeplitz operators. Here is an example:

```
# Bull head matrix
A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
ql(A - 5*I)
```

The infinite-dimensional QL factorization is a subtly thing: its defined when the operator has non-positive Fredholm index, and if the Fredholm index is not zero, it may not be unique. For the Bull head matrix `A`

, here are plots of `ql(A-λ*I).L[1,1]`

alongside the image of the symbol `A`

, which depicts the essential spectrum of `A`

and where the Fredholm index changes. Note we have two plots as the regions with negative Fredholm index have multiple QL factorizations. Where the Fredholm index is positive, the QL factorization doesn't exist and is depicted in black.