GMRES implementation does not converge for singular Hermitian problems

Question

I've just implemented the GMRES algorithm based on chapter 4 of Fundamentals of Numerical Mathematics for Physicists and Engineers using the problems in Numerical Analysis by Timothy Sauer for validation. Implementation converges for all of the worked examples in the Prof Sauer's text but does not for the singular non-textbook problem.

#!/usr/bin/python3

import math
import copy
import numpy as np

TEXTBOOK_PROBLEM = True

if TEXTBOOK_PROBLEM:
    # From Numerical Analysis by Timothy Sauer
    # https://media.pearsoncmg.com/aw/aw_sauer_num_analysis_3/solutions/sol_4_4.html
    A = np.array([
                    [ 1, 1, 0],
                    [-1, 1, 2],
                    [ 0, 0, 1]
                    ])

    b = np.array(
        [[1],
         [0],
         [0]])
else:
    A = np.array(
    [[  53,   83,  101,  119,  133,  161,    0,    0,    0],
     [  83,  130,  158,  187,  209,  253,    0,    0,    0],
     [ 101,  158,  194,  221,  247,  299,    0,    0,    0],
     [ 119,  187,  221, 1130, 1222, 1464, 1189, 1247, 1363],
     [ 133,  209,  247, 1222, 1322, 1584, 1271, 1333, 1457],
     [ 161,  253,  299, 1464, 1584, 1898, 1517, 1591, 1739],
     [   0,    0,    0, 1189, 1271, 1517, 1681, 1763, 1927],
     [   0,    0,    0, 1247, 1333, 1591, 1763, 1849, 2021],
     [   0,    0,    0, 1363, 1457, 1739, 1927, 2021, 2209]])

    b = np.array(
    [[  650],
     [ 1020],
     [ 1220],
     [ 8142],
     [ 8778],
     [10506],
     [ 9348],
     [ 9804],
     [10716]])

################################################################################
# Uses Given's rotation to convert an upper Hessenberg matrix H to an upper
# triangular matrix which is trivial to solve.
def solve(H):
    numRows, numCols = H.shape
    b = np.hstack(([1], np.zeros(numRows - 1)))
    x = np.zeros(numCols)

    # Given's rotation (converts upper Hessenberg matrix to upper triangular matrix)
    for i in range(numCols):
        d = math.sqrt(H[i][i]**2 + H[i + 1][i]**2)
        s = H[i + 1][i]/d
        c = H[i][i]/d
        e1 = b[i]
        e2 = b[i + 1]
        b[i] = c*e1 + s*e2
        b[i + 1] = c*e2 - s*e1
        for j in range(i, numCols):
            r1 = H[i][j]
            r2 = H[i + 1][j]
            H[i][j]     = c*r1 + s*r2
            H[i + 1][j] = c*r2 - s*r1

    # Solve upper triangular matrix
    for i in range(numCols - 1, -1, -1):
        accum = 0
        for j in range(numCols - 1, i, -1):
            accum = x[j]*H[i][j]
        x[i] = (b[i] - accum)/H[i][i]

    return x

################################################################################

index = 1
terminate = False

q = b/np.linalg.norm(b)
h = np.array([[0.0]])

while True:
    h = np.vstack((h, [0]*h.shape[-1])) # adds a row of zeros
    Aq = A.dot(q[:,-1])
    Q = copy.deepcopy(Aq)

    for i in range(index):
        value = q[:,-index + i]
        h[i][-1] = np.dot(value, Aq)
        Q -= h[i][-1]*value

    h[index][-1] = np.linalg.norm(Q)
    c = solve(copy.deepcopy(h))

    if h[index][-1] == 0:
        x = q.dot(c) # solution
        print(x)
        break
    else:
        q = np.column_stack((q, Q/h[index][-1]))

    h = np.column_stack((h, np.zeros(len(h))))
    index += 1

The singular (non-textbook) problem is solvable by numpy as follows:

from scipy.sparse.linalg import gmres
x, exitCode = gmres(A, b, atol=1e-5)
print(exitCode)

Answer 1

MINRES-QLP appears to be the only algorithm capable of solving singular symmetric and Hermitian Linear systems.