Move adaptive QR to InfiniteArrays.jl

[dlfivefifty]

I think it's a bad design that ApproxFun conflates operators and infinite matrices: for example, something like adaptive QR works completely on the matrix level and doesn't need to know any spaces.

I think it would be better to create a new package called InfiniteArrays.jl, to represent (by necessity lazy) infinite arrays.

We would want to be able to do things like the following, which recreates the ApproxFun \ without talking explicitly about spaces:

   B = ones(1,∞)  # creates 1 x ∞ InfiniteMatrix of all ones
   T = diagm(repeated(0.5)) + diagm(repeated(-0.5),2)  # a Toeplitz operator
   e₁ = [1; zeros(∞)]   # a PaddedInfiniteVector
   (e₁*e₁')/2 # a PaddedInfiniteMatrix: it knows to manipulate on the data
   C = T + (e₁*e₁')/2    # a PlusInfiniteMatrix, it adds the rank-1 perturbation lazily
   D = diagm(1:∞,1)   # Chebyshev derivative, represented as a DiagonalInfiniteMatrix
   A = [B; D+C]   # A ConcatenatedInfiniteMatrix
   u = A \ [2; zeros(∞)]  # Finds Chebyshev coefficients of solution to u'+u = 0, u(1) = 2, using adaptive QR

@ChrisRackauckas I think we've been dancing around this idea for a while now. Any thoughts?

[dlfivefifty]

Here's a quick sketch.

abstract type InfiniteArray{T,N} end  # maybe AbstractInfiniteArray

struct PaddedInfiniteArray{T,N} <: InfiniteArray{T,N}
   data::Array{T,N}
end

size(P::PaddedInfiniteArray, k::Int) = ∞   # this is a special infinity currently in ApproxFun
function getindex(P::PaddedInfiniteArray{T,1}, k::Int)
      k ≤ length(P.data) && return P.data[k]
      zero(T)
end


# convert an iterator to an infinite-vector lazily
struct CachedIterator{T, ITERATORTYPE} <: InfiniteVector{T}
  itr::ITTERATORTYPE
  cache::Vector{T}
end


struct DiagonalInfiniteMatrix{T,DT<:InfiniteVector} <: InfiniteMatrix{T}
   diag::DT
   k::Int
end
DiagonalInfiniteMatrix(diag::DT, k::Int) where DT<:InfiniteVector{T} where T = DiagonalInfiniteMatrix{T,DT}(diag,k)

diagm(diag::InfiniteVector, k) = DiagonalInfiniteMatrix(diag,k)
DiagonalInfiniteMatrix(itr) = DiagonalInfiniteMatrix(CachedIterator(itr))

[ChrisRackauckas]

I think this is highly desirable. Sometimes the operator form on ApproxFun Funs is desirable, but in many cases it can be good to want to use the underlying matrix form itself, and de-lazifying it to a matrix (well, BandedMatrix so it's not so bad) all the time seemed a little odd. Having the infinite matrix form be itself representable I think can have some advantages, though no "must have" feature immediately comes to mind.

[dlfivefifty]

OK, I'll try to sketch up an initial version. Let me know if you have any suggested improvements to the sketch I put above.

[dlfivefifty]

One slight issue: how to handle

[D C; 
 C D]

that is, an hvcat of infinite-matrices.

In ApproxFun this is resolved by interlacing the rows/columns, which is allowed because we spaces attached, and therefore row ordering is in some sense an undefined notion. But for infinite matrices this is definitely not allowed.

I guess the easy answer is don't allow it at all.

[ChrisRackauckas]

Yeah, I'm not exactly sure what that would mean for an infinite matrix.

[dlfivefifty]

I've started the package here

https://github.com/JuliaApproximation/InfiniteArrays.jl

1:∞ works pretty well, I think this is going to be a fun project. When this change takes place, ApproxFun will almost consist of just two types 😆:

Fun(::Space, ::InfVector)
Operator(domainspace::Space, rangespace::Space, ::InfMatrix)

And there'll be a nice duality: A Fun is an InfVector with a space attached, An Operator is an InfMatrix with domain and range space attached.

But it's pretty clear that it's better to make this a Julia v0.7 only package, as we need to mirror Base. As soon as Juno works on 0.7, I'll pick this up again.