在Python中使用Metropolis算法的蒙特卡洛算法速度极慢

Question

我正在尝试用Python实现一个简单的蒙特卡洛(我对Python来说还是个新手)。来自C语言的我可能走错了路，因为我的代码对于我所要求的来说太慢了:我有一个潜在的硬球体(参见代码中的V_pot(r) )，适用于60个3d粒子和周期边界条件(PBC)，所以我定义了以下函数

import timeit
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from numpy import inf
#
L, kb, d, eps, DIM = 100, 1, 1, 1, 3
r_c, T = L/2, eps/(.5*kb)
beta = 1/(kb*T)
#
def dist(A, B):
    d = A - B
    d -= L*np.around(d/L)
    return np.sqrt(np.sum(d**2))

#
def V_pot(r):
    V = -eps*(d**6/r**6 - d**6/r_c**6)
    if r > r_c:
        V = 0
    elif r < d:
        V = inf

    return V

#
def ener(config):
    V_jk_val, j = 0, N
    #
    while (j > 0):
        j -= 1
        i = 0
        while (i < j):
            V_jk_val += V_pot(dist(config[j,:], config[i,:]))
            i += 1

    #
    return V_jk_val

#
def acc(en_n, en_o):
    d_en = en_n-en_o
    if (d_en <= 0):
        acc_val = 1
    else:
        acc_val = np.exp(-beta*(d_en))

    return acc_val

#

然后，从配置开始(其中数组的每一行都表示3D粒子的坐标)

config = np.array([[16.24155657, 57.41672173, 94.39565792],
       [76.38121764, 55.88334066,  5.72255163],
       [38.41393783, 58.09432145,  6.26448054],
       [86.44286438, 61.37100899, 91.97737383],
       [37.7315366 , 44.52697269, 23.86320444],
       [ 0.59231801, 39.20183376, 89.63974115],
       [38.00998141,  3.84363202, 52.74021401],
       [99.53480756, 69.97688928, 21.43528924],
       [49.62030291, 93.60889503, 15.73723259],
       [54.49195524,  0.6431965 , 25.37401196],
       [33.82527814, 25.37776021, 67.4320553 ],
       [64.61952893, 46.8407798 ,  4.93960443],
       [60.47322732, 16.48140136, 33.26481306],
       [19.71667792, 46.56999616, 35.61044526],
       [ 5.33252557,  4.44393836, 60.55759256],
       [44.95897856,  7.81728046, 10.26000715],
       [86.5548395 , 49.74079452,  4.80480133],
       [52.47965686, 42.831448  , 22.03890639],
       [ 2.88752006, 59.84605062, 22.75760029],
       [ 9.49231045, 42.08653603, 40.63380097],
       [13.90093641, 74.40377984, 32.62917915],
       [97.44839233, 90.47695772, 91.60794836],
       [51.29501624, 27.03796277, 57.09525454],
       [10.30180295, 21.977336  , 69.54173272],
       [59.61327648, 14.29582325, 11.70942289],
       [89.52722796, 26.87758644, 76.34934637],
       [82.03736088, 78.5665713 , 23.23587395],
       [79.77571695, 66.140968  , 53.6784269 ],
       [82.86070472, 40.82189833, 51.48739072],
       [99.05647523, 98.63386809,  6.33888993],
       [31.02997123, 66.99709163, 95.88332332],
       [97.71654767, 59.24793618,  5.20183793],
       [ 6.79964473, 45.01258652, 48.69477807],
       [93.34977049, 55.20537774, 82.35693526],
       [17.35577815, 20.45936211, 29.27981422],
       [55.51942207, 52.22875901,  3.6616131 ],
       [61.45612224, 36.50170405, 62.89796773],
       [23.55822368,  7.09069623, 37.38274914],
       [39.57082799, 58.95457592, 48.0304924 ],
       [93.94997617, 64.34383203, 77.63346308],
       [17.47989107, 90.01113402, 81.00648645],
       [86.79068539, 66.35768515, 56.64402907],
       [98.71924121, 38.33749023, 73.4715132 ],
       [ 0.42356139, 78.32172925, 15.19883322],
       [77.75572529,  2.60088767, 56.4683935 ],
       [49.76486142,  3.01800153, 93.48019286],
       [42.54483899,  4.27174457,  4.38942325],
       [66.75777178, 41.1220603 , 19.64484167],
       [19.69520773, 41.09230171,  2.51986091],
       [73.20493772, 73.16590392, 99.19174281],
       [94.16756184, 72.77653334, 10.32128552],
       [29.95281655, 27.58596604, 85.12791195],
       [ 2.44803886, 32.82333962, 41.6654683 ],
       [23.9665915 , 49.94906612, 37.42701059],
       [30.40282934, 39.63854309, 47.16572743],
       [56.04809276, 30.19705527, 29.15729635],
       [ 2.50566522, 70.37965564, 16.78016719],
       [28.39713572,  4.04948368, 27.72615789],
       [26.11873563, 41.49557167, 14.38703697],
       [81.91731981, 12.10514972, 12.03083427]])

