Matlab 在 Research 中用得非常多，確實也是非常方便實用，只是有一個問題就是寫 Matlab 代碼的時候經(jīng)常需要用一些比較奇怪獨特的方式來思考和處理問題，否則寫出來的代碼雖然同樣能工作，但是速度上可能會差上幾百幾千倍。這里有幾個關(guān)鍵詞：向量化、緩存、稀疏性等。不過由于 Matlab 在這方面確實“問題”比較大，所以關(guān)于如何寫更高效的 Matlab 代碼的文章也已經(jīng)非常多了。但是剛巧最近 bdahz 小朋友問我一段 Matlab 代碼為什么又短又奇怪但是速度又快，所以我覺得正好可以拿這個作為一個例子，羅嗦一下 Matlab 編程時候的一些注意事項，也許會對其他人也有所幫助。

這次我們的例子是?K-medoids 算法，我在以前的 blog 中也介紹過。簡單地來說就是 K-means 算法的一個變種，只是在選取中心的時候 K-means 是直接計算所有點的平均值，而 K-medoids 則要求中心點必須是數(shù)據(jù)點中某一個，所以 K-means 的優(yōu)化如果是一個數(shù)值計算問題的話，K-medoids 應(yīng)該屬于離散優(yōu)化，通常離散優(yōu)化需要窮舉搜索來求解，所以計算上會更難一些。不過實現(xiàn)得好的話，也是可以比較高效的，比如這個版本的 Matlab K-medoids。

接下來我們就用這個作為例子分析一下寫高效的 Matlab 代碼需要注意的一些問題，先把代碼貼出來吧：

function [label, energy, index] = kmedoids(X,k)

% X: d x n data matrix

% k: number of cluster

% Written by Mo Chen (sth4nth@gamil.com)

v = dot(X,X,1);

D = bsxfun(@plus,v,v')-2*(X'*X);

n = size(X,2);

[~, label] = min(D(randsample(n,k),:));

last = 0;

while any(label ~= last)

[~, index] = min(D*sparse(1:n,label,1,n,k,n));

last = label;

[val, label] = min(D(index,:),[],1);

end

energy = sum(val);

從這段代碼里我們可以看到寫高效的 Matlab 代碼的首要注意事項是：把代碼寫得晦澀難懂……呃，開各玩笑^_^bb，不過也確實是這樣，其實這樣的問題在各種語言中都是存在的：教學(xué)或者示例用的代碼通常和真正實際項目中的代碼差別很大，實際中往往摻雜各種錯誤處理呀邊界處理之類的，變得很復(fù)雜；不過代碼清晰度最大的敵人往往還是優(yōu)化。為了讓代碼運行效率更高效所做的各種努力幾乎都會很嚴(yán)重或者非常嚴(yán)重地破壞代碼的可理解性，使得原本很清晰的算法變得面目全非。

通常人們解決這類問題的辦法就是把一些通用的優(yōu)化機制總結(jié)起來，實現(xiàn)到編譯器里面去，讓編譯器來做這些 dirty work。就 C/C++ 來說的話，現(xiàn)在的編譯器雖然離完美還差的很遠(yuǎn)，但是在優(yōu)化方面也算是已經(jīng)非常強大了。可惜的是 Matlab 在這方面似乎做得不是很好——雖然 Matlab 嚴(yán)格地來說是沒有編譯器的。比如說，Matlab 里面用?for?循環(huán)是非常慢的，導(dǎo)致大家都不太敢用?for?循環(huán)，于是 Matlab 后來說提供了?for?的 JIT 機制，據(jù)說加快速度，但是似乎結(jié)果仍然是非常慢。所以沒辦法，幸運的是 Matlab 代碼通常都是比較短的。

回到我們的例子上，首先看第 5、6 兩行，這兩行做的事情實際上就是計算所以數(shù)據(jù)點之間的 pair-wise distance，放在變量?D?里。由于 pair-wise distance 在算法中要被用到很多次，并且是不會變化的，所以一開始把它計算并存儲下來后面直接用，這是所有語言里都通用的一個加速方法，或者也可以說成是空間換時間，因為如果數(shù)據(jù)量比較大 pair-wise distance 矩陣在內(nèi)存中無法存下來的話，就沒有足夠的空間來換時間了。

