Package should be O(n) not O(n^2)

[dlfivefifty]

julia> n = 100; P = SemiclassicalJacobi(2, 1/2, 0, 1/2); @time P[0.1,n]
  0.043851 seconds (225.48 k allocations: 11.675 MiB)
-1.9278461993824771

julia> n = 200; P = SemiclassicalJacobi(2, 1/2, 0, 1/2); @time P[0.1,n]
  0.173099 seconds (771.69 k allocations: 50.952 MiB, 11.95% gc time)
-0.42808375446709973

julia> n = 400; P = SemiclassicalJacobi(2, 1/2, 0, 1/2); @time P[0.1,n]
  0.641172 seconds (2.83 M allocations: 273.142 MiB, 6.42% gc time)
0.6259343524965616

julia> n = 800; P = SemiclassicalJacobi(2, 1/2, 0, 1/2); @time P[0.1,n]
  2.965002 seconds (11.61 M allocations: 1.691 GiB, 9.94% gc time)
0.9986605877063376

Looking at the profile there's a lot of time in lanczos! but this should just be using Cholesky!

[dlfivefifty]

@DanielVandH I think this is what we were observing.

[TSGut]

Unless something was changed since I last worked on this, only integer parameter bases use decomposition methods since it can then cleanly modify from a classical basis. Modification by integers from non integer cases works as well.

It wouldn't be difficult to make it work for rational parameter starting points, probably just change a few lines that invoke lanczos in the if cases.

[MikaelSlevinsky]

Hey guys! Since it's so frustrating going through other peoples' code, I thought I'd show you just how fast life can be with the new GramMatrix architecture described in https://arxiv.org/pdf/2412.17663. The complexities here are optimal, O(bn), where b is the bandwidth of the GramMatrix, which corresponds to the number of Jacobi expansion coefficients of the semi-classical part of the weight. The code below uses the weight w(x) = (1-x)^alpha (1+x)^beta (t+x)^gamma, so an affine transformation should bring it to your package definition.

julia> using ClassicalOrthogonalPolynomials, FastTransforms, LazyArrays, LinearAlgebra

julia> function semiclassicaljacobimoments(::Type{T}, N::Int, α, β, γ, t) where T
           P = Jacobi(α, β)
           w = expand(P, x -> (t+x)^γ)
           MP = ClassicalOrthogonalPolynomials.weightedgrammatrix(P)
           wc = w.args[2].args[1]
           ν = length(wc)
           μ = PaddedVector(MP[1:ν, 1:ν]*wc, N) # See first equation, page 7
           return μ
       end
semiclassicaljacobimoments (generic function with 1 method)

julia> n = 10000
10000

julia> N = 2n-1 # we need 2n-1 moments
19999

julia> α, β, γ = 0.123, 0.456, 0.789
(0.123, 0.456, 0.789)

julia> t = 2.0
2.0

julia> @time μ = semiclassicaljacobimoments(Float64, N, α, β, γ, t);
  0.000846 seconds (140 allocations: 101.047 KiB)

julia> @time XP = jacobimatrix(Jacobi(α, β))[1:N, 1:N];
  0.000366 seconds (7 allocations: 469.125 KiB)

julia> @time W = GramMatrix(μ, XP);
  0.002420 seconds (10 allocations: 5.722 MiB)

julia> @time R = cholesky(W).U;

       # P = Q * R
  0.004740 seconds (34 allocations: 4.807 MiB)

julia> P = Jacobi(α, β)
Jacobi(0.123, 0.456)

julia> @time Qatx = R'\P[0.123, 1:n] # evaluates all the semi-classicals at the point.
  0.000296 seconds (18 allocations: 157.234 KiB)
10000-element Vector{Float64}:
  0.5590731653962269
 -0.12237755383595854
 -0.5790610626762823
 -0.012140901662257576
  0.5775595035463799
  0.15172097088437247
 -0.5411029273448453
 -0.2836817267531338
  0.4719627859801919
  0.39919522860119333
  ⋮
 -0.08480108156150551
  0.5632884860570925
  0.22337004849171394
 -0.5083394543818879
 -0.34842155369248945
  0.4226277525691457
  0.45238798034486255
 -0.3113403099177292
 -0.5289776962314543

