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):
    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.

    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
        # Check convergence
        if (rn < tol):
            # Solver terminated successfully
            status = 0
    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>
        A: m x m Matrix (numpy array)
        b: m x 1 vector (numpy array)
        N: Number of iterations
        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
            return Q,h
    return Q,h

The function for computing the error is:

def error (b, h, y, i):
    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.

    err : 

    # 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

    n : iteration counter
    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

    A = np.random.normal(1,scale=5, size=(N, N))
    determinant = np.linalg.det(A)
    if(determinant != 0):

# 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")
    print("solvers are different")

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


[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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.