然后我們來看看這個 pair-wise distance 矩陣是怎么計算出來的。首先第 5 行?dot?函數(shù)參考 Matlab 的幫助文檔就知道計算了矩陣?X?每一列和自己的點乘。然后第 6 行用了一個奇怪的函數(shù)?bsxfun。讓我們先忽略這個函數(shù)，來看一下兩個點??和??之間的距離應(yīng)該是怎么計算的，定義如下：

但是由于我們這里只需要比較距離之間的相對大小，所以可以省略一個開平方根的計算，使用“平方距離”：

當(dāng)然，我們說了，在 Matlab 里用?for?循環(huán)來計算是很慢的，所以我們要用向量化的方法來計算，可以這樣寫?sum((x-y).^2)。但是這里的問題是，我們要計算很多點之間的 pair-wise distance，雖然每一對點之間的距離可以這樣算的話，要計算所有點之間的距離，好像仍然無法避免兩重?for?循環(huán)來遍歷所有的點。但是那樣又會很慢了，所以我們需要更加深層次的向量化，首先展開距離公式

這樣把距離的計算分成了三個部分，前面兩個部分都是計算向量的 norm （的平方），而第三個部分是計算向量內(nèi)積。這樣的形式的好處是可以方便地對一堆點同時進行計算：例如，對于矩陣??的每一列的 norm 平方，就可以用我們剛才提到的?dot?函數(shù)一次算出來，也是代碼中第 5 行干的事情。接下來是內(nèi)積，這個也簡單，通過矩陣乘法的公式就可以知道，如果??的話，那么

其中??是矩陣??的第??列。所以一次矩陣乘法?X'*X?就可以把所有 pair-wise 內(nèi)積全部算出來，不用任何循環(huán)。所以接下來只要把三個部分加起來就可以了，不過這里還有一個問題：雖然?X'*X?是得到的一個形狀合適的矩陣，但是?dot(X,X)?得到的卻是一個向量。為了看得更清楚一點，我們分別用?、和??表示 pair-wise distance 計算的三個部分，按理說應(yīng)該計算得到三個形狀相同的矩陣，然后相加起來：

顯然?，而?，所以對于??來說，列坐標(biāo)是無關(guān)緊要的，如果記之前?dot?得到的結(jié)果向量為?v?的話，?應(yīng)該是向量?v?按列不斷重復(fù)而得到的矩陣；類似的，?應(yīng)該是?v?轉(zhuǎn)置之后按行重復(fù)得到的矩陣。在 Matlab 中經(jīng)常需要這樣的操作，用?repmat?即可完成，所以，下面的代碼實際上就可以計算 pair-wise distance 矩陣：

v = dot(X,X); D = repmat(v, length(v),1) + repmat(v', 1, length(v)) - 2*(X'*X);

這里又碰到了一個空間換時間的問題：由于我們希望用向量化的方式“同時”計算所有點對的距離，所以我們需要把?v?擴張成??和??這兩個矩陣，需要的存儲空間從??變到了?，并且存儲的都是重復(fù)的元素，如果用?for?循環(huán)一個一個地計算的話，這些多余的空間當(dāng)然是可以避免的，但是 Matlab 的?for又很慢。不過由于這個問題出現(xiàn)得非常多，于是 Matlab 提供了一個解決方案：bsxfun。詳細(xì)的文檔可以看 Matlab 的幫助，講得很清楚，簡單地來說，bsxfun?就是對矩陣的每個元素做同一個操作，基本等價于于寫一些?for?來對矩陣元素做計算，不同的是速度快了許多許多倍。另外還有一個特點就是傳給bsxfun?的矩陣如果某一個維度上 size 是 1 的話，在那個維度上它會根據(jù)傳進來的其他矩陣做“重復(fù)擴展”，所做的事情和我們?nèi)巳庥?repmat?是一樣的，只是實現(xiàn)方式并不是這樣，它并不會生成臨時矩陣，所以在內(nèi)存方面絕對占有。

原來代碼里其實就是用?bsxfun?做了我們剛才用?repmat?做的事情。下面的代碼對比了三種方法：

