Algorithms: Tridiagonal Hessenberg form using Krylov subspaces

In our previous discussion of Krylov subspaces (pick a vector $\bm{v}$ and build $\bm{v}_j=M^j\bm{v}$ for $j=1,2,\cdots, dim(M)-1$ ). We saw that the Ritz eigenvalues (the extreme eigenvalues) of $M$ can be constructed using a subset of the first $m$ of these vectors, the approximation improving is precision as $m$ increases.

A reminder of some of the terminology: $H$ is lower-Hessianberg if $H_{ij}=0$ for all $j>i+1$ , and is upper-Hessianberg if $H_{ij}=0$ for all $i>j+1$ . Arnoldi iteration uses Kyrlov subspaces to transform symmetric matrix $M$ into an upper-Hessianberg form $H$ .

Let $\{\mathbf{v}_1, \, M\cdot\mathbf{v}_1, \, \cdots, \, M^m\cdot\mathbf{v}_0\}$ be subjected to Gram-Schmidt orthogonalization to orthogonal set $\{\mathbf{q}_1, \, \cdots, \, \mathbf{q}_m\}$ , which is a basis for the Krylov subspace $K_m(M,\mathbf{v}_1)$ . Because $\mathbf{q}_i\cdot \sum_{k=1}^{j+1} a_k \, \mathbf{q}_k=0$ for $i>j+1$ we have conjugacy of the $\bm{q}_j$ vectors with respect to $M$ reminiscent of how normal mode vectors are conjugate with respect to the kinetic energy matrix in Physics 711

$H_{i,j}=H_{j,i}=0=\mathbf{q}_j^T\cdot M\cdot \mathbf{q}_i$

because $M\cdot \mathbf{q}_j=\sum_{k=1}^{j+1} a_k \, \mathbf{q}_k$ for some coefficients $a_k$ . Because a Hermitian matrix has such a simple Hessianberg form, Arnoldi iteration iteration for a Hermitian $M$ simplifies dramatically, because at step $m$

$M\cdot (\mathbf{q}_1, \, \mathbf{q}_2, \, \cdots, \, \mathbf{q}_m)=(\mathbf{q}_1, \, \mathbf{q}_2, \, \cdots, \, \mathbf{q}_m, \, \mathbf{q}_{m+1})\cdot \left(\begin{array}{ccccc} \alpha_1 & \beta_1 & 0 & \cdots & 0\\ \beta_1 & \alpha_2 & \beta_2 & 0 & \cdots\\ \vdots & \vdots & \vdots &\ddots & \vdots\\ 0 & \ddots & \beta_{m-2} & \alpha_{m-1} & \beta_{m-1}\\ 0 &\cdots & 0 & \beta_{m-1} & \alpha_m\\ 0 & \cdots & 0 & 0 & \beta_{m}\end{array}\right)$

the last step having only three terms

$M\cdot\mathbf{q}_m=\beta_{m-1} \, \mathbf{q}_{m-1}+\alpha_m\mathbf{q}_m+\beta_m\, \mathbf{q}_{m+1}$

$\mathbf{q}_{m+1}=(M\cdot \mathbf{q}_m-\beta_{m-1}\mathbf{q}_{m-1}-\alpha_m \mathbf{q}_m)/\beta_m$

which is a rather significant shortening of the inner loop!

So what can we use this for? Given $n\times n$ symmetric $M$ we can construct a Hessianberg $m\times m$ with $m<n$ whose extreme eigenvalues approximate those of $M$ to very good precision. We can also completely “Hessenbergerize” $M$ as a prelude to full diagonalization.

Let’s write the code in julia

hess=function(M,jmax)
n=size(M)[1]
b=rand(n,1)
r=b
v=zeros(n,jmax);alpha=zeros(jmax);beta=zeros(jmax);
beta[1]=sqrt(sum(r .* r))
for j in 1:jmax
  v[:,j]=r/beta[j]
  r=A*v[:,j]
  if j>1
    r=r-v[:,j-1]*beta[j]
  end
  alpha[j]=r'*v[:,j]
  r=r-v[:,j]*alpha[j]
  if j<jmax
  beta[j+1]=sqrt(sum(r .* r))
end
end
H=zeros(jmax,jmax)
for j in 1:jmax-1
H[j,j]=alpha[j]
H[j,j+1]=beta[j+1]
H[j+1,j]=beta[j+1]
end
H[jmax,jmax]=alpha[jmax]
return H
end

I will create a random symmetrized $8\times 8$ and approximate its extreme eigenvalues with a Hessenberg $4\times 4$ :

A=rand(8,8)
A=(A+A')/2
using LinearAlgebra
eigvals(A)
8-element Vector{Float64}:
 -0.9978484322893748
 -0.7223421236032608
 -0.311326659137159
 -0.054381916423452685
  0.525474010310458
  0.7610946030980107
  0.8889134731876
  4.0740875056699

H1=hess(A,4);
eigvals(H1)
4-element Vector{Float64}:
 -0.9921448484528467
  0.14096304088606434
  0.7822933006834016
  4.074085415567306

The agreement is excellent. $8\times 8$ is not all that big, but you can see the advantages for very large arrays, you can get very good extremal eigenvalues by Arnoldi Hessenberg reductions using an array half the dimension. Of course you can also go all the way and obtain the full tridiagonal Hessenberg similar to your original matrix

H1=hess(A,8)
8×8 Matrix{Float64}:
 3.18552  1.73103   0.0        0.0        0.0       0.0        0.0        0.0
 1.73103  0.561275  0.73341    0.0        0.0       0.0        0.0        0.0
 0.0      0.73341   0.164923   0.5864     0.0       0.0        0.0        0.0
 0.0      0.0       0.5864    -0.109781   0.306619  0.0        0.0        0.0
 0.0      0.0       0.0        0.306619  -0.415481  0.43705    0.0        0.0
 0.0      0.0       0.0        0.0        0.43705   0.302667   0.458932   0.0
 0.0      0.0       0.0        0.0        0.0       0.458932  -0.0530971  0.181784
 0.0      0.0       0.0        0.0        0.0       0.0        0.181784   0.527644

eigvals(H1)
8-element Vector{Float64}:
 -0.9978484322893754
 -0.7223421236032609
 -0.311326659137159
 -0.054381916423452525
  0.5254740103104582
  0.7610946030980107
  0.8889134731876566
  4.074087505669898

Home 2.0

Algorithms: Tridiagonal Hessenberg form using Krylov subspaces

Related Posts

Algorithms: Tridiagonal Hessenberg form using Krylov subspaces

Published by Jeff on January 8, 2022January 8, 2022

Related Posts

Algorithms: Jacobi matrix diagonalization

Pade’ approximants. Extracting critical points

Very fast recursions with memoization