如何提高M(jìn)ATLAB運(yùn)行效率

1、如何提高matlab運(yùn)行效率(1)在我們的領(lǐng)域，matlab是進(jìn)行實(shí)驗(yàn)的最主要的工具之一。我聽到過很多抱怨說matlab很慢，可是這往往不是matlab的問題，而是因?yàn)樽约撼绦驔]有寫好。matlab是我非常欣賞的一種語言，堪稱靈活易用和高效運(yùn)算的結(jié)合典范。 作為解釋語言，matlab進(jìn)行高效運(yùn)算的秘訣有幾個(gè)方面： (1) jit即時(shí)編譯技術(shù)。matlab在第一次加載運(yùn)行一個(gè)函數(shù)的時(shí)候，會(huì)在后臺(tái)對(duì)它進(jìn)行某種程度的編譯和優(yōu)化，使得后續(xù)運(yùn)行更快。因此，除了第一次運(yùn)行需要進(jìn)行語句解釋之外，后面運(yùn)行的其實(shí)是已經(jīng)放在內(nèi)存中的經(jīng)過優(yōu)化的中間碼。所以，很多時(shí)候我們可能會(huì)看到第一次運(yùn)行一個(gè)新函數(shù)，比它后面幾次的

2、運(yùn)行明顯慢一些。不過目前的jit技術(shù)還不是非常成熟，和標(biāo)準(zhǔn)的編譯語言相比還有相當(dāng)差距，僅憑這個(gè)matlab還不能稱為高效。 (2) 向量化(vectorization)。這是matlab最著名的特色。向量化配合經(jīng)過高度優(yōu)化的數(shù)值運(yùn)算引擎，是其高效運(yùn)算的最重要的基石很多matlab的使用者都明白這一點(diǎn)。能夠轉(zhuǎn)化成矩陣操作的規(guī)則運(yùn)算過程，使用matlab計(jì)算遠(yuǎn)比自己手工用c/c+實(shí)現(xiàn)高效。我自己做過很多次對(duì)比profile，matlab在關(guān)鍵的核心運(yùn)算上的實(shí)現(xiàn)可以比自己用c/c+按照標(biāo)準(zhǔn)的routine進(jìn)行實(shí)現(xiàn)快幾十倍。 其實(shí)這不完全是matlab的功勞，其實(shí)matlab是建基于blas和lapa

3、ck等的基礎(chǔ)數(shù)學(xué)運(yùn)算包的基礎(chǔ)上的。intel和amd都發(fā)布了這些東西的vendor-implementation，并且針對(duì)各自的cpu進(jìn)行了大量的優(yōu)化，這非個(gè)人之力所能企及。因此，我從來都強(qiáng)烈不建議個(gè)人用c/c+去實(shí)現(xiàn)數(shù)值運(yùn)算的關(guān)鍵代碼（學(xué)習(xí)數(shù)值分析課者除外）。 不過，對(duì)于向量化這個(gè)事情，也不宜極端化，下面的一些例子說明在哪些時(shí)候，for-loop比vectorization更適合： (a) 粗粒度的算法流程控制。比如，你要實(shí)現(xiàn)一個(gè)迭代算法，循環(huán)做一個(gè)事情直到收斂。只要循環(huán)幾次或者幾十次，但是每個(gè)循環(huán)內(nèi)部要進(jìn)行大量的復(fù)雜運(yùn)算，那么你就沒有必要費(fèi)心把這層循環(huán)給vectorize掉了。除非收斂結(jié)果

4、有某種close-form的解析解。 (b) 如果向量化可能導(dǎo)致產(chǎn)生巨型矩陣，則使用前要三思！很多情況下，向量化是一種用空間換時(shí)間的行為，就是通過把數(shù)據(jù)組織成某種方式，從而使得內(nèi)建的高效引擎能得以應(yīng)用。但是，有些時(shí)候要處理大量的數(shù)據(jù)，可能導(dǎo)致其組織過程需要耗費(fèi)額外的數(shù)百兆乃至千兆內(nèi)存空間，那么這有可能造成效率的不升反降。原因有幾個(gè)方面：第一，數(shù)據(jù)組織過程也是需要時(shí)間的，最起碼它也需要大量的操作進(jìn)行密集的內(nèi)存讀寫和用于定位的偏移量計(jì)算。這方面增加的時(shí)間 vs. 向量化節(jié)省的運(yùn)算時(shí)間，孰輕孰重，需要衡量。第二，分配額外的大塊內(nèi)存是一件非常耗時(shí)的事情，它可能導(dǎo)致虛擬內(nèi)存的使用，那么當(dāng)對(duì)這塊矩陣進(jìn)行讀

