python - 在 python 中矢量化蒙特卡洛模拟

Question

我最近一直在使用 python 编写一些代码，以使用蒙特卡洛方法模拟二维 U(1) 规范理论。本质上，我有一个 n 乘 n 乘 2 的酉复数数组(称其为 Link)(它们的大小为 1)。我随机选择我的链接数组的元素，并建议对该站点的数字进行随机更改。然后我计算由于该更改而发生的操作的最终变化。然后我接受概率等于 min(1,exp(-dS)) 的变化，其中 dS 是 Action 的变化。迭代器的代码如下

def iteration(j1,B0):
    global Link
    Staple = np.zeros((2),dtype=complex)
    for i0 in range(0,j1):
        x1 = np.random.randint(0,n)
        y1 = np.random.randint(0,n)
        u1 = np.random.randint(0,1)
        Linkrxp1 = np.roll(Link,-1, axis = 0)
        Linkrxn1 = np.roll(Link, 1, axis = 0)
        Linkrtp1 = np.roll(Link, -1, axis = 1)
        Linkrtn1 = np.roll(Link, 1, axis = 1)
        Linkrxp1tn1 = np.roll(np.roll(Link, -1, axis = 0),1, axis = 1)
        Linkrxn1tp1 = np.roll(np.roll(Link, 1, axis = 0),-1, axis = 1)


        Staple[0] = Linkrxp1[x1,y1,1]*Linkrtp1[x1,y1,0].conj()*Link[x1,y1,1].conj() + Linkrxp1tn1[x1,y1,1].conj()*Linkrtn1[x1,y1,0].conj()*Linkrtn1[x1,y1,1]
        Staple[1] = Linkrtp1[x1,y1,0]*Linkrxp1[x1,y1,1].conj()*Link[x1,y1,0].conj() + Linkrxn1tp1[x1,y1,0].conj()*Linkrxn1[x1,y1,1].conj()*Linkrxn1[x1,y1,0]


        uni = unitary()
        Linkprop = uni*Link[x1,y1,u1]
        dE3 = (Linkprop - Link[x1,y1,u1])*Staple[u1]
        dE1 = B0*np.real(dE3)
        d1 = np.random.binomial(1, np.minimum(np.exp(dE1),1))


        d = np.random.uniform(low=0,high=1)
        if d1 >= d:
            Link[x1,y1,u1] = Linkprop
        else:
            Link[x1,y1,u1] = Link[x1,y1,u1]

在程序开始时，我调用一个名为“randomize”的例程来生成 K 个具有小虚部的随机酉复数，并将它们存储在一个名为 Cnum 的长度为 K 的数组中。在同一个例程中，我还通过我的链接数组并将每个元素设置为随机酉复数。代码如下。

def randommatrix():
    global Cnum
    global Link
    for i1 in range(0,K):
        C1 = np.random.normal(0,1)
        Cnum[i1] = np.cos(C1) + 1j*np.sin(C1)
        Cnum[i1+K] = np.cos(C1) - 1j*np.sin(C1)
    for i3,i4 in itertools.product(range(0,n),range(0,n)):
        C2 = np.random.uniform(low=0, high = 2*np.pi)
        Link[i3,i4,0] = np.cos(C2) + 1j*np.sin(C2)
        C2 = np.random.uniform(low=0, high = 2*np.pi)
        Link[i3,i4,1] = np.cos(C2) + 1j*np.sin(C2)

在迭代例程中使用以下例程来获取具有较小虚部的随机复数(通过检索我们之前生成的 Cnum 数组的随机元素)。

def unitary():
    I1 = np.random.randint((0),(2*K-1))
    mat = Cnum[I1]
    return mat

这是迭代例程的用途示例。我编写了一个名为 placette 的例程，它计算给定 B0 的平均 placette(链接变量的 1 x 1 闭环的实部)。迭代例程用于生成独立于先前配置的新字段配置。在我们获得新的字段配置后，我们计算所述配置的 placette。然后我们使用 while 循环重复这个过程 j1 次，最后我们得到平均 placette。

def Plq(j1,B0):
    i5 = 0
    Lboot = np.zeros(j1)
    while i5<j1:
        iteration(25000,B0)
        Linkrxp1 = np.roll(Link,-1, axis = 0)
        Linkrtp1 = np.roll(Link, -1, axis = 1)
        c0 = np.real(Link[:,:,0]*Linkrxp1[:,:,1]*Linkrtp1[:,:,0].conj()*Link[:,:,1].conj())
        i5 = i5 + 1

我们需要在运行任何东西之前定义一些变量，所以这是我在定义任何例程之前定义的初始变量

K = 20000
n = 50
a = 1.0
Link = np.zeros((n,n,2),dtype = complex)
Cnum = np.zeros((2*K), dtype = complex)

此代码有效，但速度非常慢。有没有一种方法可以使用多处理或其他方法来加快速度？

Answer 1

你应该使用 cython和 c 数据类型。 Another cython link .它专为快速计算而构建。

你可以使用 multiprocessing ，可能，在两种情况之一。如果您有一个对象需要多个进程共享，您需要使用管理器(请参阅多进程链接)、锁和数组在进程之间共享对象。但是，不能保证这会导致速度提高，因为每个进程都需要锁定链接以保证您的预测，假设预测受到链接中所有元素的影响(如果一个进程同时修改另一个进程的元素正在对元素进行预测，则预测不会基于最新信息)。

如果您的预测没有考虑其他元素的状态，即它只关心一个元素，那么您可以将 Link 数组分成多个段并将 block 分配给进程池中的多个进程，以及何时done 将这些段组合回一个数组。这肯定会节省时间，而且您不必使用任何额外的多处理机制。