我用以下代码进行了5000个时间步长的模拟

N = 60
TIME_MC = 5000
DELTA_LIST = [d]
#d/6, d/3, d, 2*d, 3*d 
np.random.seed(19680801)
en_mc_delta = np.zeros((TIME_MC, len(DELTA_LIST)))
start = timeit.default_timer()
config_tmp = config
#
for iD, Delta in enumerate(DELTA_LIST):
    t=0

    while (t < TIME_MC):
        for k in range(N):
            RND = np.random.rand()
            config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5)
            en_o, en_n = ener(config), ener(config_tmp)
            ACC = acc(en_n, en_o)
            if (RND < ACC):
                config[k,:] = config_tmp[k,:]
                en_o = en_n


        en_mc_delta[t][iD] = en_o
        t += 1


stop = timeit.default_timer()
print('Time: ', stop-start)

遵循大都会算法的规则，用config_tmp[k,:] = config[k,:] + Delta*(np.random.random_sample((1,3))-.5)提取出接受提议的移动。

我做了一些尝试来检查代码被卡住的地方，我发现函数ener (也是因为函数dist)非常慢:它需要像~0.02s这样的东西来计算配置的能量，这意味着在~6000s附近运行完整的模拟(60个粒子，5000个建议的移动)。

它的外层只是计算不同Delta值的结果。

用TIME_MC=60运行这段代码可以让你知道这段代码(~218s)有多慢，如果用C实现，只需要几秒钟。我读到了其他一些关于如何加速Python代码的问题，但我不明白如何在这里做到这一点。

编辑：

我现在几乎可以肯定问题出在dist函数中，因为仅仅计算两个3D向量之间的PBC距离就需要大约~0.0012s，当你计算它5000*60次时，它会给出疯狂的长时间。

Answer 1

请注意，这是对原始问题的评论的部分回答。

下面是一个示例，说明“展开”numpy的函数在替换为更直接的距离计算时如何提高性能。请注意，这并未验证为等效，特别是在舍入方面。我认为，这一原则仍然适用。

import random
import time
import numpy as np

L = 100
inv_L = 0.01
vec_length = 10
repetitions = 100000

def dist_np(A, B):
    d = A - B
    d -= L*np.around(d/L)
    return np.sqrt(np.sum(d**2))

def dist_direct(A, B):
    sum = 0
    for i in range(0, len(A)):
        diff = (A[0,i] - B[0,i])
        diff -= L * int(diff * inv_L)
        sum += diff * diff
    return np.sqrt(sum)

vec1 = np.zeros((1,vec_length))
vec2 = np.zeros((1,vec_length))

for i in range(0, vec_length):
    vec1[0,i] = random.random()
    vec2[0,i] = random.random()

print("with numpy method:")
start = time.time()
for i in range(0, repetitions):
    dist_np(vec1, vec2)
print("done in {}".format(time.time() - start))

print("with direct method:")
start = time.time()
for i in range(0, repetitions):
    dist_direct(vec1, vec2)
print("done in {}".format(time.time() - start))

输出：

with numpy method:
done in 6.332799911499023
with direct method:
done in 1.0938000679016113

尝试使用平均向量长度和重复，看看最佳位置在哪里。当改变这些元参数时，我期望性能增益不是恒定的。