Inner products between functions on (un)restricted bases

[jagot]

Assume a linear operator Ĥ with a discrete and continuous
spectrum. Then any function may be expanded over its eigenfunctions:

When calculating inner products and matrix elements between such
expansions, one encounters inner products of the kinds

Usually, the eigenfunctions of the discrete spectrum have compact
support, or are at least very localized (they may have an
exponentially decaying tail, which is negligible). It would be nice to
be able to calculate such inner products, without having to represent
the discrete eigenfunctions on the entire interval. Maybe something
like this?:

B <: AbstractQuasiMatrix
B̃ = B[1:100]
i = B̃*rand(100)
k = B*rand(size(B,2))

# This is the syntax I'd like to automagically work:
i'k

As a follow-up to this, I'll need to solve the Poisson problem on the
entire interval, with source terms of the form

where the source term will be non-zero only on the interval where
is.

[dlfivefifty]

So the issue is essentially having B'(B[:,1:100]) be simplified as (B'B)[:,1:100]?

[jagot]

Yes, I think that might be correct.

[dlfivefifty]

OK I'll look into it, but first I need to finish the Applied PR in LazyArrays. This is helpful as it allows us more flexibility by overriding ApplyStyle, so that things simplify in the right way without complicating Mul for unrelated objects.

We are touching on a more complicated topic on automatically simplifying operator trees, but that's beyond my ken.

[jagot]

I discovered today that this already seems to work:

julia> t = range(0,stop=1,length=11)
0.0:0.1:1.0

julia> B = FEDVR(t, 4)
FEDVR{Float64} basis with 10 elements on 0.0..1.0

julia> B̃ = B[:,1:5]
FEDVR{Float64} basis with 10 elements on 0.0..1.0, restricted to basis functions 1..5 ⊂ 1..31

julia> u = B*ones(size(B,2))
ContinuumArrays.QuasiArrays.ApplyQuasiArray{Float64,1,LazyArrays.Applied{ContinuumArrays.QuasiArrays.LazyQuasiArrayApplyStyle,typeof(*),Tuple{FEDVR{Float64,Float64,FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}},Array{Float64,1}}}}(LazyArrays.Applied{ContinuumArrays.QuasiArrays.LazyQuasiArrayApplyStyle,typeof(*),Tuple{FEDVR{Float64,Float64,FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}},Array{Float64,1}}}(ContinuumArrays.QuasiArrays.LazyQuasiArrayApplyStyle(), *, (FEDVR{Float64} basis with 10 elements on 0.0..1.0, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])))

