<jupyter_start><jupyter_text>Convergence AccelerationThis notebook shows crude implementations of some convergence acceleration techniques. We are testing on a very simple case: Jacobi iterations for the solution of an SPD linear system.<jupyter_code>import numpy as np
from numpy import linalg as LA
import itertools
from collections import OrderedDict
from bokeh.io import push_notebook, show, output_notebook
from bokeh.layouts import row
from bokeh.plotting import figure
from bokeh.palettes import Dark2_5 as palette

output_notebook()<jupyter_output><empty_output><jupyter_text>Jacobi iterationsThe next function performs standad Jacobi iterations to solve a linear system with matrix $\mathbf{A}$ and right-hand side vector $\mathbf{b}$.Let us split the stiffness matrix into its diagonal and off-diagonal components:\begin{equation} \mathbf{A} = \mathbf{O} + \mathbf{D}\end{equation}which leads to:\begin{equation}  \mathbf{A}\mathbf{x} = (\mathbf{O} + \mathbf{D})\mathbf{x} = \mathbf{b},\end{equation}then moving the vector $\mathbf{O}\mathbf{x}$ to the right-hand side gives us the, formally trivial, diagonal system:\begin{equation}  \mathbf{D}\mathbf{x} = \mathbf{b} - \mathbf{O}\mathbf{x},\end{equation}whose solution yields the form of the fixed-point iterations we seek:\begin{equation}  \mathbf{x}^{[n+1]} = \mathbf{D}^{-1}(\mathbf{b} - \mathbf{O}\mathbf{x}^{[n]})\end{equation}<jupyter_code>def apply_jacobi(A, b, x):
    """
    Compute next solution of linear system by applying a Jacobi step.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b
      RHS vector
    x
      Current solution

    Returns
    -------
    Next iteratate in Jacobi fixed-point iteration sequence
    """
    D = np.diag(A)
    O = A - np.diagflat(D)
    return (b - np.einsum('ij,j->i', O, x)) / D


