# Kaczmarz Algorithm

Given $Ax = b$, iteratively solve with the update:
\begin{equation}
x^{(k+1)} = x^{(k)} + \frac{b_i - \langle a_i , x^{(k)}\rangle}{\lvert a_i \rvert^2}a_i^*
\end{equation}
where $i = k \mathrm{mod} m +1$


```python
import numpy as np
from numpy import random
from numpy import linalg as LA
from collections import OrderedDict


def relative_error_to_reference(x, x_ref):
    return LA.norm(x - x_ref) / LA.norm(x_ref)


def cg(A, x, b, tol=1.0e-08, max_it=200):
    # Iteration counter
    k = 0
    # Residual
    r = b - np.einsum('ij,j', A, x)
    # Search direction
    p = r
    while LA.norm(r) > tol:
        tmp = np.einsum('i,i', r, r)
        alpha = tmp / np.einsum('i,ij,j', p, A, p)
        x = x + alpha * p
        # Update residual
        r = r - alpha * np.einsum('ij,j', A, p)
        # Update search direction
        beta = np.einsum('i,i', r, r) / tmp
        p = r + beta * p
        #print('Iteration {}, Norm of residual {}'.format(k, LA.norm(r)))
        k += 1
        if k >= max_it:
            raise RuntimeError('Maximum number of iterations reached!')
            break
    return k, LA.norm(r), x

            
def kaczmarz(A, x, b, tol=1.0e-08, max_it=2000):
    m = A.shape[0]
    # Iteration counter
    k = 0
    # Residual
    r = b - np.einsum('ij,j', A, x)
    while LA.norm(r) > tol:
        i = k%m
        alpha = (b[i] - np.einsum('i,i', A[i,:], x)) / np.einsum('i,i', A[i,:], A[i,:])
        x = x + alpha * np.transpose(A[i,:])
        r = b - np.einsum('ij,j', A, x)
        #for i in range(m):
        #    j = m-i-1
        #    alpha = (b[j] - np.einsum('i,i', A[j, :], x)) / np.einsum('i,i', A[j, :], A[j, :])
        #    x = x + alpha * np.conjugate(A[j, :])
        #    r = b - np.einsum('ij,j', A, x)
        #print('Iteration {}, Norm of residual {}'.format(k, LA.norm(r)))
        k += 1
        if k >= max_it:
            raise RuntimeError('Maximum number of iterations reached!')
            break
    return k, LA.norm(r), x


def random_kaczmarz(A, x, b, tol=1.0e-08, max_it=20000):
    m = A.shape[0]
    # Iteration counter
    k = 0
    # Residual
    r = b - np.einsum('ij,j', A, x)
    # Random start row index
    i = random.randint(0, m)
    # 2-norm squared of row
    norm = np.einsum('i,i', A[i, :], A[i, :])
    while LA.norm(r) > tol:
        # Do a Kaczmarz sweep on row i
        alpha = (b[i] - np.einsum('i,i', A[i, :], x)) / np.einsum('i,i', A[i, :], A[i, :])
        x = x + alpha * np.conjugate(A[i, :])
        r = b - np.einsum('ij,j', A, x)
        #print('Iteration {}, Norm of residual {}'.format(k, LA.norm(r)))
        # Propose random index
        prop_i = random.randint(0,m)
        # 2-norm squared of row
        prop_norm = np.einsum('i,i', A[prop_i, :], A[prop_i, :])
        acceptance = min(1.0, prop_norm / norm)
        if random.uniform() < acceptance:
            i = prop_i
            norm = prop_norm
        k += 1
        if k >= max_it:
            raise RuntimeError('Maximum number of iterations reached!')
            break
    return k, LA.norm(r), x

dim = 10
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)

x_0 = np.zeros_like(b)

C = np.arange(6).reshape(2,3)
for i in C:
    print(i)

methods = OrderedDict({
    'Conjugate gradient' : cg, 
    'Kaczmarz' : kaczmarz,
    'Random Kaczmarz' : random_kaczmarz
     })
#reports = OrderedDict()
for k, v in methods.items():
    print('@ {:s} algorithm'.format(k))
    it, residual, x = v(A, x_0, 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})
```

    [0 1 2]
    [3 4 5]
    @ Conjugate gradient algorithm
    Conjugate gradient converged in 9 iterations with relative error to reference 2.10003E-09
    
    @ Kaczmarz algorithm
    Kaczmarz converged in 174 iterations with relative error to reference 7.04794E-09
    
    @ Random Kaczmarz algorithm
    Random Kaczmarz converged in 480 iterations with relative error to reference 8.09928E-09
    



```python

```

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