如何用GPU加速Matlab中Misanthrope过程的蒙特卡洛模拟并行计算

问题背景 #

我正在研究带周期性边界条件的Misanthrope过程的扩散行为：基于含M个格点的一维离散晶格，N个粒子存在排斥相互作用，每个格点可容纳任意数量粒子。

我在个人工作站（4并行线程，无GPU）上实现了Matlab仿真脚本，但速度无法满足需求。目前可使用的共享硬件配置为：Intel Core i9-7900X（10核20线程）、128GB内存、带12GB显存的Nvidia Titan V GPU，设备安装Matlab 2018a。

代码结构 #

仿真采用嵌套循环结构，其中轨迹循环（n_traj）和时间循环（t_f）是计算量最大的部分：

loop over interaction strength parameter (for loop beta)
    loop over N                          (parallelize with parfor)
        loop over trajectories           (for loop n_traj)
            loop over time               (while loop t_f)

核心疑问 #

我并非GPU专家，想咨询：是否将n_traj等变量转为Matlab的GPUarray能提升仿真速度？

测试脚本（参数建议：n_traj=100，t_f=30） #

Repulsion=[6 5 2 1 0.5];
%--- lattice size
L=2; M=10; x=1:L*M;
%--- declares hopping rates  
k0=0.1; X=2; k_r=exp(X)*k0;   k_l=exp(-X)*k0; 
%--- non-interacting single particle drift and diffusion
D=(k_l+k_r)/2; v=(k_r-k_l);
%--- pre-allocates interacting particles drift and diffusion  
driftMEx=zeros(L*M,length(Repulsion)); 
diffuMEx=zeros(L*M,length(Repulsion));
%--- gets ensemble size and time duration
prompt=sprintf('Give the value No of trajectories:\n n_traj = ');
n_traj = input(prompt);
prompt=sprintf('\nGive the value time lenght:\n t_f = ');
t_f = input(prompt);

j=1;
for beta=Repulsion
    parfor i=1:L*M 
            [driftMEx(i,j),diffuMEx(i,j)]=simulation(i,M,X,k0,beta,n_traj,t_f);
    end 
    j=j+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                        Functions used in script                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%----       Montecarlo Simulations for Misanthrope Process      ----%%%%
