Lattice su(2) gauge theory and random number generation in python

Question

I'm currently writing a simple program in python to simulate 1 + 1 dimensional SU(2) yang mills theory. For the case of SU(2) there exists a particular heatbath algorithm for updating the Link variables. In order to implement this algorithm however I need to generate a random real number X such that X is distributed according to P(x) = sqrt(1-X^2)*e^(k*X), where k is a real number from negative infinity to infinity.

Fortunately there exists an algorithm to generate X's according to said distribution. Using my limited skills in python I implemented such an algorithm. Here is the code. I'm only using numpy here.

def algorithm(k):
    count = 0
    while  1 != 0:
        r1,r2,r3,r4 = np.random.uniform(low=0,high=1),np.random.uniform(low=0,high=1),np.random.uniform(low=0,high=1),np.random.uniform(low=0,high=1)
        L1 = -1/(2*k)*(np.log(r1)+np.log(r3)*np.cos(2*np.pi*r2)**2)
        if r4**2 <= 1 - L1:
            X = 1 -2*L1
            break
        else:
            count = count +  1
            continue
    print(count)
    return X  

Basically, if we take three uniformly distributed random numbers in the intervals 0 to 1 we can generate a random variable l1 which is a function of the three random numbers.

We accept this value L1 if 1 - L1 is greater than or equal to a fourth random number squared (uniformly distributed in the interval 0 to 1). Else we loop back to the beginning and do it all over again. We do this until we accept a value of L1. After we accept L1 we compute X as being 1 - 2*L1. This algorithm ensures that X follows the required distribution.

In my program I'm going to have to generate an two dimensional array of X's. This is quite slow in my current implementation. So here's my question; is there a simpler way to do this using any preset numpy packages? If such a method doesn't exist, is there a way to vectorize this function to generate a two dimensional lattice of random X's without simply iterating it with a for loop?

Answer 1

I don't know whether a built-in function exists that returns exactly the distribution you want, but I believe vectorizing your code shouldn't be hard. Just make r1, r2, r3 and r4 vectors using the size parameter of the uniform function as mentioned by Warren and perform those operations. As mentioned by DSM you can also just use a tuple as size parameter and do everything with a single call.

You could keep the loop and repeat the operations in some way until you have N values, but I'd simply remove the loop and keep only the numbers that satisfy the condition. This yields less than N numbers but is straightforward to code.

Something like this might be what you want:

def algorithm_2(k, N):
    r1,r2,r3,r4 = np.random.uniform(low=0,high=1, size=(4,N))
    L1 = -1/(2*k)*(np.log(r1)+np.log(r3)*np.cos(2*np.pi*r2)**2)
    reduced_L1 = L1[r4**2 <= 1 - L1]
    return 1-2*reduced_L1

Running it gives:

>>> algorithm_2(1, 50)
array([-0.21110547, -0.70285195,  0.0475383 , -0.20860877, -0.07776909,
       -0.21907097,  0.70566776,  0.3207524 ,  0.71130986,  0.45789795,
        0.15865965, -0.13757757,  0.04306286,  0.46003952])

If you want a function that always returns exactly an N-ary vector you can write a wrapper that keeps calling the above function and then concatenates the arrays. Something like this:

def algorithm_3(k, N):
    total_size =0
    random_arrays = []
    while total_size < N:
        random_array = algorithm_2(k, N)
        total_size += len(random_array)
        random_arrays.append(random_array)
    return np.hstack(random_arrays)[:N]