# Convergence Acceleration

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


```python
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()
```



<div class="bk-root">
    <a href="https://bokeh.pydata.org" target="_blank" class="bk-logo bk-logo-small bk-logo-notebook"></a>
    <span id="d370c925-7561-41fc-9132-4f085eb0444d">Loading BokehJS ...</span>
</div>




## Jacobi iterations

The 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}


```python
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 
```

## Conjugate gradient and conjugate residual

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


```python
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
```

## 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$.


```python
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
```

## DIIS acceleration

Direct 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}


```python
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
```

## KAIN acceleration

The 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}$ matrix

The 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}


```python
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
```

## CROP acceleration

CHECK: CROP should only need 3 previous iterates at most for linear problems.


```python
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
```

## Regularized nonlinear acceleration

See also: https://github.com/windows7lover/RegularizedNonlinearAcceleration


```python
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
```


```python
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)
```

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


    /nix/store/i9kckiz2x7287kjn34a8k024lj6xfiwa-python3.6-bokeh0.12.14/lib/python3.6/site-packages/bokeh/models/sources.py:138: BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('x', 1473), ('y', 24)
      "Current lengths: %s" % ", ".join(sorted(str((k, len(v))) for k, v in data.items())), BokehUserWarning))
    /nix/store/i9kckiz2x7287kjn34a8k024lj6xfiwa-python3.6-bokeh0.12.14/lib/python3.6/site-packages/bokeh/models/sources.py:138: BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('x', 1473), ('y', 8)
      "Current lengths: %s" % ", ".join(sorted(str((k, len(v))) for k, v in data.items())), BokehUserWarning))
    /nix/store/i9kckiz2x7287kjn34a8k024lj6xfiwa-python3.6-bokeh0.12.14/lib/python3.6/site-packages/bokeh/models/sources.py:138: BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('x', 1473), ('y', 9)
      "Current lengths: %s" % ", ".join(sorted(str((k, len(v))) for k, v in data.items())), BokehUserWarning))
    /nix/store/i9kckiz2x7287kjn34a8k024lj6xfiwa-python3.6-bokeh0.12.14/lib/python3.6/site-packages/bokeh/models/sources.py:138: BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('x', 1473), ('y', 9)
      "Current lengths: %s" % ", ".join(sorted(str((k, len(v))) for k, v in data.items())), BokehUserWarning))
    /nix/store/i9kckiz2x7287kjn34a8k024lj6xfiwa-python3.6-bokeh0.12.14/lib/python3.6/site-packages/bokeh/models/sources.py:138: BokehUserWarning: ColumnDataSource's columns must be of the same length. Current lengths: ('x', 1473), ('y', 9)
      "Current lengths: %s" % ", ".join(sorted(str((k, len(v))) for k, v in data.items())), BokehUserWarning))




<div class="bk-root">
    <div class="bk-plotdiv" id="c05d3537-02fe-4b4e-a45b-35e225cf91aa"></div>
</div>





```python

```

10a004c4b179b55a15e3c96efba2b70a236f1736
ipynb
Jupyter Notebook
convergence-acceleration.ipynb
robertodr/solver
dced957cb3f2aa8c1f3ab085d8445c6dd7d3d284
MIT
2018-11-09T17:03:09.000Z
2018-11-19T10:30:04.000Z
convergence-acceleration.ipynb
robertodr/solver
dced957cb3f2aa8c1f3ab085d8445c6dd7d3d284
MIT
convergence-acceleration.ipynb
robertodr/solver
dced957cb3f2aa8c1f3ab085d8445c6dd7d3d284
MIT
Qwen/Qwen-72B
1. YES
2. YES