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

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

接下來(lái)我們就用這個(gè)作為例子分析一下寫(xiě)高效的 Matlab 代碼需要注意的一些問(wèn)題，先把代碼貼出來(lái)吧：

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);

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

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

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

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

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

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

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

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

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

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

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

原來(lái)代碼里其實(shí)就是用?bsxfun?做了我們剛才用?repmat?做的事情。下面的代碼對(duì)比了三種方法：

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

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

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

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

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

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

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

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

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

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

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

總結(jié)

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

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