5、寫和計(jì)算時(shí)可能涉及頻繁的內(nèi)存與外存交換區(qū)的i/o，這回造成效率的嚴(yán)重下降。我一直不贊同在out of memory的情況下，通過增大虛存來解決問題，這樣，即使你勉強(qiáng)讓程序能繼續(xù)運(yùn)行下去，運(yùn)行時(shí)間也會(huì)變得極為緩慢這時(shí)應(yīng)該做的是對(duì)程序進(jìn)行重新思考和設(shè)計(jì)，降低其對(duì)內(nèi)存的耗用。第三，vectorization有些時(shí)候還可能增加實(shí)際運(yùn)算次數(shù)，這往往出現(xiàn)在不適合的向量化過程中。這樣，即使你通過生搬硬套的向量化過程讓操作變成矩陣運(yùn)算，但是增加的無用計(jì)算使得即使是更高效的引擎也無法挽回你的損失了。說了這么些，不是想勸說從向量化的道路退回去，而是提醒在做一些事情的時(shí)候，要考慮得周全一些。 (3) inplace

6、 operation。這是一個(gè)很有趣的事情，matlab的對(duì)象管理策略是copy-on-write，就是說你把一個(gè)矩陣傳給一個(gè)函數(shù)的時(shí)候，會(huì)先傳引用，而不產(chǎn)生副本，而當(dāng)函數(shù)要對(duì)這個(gè)矩陣修改的時(shí)候，才會(huì)制造出它的副本出來，讓函數(shù)去修改。這樣聽上去很完美，既保護(hù)了輸入矩陣不被改變，又避免了大量無謂的復(fù)制。不過，這個(gè)策略真的沒有缺陷么？可以看看下面的例子 a = rand(2000, 2000);for i = 1 : 1000 a = f(a);endfunction a = f(a)a(1) = a(1) + 1; % a temporary copy of a is produced for

7、modificationreturn; % the modified temporary copy is assigned to a 在上面的代碼中，只是想對(duì)a的第一個(gè)元素做1000次加法，結(jié)果卻導(dǎo)致了對(duì)整個(gè)有四百萬元素的大矩陣做了1000次副本復(fù)制！而且這些副本其實(shí)沒有用，只要f(a)直接對(duì)a的原矩陣進(jìn)行修改，這些巨大的浪費(fèi)就能避免。你可能跟我argue說，這里完全可以寫成a(1) = a(1) + 1000，不用這么搞。我想說明的是，我只是想舉一個(gè)簡單例子說明，copy-on-write的策略為什么可能造成巨大的效率浪費(fèi)。其實(shí)，很多小修改沒法像上面那么容易合并。 以前，為了解決這個(gè)問題，f

8、函數(shù)的作者需要自己寫一個(gè)mex function來避開matlab的保護(hù)機(jī)制，直接改寫a。后來matlab自己也意識(shí)到這種策略在實(shí)際中的問題，于是他們改寫了解釋器，采用了一種辦法：a = f(a, .) 當(dāng)它發(fā)現(xiàn)這樣的定義和調(diào)用的時(shí)候，它會(huì)很智能地了解到你其實(shí)是想改寫a這個(gè)輸入，于是它就把操作改成inplace，直接對(duì)a進(jìn)行，避免拷貝。這種修改的策略，在2006b后才比較正式的運(yùn)用，在舊版matlab中仍舊不同程度地存在上述問題。我在2007a中，測試了兩種f的寫法： function y = f0(x) y = x + 1; end function x = f1(x) x = x + 1;

9、 end確實(shí)發(fā)現(xiàn)大量調(diào)用f1比f0快了好多倍。所以，如果確實(shí)是想改變input的話，函數(shù)定義中要讓input parameter和output parameter同名，以觸發(fā)這種智能的解釋機(jī)制。并且建議升級(jí)到最新的2007a版本，這個(gè)版本在這方面進(jìn)行了重要改善讓matlab更快(2)matlab是一門在數(shù)值運(yùn)算方面非常高效的語言，對(duì)于提高matlab的效率，有一些大家熟知的方法，比如避免使用for循環(huán)，向量化，使用c-mex寫核心，等等。這些原則大體是沒錯(cuò)的，但是教條地運(yùn)用有時(shí)反而會(huì)降低效率。 根據(jù)問題的規(guī)模，適當(dāng)選擇向量化的策略 向量化的基本思想是以空間換時(shí)間，通過把要處理的對(duì)象轉(zhuǎn)化成一個(gè)矩

