3
$\begingroup$

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.

$\endgroup$

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.