ClassicalOrthogonalPolynomials for Position/Fock transform

[akirakyle]

I recently encountered some numerical instability in the transform function between position and fock bases. This should reproduce it:

using QuantumOptics, Plots

b = FockBasis(99)
bp = PositionBasis(-20, 20, 20000)
ket = fockstate(b, 75);
Tpf = transform(bp, b)
plot(range(bp.xmin, bp.xmax, bp.N), real.((Tpf*ket).data))

This produces the following plot:

This just supposed to be the 76th hermite polynomial rescaled by some factors however one can see the significant numerical artifacts. They're also present for any larger fock states.

This prompted me to investigate, and I think I ruled out the horner function in QuantumOpticsBase's polynomials.jl file (by using the more numerically precise polynomial evaluation function in JuliaMath/Polynomials.jl). I couldn't figure out where exactly in the hermite.A function this numerical instability was coming from but I did find I got similar numerical instability from implementing this using Polynomials.jl with the following

import QuantumOpticsBase: transform
import Polynomials: Polynomial

function hermite_fn(N::T) where T
    A = Vector{Polynomial}(undef, N+1)
    x = Polynomial([0., 1.])
    A[1] = Polynomial([1.])
    A[2] = sqrt(2)*x
    for n=2:N
        A[n+1] = sqrt(2/n)*x*A[n] - sqrt((n-1)/n)*A[n-1]
    end
    A
end

function transform(::Type{S}, b1::PositionBasis, b2::FockBasis) where S
    T = Matrix{S}(undef, length(b1), length(b2))
    xvec = samplepoints(b1)
    A = hermite_fn(b2.N)
    delta_x = spacing(b1)
    c = 1.0/sqrt(sqrt(pi))*sqrt(delta_x)
    for i in 1:length(b1)
        u = xvec[i]
        C = c*exp(-u^2/2)
        for n=b2.offset:b2.N
            idx = n-b2.offset+1
            T[i,idx] = C*A[n+1](u)
        end
    end
    DenseOperator(b1, b2, T)
end

I also tried using Hermite from jverzani/SpecialPolynomials.jl however I still saw this numerical instability.

It seems that ClassicalOrthogonalPolynomials.jl is somehow able to avoid this numerical instability, although I haven't figured out exactly how since the code in this package is pretty advanced. Anyways, this PR does fix this numerical instability, although it does add in quite a few dependencies from ClassicalOrthogonalPolynomials.jl. Also the change in ptranspose is a workaround for a bug I found in QuasiArrays.jl: (JuliaApproximation/QuasiArrays.jl#97). Hence why I'm opening this as a draft since I'm not sure if this is the right approach to fixing this here.

[Krastanov]

This is some great investigative work!

The only worry I have is that the new dependency is a bit young and does not have many dependents yet: https://juliahub.com/ui/Packages/General/ClassicalOrthogonalPolynomials

The new dependency is quite healthy though, belonging to the JuliaApproximations organization, with the same developers as ApproxFun, so it is probably high quality. As long as it does not have long TTFX ;)

I am a bit on the fence because of the dependency, but leaning towards accepting this. I will wait on a final verdict by @david-pl

PS: This is always nice to see:

It might be a new dependency, but it does remove a ton of bespoke code that we do not need to support anymore. Instead, we now use code written by domain experts.

[Krastanov]

The test failures come from Aqua detecting a very large number of method ambiguities in the newly imported dependencies. That is not necessarily a show stopper for us, we already exclude some packages from ambiguity detection see here, but it is something that should at least be on the radar of the developers of the dependency.

@akirakyle , could you

• blacklist the packages so that aqua does not complain about them anymore (see the link above for an example -- you should be able to blacklist whole packages, no need to blacklist separate functions (that would be very teddious and fragile))
• reach out to the issue trackers of the various blacklisted packages and check whether they are aware or plan to do anything about it

Method ambiguities can go from harmless to causing significant TTFX on a small change, so we should be careful.

Other tests seem to run fine, so nothing is broken here.

[Krastanov]

Could you elaborate on why these are deleted?

[akirakyle]

Ah yes I forgot to mention in the PR that I also changed the interface for transform. It would be easy to add this back in but given this optional x0 argument wasn't documented and it wasn't at all clear how it would be used, it was easier for me to remove it as I was working on this.

[akirakyle]

Thanks! That's good to hear, I'm still never quite sure how much to trust various julia packages given both that I'm new to the ecosystem and that the ecosystem is still relatively new, but I do agree that in long term it will be good to move these bespoke functions to domain specific packages like this one.

Also I forgot to add that I did some timings. For

@time Tpf = transform(bp, b)

Before this PR I get for the first run

 0.213219 seconds (128.41 k allocations: 39.285 MiB, 7.49% gc time, 42.21% compilation time)

Then subsequent runs I get

0.151279 seconds (110 allocations: 30.717 MiB, 8.53% gc time)