function test_dist(m, n)

X = rand(m,n);

D = use_bsxfun(X);

D = use_repmat(X);

D = use_for(X);

end

function D=use_bsxfun(X)

v = dot(X,X);

D = bsxfun(@plus,v,v')-2*(X'*X);

end

function D=use_repmat(X)

v = dot(X,X);

D = repmat(v,length(v),1) + repmat(v',1,length(v)) - 2*(X'*X);

end

function D=use_for(X)

D = zeros(size(X,2));

for i=1:size(D,1)

for j=i+1:size(D,2)

D(i,j) = sum((X(:,i)-X(:,j)).^2);

end

end

D = max(D,D');

end

用 Matlab 的 Profiler 運行一下（順便說一下，Matlab 的 Profiler 是非常好用的工具，也是提升代碼性能的重要工具，用善用），在我這里，bsxfun、repmat和用循環(huán)的方式的運行時間（m=1000,n=1000）分別是 0.26、0.18 和 8.22。循環(huán)比?repmat?慢了近 50 倍，bsxfun?和?repmat?速度差不多，但是內(nèi)存更省一些，一般推薦使用?bsxfun。

然后是第 8 行，先是用?randsample?隨機選出?k?個點作為初始 center，然后為每個數(shù)據(jù)點計算 label：也就是找出它們與?k?個 center 距離最近的那個所對應(yīng)的 index。這也是用向量化的方法一次性計算的，因為 Matlab 的?min?函數(shù)能夠支持向量化操作，事實上 Matlab 的大多數(shù)基本函數(shù)都支持向量化操作，多看一下文檔會有好處。

然后是第 11 行，這一行的目的是根據(jù)每一類的數(shù)據(jù)點重新選點每類的中心點，這一步中就是 K-means 和 K-medoids 不同的地方：K-medoids 由于要求類中心必須是數(shù)據(jù)點中的某一個，所以這里需要用遍歷搜索的方法：遍歷該類中的所有數(shù)據(jù)點，選中最優(yōu)的中心。這里最優(yōu)的定義是：該中心到該類的其他點的距離之和最小，這個是和 K-means 的定義一致的。不過從代碼里來看，這里顯然又用了向量化的方法而不是循環(huán)來處理了搜索。

讓我們來看一下代碼里是怎么做的：代碼里的?sparse?函數(shù)（具體用法請參考 Matlab 幫助）構(gòu)造了一個??的稀疏矩陣，不妨?xí)簳r記為?，如果第??個數(shù)據(jù)點屬于第??類的話，那么?，否則等于?。然后用 pair-wise distance 矩陣??去乘上?，得到一個??的矩陣暫時記為?。來看一下?，它是??的第??行和??的第??列內(nèi)積的結(jié)果。?的第??列標(biāo)記了所有屬于第??類的點，其他位置全部是零，因此這樣內(nèi)積的結(jié)果就是所有第??類中的數(shù)據(jù)點到數(shù)據(jù)點??的距離之和。因此，對于第??類來說，只要求得??的第??列中數(shù)值最小的那個下標(biāo)對應(yīng)的數(shù)據(jù)點，即是最優(yōu)的中心點，而?min?函數(shù)是可以對于一個矩陣所有列同時求最小的，也就是代碼中該行達到的目的。

這里除了向量化之外還有一個注意事項就是稀疏矩陣。稀疏矩陣并不一定是很高效的，比如對里面的元素進行下標(biāo)隨機訪問就會很慢，但是有許多其他操作則可以很快（如果用了合適的函數(shù)的話），比如矩陣相乘、矩陣遍歷（尋找非零元素或者尋找最大、最小值等）、解方程、求特征向量和特征值。比如說求特征向量，如果矩陣是稀疏的，那么可以用?eigs?來進行求解，它的一個優(yōu)點是可以只求想要的幾個解，而不像?eig?那樣必須把所有解全部求出來，并且由于它是用迭代法，其中主要涉及到一些矩陣向量乘積之類的，用稀疏矩陣進行運算也會很快。當(dāng)然迭代法的缺點就是可能誤差比直接求解更大一些，數(shù)值穩(wěn)定性也更差一些。另外就是當(dāng)數(shù)據(jù)矩陣本身維度非常大但是又非常稀疏的時候，用稀疏矩陣非常節(jié)約內(nèi)存。下面是一個簡單的測試?yán)?#xff08;運行時間在注釋里）：

