0
import math
import numpy as np

#given ordered set, e.g. [A,B,C], and some value D, return 0 if D ≤ A,
#return 1 if A ≤ D < B, return 2 if B ≤ D < C, and 3 if C < D
#O(log(n)) time complexity

def findID(partition,val):   
    ub = len(partition)

    if ub == 0:
        return 0

    midx = math.floor(ub/2)
    mval = partition[midx]

    if midx == 0 and val > mval:
        return 1

    if val <= mval:
        return findID(partition[:midx],val)

    return findID(partition[midx:],val)+midx



def rKaczmarz(A,b,K):
    (m,n) = A.shape
    d = len(b)

    if m != d:
        return -1

    #normRowSq stores 2norm of all rows of A
    normRowSq = []
    for i in range(m):
        normRowSq.append(np.linalg.norm(A[i])**2)

    #frobenius norm squared
    FnormSq = sum(normRowSq)

    partition = [normRowSq[0]/FnormSq]
    for i in range(m-2):
        partition.append(normRowSq[i+1]/FnormSq + partition[i])

    x = np.random.rand(n)
    for j in range(K):
        p = np.random.uniform(0,1)
        i = findID(partition,p)

        x += (1/normRowSq[i])*(b[i]-np.dot(A[i],x))*A[i]

    return x

#examples:
X = np.array([[1,2,3],
              [4,5,6],
              [7,8,1],
              [1,1,1]], float)

Y = [-1,-2,0,3]

x = rKaczmarz(X, Y, 1000)
print(x)

I wrote the code above which implements randomized kaczmarz and according to a theorem, the expected convergence rate of randomized kaczmarz to least square is given by: $$e^{(k)} = (1-\frac 1 R)e^{(0)}, \text{ where } R = \frac{\|A\|_F^2}{\sigma_{min}(A)}$$ according to page 6 of this paper: https://arxiv.org/pdf/1205.5770.pdf, which also has the algorithm implementation

and where here $e^{(k)}$ refers to the 2-norm error at the k-th iteration. So for this particular X and Y I would expect the 2-norm error to be as small as $4.1\cdot 10^{-9}$ after 10000 iteration.

However this is not what I am seeing in reality, where there is still decently significant fluctuations after 10000 iterations. Is there something fundamentally wrong with my implementation?

(It does not seem like my implementation is significantly flawed, as replacing i = findID(....) randomization with the regular Kaczmarz i = j%m would result in fast convergence to a value very far from least square solution)

nabu1227
  • 879
  • Actually I am unsure, in the paper it writes $R = \frac{|A|F^2}{\sigma{min(m,n)}(A^\top A)}$, does that refer to singular value of $A$ or $A^\top A$? – nabu1227 Apr 15 '21 at 09:54
  • i have found the same in my implementations. Unless condition number is low, Randomized Kaczmarz convergence is slow. That said, there are ways to speed it up. – Rahul Madhavan Apr 15 '21 at 10:42
  • I think my issue is even if the condition number is low, the method hardly even converges for larger matrices. – nabu1227 Apr 15 '21 at 23:36

0 Answers0