Improved & early locking of converged Ritz vectors

[haampie]

Current implementation

Current implementation of locking is to wait until a subdiagonal element of the Hessenberg matrix becomes small, and then do a local schur decomposition of the left upper block of the Hessenberg matrix.

Problem with this implementation

Very often Ritz vectors converge without an offdiagonal element getting small. Also we know which Ritz vectors have converged at the end of an outer iteration of Arnoldi. We should use that information better!

Workaround

Use this locking technique by Sorensen [1]. Basic idea / my notes:

1. For each converged pre-Ritz pair (y, θ) determine a Householder reflection W such that Wy = e₁τ with |τ| = 1.
   1. W'HWe₁ = e₁θ, so W'HW = [θ a'; 0 B], where B is dense (Hessenberg form is lost).
   2. Further eₖ'W = eₖ' + w' where ‖w‖₂ ≤ √2 |yₖ|.
2. Post-multiply the Arnoldi relation AV = VH + hₖ₊₁,ₖvₖ₊₁eₖ' with this reflector
   1. We get AVW = VWW'HW + hₖ₊₁,ₖvₖ₊₁eₖ'W
   2. Using the above: A(VW) = (VW)[θ a'; 0 B] + hₖ₊₁,ₖvₖ₊₁eₖ' + hₖ₊₁,ₖvₖ₊₁w'
   3. We can drop the term hₖ₊₁,ₖvₖ₊₁w' because it's small according to the convergence criterion.
   4. We have to find series of Householder reflections Z s.t. Z'[θ; a'; 0 B]Z is Hessenberg again.
   5. The new Arnoldi relation is A(VWY) = (VWY)(Y'W'HWY) + hₖ₊₁,ₖvₖ₊₁eₖ'

Now I think that this step should be combined with the implicit shifts. So basically:

1. Make a maxdim Arnoldi relation
2. Compute the eigenvalues of H + eigenvectors.
3. Partition into converged & not converged
4. Create a matrix Q = I
5. For each converged eigenvalue, construct W and Y and update Q ← Q * W * Y, H ← Y'W'HWY
6. Do implicit restart with a subset of unconverged eigenvalues, where Q gets updated with the Given's rotations.
7. Only now compute V ← V * Q.

Gotcha's

Still have to find out how this is supposed to work with real arithmetic and complex conjugate eigenpairs.

---

Actually the above looks a lot like standard bulge chasing, but then with large householder reflectors... There's one "special" reflector W, and then the rest Y is about restoring the Hessenberg structure. The difference being that we place the eigenvalue in entry H[1,1] rather than H[end,end].

[1] Deflation techniques for an implicitly restarted Arnoldi iteration.

[haampie]

I'm still quite sure I could just reorder the Schur form HQ=QR with Given's rotations, rather than computing eigenvectors and applying Householder reflections.