10、陣，集中處理。這比起非向量化的處理方式往往需要消耗更多的內(nèi)存。在內(nèi)存緊張的時(shí)候，內(nèi)存分配的開銷可能是很大的。 比如，給定一個(gè)列向量v1，和一個(gè)行向量v2, 計(jì)算矩陣m，使得m(i, j) = v1(i) + v2(j)。一種常用的向量化實(shí)現(xiàn)：m = length(v1); n = length(v2);m = repmat(v1, 1, n) + repmat(v2, m, 1);這個(gè)過程中，除了m以外，還要?jiǎng)?chuàng)建兩個(gè)大小為m x n的臨時(shí)矩陣。repmat創(chuàng)建臨時(shí)矩陣有一定的overhead，主要花費(fèi)在元素索引上。由于元素加法的向量化收益本身并不顯著，扣除repmat的額外開銷后，得益就更小了

11、。另外，當(dāng)這些矩陣很大，接近物理內(nèi)存容納能力的時(shí)候，必然導(dǎo)致和虛擬內(nèi)存進(jìn)行頻繁的交換，耗費(fèi)大量的io時(shí)間，向量化可能得不償失。這種情況下，下面這種向量化程度略低的方式可能是更好的選擇： m = zeros(m, n);if m 0);ns = ns(ns 0);sp = 1, cumsum(ns(1:end-1) + 1;result = vs(cumsum(sp);對(duì)于圖像處理或者二維信號(hào)問題，善用filter。 比如，要對(duì)每個(gè)像素，基于其局部鄰域進(jìn)行一些計(jì)算。即使整個(gè)計(jì)算本身不是線性的，但是只要能分解成線性局部運(yùn)算的組合，就可以利用filter。 比如，對(duì)于圖像i, 產(chǎn)生v，使得v(i,

12、j)是像素i(i, j)周圍w x w鄰域的像素的方差。i = im2double(i);h = ones(w, w) / (w * w);m = imfilter(i, h, symmetric);m2 = imfilter(i.2, h, symmetric);v = m2 - m.2;再談matlab效率(3)關(guān)于matlab的效率問題，很多文章，包括我之前寫的一些，主要集中在使用向量化以及相關(guān)的問題上。但是，最近我在實(shí)驗(yàn)時(shí)對(duì)代碼進(jìn)行profile的過程中，發(fā)現(xiàn)在新版本的matlab下，for-loop已經(jīng)得到了極大優(yōu)化，而效率的瓶頸更多是在函數(shù)調(diào)用和索引訪問的過程中。 由于matlab

13、特有的解釋過程，不同方式的函數(shù)調(diào)用和元素索引，其效率差別巨大。不恰當(dāng)?shù)氖褂梅绞娇赡茉诒緛聿黄鹧鄣牡胤綆韲?yán)重的開銷，甚至可能使你的代碼的運(yùn)行時(shí)間增加上千倍（這就不是多買幾臺(tái)服務(wù)器或者增加計(jì)算節(jié)點(diǎn)能解決的了，呵呵）。 下面通過一些簡單例子說明問題。（實(shí)驗(yàn)選在裝有windows vista的一臺(tái)普通的臺(tái)式機(jī)(core2 duo 2.33ghz + 4gb ram)進(jìn)行，相比于計(jì)算集群, 這可能和大部分朋友的環(huán)境更相似一些。實(shí)驗(yàn)過程是對(duì)某一個(gè)過程實(shí)施多次的整體進(jìn)行計(jì)時(shí)，然后得到每次過程的平均時(shí)間，以減少計(jì)時(shí)誤差帶來的影響。多次實(shí)驗(yàn)，在均值附近正負(fù)20%的范圍內(nèi)的置信度高于95%。為了避免算上首次運(yùn)行

14、時(shí)花在預(yù)編譯上的時(shí)間，在開始計(jì)時(shí)前都進(jìn)行充分的“熱身”運(yùn)行。） 函數(shù)調(diào)用的效率 一個(gè)非常簡單的例子，把向量中每個(gè)元素加1。（當(dāng)然這個(gè)例子根本不需要調(diào)函數(shù)，但是，用它主要是為了減少函數(shù)執(zhí)行本身的時(shí)間，突出函數(shù)解析和調(diào)用的過程。） 作為baseline，先看看最直接的實(shí)現(xiàn)% input u: u is a 1000000 x 1 vectorv = u + 1; 這個(gè)過程平均需要0.00105 sec。而使用長期被要求盡量避免的for-loop n = numel(u);% v = zeros(n, 1) has been pre-allocated.for i = 1 : n v(i) = u(

15、i) + 1;end所需的平均時(shí)間大概是0.00110 sec。從統(tǒng)計(jì)意義上說，和vectorization已經(jīng)沒有顯著差別。無論是for-loop或者vectorization，每秒平均進(jìn)行約10億次“索引-讀取-加法-寫入”的過程，計(jì)算資源應(yīng)該得到了比較充分的利用。 要是這個(gè)過程使用了函數(shù)調(diào)用呢？matlab里面支持很多種函數(shù)調(diào)用方式，主要的有m-function, function handle, anonymous function, inline, 和feval，而feval的主參數(shù)可以是字符串名字，function handle, anonymous function或者inlin

