Harmonic Ritz values visualized

I thought it would be nice to clarify the previous post with some plots. Take a simple Poisson matrix $A$ of size $n \times n$ with entries $\begin{bmatrix}-1 & 2 & -1\end{bmatrix}$ on the diagonals. We know its eigenvalues are all concentrated in the interval $[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 $v_0$:

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

The well-known Arnoldi decomposition $AV_k = V_{k+1}H_{k+1,k}$ gives us an orthonormal matrix $V_k$ whose $k$ columns span the search space, and a Hessenberg matrix of size $(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
end

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

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])
    end
    H[k + 1, k] = norm(@view(V[:, k + 1]))
    V[:, k + 1] /= H[k + 1, k]
end

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)
    end

    V, H
end

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

$V_k^TAV_ky = H_{k,k}y = \theta y.$

Here the square matrix $H_{k,k}$ is $H_{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))
    end
    p
end

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 $\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 - \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 - \tau I)V_k = V_{k+1}H_{k+1,k}^\tau$ for $H_{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

$(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τ = H - UniformScaling(τ)
    A = Hτ' * Hτ
    A[end, end] += conj(f) * f
    eigvals(A, Hτ') + τ
end

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))
    end
    p
end

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.