In my previous post I presented a way to convert a nonsymmetric matrix into a shifted-skew matrix, and then exploited the simple structure of the spectrum of this matrix to derive Chebyshev iteration parameters. The purpose of this instead of say GMRES is that Chebyshev iteration has bounded memory cost even while delivering optimal convergence. It appears however that when you Arnoldi factorize a shifted-skew system \( A \) as \( AV = VH \) then the matrix \( H \) is tridiagonal. This means it should be possible to derive GMRES-like algorithm that achieves optimal convergence rates with a three-term-recurrence similar to Chebyshev iteration. This may be preferable to Chebyshev iteration because it doesn’t require estimating spectral radius and may be more robust.
I illustrate below
Arnoldi factorization arnoldi.py
import numpy as np
def arnoldi(A,v,k):
norm=np.linalg.norm
dot=np.dot
m,_=A.shape
V=np.zeros((m,k+1))
H=np.zeros((k+1,k))
V[:,0]=v/norm(v)
for j in range(0,k):
w=A@V[:,j]
for i in range(0,j+1):
H[i,j]=dot(w,V[:,i])
w=w-H[i,j]*V[:,i]
H[j+1,j]=norm(w)
V[:,j+1]=w/H[j+1,j]
return H,V
Numerical experiment
import numpy as np
import scipy.sparse as sp
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import arnoldi
seed=23949823
m = 1024
rng = np.random.default_rng(seed)
#Make a random sparse banded matrixk
A = sp.diags( [rng.uniform(-1,1,m) for _ in range(0,7)], [-100,-20,-1,0,1,20,100], shape = (m,m))
I = sp.eye(A.shape[0])
#Make into shifted-skew sparse matrix
A = I + 0.5*(A - A.T)
#Compute Arnoldi factorization
#uses random input vector v
#but for a linear solver v would be normalized initial residual
v = rng.uniform(-1,1,m)
H,V = arnoldi.arnoldi(A,v,10)
#Output to see tridiagonal structure of H
plt.imshow(H)
plt.colorbar()
plt.savefig("H.svg")
I will try to derive this apparent three-term recurrence and see if it improves on the Chebyshev approach in a later post