Rewrite operators as arrays via InfiniteArrays.jl

[dlfivefifty]

@marcusdavidwebb any objections to this change? I think this package really treats things as (infinite) matrices, so ApproxFun's domain and range space support just gets in the way.

[marcusdavidwebb]

Indeed!

[dlfivefifty]

SymTridiagonal already works!

julia> J = SymTridiagonal(Vcat(randn(5), Fill(0.0,∞)), Vcat(randn(3), Fill(0.5,∞)))
SymTridiagonal{Float64,Vcat{Float64,1,Tuple{Array{Float64,1},Fill{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
 -0.260267  -0.202958    ⋅         …  
 -0.202958   1.53811   -0.315644      
   ⋅        -0.315644  -1.34542       
   ⋅          ⋅         0.0366981     
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅         …  
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅         …  
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅         …  
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅         …  
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅         …  
   ⋅          ⋅          ⋅            
   ⋅          ⋅          ⋅            
  ⋮                                ⋱  

(One catch is the type of dv and ev need to match, which is fine for now.)

[dlfivefifty]

I looked through the operators and I think it's going to be beautiful:

julia> SymTridiagonal(Vcat([1,2], Zeros(∞)), Vcat([3], Zeros(∞))) # currently SymTriOperator
SymTridiagonal{Float64,Vcat{Float64,1,Tuple{Array{Int64,1},Zeros{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
 1.0  3.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
 3.0  2.0  0.0   ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅   0.0  0.0  0.0   ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅   0.0  0.0  0.0   ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅   0.0  0.0  0.0   ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅   0.0  0.0  0.0   ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅   0.0  0.0  0.0     
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.0  0.0     
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.0     
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
 ⋮                        ⋮              ⋱  

julia> SymTridiagonal(Vcat([1,2], Fill(0.0,∞)), Vcat([3], Fill(0.5,∞))) # currently SymTriPertToeplitz
SymTridiagonal{Float64,Vcat{Float64,1,Tuple{Array{Int64,1},Fill{Float64,1,Tuple{InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
 1.0  3.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
 3.0  2.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅   0.5  0.0  0.5   ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅   0.5  0.0  0.5   ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅   0.5  0.0  0.5   ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅   0.5  0.0  0.5   ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅   0.5  0.0  0.5     
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.5  0.0     
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.5     
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   …  
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅      
 ⋮                        ⋮              ⋱  

Some will take some work: PertToeplitz could be done via a lazy Plus, or as an infinite UpperTriangular(::BandedMatrix). HessenbergUnitary should actually be a finite dimensional matrix-type, where the infinite case is handled by c and s being infinite. BandedUnitary can be a lazy Mul.

[dlfivefifty]

This will probably have to wait until banded matrices supports infinite rows (which will be made possible by supporting offset indexing).

A minor question is how to store a 3 x ∞ matrix which is constant on each rows, because when L is stored as a banded matrix we really store the data. Probably the simplest is in low rank form:

Mul([1,2,3], Ones(1,∞))

so that

L_data = ... # finite entries
L = _BandedMatrix(Hcat(L_data, Mul([l⁰,l¹,l²], Ones(1,∞))), Base.OneTo(∞), 2, 0)

(it just occurred to me that LowRankMatrix is really just a lazy multiplication.)