I present a new-to-me algorithm for solving a sparse system of linear equations where the sparse matrix is positive definite (but not necessarily symmetric). This solver has some interesting properties:
Achieving (1) through (4) doesn’t guarantee a practical solver but they are interesting to me because nonsymmetric solvers frequently require the following for robust and fast convergence:
I take a simple approach by splitting a matrix into its symmetric and anti-symmetric parts. I form the Cholesky factorization of the symmetric part and apply the inverse of its factors to the left and right. I prove that the system preconditioned in this way has a statically knowable spectral enclosure in the complex plane, which means we may choose minimax parameters for a Chebyshev iteration targeting this enclosure. If the enclosure is valid, the resulting polynomial gives a convergence bound. It is optimal for the worst case over that segment, not necessarily for the particular discrete spectrum and right-hand side in hand; GMRES can choose a better residual polynomial for an individual problem. In this way I rely on the pivoting-free nature of the Cholesky factorization to factorize the symmetric part of a matrix, and the optimality properties of Chebyshev iteration to deliver convergence.
I first describe the mathematical derivation of this approach including proof of the spectral properties of the preconditioned operator. I derive the optimal Chebyshev parameters for this preconditioned system. I show numerical results confirming this does work. I then finish with some conclusions about the utility of this approach. I think this method largely has theoretical rather than practical value.
\(\newcommand{\norm}[1]{\left\lVert#1\right\rVert}\)
I wish to solve a sparse linear system of the form
where \(A\) is a nonsymmetric sparse matrix that is positive definite in the quadratic-form sense. For a real matrix, this is equivalent to requiring its symmetric part \( H=\frac{1}{2}(A + A^T) \) to be symmetric positive definite. I also assume \(A\) is real and do not consider the complex number case here.
One can decompose any real matrix into a symmetric and a skew-symmetric matrix as follows:
Here \(H\) and \(S\) denote the symmetric and skew-symmetric parts, respectively.
Thus since \(H\) is positive definite it has a Cholesky Factorization
One can compute this with sparse \(H\) using e.g. CHOLMOD.
If we left and right precondition \(A\) with this factor we get
this transforms \(A\) into a shifted skew form i.e. since \(S\) is skew-symmetric, so is \( L^{-1} S L^{-T} \).
A shifted-skew matrix \(I+S\) is normal. For normal matrices, polynomial residual bounds can be expressed directly in terms of the spectrum, without the additional eigenvector-conditioning and pseudospectral effects that occur for nonnormal matrices. The actual residual still depends on the initial residual’s components in the eigenvectors. We also know the eigenvalues of \(I+S\) take the form \(1 + \alpha i\). Furthermore since the matrix \(A\) is real we know that the eigenvalues come in conjugate pairs, so that if \(1 + \alpha i\) is an eigenvalue of \(I+S\) then so is \(1 - \alpha i\). In particular if we know the largest \( |\alpha| \), then we know the vertical segment containing all eigenvalues of \(I+S\). A reliable way to estimate this endpoint is to estimate the largest eigenvalue of \(-S^2=S^TS\) or the largest singular value of \(S\). Plain power iteration on \(I+S\) can be awkward because the two conjugate endpoint eigenvalues have the same magnitude.
The Chebyshev method for solving a system of linear equations is a Krylov method like GMRES or Conjugate Gradients but instead of building up a basis through an explicit or implicit orthogonalization procedure it uses prior knowledge of the input matrix spectrum. This knowledge generally must define a valid enclosure or else the procedure may converge slowly or even diverge, and since accurately estimating spectral bounds can itself be expensive this limits the usefulness of Chebyshev iteration in general. However in light of above comments on the spectrum of \(I + S\) we can achieve very accurate spectral estimates using simply power iteration.
I give simple Python code below for Chebyshev iteration
def chebyshev(m,A,b,alpha,c,maxiter=100,verbose=True,tol=1e-10):
nu=-alpha/c
beta=0.5*c * (1.0/nu)
gamma=-alpha
x=np.zeros(m)
xm1=np.zeros(m)
xp1=np.zeros(m)
r=b.copy()
rm1=np.zeros(m)
rp1=np.zeros(m)
for n in range(0,maxiter):
if n>=2:
beta=((0.5*c)**2)*(1.0/gamma)
if n>=1:
gamma=-(alpha+beta)
if n==0:
xp1=-(r+alpha*x)/gamma
rp1=(A(r)-alpha*r)/gamma
else:
xp1=-(r+alpha*x+beta*xm1)/gamma
rp1=(A(r)-alpha*r-beta*rm1)/gamma
if verbose:
print(f"iter={n}, residual={np.linalg.norm(rp1)}")
xm1[:]=x
rm1[:]=r
x[:]=xp1
r[:]=rp1
if np.linalg.norm(r)<tol:
break
return x
The meaning of the parameters alpha,c is that the eigenvalues of the input matrix \(A\) must all lie in the ellipse defined by complex numbers alpha,c. Smaller ellipses result in faster Chebyshev iteration. If the ellipse degenerates to the smallest valid segment, the corresponding Chebyshev polynomial is minimax over that enclosure. This is an optimal worst-case polynomial statement, rather than a claim that it matches GMRES for every finite spectrum and starting residual.
Note that for simplicity I kept the parameters alpha,c complex but in the shifted-skew case they are either pure-real or pure-imaginary and this makes it possible to deal with only pure-real arrays (to avoid duplicating memory costs for storing residual and successive iteration vectors).
I give Python code below for solving a system \(Ax=b\) using the above ideas.
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
from sksparse.cholmod import cholesky
def symm_chebyshev(A,b):
m=A.shape[0]
#Split into hermitian and skew-hermitian parts
H=0.5*(A+A.T)
S=0.5*(A-A.T)
#Compute cholesky factorization of hermitian part
L=cholesky(H)
#Reorder original matrices according to fill-reduce ordering from cholesky
P=L.P()
b=b[P]
#A=A[P[:, np.newaxis], P[np.newaxis, :]]
H=H[P[:, np.newaxis], P[np.newaxis, :]]
S=S[P[:, np.newaxis], P[np.newaxis, :]]
#Convenience functions
def invL(x):
return L.solve_L(x,use_LDLt_decomposition=False)
def invLT(x):
return L.solve_Lt(x,use_LDLt_decomposition=False)
#Get maximum eigenvalue for the skew-symmetric matrix inv(L)*S*inv(L^T)
#"neval" below to keep track of how many times we evaluate this
#skew symmetric part as it's expensive
neval=0
def skew_factor(x):
nonlocal neval
neval+=1
y=invLT(x)
return invL(S@y)
es,_=spla.eigs(spla.LinearOperator((m,m),matvec=skew_factor),k=6,which="LM",tol=1e-10)
eig_neval=neval
#Define ellipse parameters for Chebyshev iteration
alpha=1.0+0.0j
print(es)
c=max(np.abs(np.imag(es)))*1.0j
#Left-precondition b
b=invL(b)
#Now solve by chebyshev iteration
neval=0
x=chebyshev.chebyshev_pureimagc(m,lambda x: x + skew_factor(x),b,alpha,c,maxiter=5000,tol=1e-14,verbose=False)
#Apply right-preconditioner
x[P]=invLT(x)
cheb_neval=neval
return x
I did not yet extensively test this algorithm or compare against alternatives but I did confirm it works on some challenging nonsymmetric systems which I generated randomly. Convergence of the Chebyshev method is generally rapid when the skew-symmetric part is not very large compared to the symmetric part - like below:
iter=0, residual=0.0001340813560623113
iter=1, residual=9.845853216344927e-06
iter=2, residual=1.013103495865824e-11
iter=3, residual=2.5020940334707195e-13
iter=4, residual=4.275895432223994e-19
convergence can take a little longer when the skew-symmetric part is large compared to the symmetric part:
iter=0, residual=4.186754323212551
iter=1, residual=210.66587568527942
iter=2, residual=3.474588981192882
iter=3, residual=167.0642639225858
iter=4, residual=3.04771078119865
iter=5, residual=120.54588772351465
iter=6, residual=2.581557296578431
iter=7, residual=82.71211926014004
...
...
iter=164, residual=1.6864217663823996e-13
iter=165, residual=3.096346201216635e-13
iter=166, residual=1.1169909961323464e-13
iter=167, residual=2.0317413138929607e-13
iter=168, residual=7.397292905249339e-14
iter=169, residual=1.3331742365813827e-13
iter=170, residual=4.898175152676594e-14
iter=171, residual=8.747942227362621e-14
iter=172, residual=3.24292721522201e-14
iter=173, residual=5.740166656201045e-14
iter=174, residual=2.146760662263503e-14
iter=175, residual=3.766550578885808e-14
iter=176, residual=1.4209256729454957e-14
iter=177, residual=2.4715157991933956e-14
iter=178, residual=9.403862049980998e-15
This algorithm is interesting because it removes nonnormality from the transformed operator. The Cholesky factorization and two-sided preconditioning by the Cholesky factor transform the system into shifted-skew form even if the underlying matrix \(A\) is highly nonsymmetric and nonnormal. The predictable spectral enclosure of a shifted-skew matrix then gives a low-memory Chebyshev algorithm with an explicit worst-case convergence bound.
I am interested to see if there could be utility in using incomplete cholesky factorizations to construct preconditioners to nonsymmetric systems using the above ideas. I have not tested this idea yet however.