Orthogonal connection matrix

[dlfivefifty]

There are issues with the current connection coefficient matrix introducing bad conditioning:

V = SymTriOperator([zeros(5); -ones(6); zeros(10); -ones(7)],zeros(4))
J = -Δ + V
Λ, U = eig(J)
norm(Matrix(U[1:100,1:100])) # 2.3265722544996794e28

However, from the spectral measure we know we can construct an orthogonal connection coefficients Q:

Δ = freejacobioperator()
K = SymTriOperator(0.01ones(3),zeros(4))
J = Δ + K
μ = spectralmeasure(J)


Λ, U = eig(J)

M = Multiplication(1./sqrt(μ.q), rangespace(U))
Q = M*U
norm((Q'*Q - I)[1:10,1:10])  # ≈ 0

# show it's still a conversion
C = connectioncoeffsoperator(J)
X = Multiplication(Fun(), Ultraspherical(1)) # Δ with right spaces
(X*Q - Q*J)[1:10,1:10] |> norm   # ≈ 0

By definition, Q is well conditioned. Can we recover Q in a reasonable form without calculating the badly conditioned U?

[dlfivefifty]

FYI we know it can't be Hessenberg since the lower bandwidth is equal to the lowerbandwidth of M

Note also that while M is not technically banded, it of course can be approximated by a polynomial. But this isn't great since we want complexity given by the perturbation size, nothing else. However, that said we know that f^(-1/2) is a solution to a banded system

f_nsqrt = [ldirichlet() ; f*D + 1/2*I] \ [sqrt(f(-1)); 0]

Since in this case f = μ.q is a polynomial of degree dictated by the perturbation, we know that Q can be encoded by a fixed number of entries.