向量化代碼 - 如何減少MATLAB計算時間

Question

我有這樣一段代碼

N=10^4; 
for i = 1:N 
    [E,X,T] = fffun(); % Stochastic simulation. Returns every time three different vectors (whose length is 10^3). 
    X_(i,:)=X; 
    T_(i,:)=T; 
    GRID=[GRID T]; 
end 
GRID=unique(GRID); 
% Second part 
for i=1:N 
for j=1:(kmax) 
    f=find(GRID==T_(i,j) | GRID==T_(i,j+1)); 
    s=f(1); 
    e=f(2)-1; 

counter(X_(i,j), s:e)=counter(X_(i,j), s:e)+1; 
end 
end 

代碼執行一個隨機過程（它由10^3點的事件，在離散的時刻發生的N個不同的模擬（T （第二部分）作爲時間的函數，我想知道在一個特定的狀態下有多少個模擬（X取值在1到10之間）。我有這樣的想法：創建一個具有所有時刻的網格矢量，然後，在模擬上循環，循環發生某些事情的時間步並遞增與此p相對應的所有計數器餘數時間的關鍵片段。

然而，這第二部分是非常重的（我的意思是在一個標準的四核CPU上處理幾天）。它不應該。 是否有任何想法（可能是以更高效的方式比較向量）來縮短CPU時間？

這是一個獨立 'SECOND_PART'

N=5000; 
counter=zeros(11,length(GRID)); 

for i=1:N 
    disp(['Counting sim #' num2str(i)]); 
    for j=1:(kmax) 
     f=find(GRID==T_(i,j) | GRID==T_(i,j+1),2); 
     s=f(1); 
     e=f(2)-1; 

     counter(X_(i,j), s:e)=counter(X_(i,j), s:e)+1; 

    end 
end 

counter=counter/N; 
stop=find(GRID==Tmin); 
stop=stop-1; 
plot(counter(:,(stop-500):stop)') 

具有相關聯的僞數據（filedropper.com/data_38）。在真實的情況下，矩陣有2x行和10x列。

Answer 1

這裏是我的理解：

T_是從N個模擬時間步的矩陣。
X_是在那些模擬中的模擬狀態矩陣T_。

所以，如果你這樣做：

[ut,~,ic]= unique(T_(:)); 

你得到ic這是指數在T_所有獨特元素的矢量。那麼你可以寫：

counter = accumarray([ic X_(:)],1); 

並得到counter沒有。的行作爲你獨特的時間步調，沒有。作爲X_（它們都是，必須是整數）中的唯一狀態列。現在你可以說，對於每個時間步ut(k)模擬處於狀態m的時間數是counter(k,m)。

在您的數據中，m和k的值大於1的唯一組合是(1,1)。

---

編輯：

從下面的評論，我知道你記錄所有的狀態變化，以及何時發生的時間步驟。然後每次模擬改變一個狀態時，你想從所有模擬中收集所有狀態，並計算每種類型有多少個狀態。

這裏的主要問題是，你的時間是連續的，所以基本上在T_每個元素都是獨特，你有超過一百萬的時間步長遍歷。完全矢量化這樣的過程將需要大約80GB的內存，這可能會卡住你的電腦。

所以我尋找矢量化和循環的時間步驟的組合。首先，我們要找出所有獨特的間隔和預分配counter：

ut = unique(T_(:)); 
stt = 11; % no. of states 
counter = zeros(stt,numel(ut));r = 1:size(T_,1); 
r = 1:size(T_,1); % we will need that also later 

然後我們循環遍歷所有ut元素，每一次找相關的時間步長T_在所有模擬的量化方法。最後，我們使用histcounts計算所有的狀態：

for k = 1:numel(ut) 
    temp = T_<=ut(k); % mark all time steps before ut(k) 
    s = cumsum(temp,2); % count the columns 
    col_ind = s(:,end); % fins the column index for each simulation 
    % convert the coulmns to linear indices: 
    linind = sub2ind(size(T_),r,col_ind.'); 
    % count the states: 
    counter(:,k) = histcounts(X_(linind),1:stt+1); 
end 

這在我的電腦模擬1000大約需要4秒鐘，因此將其添加到略多於一個小時的整個過程。不是很快......

您也可以嘗試一個或兩個以下調整的擠運行時間多一點點：

1. 正如你can read here，accumarray似乎在小數組更快的工作，然後histcouns。所以可能想切換到它。

2. 此外，直接計算線性索引是比sub2ind更快的方法，因此您可能需要嘗試。

實施上述循環這些建議，我們得到：

R = size(T_,1); 
r = (1:R).'; 
for k = 1:K 
    temp = T_<=ut(k); % mark all time steps before ut(k) 
    s = cumsum(temp,2); % count the columns 
    col_ind = s(:,end); % fins the column index for each simulation 
    % convert the coulmns to linear indices: 
    linind = R*(col_ind-1)+r; 
    % count the states: 
    counter(:,k) = accumarray(X_(linind),1,[stt 1]); 
end 

在我的電腦切換到accumarray和或刪除sub2ind獲得略有改善，但它不是一致的（使用timeit，以測試在ut中有100個或1K個元素），所以你最好自己測試一下。但是，這仍然很長。你可能要考慮

---

有一件事是想你的離散時間步長，所以你必須遍歷更獨特的元素。在你的數據中，約8％的時間間隔小於1.如果你可以認爲這個時間足夠短，那麼你可以圍繞你的T_，只得到約12.5K的獨特元素，一分鐘循環。你可以在0.1的時間間隔（小於1％的時間間隔）內完成相同的操作，並獲得122K個元素循環，大約需要8個小時...

當然，以上所有時間都是粗略估計使用相同的算法。如果你選擇做選擇循環時間可能有更好的方法來解決這個問題。