April 16, 2022

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

# Simple numerical example illustrating

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


# Thoughts

I will try to derive this apparent three-term recurrence and see if it improves on the Chebyshev approach in a later post