def jacobi(A, b, rtol=1.0e-8, max_it=50, x_0=None, print_report=False):
    """
    Solve linear system of equations with Jacobi iterations.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b 
      RHS vector
    rtol 
      Tolerance for solution norm
    max_it
      Maximum number of iterations
    x_0
      Initial guess
    print_report
      Whether to print iteration report inside the function
    
    Returns
    -------
    x 
      Linear system solution vector
    n_it 
      Number of iterations performed
    report
      Iteration number, residual norm sequence
    """

    x = np.zeros_like(b) if x_0 is None else x_0
    it = 0
    report = []
    if print_report:
        print('Iteration #    Residual norm')
    while it < max_it:
        # Compute residual
        r = apply_jacobi(A, b, x) - x
        rnorm = LA.norm(r)

        # Report
        report.append(rnorm)
        if print_report:
            print('    {:4d}        {:.5E}'.format(it, rnorm))

        # Check convergence
        if rnorm < rtol:
            break

        # Update solution vector
        x = apply_jacobi(A, b, x)
        it += 1
    else:
        raise RuntimeError(
            'Maximum number of iterations ({0:d}) exceeded, but residual norm {1:.5E} still greater than threshold {2:.5E}'.
            format(max_it, rnorm, rtol))
    return x, it, report<jupyter_output><empty_output><jupyter_text>Conjugate gradient and conjugate residualAssuming that the stiffness matrix is SPD, the conjugate gradient (CG) method would be a more robust choice of iterative method than Jacobi iterations. If the stiffness matrix is only known to be Hermitian, the conjugate residual (CR) method can be used instead.<jupyter_code>def cg(A, b, rtol=1.0e-8, max_it=25, x_0=None, print_report=False):
    """
    Solve linear system of equations with conjugate gradient.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b 
      RHS vector
    rtol 
      Tolerance for solution norm
    max_it
      Maximum number of iterations
    x_0
      Initial guess
    print_report
      Whether to print iteration report inside the function
    
    Returns
    -------
    x 
      Linear system solution vector
    n_it 
      Number of iterations performed
    report
      Iteration number, residual norm sequence
    """
       
    x = np.zeros_like(b) if x_0 is None else x_0
    it = 0
    report = []
    r = b - np.einsum('ij,j->i', A, x)
    p = r
    if print_report:
        print('Iteration #    Residual norm')
    while it < max_it:
        Ap = np.einsum('ij,j->i', A, p)
        rtr = np.einsum('i,i', r, r)
        alpha = rtr / np.einsum('i,i', p, Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rnorm = LA.norm(r)

        # Report
        report.append(rnorm)
        if print_report:
            print('    {:4d}        {:.5E}'.format(it, rnorm))

        # Check convergence
        if rnorm < rtol:
            break

        beta = np.einsum('i,i', r, r) / rtr
        p = r + beta * p
        it += 1
    else:
        raise RuntimeError(
            'Maximum number of iterations ({0:d}) exceeded, but residual norm {1:.5E} still greater than threshold {2:.5E}'.
            format(max_it, rnorm, rtol))
    return x, it, report


def cr(A, b, rtol=1.0e-8, max_it=25, x_0=None, print_report=False):
    """
    Solve linear system of equations with conjugate residual.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b 
      RHS vector
    rtol 
      Tolerance for solution norm
    max_it
      Maximum number of iterations
    x_0
      Initial guess
    print_report
      Whether to print iteration report inside the function
    
    Returns
    -------
    x 
      Linear system solution vector
    n_it 
      Number of iterations performed
    report
      Iteration number, residual norm sequence
    """  

    x = np.zeros_like(b) if x_0 is None else x_0
    it = 0
    report = []
    r = b - np.einsum('ij,j->i', A, x)
    p = r
    if print_report:
        print('Iteration #    Residual norm')
    while it < max_it:
        Ap = np.einsum('ij,j->i', A, p)
        rtAr = np.einsum('i,ij,j', r, A, r)
        alpha = rtAr / np.einsum('i,i', Ap, Ap)
        x = x + alpha * p
        r = r - alpha * Ap
        rnorm = LA.norm(r)

        # Report
        report.append(rnorm)
        if print_report:
            print('    {:4d}        {:.5E}'.format(it, rnorm))

        # Check convergence
        if rnorm < rtol:
            break

        beta = np.einsum('i,ij,j', r, A, r) / rtAr
        p = r + beta * p
        it += 1
    else:
        raise RuntimeError(
            'Maximum number of iterations ({0:d}) exceeded, but residual norm {1:.5E} still greater than threshold {2:.5E}'.
            format(max_it, rnorm, rtol))
    return x, it, report<jupyter_output><empty_output><jupyter_text>Kaczmarz iterations\begin{equation}  \mathbf{x}^{[k+1]} = \mathbf{x}^{[k]} + \frac{b_i - \mathbf{A}_{i}\cdot\mathbf{x}^{[k]}}{||\mathbf{A}_{i}||^{2}}\mathbf{A}_{i}^{t},\end{equation}with $i = k \mod m$.<jupyter_code>def kaczmarz(A, b, rtol=1.0e-8, max_it=2005, x_0=None, print_report=False):
    
    x = np.zeros_like(b) if x_0 is None else x_0
    it = 0
    report = []
    r = b - np.einsum('ij,j->i', A, x)
    # Number of rows
    m = A.shape[0]
    if print_report:
        print('Iteration #    Residual norm')
    while it < max_it:
        # Kaczmarz sweeps
        i = it%m
        alpha = (b[i] - np.einsum('i,i->', A[i,:], x)) / np.einsum('i,i->', A[i,:], A[i,:])
        x = x + alpha * A[i,:].T
        r = b - np.einsum('ij,j->i', A, x)
        rnorm = LA.norm(r)
        
        # Report
        report.append(rnorm)
        if print_report:
            print('    {:4d}        {:.5E}'.format(it, rnorm))
            
        # Check convergence
        if rnorm < rtol:
            break
            
        it +=1
    else:
        raise RuntimeError(
            'Maximum number of iterations ({0:d}) exceeded, but residual norm {1:.5E} still greater than threshold {2:.5E}'.
            format(max_it, rnorm, rtol))
    return x, it, report<jupyter_output><empty_output><jupyter_text>DIIS accelerationDirect inversion in the iterative subpace (DIIS) is an acceleration method ubiquitous in quantum chemistry. Given a queue of iterates $\lbrace \mathbf{x}^{[i]} \rbrace_{i=0}^{h}$, where $h$ is the maximum size of the queue, the next iterate in the sequence is obtained as a linear combination of the iterates in the queue:\begin{equation}  \tilde{\mathbf{x}} := \sum_{i=0}^{h} c_{i} \mathbf{x}_{i}\end{equation}with the coefficients subject to the normalization condition:\begin{equation}  \mathbf{1} \cdot \mathbf{c} = \sum_{i=0}^{h} c_{i} = 1.\end{equation}The optimal coefficients are found by minimizing the squared norm of the extrapolated iterate, under the normalization constraint. This leads to a linear system, sometimes known as _Pulay's equations_:\begin{equation}  \begin{pmatrix}    \mathbf{x}_{0}\cdot\mathbf{x}_{0} & \cdots & \mathbf{x}_{0}\cdot\mathbf{x}_{n} & -1 \\    \vdots & \cdots & \vdots \\    \mathbf{x}_{n}\cdot\mathbf{x}_{0} & \cdots & \mathbf{x}_{n}\cdot\mathbf{x}_{n} & -1 \\    -1  & \cdots & -1 & 0 \\  \end{pmatrix}  \begin{pmatrix}  c_{0} \\  \vdots \\  c_{n} \\  \lambda \\  \end{pmatrix}  =  \begin{pmatrix}  0 \\  \vdots \\  0 \\  -1 \\  \end{pmatrix}\end{equation}<jupyter_code>def jacobi_diis(A, b, rtol=1.0e-8, max_it=25, max_hist=8, x_0=None, print_report=False):
    """
    Solve linear system of equations with Jacobi iterations and DIIS acceleration.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b 
      RHS vector
    rtol 
      Tolerance for solution norm
    max_it
      Maximum number of iterations
    max_hist
      Maximum size of DIIS queue
    x_0
      Initial guess
   print_report
      Whether to print iteration report inside the function
    
    Returns
    -------
    x 
      Linear system solution vector
    n_it 
      Number of iterations performed
    report
      Iteration number, residual norm sequence      
    """

    x = np.zeros_like(b) if x_0 is None else x_0

    # Lists of iterates and residuals
    xs = []
    rs = []

    report = []
    it = 0
    if print_report:
        print('Iteration #    Residual norm')
    while it < max_it:
        # Compute residual
        r = apply_jacobi(A, b, x) - x
        rnorm = LA.norm(r)

        # Collect DIIS history
        xs.append(x)
        rs.append(r)

        # Report
        report.append(rnorm)
        if print_report:
            print('    {:4d}        {:.5E}'.format(it, rnorm))

        # Check convergence
        if rnorm < rtol:
            break

        if it >= 2:
            # Prune DIIS history
            diis_hist = len(xs)
            if diis_hist > max_hist:
                xs.pop(0)
                rs.pop(0)
                diis_hist -= 1

            # Build error matrix B
            B = np.empty((diis_hist + 1, diis_hist + 1))
            B[-1, :] = -1
            B[:, -1] = -1
            B[-1, -1] = 0
            for i, ri in enumerate(rs):
                for j, rj in enumerate(rs):
                    if j > i: continue
                    val = np.einsum('i,i', ri, rj)
                    B[i, j] = val
                    B[j, i] = val

            # Normalize B
            B[:-1, :-1] /= np.abs(B[:-1, :-1]).max()

            # Build rhs vector
            rhs = np.zeros(diis_hist + 1)
            rhs[-1] = -1

            # Solve Pulay equations
            cs = LA.solve(B, rhs)

            # Calculate new solution as linear
            # combination of previous solutions
            x = np.zeros_like(x)
            for i, c in enumerate(cs[:-1]):
                x += c * xs[i]

        # Update solution vector
        x = apply_jacobi(A, b, x)
        it += 1
    else:
        raise RuntimeError(
            'Maximum number of iterations ({0:d}) exceeded, but residual norm {1:.5E} still greater than threshold {2:.5E}'.
            format(max_it, rnorm, rtol))
    return x, it, report<jupyter_output><empty_output><jupyter_text>KAIN accelerationThe Krylov accelerated inexact Newton (KAIN) method is another example of convergence acceleration algorithm similar to DIIS.Given the nonlinear equation $\mathbf{\omega}(\mathbf{x}) = 0$, with $\mathbf{\omega}: \mathbb{R}^{N} \rightarrow \mathbb{R}^{N}$,the Newton method formulates the solution as iterative updates $\mathbf{x}^{[n+1]} = \mathbf{x}^{[n]} + \mathbf{\delta}^{[n]}$ where $\mathbf{\delta}^{[n]}$ is of the form:\begin{equation}  \mathbf{A}^{[n]}\mathbf{\delta}^{[n]} = -\mathbf{\omega}^{[n]},\end{equation}with $\mathbf{A}^{[n]}$ the Jacobian computed at the $n$-th iteration, _i.e._ from the $\mathbf{x}^{[n]}$ iterate.We assume that we have generated a queue of iterates $\lbrace \mathbf{x}^{[i]}\rbrace_{i=0}^{h}$ and corresponding queue of function iterates $\lbrace \mathbf{\omega}^{[i]} \equiv \mathbf{\omega}(\mathbf{x}^{[i]})\rbrace_{i=0}^{h}$. As before, here $h$ is the maximum size of the queues. Similarly to DIIS we seek to construct the update _without_ ever forming the full Jacobian. However, at variance with DIIS, we would like to include information from _outside_ the iterative subspace in our update.An approximation of the Jacobian _within the subspace of iterates differences_ can be built as:\begin{equation} \omega^{[i]} = \omega(\mathbf{x}^{[i]} + \mathbf{x}^{[h]} - \mathbf{x}^{[h]})  \simeq  \omega^{[h]} + \mathbf{A}^{[h]}(\mathbf{x}^{[i]} - \mathbf{x}^{[h]})\end{equation}such that _within the subspace_ $\mathcal{D} = \lbrace \mathbf{y}^{[k]} | \mathbf{y}^{[k]} = \mathbf{x}^{[i]} - \mathbf{x}^{[h]},\, \forall k = 0,\ldots, h-1 \rbrace$:\begin{equation} \mathbf{A}^{[h]}(\mathbf{x}^{[i]} - \mathbf{x}^{[h]}) \simeq \omega^{[i]} - \omega^{[h]}.\end{equation}We introduce the projector into the $\mathcal{D}$ subspace $\hat{P}$ and its complement $\hat{Q}$ and we partition the Newton equation accordingly:\begin{equation}  \mathbf{A}^{[h]}\hat{P}\delta^{[h]} + \mathbf{D}^{[h]}\hat{Q}\delta^{[h]} = - \omega^{[h]},\end{equation}where $\mathbf{D}^{[h]}$ is a _readily invertible approximation_ to the full Jacobian _outside_ the the $\mathcal{D}$ subspace.Now, the update _within_ the subspace is a clearly a linear combination of differences of previous iterates:\begin{equation}  \hat{P}\delta^{[h]} = \sum_{i=0}^{h-1} c_{i}(\mathbf{x}^{[i]} - \mathbf{x}^{[h]})\end{equation}with coefficients yet to be determined. The update _outside_ the subspace is then:\begin{equation} \hat{Q}\delta^{[h]} = - (\mathbf{D}^{[h]})^{-1}\left( \mathbf{A}^{[h]}\hat{P}\mathbf{\delta}^{[h]} + \omega^{[h]}\right).\end{equation}Projection onto the the $\mathcal{D}$ subspace and insertion of the equation for $\hat{P}\delta^{[h]}$ yields the subspace equations for the mixing coefficients $\mathbf{c}$:\begin{equation} \sum_{j=0}^{h-1}  \langle \mathbf{x}^{[i]} - \mathbf{x}^{[h]} | (\mathbf{D}^{[h]})^{-1}(\omega^{[j]} - \omega^{[h]}) \rangle c_{j}  = - \langle \mathbf{x}^{[i]} - \mathbf{x}^{[h]} | (\mathbf{D}^{[h]})^{-1}\omega^{[h]} \rangle\end{equation}The final form for the inexact Newton update is then:\begin{equation} \delta^{[h]} = \sum_{i=0}^{h-1} c_{i}(\mathbf{x}^{[i]} - \mathbf{x}^{[h]})  - (\mathbf{D}^{[h]})^{-1}\left(\sum_{i=0}^{h-1} c_{i}(\omega^{[i]} - \omega^{[h]}) + \omega^{[h]} \right)\end{equation} The $\mathbf{Q}$ matrixThe KAIN subspace equations can be manipulated a bit further:\begin{equation} \sum_{j=0}^{h-1}    (Q^{[h]}_{ij} - Q^{[h]}_{ih} - Q^{[h]}_{hj} + Q^{[h]}_{hh})c_{j} =    - Q^{[h]}_{ih} + Q^{[h]}_{hh}\end{equation}with the $\mathbf{Q}$ matrix being:\begin{equation} Q^{[h]}_{ij} = \langle \mathbf{x}^{[i]} | (\mathbf{D}^{[h]})^{-1}\omega^{[j]} \rangle\end{equation}Some further manipulation yields:\begin{equation} \sum_{j=0}^{h-1}    (Q^{[h]}_{ij} - Q^{[h]}_{hj})c_{j} + (Q^{[h]}_{hh} - Q^{[h]}_{ih})\left(\sum_{j=0}^{h-1} c_{j} \right) =    - Q^{[h]}_{ih} + Q^{[h]}_{hh}.\end{equation}We now impose the normalization condition $\sum_{j=0}^{h} c_{j} = 1$ such that the subspace equation becomes:\begin{equation} \sum_{j=0}^{h}    (Q^{[h]}_{ij} - Q^{[h]}_{hj})c_{j} =    0,\end{equation}while the update now reads:\begin{equation}\tilde{\mathbf{x}} = \sum_{j=0}^{h} c_{j}\mathbf{x}^{[j]} - (\mathbf{D}^{[h]})^{-1}\left(\sum_{j=0}^{h} c_{j}\omega^{[j]}\right)\end{equation}<jupyter_code>def jacobi_kain(A, b, rtol=1.0e-8, max_it=25, max_hist=8, x_0=None, print_report=False):
    """
    Solve linear system of equations with Jacobi iterations and KAIN acceleration.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b 
      RHS vector
    rtol 
      Tolerance for solution norm
    max_it
      Maximum number of iterations
    max_hist
      Maximum size of KAIN queue
    x_0
      Initial guess
   print_report
      Whether to print iteration report inside the function
    
    Returns
    -------
    x 
      Linear system solution vector
    n_it 
      Number of iterations performed
    report
      Iteration number, residual norm sequence      
    """
        
    x = np.zeros_like(b) if x_0 is None else x_0

    # Lists of iterates and residuals
    xs = []
    rs = []

    report = []
    it = 0
    if print_report:
        print('Iteration #    Residual norm')
    while it < max_it:
        # Compute residual
        r = apply_jacobi(A, b, x) - x
        rnorm = LA.norm(r)

        # Collect DIIS history
        xs.append(x)
        rs.append(r)

        # Report
        report.append(rnorm)
        if print_report:
            print('    {:4d}        {:.5E}'.format(it, rnorm))

        # Check convergence
        if rnorm < rtol:
            break

        if it >= 2:
            # Prune KAIN history
            kain_hist = len(xs)
            if kain_hist > max_hist:
                xs.pop(0)
                rs.pop(0)
                kain_hist -= 1

            # Build Q matrix
            Q = np.empty((kain_hist, kain_hist))
            for i, x in enumerate(xs):
                for j, r in enumerate(rs):
                    Q[i, j] = np.einsum('i,i', x, r)

            # Build subspace matrix B and vector rhs
            B = np.empty((kain_hist - 1, kain_hist - 1))
            rhs = np.zeros(kain_hist - 1)
            for i in range(kain_hist - 1):
                rhs[i] = -Q[i, -1] + Q[-1, -1]
                for j in range(kain_hist - 1):
                    B[i, j] = Q[i, j] - Q[i, -1] - Q[-1, j] + Q[-1, -1]

            # Solve KAIN equations
            cs = LA.solve(B, rhs)
            cs = np.append(cs, 1.0 - np.sum(cs))

            # Calculate new solution as linear
            # combination of previous solutions
            x = np.zeros_like(x)
            for i, c in enumerate(cs):
                x += c * xs[i] - c * rs[i]

        # Update solution vector
        x = apply_jacobi(A, b, x)
        it += 1
    else:
        raise RuntimeError(
            'Maximum number of iterations ({0:d}) exceeded, but residual norm {1:.5E} still greater than threshold {2:.5E}'.
            format(max_it, rnorm, rtol))
    return x, it, report<jupyter_output><empty_output><jupyter_text>CROP accelerationCHECK: CROP should only need 3 previous iterates at most for linear problems.<jupyter_code>def jacobi_crop(A, b, rtol=1.0e-8, max_it=25, max_hist=8, x_0=None):
    """
    Solve linear system of equations with Jacobi iterations and CROP acceleration.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b 
      RHS vector
    rtol 
      Tolerance for solution norm
    max_it
      Maximum number of iterations
    max_hist
      Maximum size of CROP queue
    x_0
      Initial guess
    
    Returns
    -------
    x 
      Linear system solution vector
    """
    
    pass<jupyter_output><empty_output><jupyter_text>Regularized nonlinear accelerationSee also: https://github.com/windows7lover/RegularizedNonlinearAcceleration<jupyter_code>def jacobi_rna(A, b, rtol=1.0e-8, max_it=25, max_hist=8, x_0=None, reg=0.1, print_report=False):
    """
    Solve linear system of equations with Jacobi iterations and RNA.
    
    Parameters
    ----------
    A
      Stiffness matrix
    b 
      RHS vector
    rtol 
      Tolerance for solution norm
    max_it
      Maximum number of iterations
    max_hist
      Maximum size of RNA queue
    x_0
      Initial guess
    reg
      Regularization parameter
   print_report
      Whether to print iteration report inside the function
    
    Returns
    -------
    x 
      Linear system solution vector
    n_it 
      Number of iterations performed
    report
      Iteration number, residual norm sequence      
    """

    x = np.zeros_like(b) if x_0 is None else x_0

    # Lists of iterates
    xs = []

    report = []
    it = 0
    if print_report:
        print('Iteration #    Residual norm')
    while it < max_it:
        # Compute residual
        r = apply_jacobi(A, b, x) - x
        rnorm = LA.norm(r)
        
        # Collect history
        xs.append(x)

        # Report
        report.append(rnorm)
        if print_report:
            print('    {:4d}        {:.5E}'.format(it, rnorm))

        # Check convergence
        if rnorm < rtol:
            break

        if it >= 2:
            # Prune history
            hist = len(xs)
            if hist > max_hist:
                xs.pop(0)
                hist -= 1

            # Build R matrix of residuals
            X = np.asarray(xs)
            R = np.diff(X, axis=0)
            
            # Build R^{t}R and normalize it
            RR = np.einsum('ik,jk->ij', R, R)
            print('RR is ', RR)
            my_RR = RR = np.dot(np.transpose(R),R)
            print('my_RR is ', my_RR)
            RR /= np.sqrt(np.einsum('ii', RR))

            # Solve RNA equations
            zs = LA.solve(RR + reg*np.eye(RR.shape[0]), np.ones(RR.shape[0]))
            
            # Recover mixing weights
            cs = zs / np.sum(zs)
            # Calculate new solution as linear
            # combination of previous solutions
            x = np.zeros_like(b)
            for i, c in enumerate(cs[:-1]):
                x += c * xs[i]

        # Update solution vector
        x = apply_jacobi(A, b, x)

        # Update iteration counter
        it += 1
    else:
        raise RuntimeError(
            'Maximum number of iterations ({0:d}) exceeded, but residual norm {1:.5E} still greater than threshold {2:.5E}'.
            format(max_it, rnorm, rtol))
    return x, it, report
def relative_error_to_reference(x, x_ref):
    return LA.norm(x - x_ref) / LA.norm(x_ref)


print('Experiments with linear solvers')
dim = 100
M = np.random.randn(dim, dim)
# Make sure our matrix is SPD
A = 0.5 * (M + M.transpose())
A = A * A.transpose()
A += dim * np.eye(dim)
b = np.random.rand(dim)
x_ref = LA.solve(A, b)

methods = OrderedDict({
    'Jacobi' : jacobi, 
    'Jacobi-DIIS' : jacobi_diis, 
    'Jacobi-KAIN' : jacobi_kain,
    #'Jacobi-RNA' : jacobi_rna, # NOT working
    'Kaczmarz' : kaczmarz,
    'Conjugate gradient' : cg, 
    'Conjugate residual': cr
     })
reports = OrderedDict()
for k, v in methods.items():
    print('@ {:s} algorithm'.format(k))
    x, it, report = v(A, b)
    rel_err = relative_error_to_reference(x, x_ref)
    print('{:s} converged in {:d} iterations with relative error to reference {:.5E}\n'.format(k, it, rel_err))
    reports.update({k: report})

print('--- Final report ---')
# TODO generate table
#header = ' '.join(reports.keys())

colors = itertools.cycle(palette)    
p = figure(plot_width=800, plot_height=400, title='Convergence', y_axis_type='log')
for method, report, color in zip(reports.keys(), reports.values(), colors):
    p.line(np.arange(len(reports['Kaczmarz'])), report, line_width=2, legend=method, color=color)
show(p)<jupyter_output>Experiments with linear solvers
@ Jacobi algorithm
Jacobi converged in 23 iterations with relative error to reference 1.33673E-07

@ Jacobi-DIIS algorithm
Jacobi-DIIS converged in 7 iterations with relative error to reference 3.99412E-08

@ Jacobi-KAIN algorithm
Jacobi-KAIN converged in 8 iterations with relative error to reference 1.89029E-07

@ Kaczmarz algorithm
Kaczmarz converged in 1472 iterations with relative error to reference 1.92513E-09

@ Conjugate gradient algorithm
Conjugate gradient converged in 8 iterations with relative error to reference 1.55673E-10

@ Conjugate residual algorithm
Conjugate residual converged in 8 iterations with relative error to reference 1.56918E-10

--- Final report ---
{'fraction_non_alphanumeric': 0.13924464994476043, 'fraction_numerical': 0.012930152624902819, 'mean_word_length': 2.295132803020089, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 49, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '38891', 'license': 'MIT', 'n_tokens_mistral': 8357, 'n_tokens_neox': 7850, 'n_words': 2299, 'path': 'convergence-acceleration.ipynb', 'repo_name': 'robertodr/solver', 'size': '24439'}