Monte Carlo模拟在PC上使用Numba更快,但在另一台电脑上不使用

2024-04-26 00:28:13 发布

您现在位置:Python中文网/ 问答频道 /正文

我正在研究一种算法(用于某种PCA),并根据一个指标测试其性能。通过蒙特卡罗方法对其进行评估,即对多个样本进行平均,以选择参数。 初始,我的代码由参数的简单循环和算法的Numpy函数样本组成

我听说Numba软件包加快了我的代码速度。我发现我在速度上得到了5倍的提高,只需要重新编写Numpy部分以适应Numba standarts并添加@njit装饰器。 然而,我在另一台电脑上测试了相同的改进版本(两台都是采用英特尔i7的笔记本电脑),发现在这第二台电脑上没有速度提升

我很惊讶为什么会这样

编辑:

关于这两台电脑的更多细节

PC 1:

  • Python 3.7.6
  • 麻木0.48.0
  • 英特尔酷睿i7-8565U处理器@1.80GHz(4核)
  • RAM 8开始

PC 2:

  • Python 3.7.4
  • 麻木0.49.0
  • 英特尔酷睿i7-8750H处理器@2.20GHz(6核)
  • RAM 16,0 Go

这是一个最小的可重复的例子

import numpy as np
from numba import jit
import matplotlib.pyplot as plt
import time

@jit(nopython=True)
def multivariate_normal(m,C):
    L = np.linalg.cholesky(C)
    X = np.random.randn(m.shape[0])
    return(L @ X + m)


@jit(nopython=True)
def gendata(p, d, k, n, sig_s):

    H = np.random.randn(d,k)
    u, s, vh = np.linalg.svd(H, full_matrices=False)
    H = u @ vh
    rest = np.zeros((p-d, k))
    U0 = np.concatenate((H,rest),axis=0)

    X = np.empty((p,n))
    for i in range(n):
        b = multivariate_normal(np.zeros(p),np.identity(p))
        s = multivariate_normal(np.zeros(k),(sig_s**2) * np.identity(k)) 
        x = U0 @ s + b
        X[:,i] = x

    return((X,U0))

@jit(nopython=True)
def geomspace(a,b,n):
    seq = np.linspace(np.log(a),np.log(b),n)
    return(np.exp(seq))


@jit(nopython=True)
def Afe(U,U0):
    '''
    Average Fraction of Energy of subspace U based on original subspace U0.
    '''
    Oorth,r = np.linalg.qr(U)
    afe = np.trace(Oorth.T @ U0 @ U0.T @ Oorth) / np.trace(U0 @ U0.T) 
    return(afe)
    
@jit(nopython=True)
def sim2a():
    p = 50
    d = int(0.99 * p)
    k = 5
    sig_s = np.sqrt(10)

    
    ns = geomspace(3*k,3*p,10).astype(np.int32)
    echs = 50
    
    mc_afe_pca = np.empty(ns.shape[0])
    mc_afe_sphpca = np.empty(ns.shape[0])
    
    for n in range(ns.shape[0]):
        res_pca = np.empty(echs)
        res_sphpca = np.empty(echs)
    
    
        for e in range(echs):
            
            #print(n,e)
            X,U0 = gendata(p,d,k,ns[n],sig_s)
            sX = np.diag(X.T @ X)
            
            #PCA
            u, sig, vh = np.linalg.svd(X.T, full_matrices=False)
            Upca = vh.T[:,:k]
            res_pca[e]= Afe(Upca,U0)
            
            #Spherical PCA
            Xn = X / np.sqrt(sX)
            u, sig, vh = np.linalg.svd(Xn.T, full_matrices=False)
            Usph = vh.T[:,:k]
            res_sphpca[e] = Afe(Usph,U0)
            
        mc_afe_pca[n] = np.mean(res_pca)
        mc_afe_sphpca[n] = np.mean(res_sphpca)
        
    return((mc_afe_pca,mc_afe_sphpca))

#first run for compilation before measurements
mc_afe_pca,mc_afe_sphpca = sim2a()

#Graph 
p = 50
k = 5
ns = geomspace(3*k,3*p,10).astype(np.int32)

plt.figure(figsize=(8,4))
plt.plot(ns,mc_afe_pca,'+-',label="PCA")
plt.plot(ns,mc_afe_sphpca,'x-',label="Spherical PCA")
plt.xlabel("n")
plt.ylabel("AFE")   
plt.legend()
plt.plot()

#Measure
def sim2atimer():
    start = time.time()  
    a,b= sim2a()
    end = time.time()
    print(end - start)
    
sim2atimer()

#Now, replace @jit by #@jit and reexecute the code to compare.
# The option Parallel=True was erased from @jit()

事实上,在处理这个示例时,我想删除@jit装饰器中的Parallel=True选项,这会提示未使用该选项,并修复了差异问题。

现在,我仍然不知道为什么只有一台电脑受到影响


Tags: truedefnprespltmcjitsig