julia> v = B̃*ones(size(B̃,2))
ContinuumArrays.QuasiArrays.ApplyQuasiArray{Float64,1,LazyArrays.Applied{ContinuumArrays.QuasiArrays.LazyQuasiArrayApplyStyle,typeof(*),Tuple{FEDVR{Float64,Float64,FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}},Array{Float64,1}}}}(LazyArrays.Applied{ContinuumArrays.QuasiArrays.LazyQuasiArrayApplyStyle,typeof(*),Tuple{FEDVR{Float64,Float64,FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}},Array{Float64,1}}}(ContinuumArrays.QuasiArrays.LazyQuasiArrayApplyStyle(), *, (FEDVR{Float64} basis with 10 elements on 0.0..1.0, [1.0, 1.0, 1.0, 1.0, 1.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> u'v
5.0

julia> B'B̃
31×5 BandedMatrices.BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 1.0   ⋅    ⋅    ⋅    ⋅
  ⋅   1.0   ⋅    ⋅    ⋅
  ⋅    ⋅   1.0   ⋅    ⋅
  ⋅    ⋅    ⋅   1.0   ⋅
  ⋅    ⋅    ⋅    ⋅   1.0
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅

To make it more efficient, I basically need to dispatch copyto! on the restricted Applied, to avoid expanding v into the full vector (with trailing zeros) before the dot product.

[dlfivefifty]

Cool! Does the “extra” allocation actually affect much?

[jagot]

I realized that will not work :( since when forming v above, the Vector is already then expanded to the full size with trailing zeros:

julia> v.applied.args[2]
31-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
 1.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
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

[jagot]

Yes, it will affect me in two ways:

1. Calculating scalar products between bound states (which are fairly compact, maybe 100 points or so) and continuum states (which are of infinite extents, i.e. thousands of points) will be more expensive.
2. For solving Poisson's problem, I need to preserve the restrictedness of the basis to compute the right Laplacian (e.g. excluding r=0).

[dlfivefifty]

Is the goal to have it be allocation free:

Vcat(Ones(5), Zeros(26))

?

In the ∞-Jacobi case I do achieve this allocation-free form so it should be possible, though may require a few "hacks" for special forms of Vcat, which is less than ideal. (I'd like LazyArrays.jl to become less hacky, not more...)

[jagot]

Yes and no, as few allocations/unnecessary computations as possible is of course nice, but a bigger problem is the Laplacian; if you look at the operator in square brackets in the top post of this issue, you'll notice that it is undefined for r=0, hence I must exclude the first basis function (my boundary conditions anyway require y(0)=0).

When I created the issue, I was mainly thinking about efficiency, but today when working on getting the Poisson problem to work with FEDVRQuasi, I just got a lot of NaNs back, due to the first element of the Laplacian being Inf.

I have a solution that is semi-hacky, I'll post a link with code for you to critique :) when I've put together a test.

[dlfivefifty]

Special casing r = 0 seems like exactly what Vcat is good for... unless I'm missing something.

[jagot]

It may be I simply misunderstood how to achieve what you meant. Anyway, here's what I do:

I first define a set of type aliases, essentially unions of (un)restricted quasiarrays and their adjoints:
https://github.com/jagot/FEDVRQuasi.jl/blob/3be9d880275fb4e7629a709d63408bed982f592d/src/FEDVRQuasi.jl#L25
https://github.com/jagot/FEDVRQuasi.jl/blob/3be9d880275fb4e7629a709d63408bed982f592d/src/FEDVRQuasi.jl#L108

To calculate the inner products, I do the following:
https://github.com/jagot/FEDVRQuasi.jl/blob/3be9d880275fb4e7629a709d63408bed982f592d/src/FEDVRQuasi.jl#L317

The same example as above:

t = range(0,stop=1,length=11)
B = FEDVR(t, 4)
B̃ = B[:,1:5]

uv = ones(size(B,2))
u = B*uv
lu = B ⋆ uv

vv = ones(size(B̃,2))
v = B̃*vv
lv = B̃ ⋆ vv

julia> @test u'v == 5
Test Passed

julia> lazyip = lu' ⋆ lv
[1.0 1.0 … 1.0 1.0]⋆ContinuumArrays.QuasiArrays.QuasiAdjoint{Float64,FEDVR{Float64,Float64,FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}}}(FEDVR{Float64} basis with 10 elements on 0.0..1.0)⋆FEDVR{Float64} basis with 10 elements on 0.0..1.0, restricted to basis functions 1..5 ⊂ 1..31⋆[1.0, 1.0, 1.0, 1.0, 1.0]

julia> @test lazyip isa FEDVRQuasi.LazyFEDVRInnerProduct
Test Passed

julia> @test materialize(lazyip) == 5
Test Passed

Both the lazy and the non-lazy versions will only sum those elements that actually are non-zero, but the non-lazy versions will first have expanded the vectors to the full lengths, with trailing zeros.

[jagot]

I find my returning to this issue every so often :)

Something changed in the materialization algorithm, I think the current approach is to flatten as much as possible. I found some discrepancies:

julia> R = FEDVR(range(0,stop=1,length=3),3)[:,2:end-1]
FEDVR{Float64} basis with 2 elements on 0.0..1.0, restricted to basis functions 2..4 ⊂ 1..5

julia> c = rand(size(R,2))
3-element Array{Float64,1}:
 0.6897536131925821
 0.29789024566149047
 0.043702178352102994

julia> R*c
QuasiArrays.ApplyQuasiArray{Float64,1,typeof(*),Tuple{FEDVR{Float64,Float64,Fill{Int64,1,Tuple{Base.OneTo{Int64}}}},Array{Float64,1}}}(*, (FEDVR{Float64} basis with 2 elements on 0.0..1.0, [0.0, 0.6897536131925821, 0.29789024566149047, 0.043702178352102994, 0.0]))

julia> (R*c).args[2]
5-element Array{Float64,1}:
 0.0
 0.6897536131925821
 0.29789024566149047
 0.043702178352102994
 0.0

julia> materialize(applied(*, R, c))
QuasiArrays.ApplyQuasiArray{Float64,1,typeof(*),Tuple{FEDVR{Float64,Float64,Fill{Int64,1,Tuple{Base.OneTo{Int64}}}},BandedMatrices.BandedMatrix{Int64,Ones{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}},Array{Float64,1}}}(*, (FEDVR{Float64} basis with 2 elements on 0.0..1.0, [0 0 0; 1 0 0; … ; 0 0 1; 0 0 0], [0.6897536131925821, 0.29789024566149047, 0.043702178352102994]))

julia> materialize(applied(*, R, c)).args[3]
3-element Array{Float64,1}:
 0.6897536131925821
 0.29789024566149047
 0.043702178352102994

R*c has the effect of dropping the restriction matrix and padding with zeros (which I don't want, because I lose the restriction information), whereas materialize(applied(*, R, c)) now produces a flat structure with 3 arguments: the basis, the restriction matrix, and the coefficients. Previously, this was a nested structure: ((basis, restriction), coefficients), which simplified dispatching on the basis type, i.e. (basis, restriction) or simply basis. I guess I should stick to applied(*, R, c), but the R*c syntax is very nice. Is this where I should start writing @~ R*c instead?

[dlfivefifty]

Sorry for the constant changes, getting the ∞-dimensional case of OPs working is quite finicky, though I think its there now.

Note you can use apply(*, R, c) instead of materialize(applied(*, R, c)).

I see the problem. I think representing restrictions as a lazy multiplication was probably a bad idea, and we should leave restrictions un-simplified, that is, as a view. I'll experiment on the linear splines in ContinuumArrays.jl

[jagot]

No problem, it's good to get a long-term solution that's nice.

If I use LazyArrays' approach, I get this

julia> v = @~ R*c
FEDVR{Float64} basis with 2 elements on 0.0..1.0, restricted to basis functions 2..4 ⊂ 1..5⋆[0.6897536131925821, 0.29789024566149047, 0.043702178352102994]

julia> materialize(@~ v'v)
0.5664085257652043

but I can't do this:

julia> v'v
ERROR: MethodError: no method matching *(::Applied{QuasiArrays.QuasiArrayApplyStyle,typeof(*),Tuple{Adjoint{Float64,Array{Float64,1}},QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{Adjoint{Int64,BandedMatrices.BandedMatrix{Int64,Ones{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}}},QuasiArrays.QuasiAdjoint{Float64,FEDVR{Float64,Float64,Fill{Int64,1,Tuple{Base.OneTo{Int64}}}}}}}}}, ::Applied{LazyArrays.FlattenMulStyle,typeof(*),Tuple{QuasiArrays.ApplyQuasiArray{Float64,2,typeof(*),Tuple{FEDVR{Float64,Float64,Fill{Int64,1,Tuple{Base.OneTo{Int64}}}},BandedMatrices.BandedMatrix{Int64,Ones{Int64,2,Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},Base.OneTo{Int64}}}},Array{Float64,1}}})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:529
  *(::PyCall.PyObject, ::Any) at /Users/jagot/.julia/packages/PyCall/ttONZ/src/pyoperators.jl:13
  *(::AlgebraicMultigrid.Preconditioner, ::Any) at /Users/jagot/.julia/packages/AlgebraicMultigrid/eCNmE/src/preconditioner.jl:37
  ...
Stacktrace:
 [1] top-level scope at REPL[121]:1

which is maybe a bit surprising.

EDIT:

julia> apply(*, v', v)
0.5664085257652043

works.

[dlfivefifty]

That's an easy fix I think: just need to add *(A::Applied, B::Applied) = apply(*, A, B) in LazyArrays.jl