JuliaLinearAlgebra / ArnoldiMethod.jl

The Arnoldi Method with Krylov-Schur restart, natively in Julia.
https://julialinearalgebra.github.io/ArnoldiMethod.jl/dev
MIT License
102 stars 18 forks source link

Use a Krylov-Schur implementation to make life easier #82

Closed haampie closed 6 years ago

haampie commented 6 years ago

Suppose that we have an Arnoldi relation

AV = VH + hveₖ' = [V v]B   where B = [H; heₖ']

and some unitary matrix Q, then a similar Krylov-relation would be

AVQ = VHQ + hveₖ'Q

If you update B ← [Q 0; 0 1]'BQ and V ← VQ

AV = [V v]B

But B is no longer Hessenberg. Fortunately it's easy to restore via householder reflections like this:

xxxx    ****    ****    ***x    ****    **xx    ****
xxxx    ****    ****    ***x    ****    **xx    ****
xxxx -> **** -> **** -> ***x -> **** ->  *xx ->  xxx
xxxx    ****    ****      *x      xx      xx      xx
xxxx       *       x       x       x       x       x
  B     B←BW₁   B←W₁B   B←BW₂   B←W₂B   B←BW₃   B←W₃B

The proces above is just B ← [W 0; 0 1]'BW repeatedly. And then in fact the unitary matrix Z = QW₁⋯Wₙ would result in another Arnoldi relation with V←VZ and B←[Z 0; 0 1]'BZ.

Now the interesting bit is that we can pick Q as the Schur matrix in HQ = QR where R is upper triangular. We can freely reorder the diagonal blocks of R s.t. the Ritz values we don't want appear bottom right. Now we can just drop the last columns of the new R and only then restore Hessenberg form.

This seems a lot easier than the convoluted purging strategies etc etc. Worth a shot.

haampie commented 6 years ago

Closed via #81