How to do a hierarchy of OPs?

[dlfivefifty]

We need to decide how to redesign the implemention with a hierarchy of OPs in mind. The cases we have are

1. Annuli: SemiclassicalJacobi.(t, 1, 1, -1:M) needed for annuli.
2. Spherical cap: SemiclassicalJacobi.(t, (0:M) .+ 1/2, 0, (0:M) .+ 1/2)

Two options:

1. Have OPs link to a lower order one, e.g., SemiclassicalJacobi(t, 1, 1, M) will contain a pointer to SemiclassicalJacobi(t, 1, 1, M-1). Each OP would cache its Jacobi operators which contain all the necessary information to build other operators. To extend this to the other case we would always lower the higher parameter, so SemiclassicalJacobi(t, M + 1/2, 0, M + 1/2) would have a pointer to SemiclassicalJacobi(t, M - 1/2, 0, M + 1/2) which would have a pointer to SemiclassicalJacobi(t, M - 1/2, 0, M - 1/2).
2. Store a dictionary for the whole family. This is what @snowball13 does in his ApproxFun-based code. Then each OP also has a pointer back to the dictionary. This is arguably more versatile and allows more to be cached.

In terms of interface I'm picturing something like:

M = 10 # Total degree
Ps = SemiclassicalJacobi.(t, 1, 1, -1:M) # overloaded to either return a vector of OPs as in (1), or a special type `SemiclassicalJacobiFamily` which caches the info.

We would probably want to support growing the cache for adaptivity. Either approach would support this naturally.

[dlfivefifty]

Er we probably need to also consider that we actually want sqrt-weighted OPs, so for

sqrtsemijacobi(t, a, b, c) = SemiclassicalJacobiWeight(t, a/2, b/2, c/2) .* SemiclassicalJacobi(t, a, b, c)

we actually want to work with

1. Annuli: sqrtsemijacobi.(t, 1, 1, -1:M)
2. Disk: sqrtsemijacobi.(t, (0:M) .+ 1/2, 0, (0:M) .+ 1/2)

Arguably this should just be a special type SqrtWeightedSemiclassicalJacobi. So the interface would be:

M = 10 # Total degree
Ps = SqrtWeightedSemiclassicalJacobi.(t, 1, 1, -1:M) # overloaded to either return a vector of SqrtWeightedSemiclassicalJacobi as in (1), or a special type `SqrtWeightedSemiclassicalJacobiFamily` which caches the info.

[dlfivefifty]

Actually we probably want to add Weighted(...) and SqrtWeight(...) for any orthogonal polynomial basis.

[MikaelSlevinsky]

For 2-normalized SqrtWeightedSemiclassicalJacobi, the connection coefficients within the hierarchy can be described by Givens rotations, either using Theorem 4.6 or 4.9 of our review, which relate the sines and cosines of the Givens rotations to leading-degree coefficients (which can be derived from the recurrence coefficients / Jacobi operator).

[dlfivefifty]

Now SemiclassicalJacobi are normalized, though with SemiclassicalJacobi(t,a,b,c)[x,1] == 1.0 Thanks for the reminder about Theorem 4.6. Maybe a starting point is in ClassicalOrthogonalPolynomials.jl to implement the orthogonal connection matrices for SqrtWeighted(Normalized(Jacobi(a,b))).

[dlfivefifty]

Probably SqrtWeighted is to ambiguous: e.g. for Hermite are we weighting by exp(-x^2/2) or exp(-x^2/2) * const? OrthonormalWeighted might be a better term since we are "weighted so that they are orthonormal with respect to L^2`.

The same question exists with Weighted, though I guess there we can take it to be equivalent to orthogonalityweight(P) .* P. But perhaps we don't need Weighted for the time being.

Let me cook up an example for Hermite of how I see this being implemented.

[MikaelSlevinsky]

I'd say we must define a given weight in the docs. For the classical ones, it only makes sense to me to use the DLMF standards since otherwise pointwise evaluation of the weight is more expensive: no need to say w(x) = const * exp(-x^2), when it's faster to say const == 1. Same thing with Jacobi, where the constant to make its integral one involves the gamma/beta functions, which is more complicated and expensive than the x-dependence.

When 2-normalized, there's always gonna be ambiguity with respect to phase, so as long as we're consistent with our own docs it's gonna be fine.

[dlfivefifty]

The constant only impacts evaluation, it doesn't impact operators, and so I don't think it will actually make a difference performance wise.

