如何加速Matlab嵌套for循環，當我無法矢量化計算？

Question

我有三個大小相同的大型3D數組[41 * 141 * 12403]，在alpha，beta和ni下面的Matlab代碼中命名。從他們我需要計算另一個相同大小的3D數組，通過使用每個元素的值，通過結合無限和和定積分計算的計算，從原始矩陣獲得元素。因此，不得不使用幾個嵌套循環來進行計算。該代碼已經運行了幾個小時（！），並且仍然在外循環的第一次迭代中（這需要執行41次！根據我的計算，這樣程序將不得不運行超過兩年!!!）。我不知道如何優化代碼。請幫幫我 ！！

的代碼我使用：

z_len=size(KELDYSH_PARAM_r_z_t,1); % 41 rows 
    r_len=size(KELDYSH_PARAM_r_z_t,2); % 141 columns 
    t_len=size(KELDYSH_PARAM_r_z_t,3); % 12403 slices 

    sumRes=zeros(z_len,r_len,t_len); 

    for z_ind=1:z_len 
     z_ind  % in order to track the advancement of the calculation 
     for r_ind=1:r_len 
      for t_ind=1:t_len 
       sumCurrent=0; 
       sumPrevious=inf; 
       s=0; 

       while abs(sumPrevious-sumCurrent)>1e-6 
        kapa=kapa_0+s; %some scalar 
        x_of_w=(beta(z_ind,r_ind,t_ind).*(kapa-ni... 
         (z_ind,r_ind,t_ind))).^0.5;    
        sumPrevious=sumCurrent; 
        sumCurrent=sumCurrent+exp(-alpha(z_ind,r_ind,t_ind).* ... 
         (kapa-ni(z_ind,r_ind,t_ind))).*(x_of_w.^(2*abs(m)+1)/2).* ... 
          w_m_integral(x_of_w,m); 
        s=s+1; 
       end 

       sumRes(z_ind,r_ind,t_ind)=sumCurrent; 
      end 
     end 
    end 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function res=w_m_integral(x_of_w,m) 

    res=quad(@integrandFun,0,1,1e-6); 

    function y=integrandFun(t) 
      y=exp(-x_of_w^2*t).*t.^(abs(m))./((1-t).^0.5); 
    end 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Answer 1

選項1 - 更多的向量化

這是你正在使用一個相當複雜的模型，而不是所有的術語解釋，但一些地方仍然可以進一步被矢量化。您的alpha,beta和ni矩陣大概是靜態的並且是預先計算好的？您的s值是一個標量，kapa也可以是，所以您可以預先計算x_of_w矩陣。這會給你一個非常輕微的加速，儘管你會花費內存來獲得它 - 現在有7千1百萬點是可行的，但會需要大量的硬件。爲41行中的每一行做一次就可以減輕整體的負擔。

這留下了積分本身。 quad函數不接受向量輸入 - 這將是一場噩夢不是嗎？ - integral，Mathworks都建議您改用它。但是，如果你的整合限制在每種情況下都是一樣的，那麼爲什麼不按照老式的方式去做積分呢？計算被積函數值爲1的矩陣，計算被積函數值爲0的另一個矩陣，然後取差值。

然後，您可以編寫一個計算整個輸入空間積分的單個循環，然後測試所有矩陣元素的收斂性。製作一個面具，記下那些未收斂的面具，並重新計算增加的​​s。重複，直到所有人都收斂（或你達到迭代閾值）。

選擇2 - parallelise它

它曾經是MATLAB快得多比循環向量化操作的情況。我現在找不到它的源代碼，但我想我已經讀過，最近for循環也變得快了很多，所以根據您可用的資源，通過並行處理當前的代碼可能會獲得更好的結果。這也需要進行一些重構 - 大問題是將數據複製到工作人員時的開銷（可以通過將輸入切分爲塊並只輸入相關的數據）和循環不允許你使用某些變量，通常是覆蓋整個空間的變量。再次砍他們幫助。

但是如果你有一個2年的運行時間，你將需要一個至少100的因子我猜，這意味着一個集羣！如果你在大學或某個地方，你可能可以在500個核心羣上獲得幾天，那就去...

如果你可以在封閉形式下編寫積分，那麼它可能是適合GPU計算。這些東西可以非常快地完成某些計算類，但是您必須能夠將作業並行化並將實際計算減少到主要由加法和乘法組成的基本計算。 CUDA庫已經完成了很多工作，並且matlab has an interface to them因此閱讀了這些內容。

選項3 - 縮小範圍

最後，如果沒有上述兩個結果足夠的加速，那麼你可能不得不減少你的計算範圍。儘可能多地修剪輸入空間，並可能接受較低的收斂閾值。如果您知道在最內層的while循環（其中s計數器在內）中需要多少次迭代，那麼可能會發現減少收斂準則會減少您需要的迭代次數，從而加快迭代速度。分析器可以幫助您查看您的時間。

雖然底線是7100萬點需要一些時間來計算。只有到目前爲止您可以優化計算，但出現這種大小問題的可能性是，您將不得不拋硬件。