JuliaApproximation
Rewrite operators as arrays via InfiniteArrays.jl
@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.
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.)
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.
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.)