julia> 

[dlfivefifty]

So I guess we'd need to add support for ∞ Jacobi matrices to GramMatrix to support adaptivity?

[dlfivefifty]

Though it's so fast it can probably just be used on a bigger truncation each time and still be way faster

[MikaelSlevinsky]

When the moments are finite (a PaddedVector), the algorithm can be written adaptively (https://github.com/JuliaApproximation/FastTransforms.jl/blob/a40b3dc339f8cceaff5a956814079c44349c8124/src/GramMatrix.jl#L185-L217 no need to compute the skew generators). The savings would be O(bn) on O(bn), but this would have to be weighed against the resizing/re-allocating strategy. (I wouldn't bother.)

[dlfivefifty]

So let's say t is close to 1, i.e., t = 1+1/m.

Then probably the bandwidth is $O(m)$? And the complexity is $O(nm)$?

For a hierarchy this will be a one-off.

That's probably about as good as we need!

[dlfivefifty]

What is fastcholesky! for?

[MikaelSlevinsky]

That's the code that makes Cholesky ... fast

[TSGut]

Why does it live in FastTransforms then?

[MikaelSlevinsky]

Because it's specific to factoring a matrix that satisfies a certain displacement equation, such as the Gram matrix.

[dlfivefifty]

Probably the function name grammatrix_cholesky! would be clearer...

[dlfivefifty]

I think your documentation should be using unicode not latex as it's currently unreadable:

help?> GramMatrix
search: GramMatrix BitMatrix Matrix AbstractMatrix DenseMatrix StridedMatrix

  GramMatrix(W::AbstractMatrix, X::AbstractMatrix)

  Construct a symmetric positive-definite Gram matrix with data stored in W.
  Given a family of orthogonal polynomials {\bf P}(x) = \{p_0(x),
  p_1(x),\ldots\} and a continuous inner product \langle f, g\rangle, the Gram
  matrix is defined by:

  W_{i,j} = \langle p_{i-1}, p_{j-1}\rangle.

  Moreover, given X, the transposed Jacobi matrix that satisfies x {\bf P}(x)
  = {\bf P}(x) X, the Gram matrix satisfies the skew-symmetric rank-2
  displacement equation (X = X_{1:n, 1:n}):

  X^\top W - WX = GJG^\top,

  where J = \begin{pmatrix} 0 & 1\\ -1 & 0\end{pmatrix} and where:

  G_{:, 1} = e_n,\quad{\rm and}\quad G_{:, 2} = W_{n-1, :}X_{n-1, n} - X^\top W_{:, n}.

  Fast (O(n^2)) Cholesky factorization of the Gram matrix returns the
  connection coefficients between {\bf P}(x) and the polynomials {\bf Q}(x)
  orthogonal in the modified inner product, {\bf P}(x) = {\bf Q}(x) R.

  ────────────────────────────────────────────────────────────────────────────

  GramMatrix(μ::AbstractVector, X::AbstractMatrix)

  Construct a GramMatrix from modified orthogonal polynomial moments and the
  multiplication operator. In the standard (classical) normalization, p_0(x) =
  1, so that the moments \mu_n = \langle p_{n-1}, 1\rangle are in fact the
  first column of the Gram matrix. The recurrence is built from X^\top W = WX.

[dlfivefifty]

JuliaApproximation/FastTransforms.jl#258

[dlfivefifty]

I'm not sure I understand the motivation for the GramMatrix type since it appears to be identical to the stored matrix W.

I guess it's to flag that W has the extra structure of a gram matrix for the cholesky? I.e. it's morally equivalent to wrapping a matrix whose entries are symmetric in the type Symmetric?

I'm wondering if it makes more sense to leave GramMatrix lazy and only construct W when cholesky is called.