JuliaApproximation / JuliaApproximation/ApproxFunOrthogonalPolynomials.jl

ApproxFun+IntervalArithmetic: Support evaluation with infinite vectors

Open
#3 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
Julia
Stars
4
Forks
11
Avg merge
5h 31m
Merged PRs (30d)
4

Description

It's possible to represent functions with infinite coefficients using InfiniteArrays.jl, however, evaluation is broken:

julia> f = Fun(Chebyshev(), -2.0 .^(-(1:∞)))
Fun(Chebyshev(),[-0.5, -0.25, -0.125, -0.0625, -0.03125, -0.015625, -0.0078125, -0.00390625, -0.00195313, -0.000976563  …  ])

julia> f(0.1)
ERROR: ArgumentError: Cannot create range starting at infinity
Stacktrace:
 [1] (::Colon)(::InfiniteArrays.Infinity, ::Int64, ::Int64) at /Users/solver/Projects/InfiniteArrays.jl/src/infrange.jl:10
 [2] chebyshev_clenshaw(::BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/LinearAlgebra/clenshaw.jl:243
 [3] clenshaw(::Chebyshev{ChebyshevInterval{Float64},Float64}, ::BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:73
 [4] evaluate at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/PolynomialSpace.jl:23 [inlined]
 [5] evaluate at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:216 [inlined]
 [6] (::Fun{Chebyshev{ChebyshevInterval{Float64},Float64},Float64,BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}})(::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:220
 [7] top-level scope at none:0

It's actually a bit simpler to think about in the interval arithmetic setting where we only need a bound on evaluation (otherwise we would need to talk about tolerances). That is, in the format that #2 would return:

julia> f = Fun(Chebyshev(), BroadcastArray(IntervalArithmetic.Interval, -2.0 .^(-(1:∞)), 2.0 .^(-(1:∞))))
Fun(Chebyshev(),IntervalArithmetic.Interval{Float64}[[-0.5, 0.5], [-0.25, 0.25], [-0.125, 0.125], [-0.0625, 0.0625], [-0.03125, 0.03125], [-0.015625, 0.015625], [-0.0078125, 0.0078125], [-0.00390625, 0.00390625], [-0.00195313, 0.00195313], [-0.000976563, 0.000976563]  …  ])

julia> f(0.1)
ERROR: ArgumentError: Cannot create range starting at infinity
Stacktrace:
 [1] (::Colon)(::InfiniteArrays.Infinity, ::Int64, ::Int64) at /Users/solver/Projects/InfiniteArrays.jl/src/infrange.jl:10
 [2] chebyshev_clenshaw(::BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/LinearAlgebra/clenshaw.jl:243
 [3] clenshaw(::Chebyshev{ChebyshevInterval{Float64},Float64}, ::BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}, ::Float64) at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/Chebyshev/Chebyshev.jl:73
 [4] evaluate at /Users/solver/Projects/ApproxFunOrthogonalPolynomials.jl/src/Spaces/PolynomialSpace.jl:23 [inlined]
 [5] evaluate at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:216 [inlined]
 [6] (::Fun{Chebyshev{ChebyshevInterval{Float64},Float64},IntervalArithmetic.Interval{Float64},BroadcastArray{IntervalArithmetic.Interval{Float64},1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},Type{IntervalArithmetic.Interval},Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(-),Tuple{BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}},BroadcastArray{Float64,1,Base.Broadcast.Broadcasted{LazyArrays.LazyArrayStyle{1},Tuple{InfiniteArrays.OneToInf{Int64}},typeof(^),Tuple{Float64,InfiniteArrays.InfStepRange{Int64,Int64}}}}}}}})(::Float64) at /Users/solver/Projects/ApproxFunBase.jl/src/Fun.jl:220
 [7] top-level scope at none:0

For Chebyshev it should be easy to work out bounds and override chebyshev_clenshaw(::BroadcastArray{...}, x) to return these bounds. We could also then override chebyshev_clenshaw(a::Vcat, x) to do Clenshaw recursively, using the result of chebyshev_clenshaw(a.arrays[end], x) to determine the initial condition for the rest of Clenshaw. This would give rigorous evaluation of functions.

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start in src/LinearAlgebra/clenshaw.jl at the failing range construction and trace the Chebyshev dispatch through src/Spaces/Chebyshev/Chebyshev.jl. Determine how infinite BroadcastArray coefficients and Vcat inputs should produce rigorous evaluation bounds, then verify that evaluating the examples no longer raises the infinite-range error.

Written by the indexing model from the issue text.

Assessment

Tech stack
julia
Domain
backend
Issue type
Feature
Difficulty
5/5
Estimated time
Over a week
Activity status
Stale
Clarity
Mostly clear
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.