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 attraction compared with unrestarted GMRES is that Chebyshev iteration has bounded memory cost while delivering a minimax polynomial over the chosen spectral enclosure. It appears, however, that when you Arnoldi factorize a shifted-skew system \( A \) as \( AV_k = V_{k+1} \bar H_k \), the projected matrix is tridiagonal. This means it should be possible to derive a GMRES-like minimal-residual algorithm 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.
The tridiagonal structure follows directly in exact arithmetic. Write \(A=\sigma I+S\), where \(S^T=-S\), and let \(H_k=V_k^TAV_k\). Then \(H_k-\sigma I=V_k^TSV_k\) is skew-symmetric. Arnoldi also makes \(H_k\) upper Hessenberg. A matrix that is both shifted skew-symmetric and upper Hessenberg cannot have entries more than one position away from the diagonal: an entry below the first subdiagonal is zero by Hessenberg structure, and skew symmetry forces the mirrored entry above the first superdiagonal to be zero as well. Thus \(H_k\) is tridiagonal, with diagonal \(\sigma\) and opposite-signed off-diagonal entries in the real case.
In floating-point arithmetic, loss of orthogonality can introduce small entries outside this band, so the exact three-term structure may slowly degrade.
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)
if H[j+1,j] == 0.0:
return H[:j+2,:j+1],V[:,:j+1]
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 matrix
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")
The three-term recurrence is therefore structural rather than accidental. I will investigate whether a stable minimal-residual implementation improves on the Chebyshev approach in a later post.