Lanczos via Givens

[dlfivefifty]

I'm trying to rewrite lanczos (i.e. Stieltjes) in ∞-arrays which made me realise that it essentially boils down to tridiagonalising W, a symmetric positive definite matrix such that

w .* P = P * W

where P are normalized Legendre polynomials. In practice W is banded.

@MikaelSlevinsky I assume https://github.com/JuliaMatrices/BandedMatrices.jl/blob/master/src/symbanded/tridiagonalize.jl will work here? We'll just have to modify it to make it adaptive.

But you aren't storing the Q? Should I instead call sbtrd!? Will need to think how to do this adaptively...

[dlfivefifty]

PS Conjecture: Steiltjes sucks, Given's tridiagonalisation doesn't. Do you know if anyone has written this up?

On second thought one doesn't actually need the Q, which I believe encodes the conversion matrix from legendre. But it would be cool to have.

[MikaelSlevinsky]

I think you want the algorithms in section 6 https://www.sciencedirect.com/science/article/pii/S0021999120301571

Prototypes exist in a draft PR somewhere in JuliaApproximation.

[MikaelSlevinsky]

Oh I guess you only want tridiagonalization, not full spectral decomposition? Probably the idea is nearly the same.

[dlfivefifty]

Yes I think you are right. It shouldn't be necessary to perturb M in this case, perhaps they are automatically exactly 0?

Probably just need to try...

[dlfivefifty]

Is there a reason we aren't dealing with hessenberg? Is it just because we don't need to store the Q?

Can we store the Q in tridiagonalize! as a product of Given's rotations in optimal complexity? (Ref https://github.com/JuliaMatrices/InfiniteLinearAlgebra.jl/blob/409083ce26c85f20091f46bebed414a48647ab2d/src/banded/infqltoeplitz.jl#L86)

[dlfivefifty]

I'm picturing the implementation of Stieltjes looking something like this:

P = Normalized(Legendre())
x = axes(P,1)
w = exp.(x)
W = P' * (w .* P)
J = hessenberg(W).H
Q = RecurrenceOrthogonalPolynomial(J) # OPs wrt exp(x)

[dlfivefifty]

Even better, to make it 100% linear algebra:

P = qr(Legendre()).Q
Q = eigen(J).Q

[dlfivefifty]

Every time I think about this I realise I'm very confused. These are facts:

1. Given a sym. def. matrix A tridiagonalisation gives a construction for find an orthogonal basis w.r.t. A:

julia> A = Symmetric(randn(5,5)) + 10I;

julia> V,J = hessenberg(A);

julia> λ,W = eigen(J); Λ = Diagonal(λ);

julia> Q = V * W * Λ^(-1/2);

julia> Q'*A*Q
5×5 Array{Float64,2}:
  1.0           1.35955e-16  -9.218e-17     1.14058e-16  -2.49781e-17
  2.55151e-17   1.0          -1.15394e-17  -1.30977e-16  -7.15788e-17
 -7.8141e-17    6.46896e-17   1.0           4.57227e-16   3.5929e-16
  2.05968e-17  -3.4387e-17    4.28762e-16   1.0           8.70383e-17
 -1.72908e-16  -8.55082e-17   4.35528e-16   1.84408e-16   1.0

2. Lanczos also gives a construction for constructing OPs with respect to a weight w, which we can think of as a symmetric matrix W in the Legendre basis.
3. This suggest using (1) on W is another way of constructing OPs. I don't see how to show this though.
4. What's weird about this is (1) deals with the spectrum of A, so in the case of the weight this is actually the image of the weight, not the support. So the J is not the OPs since the support is wrong!
5. Actually in ∞-dimensions this touches on an issue that @marcusdavidwebb is familiar with: tridiagonalization is not necessarily unique.

So what the hell is the ∞-dimensional version of (1): tridiagonalisation followed by eigencomputation, and how does it relate to Stieltjes??

[dlfivefifty]

Hmm, we want the columns of Q to be polynomials in A, which they won't necessarily be in the above: tridiagonalization is only unique up to first entries...

[dlfivefifty]

OK we need to be careful about inner products defined by a weight W and operators A. Let's start with Householder for upper triangularising w.r.t. a different inner product, <x, y> := x'Wy where W is spd. Then the Householder reflector w.r.t. the inner product is

Q = I - 2 v v' W
Q' = I - 2 W v v'
Q† = inv(W) Q' W = Q

where Q† is the adjoint w.r.t. the inner product. Q is orthogonal because it preserves norms, that is <Qx, Qx> = <x,x> = x'Wx, which follows from the following equivalent relationships:

Q'WQ = W
Q†Q = I

Now we need to upper-triangular an operator A w.r.t. the inner product W, column by column, starting with the first column. That is, we want to choose v in the Householder reflection so that:

<e_k, QAe_1> = e_k' W  Q Ae_1 = 0

for k > 0. Let x = A*e_1 be the first column and the role of e_1 in standard Householder will be played by r = W \ e_1. Thus we choose:

e_1 = [1; zeros(size(A,1)-1)]
α = -sqrt(x'W*x)
r = W\e_1
u = x - α*r
v = u/sqrt(u'W*u)

So that

Qx = x - 2 v v' W x = α * r

For the next column I suppose we need r = W\ e_2 and so on... in other words we need access to inv(W). But how does that make sense if the weight vanished? (I suppose a vanishing weight is only positive semi-definite...)

@ajt60gaibb Have you thought about Householder with other inner products / what it means to be upper triangular?

[dlfivefifty]

Hmm the more I understand the more I'm confused. I think it's more of a similarity transform. That is if Q = P * R where P = Normalized(Legendre()) and R is the upper triangular converison matrix, Q'wQ = I, we have R'WR = I. That is, R is orthogonal from $\ell^2$ to $\ell_W^2$: it is not orthogonal w.r.t. <f,g> = f'W*g. We can define U = sqrt(W) * R (that is, consider sqrt(w) * Q = P * sqrt(W) * R) so that U'U = I.

If Δ is the Jacobi operator for P, finding the Jacobi operator J for Q consists of

 J = Q' w x Q = R' W Δ R = U' Δ U

That is, conjugating by the upper triangular matrix R tridiagonalises the banded operator W * Δ, while conjugating by the orthogonal (w.r.t. $\ell^2$) operator U maps between the two Jacobi operators.

I don't see how orthogonality w.r.t. W comes up naturally... U could be representable by a sequence of Householder reflections and if we approximate sqrt(w) by a polynomial then sqrt(W) is banded and U is banded below... as used in the paper with @marcusdavidwebb this fits in nicely with the Jacobi operator J being a finite rank perturbation of Δ. w = (1-x)^2 is a reasonable test problem.

[TSGut]

This is basically resolved via https://arxiv.org/abs/2302.08448 now, right?