[dlfivefifty]

Also we already do the constant in Normalized and it is cached, so doesn't really make any difference...

[dlfivefifty]

We also need to combine an infinite number of banded operators. For example, conversion for Zernike polynomials involves [Normalized(Jacobi(1,m)) \ Normalized(Jacobi(0,m)) for m = 0:∞]:

https://github.com/JuliaApproximation/MultivariateOrthogonalPolynomials.jl/blob/763954b5904e73c40071d98f30977fb94c90435f/src/disk.jl#L203

I guess this could be generalised and syntax could be, in the semiclassical case,

ZernikeInterlace(SemiclassicalJacobi.(t, 1, 1, 1:∞) .\ SemiclassicalJacobi.(t, 0, 0, 0:∞))

That is, SemiclassicalJacobi.(t, 1, 1, 1:∞) would return a special type representing an infinite hierarchy of OPs and then we can overload .\ to represent the infinite operators.

[dlfivefifty]

I don't think we need a special type after all: SemiclassicalJacobi.(t, 0, 0, 0:100) is quite fast already, building in order.

The only possible issue is weighted derivatives such as D * Weighted(SemiclassicalJacobi(t, 1, 1, 100)), which need to construct SemiclassicalJacobi(t, 0, 0, 99). At the moment I know how to build a SemiclassicalJacobi from lower parameters but not from higher parameters. While hidden in the Jacobi operator of say SemiclassicalJacobi(t, 1, 1, 100)) is the Jacobi operator for SemiclassicalJacobi(t,1,1,99), this doesn't help getting to SemiclassicalJacobi(t, 0, 0, 99).

There are some options for handling this:

1. Figure out the formula for "lowering", e.g. constructing SemiclassicalJacobi(t, a, b, c-1) given knowledge of SemiclassicalJacobi(t, a, b, c). @TSGut does your approach for SemiclassicalJacobi(t, 0, 0, -1) extend to this?
2. Keep a big "database" of constructed OPs, so that SemiclassicalJacobi.(t, 0, 0, 0:∞)[5] will return a view into this database, with the ability of looking up other parameters and caching on the fly. We could just cache on the fly. This is close to the design @snowball13 did for the disk OPs.
3. Leave D * Weighted(SemiclassicalJacobi(t, 1, 1, 100)) unmaterialized, which is materialized at the stage Weighted(SemiclassicalJacobi(t, 0, 0, 99)) \ SemiclassicalJacobi(t, 1, 1, 100)). Since this is passing in an already construct SemiclassicalJacobi there is no issue.

[TSGut]

Figure out the formula for "lowering", e.g. constructing SemiclassicalJacobi(t, a, b, c-1) given knowledge of SemiclassicalJacobi(t, a, b, c). @TSGut does your approach for SemiclassicalJacobi(t, 0, 0, -1) extend to this?

I'll give this some thought. There definitely are generalizations of the approach but I need to check if this particular one is among them.

[TSGut]

Going back to the inspiration for my approach in Price 1979 again, https://epubs.siam.org/doi/abs/10.1137/0716073, I believe the approach does generalize in the way you say. The approach I used is equivalent to using Psi = Legendre polynomials in Example 1 (again, careful the t and x convention is inverted from our own).

But nothing says we have to use Legendre there - if we know the polynomials wrt w(x) using that approach we can get the orthogonal polynomials wrt w(x)/(t-x) which is exactly the c-1 step you seem to want?

[TSGut]

Actually, it looks like this would probably more naturally lead to raising operators. I'll see if there is a way to get a lowering operator out of it.

[dlfivefifty]

Raising operators are already done, that's who we construct P(a+1,b,c) from P(a,b,c) (etc.)

[TSGut]

The only lowering operator that I think can be obtained without too much effort from this approach would be that wrt multiplication by (t-x), i..e

(t-x) Q_n^{t,a,b,c}(x) = A*Q_n^{t,a,b,c-1}(x)+B*Q_{n+1}^{t,a,b,c-1}(x)

I implicitly already have this for the Legendre / c=-1 case and I think it could be generalized. Useful for what you are doing?

[dlfivefifty]

That's equivalent to the raising operator from a,b,c-1 to a,b,c...

it's not actually the raising operators we need, just the Jacobi operators. That is, given the Jacobi operator for a,b,c I want to derive the Jacobi operator for a,b,c-1. I thought this is exactly what you do just for c = 0 case?

[TSGut]