After this PR I get first run

1.012682 seconds (12.94 M allocations: 450.505 MiB, 7.84% gc time, 78.97% compilation time)

Then subsequent runs I get

0.232888 seconds (9.80 M allocations: 240.023 MiB, 18.45% gc time)

Which I honestly find quite impressive as it seems ClassicalOrthogonalPolynomials.jl has a lot of useful abstractions and fixes numerical instability with seemingly very little cost!

[codecov]

Codecov Report

Merging #137 (b478f18) into master (032dd82) will decrease coverage by 0.15%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #137      +/-   ##
==========================================
- Coverage   93.66%   93.51%   -0.15%     
==========================================
  Files          25       24       -1     
  Lines        3125     3054      -71     
==========================================
- Hits         2927     2856      -71     
  Misses        198      198              

Files Changed	Coverage Δ
src/QuantumOpticsBase.jl	100.00% <ø> (ø)
src/transformations.jl	100.00% <100.00%> (ø)

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[akirakyle]

• blacklist the packages so that aqua does not complain about them anymore (see the link above for an example -- you should be able to blacklist whole packages, no need to blacklist separate functions (that would be very teddious and fragile))

Unfortunately it doesn't appear there's a way to remove entire packages from the ambiguity test (https://juliatesting.github.io/Aqua.jl/dev/#Aqua.test_ambiguities-Tuple{Any}). Instead I just disabled recursive checking of ambiguities. Another possibility would be to keep recursive ambiguity testing but remove checks for ambiguities with Base via

Aqua.test_ambiguities([QuantumOpticsBase, Core]; exclude=[Infinities.ComplexInfinity])

There's a relevant issue in the Aqua.jl repo discussing this here: JuliaTesting/Aqua.jl#77

• reach out to the issue trackers of the various blacklisted packages and check whether they are aware or plan to do anything about it

There's 9 direct dependencies of ClassicalOrthogonalPolynomials and two transitive dependencies (Infinities and SemiseparableMatrices) with ambiguity issues. Given they're all under various Julia* github organizations and it appears they share developers in common I went ahead and just opened an issue ClassicalOrthogonalPolynomials explaining the situation (JuliaApproximation/ClassicalOrthogonalPolynomials.jl#154).

[dlfivefifty]

Note OrthonormalWeighted(Normalized(Hermite())) still has issues with underflow:

julia> norm(H[40,1:100])
0.0

There's a much better way to do this: take the n-th root w_n = exp(-x^2/(2n)) and incorporate it as you build the recurrence. This is something I haven't got around to but there is some version of this method here:

https://github.com/JuliaApproximation/FastGaussQuadrature.jl/blob/a5932947c43186d1e13730015da705b5052012f5/src/gausshermite.jl#L130

Given your concerns with ambiguities and since you only need evaluation I'm wondering if a separate package (OrthogonalPolynomialRecurrences.jl) would be worthwhile? The design of ClassicalOrthogonalPolynomials.jl depends on InfiniteArrays.jl to avoid rewriting code for each family but this can be avoided by just rewriting the code for each family.

[dlfivefifty]

Or just use FastGaussQuadrature.jl? This is the relationship:

julia> OrthonormalWeighted(Hermite())[x, 1:n]  ≈ FastGaussQuadrature.hermpoly_rec(0:n-1, sqrt(2)* x) * π^(-1/4)
true

(So you would also need to add the normalisation constant.) But then this avoids underflow:

julia> x = 40; n = 100; FastGaussQuadrature.hermpoly_rec(0:n-1, sqrt(2)* x) * π^(-1/4)
100-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 9.961817500563452e-262
 5.804365934206908e-261
 3.361991950596754e-260
 1.9359355316676075e-259
 1.1083157837936175e-258
 6.308702574489045e-258
 3.570634917149043e-257
 2.0095745538076928e-256
 1.1247085178318906e-255
 6.2600303429468115e-255
 3.4652646464337584e-254
 1.9078417421194034e-253

[david-pl]

@Krastanov thanks for the ping. The dependency looks solid to me, so feel free to merge as soon as you're happy with the code.

@akirakyle thanks for the deep dive you did here ;)

[akirakyle]

Or just use FastGaussQuadrature.jl?

Thanks @dlfivefifty for pointing out hermpoly_rec in FastGaussQuadrature.jl! I somehow missed that when I was looking through it and thought from the documentation that FastGaussQuadrature.jl only implemented functions that efficiently found the roots of the hermite polynomials without actually needing to evaluate them directly.

(So you would also need to add the normalisation constant.) But then this avoids underflow:

Ah I also missed this underflow in my testing. Thanks for explicitly giving the relationship between Hermite and hermpoly_rec, I've created a new PR at #142 that should be preferred to this one.

[Krastanov]

superseded by #142

[Krastanov]

Thank you @dlfivefifty for the helpful suggestions!