Harmonic Ritz values visualized

A hopefully insightful visualizion of harmonic Ritz values versus ordinary Ritz values.

June 7, 2017 - 11 minute read -
gsoc julia

I thought it would be nice to clarify the previous post with some plots. Take a simple Poisson matrix AA of size n×nn \times n with entries [121]\begin{bmatrix}-1 & 2 & -1\end{bmatrix} on the diagonals. We know its eigenvalues are all concentrated in the interval [0,4][0, 4] on the real line.

As a preliminary note: I’m not exploiting the symmetry of the matrix in my code. I want to test it on general matrices as well.

At any rate, in Julia we generate this Poisson matrix as

my_matrix(n) = SymTridiagonal(fill(2.0, n), fill(-1.0, n - 1))

Sometimes Jacobi-Davidson’s search space is initialized with the search space of Arnoldi’s method, namely a Krylov subspace with a random initial vector v0v_0:

V=Kk(A,v0)=span{v0,Av0,,Ak1v0}. \mathcal{V} = \mathcal{K}_k(A, v_0) = span\{v_0, Av_0, \cdots, A^{k-1}v_0\}.

The well-known Arnoldi decomposition AVk=Vk+1Hk+1,kAV_k = V_{k+1}H_{k+1,k} gives us an orthonormal matrix VkV_k whose kk columns span the search space, and a Hessenberg matrix of size (k+1)×k(k + 1) \times k . In Julia we can generate these using

function init(n, dimension)
    V = zeros(n, dimension + 1)
    H = zeros(dimension + 1, dimension)
    v0 = rand(n)
    V[:, 1] = v0 / norm(v0)

    V, H

function expand!(A, V, k)
    V[:, k + 1] = A * @view(V[:, k])

function orthogonalize!(V, H, k)
    for i = 1 : k
        H[i, k] = dot(@view(V[:, i]), @view(V[:, k + 1]))
        V[:, k + 1] -= H[i, k] * @view(V[:, i])
    H[k + 1, k] = norm(@view(V[:, k + 1]))
    V[:, k + 1] /= H[k + 1, k]

function arnoldi(A, m)
    n = size(A, 1)
    V, H = init(n, m)

    for k = 1 : m
        expand!(A, V, k)
        orthogonalize!(V, H, k)

    V, H

Now the ordinary Ritz values θ\theta are obtained by solving the small eigenvalue problem

VkTAVky=Hk,ky=θy. V_k^TAV_ky = H_{k,k}y = \theta y.

Here the square matrix Hk,kH_{k,k} is Hk+1,kH_{k+1,k} with its last row removed. We can compute and plot the Ritz values for each iteration as follows:

using Plots

ritz(H) = eigvals(H)

function ritz_example(; n = 60, m = 30)
    A = my_matrix(n)
    V, H = arnoldi(A, m)
    p = plot(;legend = :none)

    for k = 1 : m
        ritz_values = ritz(@view(H[1 : k, 1 : k]))
        scatter!(real(ritz_values), k * ones(ritz_values), marker = (:x, :green))

What we get is something like


On the vertical axis we have the iteration (or size of our search space), and horizontally we have the real line. The black stars are the actual eigenvalues, while the green crosses the Ritz values. Now if we were to target interior eigenvalues near τ=3\tau = 3, we can clearly see that many Ritz values start out near this target, but actually converge away from it. This is undesirable for restarts.

Enter harmonic Ritz values. As shown in the previous post harmonic Ritz values satisfy the Petrov-Galerkin condition

(AτI)zθz(AτI)Vk for z=Vky. (A - \tau I)z - \theta z \perp (A - \tau I)V_k \text{ for } z = V_k y.

Now it’s not hard to verify we can write (AτI)Vk=Vk+1Hk+1,kτ(A - \tau I)V_k = V_{k+1}H_{k+1,k}^\tau for Hk+1,kτ:=Hk+1,kτIH_{k+1,k}^\tau := H_{k+1,k} - \tau I. Note that the Hessenberg matrix is not square, but by abuse of notation I mean that τ\tau must be substracted from the diagonal of the Hessenberg matrix. Some more work shows that the Petrov-Galerkin condition is equivalent to solving the generalized eigenvalue problem

(Hk+1,kτ)Hk+1,kτy=θ(Hk,kτ)y. (H^\tau_{k+1,k})^* H^\tau_{k+1,k} y = \theta (H^\tau_{k,k})^* y.

This can be implemented using:

function harmonic_ritz(H, f, τ)
     = H - UniformScaling(τ)
    A = ' * 
    A[end, end] += conj(f) * f
    eigvals(A, ') + τ

function harmonic_ritz_example(; n = 60, m = 30)
    A = my_matrix(n)
    V, H = arnoldi(A, m)
    p = plot(;legend = :none)

    for k = 1 : m
        Hk = @view(H[1 : k, 1 : k]
        harmonic = harmonic_ritz(Hk, H[k + 1, k], τ))
        scatter!(real(harmonic), k * ones(harmonic), marker = (:x, :green))

And produces harmonic Ritz values like this


The point is that harmonic Ritz values come close to our target τ\tau only when they converge to an eigenvalue near τ\tau.

For restarts this means that retaining harmonic Ritz vectors with harmonic Ritz values close to τ\tau will be a good strategy. In sharp contrast with ordinary Ritz values, where lots of Ritz values near our target are useless to retain, since they are on their way converging to exterior eigenvalues.