Indeed, which is why I was a bit confused about what was being asked, I think I got it now though.

Since this derivation of the general case will be a pain to go through (it already was for the -1 case), just to be fully explicit so I only do this if it's actually useful: You want to find values α, β and γ such that they satisfy

x Q_n^{t,a,b,c-1}(x) = α*Q_n^{t,a,b,c-1}(x) + β*Q_{n-1}^{t,a,b,c-1}(x) + γ*Q_{n+1}^{t,a,b,c-1}(x)

depending on values A, B, C which we assume satisfy

x Q_n^{t,a,b,c}(x) = A*Q_n^{t,a,b,c}(x) + B*Q_{n-1}^{t,a,b,c}(x) + C*Q_{n+1}^{t,a,b,c}(x) ? With obvious implication for the Jacobi operators.

If you say this is exactly what we want, I'll start chipping away at the derivation.

Also, if yes: More for convenience in said derivation than anything else, it should be fine to assume the non-normalized OP non-symmetric Jacobi operators are given for a,b,c, then do the sqrt() thing in the implementation right?

[dlfivefifty]

If you say this is exactly what we want

Yes. I already do the other way around (find A,B,C from α,β,γ) by first finding the lowering operator using CD and then using that to determine the Jacobi operator, the sketch is in the draft. But I'm not sure how to find α,β,γ from A,B,C.

then do the sqrt() thing in the implementation right?

Sure.

[TSGut]

But I'm not sure how to find α,β,γ from A,B,C.

If it's like the c=-1 case I just finished, then this result should come out of the knowledge that Q_n^{t,a,b,c}(x) satisfies a three term recurrence with assumed known coefficients in addition to the explicit form of the c-1 case in terms of the c case which is Q_n^{t,a,b,c-1}(x) = α(n,t)*Q_{n-1}^{t,a,b,c}(x)+Q_n^{t,a,b,c}(x) (possibly an additional normalization constant somewhere) where α(n,t) is given by the moments computed via recurrence.

Alright, I'll get to work on deriving said relationship.

[TSGut]

I think I may have it. @dlfivefifty Is there a straightforward way to get non-symmetric Jacobi matrices for non-normalized SemiclassicalJacobi (e.g. via Lanczos) so I can test the result?

[dlfivefifty]

I don't know what you mean by non-normalized SemiclassicalJacobi...

can't you just check symtridiagonalize(X) == jacobimatrix(LanczosPolynomial(...))?

[TSGut]

Yes but there would be a simpler way to sanity check the result if we had access to the Tridiagonal instead of SymTridiagonal version. But in the absence of that I'll do it in the way you said.

[dlfivefifty]

I don't know what you mean as you need to choose a normalization for the OPs to be well-defined

[TSGut]

What I mean is that the approach I'm using makes a particular non-orthonormal semiclassical Jacobi basis natural to use for each c. In the hierarchy of this normalization I can now find the c-1 Jacobi matrix if I have the c Jacobi matrix.

Now to convert all of this into the hierarchy we are using and refer to as SemiclassicalJacobi() takes some extra work because it's not as simple as SymTridiagonal, since both of the bases involved have to be adjusted. In the Legendre case this was easy because, well, we know both Legendre() and Normalized(Legendre()), which correspond to these two in the special case c=0. The general case now is as if I only have access to Normalized(Legendre()) but want the recurrence coefficients of Legendre(). Hope that makes sense. It's doable to get it for our convention by making some adjustments and I'm working on it but just takes some extra care to get right.

[TSGut]

Okay, this works. The caveat is that just like in the c=0 case we need to be able to compute the integral of the c basis with respect to the c-1 weight at least for n=1 (if we can do it for any n it would be just as nice as Legendre case but with n=1 we can still make it work via a recurrence argument).

Basically this needs \int_0^1 w_{t,a,b,c-1}(x)*Q_1^{t,(a,b,c)}(x) dx. Note that it's c-1 and c in weight and polynomial. Any good ideas?

[dlfivefifty]

Easy:

A,B,_ = recurrencecoefficients(SemiclassicalJacobi(t,a,b,c)) # Note Q_1 == A[1]x + B[1]
A[1]*sum(SemiclassicalJacobiWeight(t, a+1, b, c-1)) + B[1]*sum(SemiclassicalJacobiWeight(t, a, b, c-1))

[TSGut]

Awesome, that's all we should need. Did you want me to implement it as well or do you want me to just send you a collection of the resulting formulas? I'm happy with either.