有關matlab循環怎麼改成矩陣運算？

01-07

矩陣里還嵌套矩陣的二重循環怎麼優化成矩陣運算，圖示為一個例子

從你給的例子中不能明確知道dmat和pop到底是矩陣還是函數，但是從你「矩陣里還嵌套矩陣」這句話中推斷dmat和pop大概都是數組，以下所有程序都只針對dmat和pop都是數組的情況，如果不是這樣你需要給出dmat函數或者pop函數的內容。因為你沒有給定size，所以假定size設置如下以便測試：

n = 3e3; popSize = 3e3; mDmat = 3e3; nDmat = 3e3; dmat = rand(mDmat,nDmat); pop = randi(min(mDmat,nDmat),popSize,n);

你原始的方法是：

totalDist = zeros(popSize,1); for p = 1:popSize d = dmat(1,pop(p,1)); for k = 2:n d = d + dmat(pop(p,k-1),pop(p,k)); end totalDist(p) = d; end % Elapsed time is 1.136523 seconds.

如果dmat和pop如假定一般都是數組，那麼這種和索引有關的向量化一般會考慮轉化為線性索引.

將二維索引轉化為一維索引（即線性索引）可以用sub2ind也可以自己根據索引的順序計算，但是sub2ind會有overhead，所以我們選擇直接計算，例如M-by-N的矩陣A中，A(i,j)的線性索引是：

i+(j-1)*M

應用到問題中，內層循環可以向量化為：

totalDist = zeros(popSize,1); for p = 1:popSize totalDist(p) = dmat(1,pop(p,1))+sum(dmat(pop(p,1:n-1)+(pop(p,2:n)-1)*mDmat)); end % Elapsed time is 0.395223 seconds.

但這樣顯然不過癮，那就把內層外層一起向量化吧：

totalDist = sum(dmat([ones(popSize,1) pop(:,1:end-1)]+(pop-1)*mDmat),2); % Elapsed time is 0.256587 seconds.

這裡有大塊的索引運算(pop(:,1:end-1))和cat運算([])，通常這樣的運算會比較消耗時間，可以考慮避免，首先將求和的兩部分分開來寫避免cat運算：

totalDist = dmat((pop(:,1)-1)*mDmat+1)+... sum(dmat(pop(:,1:end-1)+(pop(:,2:end)-1)*mDmat),2); % Elapsed time is 0.254859 seconds.

這樣雖然避免了cat運算，但是又引入了一個大塊的索引運算(pop(:,2:end)),所以可以看到性能並沒有什麼變化，可以進一步考慮消除索引運算：

totalDist2 = dmat((pop(:,1)-1)*mDmat+1)+... sum(dmat(conv2(pop,[mDmat 1],"valid")-mDmat),2); % Elapsed time is 0.227807 seconds.

可以看到性能提升了10%，另外需要注意上邊兩種分塊求和再相加的運算與放在一起求和的計算結果略有不同，浮點數計算誤差導致的，這種相對誤差很小，一般不影響運算

其實索引部分的運算還可以寫的更緊湊一些：

totalDist = sum(dmat(filter([mDmat 1],1,pop,1,2)-mDmat),2); % Elapsed time is 0.298399 seconds.

但是可以看到這種方法是以上幾種向量化方法中性能最低的，但是還是明顯要快於原始的循環方法。

不過在這個問題中向量化有個問題，就是要生成一個和pop差不多大的矩陣用於求和，當pop足夠大的時候開闢這麼大的空間也是很消耗時間的，所以重新考慮下循環的方法。由於MATLAB的矩陣存儲順序是列優先，所以沿著列方向運算通常要快一些，自然想到下邊的方法：

totalDist = dmat(1,pop(:,1))."; for k = 2:n totalDist = totalDist + dmat(pop(:,k-1)+mDmat*(pop(:,k)-1)); end % Elapsed time is 0.186153 seconds.

可以看到這樣的循環比上邊任何一種向量化方法都要快。不過你可以發現這裡的循環次數為n-1，於是可以預見，如果n遠大於popSize時上述循環性能就會不如其他方法，在我這裡測試，上述例子中其他參數不變，將n改為3e4，popSize取300，兩者相差100倍時，該方法性能與第三種向量化演算法性能差不多，而當n取3e5,popsize取30的時候這種方法就明顯比向量化方法慢了

所以說優化MATLAB程序時向量化只是一種手段，並不是目標。上邊的例子說明循環並不一定比向量化慢，當然你有可能找到一種新的向量化方法比上邊最快的那個循環在任何情況下都快，但是你不能保證每次遇到問題你都能找到這樣的向量化方法。所以優化的目標歸根結底是提升程序的性能，對於性能敏感的程序最好還是具體情況具體分析（例如這裡的pop的size），而不是一味的追求向量化。

前提是保證你得dmat支持參數為向量情況，也就是輸入標量和矢量或矢量與矢量，並且返回一個矢量結果

% 標量與矢量 d = dmat(1, pop(:,1)); % 返回d為列向量 % 矢量與矢量 d = dmat(pop(:,1), pop(:,2)); % 返回d為列向量

改矩陣運算並非沒有缺點，缺點就是得多生成一個數據，計算使用內存加倍。現在電腦內存都夠大，只要給Matlab足夠內存應該就行

n = size(pop, 1); pop_m1 = [ones(n,1), pop(:, 1:end-1)]; d = dmat(pop_m1, pop); totalDist = sum(d, 2);

憑經驗寫的，沒有運行，有哪些不對的歡迎指正