16、e。 用m-function，就是專門定義一個(gè)函數(shù)function y = fm(x) y = x + 1;在調(diào)用時(shí) for i = 1 : n v(i) = fm(u(i);endfunction handle就是用來把一個(gè)function賦到一個(gè)變量上，類似于c/c+的函數(shù)指針，或者c#里面的delegate的作用fh = fm;for i = 1 : n v(i) = fh(u(i);endanonymous function是一種便捷的語法來構(gòu)造簡單的函數(shù)，類似于lisp, python的lambda表達(dá)式fa = (x) x + 1;for i = 1 : n v(i) = fa(u

17、(i);endinline function是一種傳統(tǒng)的通過表達(dá)式字符串構(gòu)造函數(shù)的過程fi = inline(x + 1, x);for i = 1 : n v(i) = fi(u(i);endfeval的好處在于可以以字符串方式指定名字來調(diào)用函數(shù)，當(dāng)然它也可以接受別的參數(shù)。v(i) = feval(fm, u(i);v(i) = feval(fh, u(i);v(i) = feval(fa, u(i); 對(duì)于100萬次調(diào)用（包含for-loop本身的開銷，函數(shù)解析(resolution)，壓棧，執(zhí)行加法，退棧，把返回值賦給接收變量），不同的方式的時(shí)間差別很大： m-function 0.38

18、5 sec function handle 0.615 sec anonymous function 0.635 sec inline function 166.00 sec feval(fm, u(i) 8.328 sec feval(fh, u(i) 0.618 sec feval(fa, u(i) 0.652 sec feval(fm, u(i) 2.788 sec feval(fa, u(i) 34.689 sec 從這里面，我們可以看到幾個(gè)有意思的現(xiàn)象： 首先，調(diào)用自定義函數(shù)的開銷遠(yuǎn)遠(yuǎn)高于一個(gè)簡單的計(jì)算過程。直接寫 u(i) = v(i) + 1 只需要 0.0011 秒左右，而寫u

19、(i) = fm(v(i) 則需要0.385秒，時(shí)間多了幾百倍，而它們干的是同一件事情。這說明了，函數(shù)調(diào)用的開銷遠(yuǎn)遠(yuǎn)大于for-loop自己的開銷和簡單計(jì)算過程在不同情況可能有點(diǎn)差別，一般而言，一次普通的函數(shù)調(diào)用花費(fèi)的時(shí)間相當(dāng)于進(jìn)行了幾百次到一兩千次雙精度浮點(diǎn)加法。 使用function handle和anonymous function的開銷比使用普通的m-函數(shù)要高一些。這可能是由于涉及到句柄解析的時(shí)間，而普通的函數(shù)在第一次運(yùn)行前已經(jīng)在內(nèi)存中進(jìn)行預(yù)編譯。 inline function的運(yùn)行時(shí)間則令人難以接受了，竟然需要一百多秒（是普通函數(shù)調(diào)用的四百多倍，是直接計(jì)算的十幾萬倍）。這是因?yàn)閙a

20、tlab是在每次運(yùn)行時(shí)臨時(shí)對(duì)字符串表達(dá)式（或者它的某種不太高效的內(nèi)部表達(dá)）進(jìn)行parse。 feval(fh, u(i)和fh(u(i)，feval(fa, u(i)和fa(u(i)的運(yùn)行時(shí)間很接近，表明feval在接受句柄為主參數(shù)時(shí)本身開銷很小。但是，surprising的是 for i = 1 : n v(i) = feval(fm, u(i);end比起fh = fm;for i = 1 : n v(i) = feval(fh, u(i);end慢了4.5倍 (前者0.618秒，后者2.788秒)。for i = 1 : n v(i) = feval(x) x + 1, u(i);end

21、比起fa = (x) x + 1;for i = 1 : n v(i) = feval(fa, u(i);end竟然慢了53倍(前者0.652秒，后者34.689秒）。由于在matlab的內(nèi)部實(shí)現(xiàn)中，function handle的解析是在賦值過程中進(jìn)行的，所以預(yù)先用一個(gè)變量把句柄接下，其效果就是預(yù)先完成了句柄解析，而如果直接把fm或者(x) x + 1寫在參數(shù)列上，雖然寫法簡潔一些，但是解析過程是把參數(shù)被賦值到所調(diào)函數(shù)的局部變量時(shí)才進(jìn)行，每調(diào)用一次解析一次，造成了巨大的開銷。 feval使用字符串作為函數(shù)名字時(shí)，所耗時(shí)間比傳入句柄大，因?yàn)檫@涉及到對(duì)函數(shù)進(jìn)行搜索的時(shí)間（當(dāng)然這個(gè)搜索是在一個(gè)索引

22、好的cache里面進(jìn)行（除了第一次），而不是在所有path指定的路徑中搜索。）在2007年以后，matlab推出了arrayfun函數(shù)，上面的for-loop可以寫為 v = arrayfun(fh, u)這平均需要4.48 sec，這比起for-loop(需時(shí)0.615 sec)還慢了7倍多。這個(gè)看上去“消除了for-loop的函數(shù)，由于其內(nèi)部設(shè)計(jì)的原因，未必能帶來效率上的正效果。 元素和域的訪問 除了函數(shù)調(diào)用，數(shù)據(jù)的訪問方式對(duì)于效率也有很大影響。matlab主要支持下面一些形式的訪問： array-index a(i): cell-index: ci; struct field: s.fi

23、eldname struct field (dynamic): s.(fieldname)這里主要探索單個(gè)元素或者域的訪問（當(dāng)然，matlab也支持對(duì)于子數(shù)組的非常靈活整體索引）。 對(duì)于一百萬次訪問的平均時(shí)間 a(i) for a numeric array 0.0052 sec ci for a cell array 0.2568 sec struct field 0.0045 sec struct field (with dynamic name) 1.0394 sec 我們可以看到matlab對(duì)于單個(gè)數(shù)組元素或者靜態(tài)的struct field的訪問，可以達(dá)到不錯(cuò)的速度，在主流臺(tái)式機(jī)約每秒

24、2億次（連同for-loop的開銷）。而cell array的訪問則明顯緩慢，約每秒400萬次（慢了50倍）。matlab還支持靈活的使用字符串來指定要訪問域的語法（動(dòng)態(tài)名字），但是，是以巨大的開銷為代價(jià)的，比起靜態(tài)的訪問慢了200倍以上。 關(guān)于object-oriented programming matlab在新的版本中（尤其是2008版），對(duì)于面向?qū)ο蟮木幊烫峁┝藦?qiáng)大的支持。在2008a中，它對(duì)于oo的支持已經(jīng)不亞于python等的高級(jí)腳本語言。但是，我在實(shí)驗(yàn)中看到，雖然在語法上提供了全面的支持，但是matlab里面面向?qū)ο蟮男屎艿停_銷巨大。這里僅舉幾個(gè)例子。 object中的pro

25、perty的訪問速度是3500萬次，比struct field慢了6-8倍。matlab提供了一種叫做dependent property的屬性，非常靈活，但是，效率更低，平均每秒訪問速度竟然低至2.6萬次（這種速度基本甚至難以用于中小規(guī)模的應(yīng)用中）。 object中method調(diào)用的效率也明顯低于普通函數(shù)的調(diào)用，對(duì)于instance method，每百萬次調(diào)用，平均需時(shí)5.8秒，而對(duì)于static method，每百萬次調(diào)用需時(shí)25.8秒。這相當(dāng)于每次調(diào)用都需要臨時(shí)解析的速度，而matlab的類方法解析的效率目前還明顯偏低。 matlab中可以通過改寫subsref和subsasgn的方法，對(duì)于對(duì)象的索引和域的訪問進(jìn)行非常靈活的改造，可以通過它創(chuàng)建類似于數(shù)組的對(duì)象。但是，一個(gè)符合要求的subsref的行為比較復(fù)雜。在一個(gè)提供了subsref的對(duì)象中，大部分行為都需要subsref進(jìn)行調(diào)度，而默認(rèn)的較優(yōu)的調(diào)度方式將被掩蓋。在一個(gè)提供了subsref的類中（使用一種最快捷的實(shí)現(xiàn)），object property的平均訪問速度竟然降到1萬5千次每秒。建議 根據(jù)上面的分析，對(duì)于撰寫高效matlab代碼，我有下面一些建議： 雖然for-loop的速度有了很大改善，vectorization（向量化）仍舊是改善效率的重要途徑，尤其是在能把運(yùn)算