<jupyter_start><jupyter_text>Randomized Kaczmarz AlgoritihmDemonstration of the basic Randomized Kaczmarz Algorithm  (cf. ["A Randomized Kaczmarz Algorithm with Exponential Convergence"](http://www.springerlink.com/index/10.1007/s00041-008-9030-4)   by Strohmer and Vershynin, 2008 )This solves **overdetermined**, **consistent** systems of linear equations$$ Ax=b $$where $A$ is $M\times N$ (so $M\ge N$).For fancier versions, see:- "Acceleration of randomized Kaczmarz method via theJohnson?Lindenstrauss Lemma" by <NAME> and <NAME> 2011- "Paved with good intentions: analysis of a randomized block kaczmarzmethod" by <NAME> and <NAME> 2012- "Stochastic gradient descent, weighted sampling, and the randomizedKaczmarz algorithm" by <NAME>, <NAME> and <NAME> 2016Code: <NAME> 2019 (Matlab), 2021 (Python)<jupyter_code>import numpy as np
import numpy.linalg
from numpy.random import default_rng
rng = default_rng()
import scipy.linalg
from matplotlib import pyplot as plt
M, N = int(3e5), int(1e2)

# A   = rng.standard_normal( (M,N)) # this is well-conditioned

# This case is much harder...
A   = rng.standard_normal( (M,N) )@np.diag(np.logspace(0,2,N))@(
    rng.standard_normal((N,N) ) + 0.9*np.eye(N) )


x0  = rng.standard_normal( (N,1))
b   = A@x0  # important: no noise, we assume the system is consistent

errFcn    = lambda x : scipy.linalg.norm( x.ravel() - x0.ravel() )
# Solve using classical dense method
#  [ scipy.linalg.solve(A,b) won't work since it needs a square matrix
#    so we'll ask for a least-squares solution ]
print("Solving via classical dense method")
%time xClassical, residLS, rank, singVals = np.linalg.lstsq(A,b,rcond=None)
print(f'  Condition number of A is {singVals[0]/singVals[-1]:.3e}')
print(f'  Error is {errFcn(xClassical):.3e}')<jupyter_output>Solving via classical dense method
CPU times: user 5.66 s, sys: 207 ms, total: 5.86 s
Wall time: 3.24 s
  Condition number of A is 7.473e+03
  Error is 2.911e-13<jupyter_text>Try the randomized Kaczmarz approachClassical Kaczmarz cycles through indices in orderThe Strohmer/Vershynin randomized version chooses indices based on iid samples from the weighted distribution determined by the squared Euclidean norms of the rows<jupyter_code>rowNorms = scipy.linalg.norm( A, ord=2, axis=1)**2
# NOTE: run this cell once without the %%timeit statement,
#   since if you run with that statement, it doesn't save
#   the value of the variables into the workspace.
%%timeit -o -r 1 -n 1
rowNorms = scipy.linalg.norm( A, ord=2, axis=1)**2
tmStructure = _
tmPreprocess = tmStructure.best
plt.hist( rowNorms )
plt.show()<jupyter_output><empty_output><jupyter_text>Now the main iteration<jupyter_code>def RK(N,wts,maxIter,blockSize=1000,tol=1e-13):

  x         = np.zeros(N) # have 1D not 2D column vector
  errList   = []

  #  We will really be taking maxIter * blockSize iterations
  #   but chose blockSize > 1 to reduce overhead costs
  for k in range(int(maxIter)):

    indList = rng.choice(M,int(blockSize),replace=True,p=wts,shuffle=True)

    for i in indList:
      res = b[i] - A[i,:].dot(x)
      x += res/rowNorms[i] * A[i,:]

    errList.append( errFcn(x) )

    if np.abs(res) > 1e100:
      print('Looks like divergence, stopping early')
      break
    if errList[-1] < tol:
      print(f'Reached tolerance, stopping early: err is {errList[-1]}')
      break

  return x, errList
wts = rowNorms.copy()/rowNorms.sum()
maxIter = 5e2
xRK, errList = RK( N, wts, maxIter=maxIter )
%%timeit -o -r 1 -n 1
# Re-run it and get the time (but the outputs are not saved)
xRK, errList = RK( N, wts, maxIter=maxIter )
tmStructure = _
tmIterative = tmStructure.best
plt.semilogy( np.arange(0,len(errList)*blockSize,blockSize)/M, errList,'o-' )
plt.xlabel('Iterations, divided by # of rows')
plt.ylabel('Error in "x" variable')
plt.show()
plt.semilogy( tmPreprocess + np.arange(len(errList))*tmIterative/len(errList), errList,'o-' )
plt.xlabel('Time (seconds)')
plt.ylabel('Error in "x" variable')
plt.xlim( left=0)
plt.show()<jupyter_output><empty_output><jupyter_text>And try with an iterative solver<jupyter_code>%time xLSQR, flag, iter, nrm = scipy.sparse.linalg.lsqr( A, b, show=True, iter_lim=int(1e2))[:4]
errFcn(xLSQR)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1268004722550177, 'fraction_numerical': 0.027390791027154665, 'mean_word_length': 4.739837398373984, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '119085', 'license': 'MIT', 'n_tokens_mistral': 1672, 'n_tokens_neox': 1542, 'n_words': 449, 'path': 'Demos/demo09_RandomizedKaczmarz.ipynb', 'repo_name': 'stephenbeckr/randomized-algorithm-class', 'size': '4263'}