# Help with restart functionnality in sef-made GMRES solver in python

I am new to this forum and to computational science in general.

I started to learn numerical liner algebra on my own and would like to code a GMRES solver in python (no preconditioner for the time being) with the ability to restart. My end goal is to code a Flexible GMRES [2]. However, I am facing some difficulties in understanding the restart part of GMRES(m).

From [1], the algorithm for the GMRES(m) is:

1. Start: Choose $$x_0$$ and compute $$r_0 = f - Ax_0$$ and $$v_1 = \frac{r_0}{\left \| r_0 \right \|}$$

2. Iterate: for $$j= 1,2,3,...m$$ do: [note: $$m$$ should be the number of restart iterations]

• $$h_{i,j} = (Av_j, v_i), i=1,2,...,j$$,
• $$\hat{v}_{j+1}=Av_j-\sum_{i=1}^j h_{i,j}v_i$$,
• $$h_{j+1,j} =\left \| \hat{v}_{j+1} \right \|$$,
• $$v_{j+1} = \frac{\hat{v}_{j+1}}{h_{j+1,j}}$$
3. Form the approximate solution:

• $$x_m = x_0 + V_m y_m$$, where $$y_m$$ minimizes $$\left \| \beta e_1 - \bar{H}_m y \right \|$$, $$y \in R^m$$
4. Restart:

• Compute $$r_m = f-A x_m$$;
• if satisfied then stop, else compute $$x_0:=x_m$$ and $$v_1:=\frac{r_m}{\left \| r_m \right \|}$$ and go to point 2.

Note: Point 2 is known as Arnoldi iterations.

I have the following code for a GMRES solver without restart.

def GMRES (A, b, max_iterations, x0=None, tol = 1e-8):
"""
Parameters
----------
A : np.array
(m x m) numpy array. Coefficient matrix
b : np.array
(m x 1) numpy array. Source term vector
max_iterations : int
maxium number of iterations to perform
x0 : np.array
(m x 1) numpy array. Initial guesss for the iterative procedure
tol : double
tolerance. The default is 1e-8.

Returns
-------
x : np.array
(m x 1) numpy array.
residualVector : list
values of residual throughout the iterative procedure

"""
# Algorithm status variable
status = -1

# Size of solution
m = len(b)

if np.all(x0==None):
x0 = np.zeros((m, 1))

# Compute the residual
r0 = b - (A @ x0)

# Compute beta
beta = norm(r0)

# List to store the residual values
residualVector = [beta]

for i in np.arange(1, (max_iterations + 1)):

# Compute Arnoldi iterations
Q, h = arnoldi_iter(A, r0, i)

# Solve least square problem
y = np.linalg.lstsq(h, beta*e1(i), rcond=None)[0]

# Computes the error of the iteration
rn = error (r0, h, y, i)

# Appends the value of error to the list
residualVector.append(rn)

# Check convergence
if (rn < tol):
# Solver terminated successfully
status = 0
break

if (i == max_iterations):
# Solver did not converge
status = max_iterations

# Form solution
x = x0 + np.matmul(Q[:,:-1], y )

return x, residualVector, status


The code for the Arnoldi Iterations is:

def arnoldi_iter(A, b, N):
"""
Computes a basis of the (n+1)-Krylov Subspace of A

Krylov subspace A = <b, bA, A²b, A³b, ... A^nb>

Parameters:
A: m x m Matrix (numpy array)
b: m x 1 vector (numpy array)
N: Number of iterations

returns:
Q: The columns of an orthogonal basis of A's Krylov subspace. m x (n+1) Matrix
h: A on basis Q, It is an upper-Hessenber Matrix of size (n+1) x n
"""

# Set m size of A
m = len(b)

# Initializes the Matrices to the required dimensions
Q = np.zeros((m, N+1))

h = np.zeros((N+1, N))

# Computes the first q-vector
q = b/norm(b)

# Sets the first column of Q to the transpose of the q-vector
Q[:, [0]] = q

# Loop
for n in np.arange(N):

# Calculate v-vector (New candidate vector)
v = A @ q

for j in np.arange(n + 1):
qj= Q[:, j].reshape(m, 1) # make vector (m, 1) 2-D array instead of (m,). Causes problems with subtraction when using (m,1) - (m,) vectors

h[j, n] = ( qj.conj() ).T @ v

v = v - ( h[j, n] * qj)

h[n + 1, n] = norm(v)

tol = 1e-15 # avoid division by 0

if (h[n+1, n] > tol):
q = v/h[n + 1, n]

Q[:, [n + 1]] = q

else:
return Q,h
return Q,h


The function for computing the error is:

def error (b, h, y, i):
"""
Parameters
----------
b : np.array
residual vector.
h : np.array
Hessenberg matrix at the i'th iteration
y : np.array
Solution value at i'th iteration.
i : int
iteration number.

Returns
-------
err :
DESCRIPTION.

"""
# Compute norm of the b vector
beta = norm(b)

# Compute the error

err = norm(np.matmul(h, y) - ( beta*e1(i) ))

return err


And finally, the function to return the first column vector of the identity matrix of size n+1 is:

def e1(n):
"""
Compute the first column vector of the identity matrix with size (n+1) x 1

Parameters
----------
n : iteration counter

Returns
-------
tmp : np.array
first column vector of the identity matrix with size (n+1) x 1

"""
tmp = np.zeros((n + 1 , 1))
tmp[0] = 1
return tmp


I have been trying out the above algorithm as:

N = 30

b = np.random.normal(1, scale=5, size=(N, 1))

# Create non-singular matrix. Singular matrices are not inverteble ...
matrixQualityCheck = False

while(not(matrixQualityCheck)):
A = np.random.normal(1,scale=5, size=(N, N))
determinant = np.linalg.det(A)

if(determinant != 0):
matrixQualityCheck=True

# Solve system exactly
x = np.linalg.solve(A,b)

# Solve with GMRES
x2, RESVEC, status = GMRES(A, b, N, tol=1e-12)

# Compute the norm of the difference between exact and GMRES solution
normDiff = np.abs(norm(x - x2))

if (  normDiff < 1e-10  ):
print("Same result")
else:
print("solvers are different")
print(normDiff)


The above code appears to work correctly.

Now, for the restart part of the code I have:

def restart_GMRES(A, b, max_iterations, restart_iterations, x0=None, tol = 1e-8):

for i in np.arange(restart_iterations):
if (i == 0):
# First iteration computes normal GMRES with x0 by the user
x, RESVEC, status = GMRES(A, b, max_iterations, x0=x0, tol=tol)
else:
x, RESVEC, status = GMRES(A, b, max_iterations, x0=x, tol=tol)

if (status == 0):
# solver finished with success
return x, RESVEC, [status, i]
return x, RESVEC, [status, i]


Since I have a status variable inside the code that is telling me that the solution convergence to the correct value I am using it to exist the restart algorithm. However, the results do not appear to be Ok and I do not know why...

I would really appreciate some help!

Sorry for the long post.

References

[1] - Saad, Y., & Schultz, M. H. (1986). GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7(3), 856–869. https://doi.org/10.1137/0907058

[2] - SAAD, Y. (1993). A FLEXIBLE INNER-OUTER PRECONDITIONED GMRES ALGORITHM. SIAM Journal on Scientific Computing, 14(2), 461–469.