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