function [drift,diffusion] = simulation(N,M,X,k_o,beta,n_traj,t_f)
%%%%%%%%%%%%%%%%%%%%%%% parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k_r=exp(X)*k_o; k_l=exp(-X)*k_o; 
s_min=1; % lowest state in boundary conditions
%%%%%%%%%%%%%%%%%%%%%%%  Pre-allocation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_final=zeros(n_traj,1);
P=zeros(2*N+1,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for jj=1:n_traj
    tc=0; % current time starts at 0
    state0=1:N;
    state0=mod(state0,M);
    state0(state0==0)=M;
    state=histcounts(state0,0.5+(0:M));
    bc_counter=0;
    hop=1;
    while tc<t_f  % loop over time
        % findMisanthrope finds index of where motors are this state
        N_pos=findMisanthrope(state,N);  
        T=zeros(N,M); % reset transition matrix
        % temporarily assign state of hopping fowards and backwards
        N_plus=N_pos+1;
        N_minus=N_pos-1;
        nn=histcounts(N_pos,0.5+(0:M));
        % finds possible moves
        for qq=1:N     
            % findMisanthrope finds index of where motors are this state
            N_pos=findMisanthrope(state,N);              
            if N_pos(qq)==1
                N_minus(qq)=M;
            end            
            if N_pos(qq)==M
                N_plus(qq)=1;
            end
            T(qq, N_minus(qq))=k_l; % rate of hopping left
            T(qq,N_plus(qq))=k_r;   % rate of hopping right
            
            alpha=0;
            bb=exp(-beta.*nn);
            aa=exp(alpha.*(nn-1));

            temp_qq=N_pos(qq);
            T(qq,setdiff(1:end,temp_qq))=T(qq,setdiff(1:end,temp_qq)).*bb(setdiff(1:end,temp_qq)).*aa(temp_qq);           
        end

        for mm=2:2:2*N
            P(mm)=T(mm/2,N_minus(mm/2))/sum(T(:)) + P(mm-1);
            P(mm+1)=T(mm/2,N_plus(mm/2))/sum(T(:)) + P(mm);
        end
        tau=(-log(rand))/(sum(T(:)));  % find time until next hop
        tc=tc+tau; % update current time
        R=rand; % pick another random number from 0 to 1 (r_2)
        iii=find(P<=R,1,'last'); % find which entry in vector correspons to r_2
        N_hopp=round(iii/2);  % which motor changes position
            state(N_pos(N_hopp))=state(N_pos(N_hopp))-1; % annihilate motor to hop
            N_plus=N_pos(N_hopp)+1;
            N_minus=N_pos(N_hopp)-1;
            
            if mod(iii/2,1)~=0   % if odd then motor hops left
                if N_pos(N_hopp)==s_min
                    N_minus=M;
                    bc_counter=bc_counter-1; % take note of how many times a motor crosses the boudary
                end
                state(N_minus)=1+state(N_minus); % create motor where it hops to
            else   % even then hops right
                if N_plus>M
                    N_plus=s_min;
                    bc_counter=bc_counter+1;
                end               
                state(N_plus)=1+state(N_plus);       
            end
        hop=hop+1; % advance hop index
    end
    N_f=findMisanthrope(state,N); % final position of motors
    delta=N_f-state0; % find change in state for drift calculation
    x_final(jj)=(sum(delta)+M*bc_counter); % add boundary condition info on. To find total distance travelled
end
% calculates diffusion
ave=mean(x_final);
diff_n(1:n_traj)=(x_final(1:n_traj)-ave).^2;
diffusion=sum(diff_n)/(2*n_traj*t_f);
% calculates drift
drift=sum(x_final(1:n_traj)/t_f)/n_traj;
end


%%%%----     Find indices of particles for Misanthrope Process     ----%%%%
function [rr]=findMisanthrope(tt,m)

    if 1*isempty(tt(tt>=2))==0    
        ttt=tt(tt~=0);
        rrr=find(tt);
        rr=zeros(1,m);
        rrt=cell(1,length(rrr));    
        for i=1:length(rrr)
            rrt{i}=rrr(i)*ones(1,ttt(i));
        end
        rr(1:m)=cell2mat(rrt);    
    else
        rr=find(tt);    
    end

end

优化建议 #

1. GPU加速的可行性 #

将n_traj相关的循环转为GPUarray是可行的，核心原因是轨迹之间完全独立，符合GPU擅长的批量并行计算场景。但需要注意：

• 单轨迹内的时间循环（t_f）是串行过程，GPU对这类串行逻辑加速有限，但可以将多条轨迹的时间循环同时放在GPU上运行，实现批量并行。
• Matlab 2018a的GPU支持存在局限性，部分函数（如cell操作、部分版本的histcounts）无法直接在GPU上运行，需要适配改造。

2. 具体GPU改造步骤 #

• 变量GPU化：将x_final、state等轨迹相关变量初始化为gpuArray，例如x_final = gpuArray.zeros(n_traj,1);。
• 替换GPU不兼容函数：
  • findMisanthrope中的cell操作无法在GPU上运行，重写为向量化实现，比如用repelem(find(tt), tt(tt>0))替代cell拼接逻辑。
  • 随机数生成替换为GPU版本：rand()改为gpuArray.rand()，log(rand())改为log(gpuArray.rand())。
  • 若histcounts在GPU上不支持，手动实现基于GPU的直方图统计（比如用accumarray的GPU版本）。
• 减少数据传输：尽量让所有计算在GPU上完成，仅在最终计算漂移和扩散时将结果传回CPU，避免频繁的gather()操作。

3. CPU并行的前置优化 #

在尝试GPU之前，可先优化CPU并行效率，充分利用i9的20线程：

• 将parfor的层级调整到n_traj循环，或者修改simulation函数，让它接受批量轨迹参数，让parfor处理更多并行任务。
• 优化findMisanthrope函数，当前的cell循环效率极低，改为向量化实现（如repelem），即使在CPU上也能大幅提速。

内容的提问来源于stack exchange，提问作者Jared Lo