X = sprand(300, 2000, 0.01);

S = X'*X;

[A B] = eigs(S); % time: 0.16

[A B] = eigs(full(S)); % time: 0.91

[A B] = eig(full(S)); % time: 13.76

運行時間差異還是比較清楚的，也就不用我多解釋了。不過有一點需要注意的是，eigs?在沒有指定個數(shù)的情況下默認(rèn)是只求 6 個特征值和特征向量的，所以和?eig?把所有的特征值特征向量全部求出來其實也并不是可以直接比較的。但是許多時候我們需要的都不是所有的特征向量和特征值，如果真的需要全部的話，用?eigs?來計算可能并不是一個合適的選擇，屆時可以自己嘗試和比較一下。由于求特征值和特征向量在各種基于 Graph 的方法（像 Laplacian Eigenmaps、Laplacian Regularized Least Square 等）中用得非常多，所以這些還是很有用的。這里可以順便簡單說一下稀疏矩陣的存儲。當(dāng)然是有各種存儲方式的，比較基本的比如按行存儲：矩陣是一個鏈表，把矩陣的每一行鏈起來，而每一行也是一個鏈表，把該行的非零元素鏈起來；類似的有按列存儲；此外可能還有完全按元素存儲，可以看成一個表格，如果(i,j)?位置有非零值的話，就索引到該值。根據(jù)不同的存儲方式，計算效率也會不一樣，比如一個按行存儲的矩陣乘以一個按列存儲的矩陣，就可以很快，因為矩陣相乘的計算方式就是左邊的行和右邊的列做內(nèi)積；但是如果反過來，一個按列存儲的矩陣乘以一個按行存儲的矩陣的話，就會比較麻煩了。Matlab 里做得比較好的是把稀疏矩陣搞得很透明，你不用關(guān)心它底層到底是怎么存儲的，大多數(shù)時候就像使用普通矩陣一樣用就 OK 了，并且性能也挺不錯。

回到我們原來的代碼，第 13 行和第 8 行是一樣的，只是現(xiàn)在使用計算出來的中心而不是隨機選出來的。剩下的就沒有什么好解釋的了。這里的 stop condition 是中心點不再變動，intuitively 想一想對于 K-medoids 來說這樣的 stop condition 似乎在一些比較特殊的情況下可能會出現(xiàn)來回振蕩不停下來的結(jié)果，不過那個不是我們今天要關(guān)注的問題了。

最后總結(jié)一下，要寫出高效的 Matlab 代碼的一些注意事項：

• Profiler: Matlab 的 Profiler 是非常好用的，要善于利用這個工具，同所有其他編程語言一樣，找準(zhǔn) bottleneck 是進行優(yōu)化的最重要的一步，如果只是想當(dāng)然地去搞的話可能浪費了大把的精力又沒有把性能改善多少而且還把代碼搞得一團糟。
• Sparsity: 如果問題是有稀疏性質(zhì)的，那么可以嘗試一下用稀疏矩陣和配套的那些操作。
• Vectorization: 向量化可以說是 Matlab 編程的一個特點，就好像函數(shù)式編程總是一堆?map?呀filter?呀?reduce?呀之類的一樣。用好向量化是改善 Matlab 性能的關(guān)鍵。要多嘗試和練習(xí)，逐漸習(xí)慣向量化的思維方式。特別是矩陣相乘呀、分塊之類的要熟練，例如我們在介紹代碼第 11 行的時候構(gòu)造的那個矩陣?，通常稱作 indicator matrix，元素只有 0 和 1，一般用于表示哪些元素被選出來了。這個矩陣不論是在計算上還是在公式推導(dǎo)上都經(jīng)常被用到。

from:?http://freemind.pluskid.org/programming/recipes-for-faster-matlab-code/

總結(jié)

以上是生活随笔為你收集整理的让书写的Matlab代码运行更快 Recipes for Faster Matlab Code的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。