2010-11-08  16:23:50 zz:  http://dahua.spaces.live.com/blog/cns!28AF4251DF30CA42!2459.entry June 02

MATLAB 效率再議

關于MATLAB的效率問題，很多文章，包括我之前寫的一些，主要集中在使用向量化以及相關的問題上。但是，最近我在實驗時對代碼進行profile的過程中，發現在新版本的MATLAB下，for-loop已經得到了極大優化，而效率的瓶頸更多是在函數調用和索引通路的過程中。

由于MATLAB特有的解釋過程，不同方式的函數調用和元素索引，其效率差别巨大。不恰當的使用方式可能在本來不起眼的地方帶來嚴重的開銷，甚至可能使你的代碼的運作時間增加上千倍（這就不是多買幾台伺服器或者增加計算節點能解決的了，呵呵）。

下面通過一些簡單例子說明問題。（實驗選在裝有Windows Vista的一台普通的桌上型電腦(Core2 Duo 2.33GHz + 4GB Ram)進行，相比于計算叢集, 這可能和大部分朋友的環境更相似一些。實驗過程是對某一個過程實施多次的整體進行計時，然後得到每次過程的平均時間，以減少計時誤差帶來的影響。多次實驗，在均值附近正負20%的範圍内的置信度高于95%。為了避免算上首次運作時花在預編譯上的時間，在開始計時前都進行充分的“熱身”運作。）

函數調用的效率

一個非常簡單的例子，把向量中每個元素加1。（當然這個例子根本不需要調函數，但是，用它主要是為了減少函數執行本身的時間，突出函數解析和調用的過程。）

作為baseline，先看看最直接的實作

% input u: u is a 1000000 x 1 vector
v = u + 1; 
      

這個過程平均需要0.00105 sec。

而使用長期被要求盡量避免的for-loop

n = numel(u);
% v = zeros(n, 1) has been pre-allocated.
for i = 1 : n
    v(i) = u(i) + 1;
end
      

所需的平均時間大概是0.00110 sec。從統計意義上說，和vectorization已經沒有顯著差别。無論是for-loop或者vectorization，每秒平均進行約10億次“索引-讀取-加法-寫入”的過程，計算資源應該得到了比較充分的利用。

要是這個過程使用了函數調用呢？MATLAB裡面支援很多種函數調用方式，主要的有m-function, function handle, anonymous function, inline, 和feval，而feval的主參數可以是字元串名字，function handle, anonymous function或者inline。

用m-function，就是專門定義一個函數

function y = fm(x)
    y = x + 1;
      

在調用時

for i = 1 : n
    v(i) = fm(u(i));
end
      

function handle就是用@來把一個function賦到一個變量上，類似于C/C++的函數指針，或者C#裡面的delegate的作用

fh = @fm;
for i = 1 : n
    v(i) = fh(u(i));
end
      

anonymous function是一種便捷的文法來構造簡單的函數，類似于LISP, Python的lambda表達式

fa = @(x) x + 1;
for i = 1 : n
    v(i) = fa(u(i));
end
      

inline function是一種傳統的通過表達式字元串構造函數的過程

fi = inline('x + 1', 'x');
for i = 1 : n
    v(i) = fi(u(i));
end
      

feval的好處在于可以以字元串方式指定名字來調用函數，當然它也可以接受别的參數。

v(i) = feval('fm', u(i));
v(i) = feval(fh, u(i));
v(i) = feval(fa, u(i)); 
      

對于100萬次調用（包含for-loop本身的開銷，函數解析(resolution)，壓棧，執行加法，退棧，把傳回值賦給接收變量），不同的方式的時間差别很大：

m-function	0.385 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

從這裡面，我們可以看到幾個有意思的現象：

• 首先，調用自定義函數的開銷遠遠高于一個簡單的計算過程。直接寫 u(i) = v(i) + 1 隻需要  0.0011 秒左右，而寫u(i) = fm(v(i)) 則需要0.385秒，時間多了幾百倍，而它們幹的是同一件事情。這說明了，函數調用的開銷遠遠大于for-loop自己的開銷和簡單計算過程——在不同情況可能有點差别，一般而言，一次普通的函數調用花費的時間相當于進行了幾百次到一兩千次雙精度浮點加法。
• 使用function handle和anonymous function的開銷比使用普通的m-函數要高一些。這可能是由于涉及到句柄解析的時間，而普通的函數在第一次運作前已經在記憶體中進行預編譯。
• inline function的運作時間則令人難以接受了，竟然需要一百多秒（是普通函數調用的四百多倍，是直接計算的十幾萬倍）。這是因為matlab是在每次運作時臨時對字元串表達式（或者它的某種不太高效的内部表達）進行parse。
• feval(fh, u(i))和fh(u(i))，feval(fa, u(i))和fa(u(i))的運作時間很接近，表明feval在接受句柄為主參數時本身開銷很小。但是，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
  比起
  fa = @(x) x + 1;
  for i = 1 : n
      v(i) = feval(fa, u(i));
  end
  竟然慢了53倍(前者0.652秒，後者34.689秒）。
        

  由于在MATLAB的内部實作中，function handle的解析是在指派過程中進行的，是以預先用一個變量把句柄接下，其效果就是預先完成了句柄解析，而如果直接把@fm或者@(x) x + 1寫在參數列上，雖然寫法簡潔一些，但是解析過程是把參數被指派到所調函數的局部變量時才進行，每調用一次解析一次，造成了巨大的開銷。
• feval使用字元串作為函數名字時，所耗時間比傳入句柄大，因為這涉及到對函數進行搜尋的時間（當然這個搜尋是在一個索引好的cache裡面進行（除了第一次），而不是在所有path指定的路徑中搜尋。）

在2007年以後，MATLAB推出了arrayfun函數，上面的for-loop可以寫為

v = arrayfun(fh, u)
      

這平均需要4.48 sec，這比起for-loop(需時0.615 sec)還慢了7倍多。這個看上去“消除了for-loop"的函數，由于其内部設計的原因，未必能帶來效率上的正效果。

元素和域的通路

除了函數調用，資料的通路方式對于效率也有很大影響。MATLAB主要支援下面一些形式的通路：

• array-index  A(i): 
• cell-index:  C{i}; 
• struct field:  S.fieldname
• struct field (dynamic):  S.('fieldname')

這裡主要探索單個元素或者域的通路（當然，MATLAB也支援對于子數組的非常靈活整體索引）。

對于一百萬次通路的平均時間

A(i) for a numeric array	0.0052 sec
C{i} for a cell array	0.2568 sec
struct field	0.0045 sec
struct field (with dynamic name)	1.0394 sec

我們可以看到MATLAB對于單個數組元素或者靜态的struct field的通路，可以達到不錯的速度，在主流桌上型電腦約每秒2億次（連同for-loop的開銷）。而cell array的通路則明顯緩慢，約每秒400萬次（慢了50倍）。MATLAB還支援靈活的使用字元串來指定要通路域的文法（動态名字），但是，是以巨大的開銷為代價的，比起靜态的通路慢了200倍以上。

關于Object-oriented Programming

MATLAB在新的版本中（尤其是2008版），對于面向對象的程式設計提供了強大的支援。在2008a中，它對于OO的支援已經不亞于python等的進階腳本語言。但是，我在實驗中看到，雖然在文法上提供了全面的支援，但是matlab裡面面向對象的效率很低，開銷巨大。這裡僅舉幾個例子。

• object中的property的通路速度是3500萬次，比struct field慢了6-8倍。MATLAB提供了一種叫做dependent property的屬性，非常靈活，但是，效率更低，平均每秒通路速度竟然低至2.6萬次（這種速度基本甚至難以用于中小規模的應用中）。
• object中method調用的效率也明顯低于普通函數的調用，對于instance method，每百萬次調用，平均需時5.8秒，而對于static method，每百萬次調用需時25.8秒。這相當于每次調用都需要臨時解析的速度，而matlab的類方法解析的效率目前還明顯偏低。
• MATLAB中可以通過改寫subsref和subsasgn的方法，對于對象的索引和域的通路進行非常靈活的改造，可以通過它建立類似于數組的對象。但是，一個符合要求的subsref的行為比較複雜。在一個提供了subsref的對象中，大部分行為都需要subsref進行排程，而預設的較優的排程方式将被掩蓋。在一個提供了subsref的類中（使用一種最快捷的實作），object property的平均通路速度竟然降到1萬5千次每秒。

建議

根據上面的分析，對于撰寫高效MATLAB代碼，我有下面一些建議：

1. 雖然for-loop的速度有了很大改善，vectorization（向量化）仍舊是改善效率的重要途徑，尤其是在能把運算改寫成矩陣乘法的情況下，改善尤為顯著。
2. 在不少情況下，for-loop本身已經不構成太大問題，尤其是當循環體本身需要較多的計算的時候。這個時候，改善機率的關鍵在于改善循環體本身而不是去掉for-loop。
3. MATLAB的函數調用過程（非built-in function）有顯著開銷，是以，在效率要求較高的代碼中，應該盡可能采用扁平的調用結構，也就是在保持代碼清晰和可維護的情況下，盡量直接寫表達式和利用built-in function，避免不必要的自定義函數調用過程。在次數很多的循環體内（包括在cellfun, arrayfun等實際上蘊含循環的函數）形成長調用鍊，會帶來很大的開銷。
4. 在調用函數時，首選built-in function，然後是普通的m-file函數，然後才是function handle或者anonymous function。在使用function handle或者anonymous function作為參數傳遞時，如果該函數被調用多次，最好先用一個變量接住，再傳入該變量。這樣，可以有效避免重複的解析過程。
5. 在可能的情況下，使用numeric array或者struct array，它們的效率大幅度高于cell array（幾十倍甚至更多）。對于struct，盡可能使用普通的域（字段，field）通路方式，在非效率關鍵，執行次數較少，而靈活性要求較高的代碼中，可以考慮使用動态名稱的域通路。
6. 雖然object-oriented從軟體工程的角度更為優勝，而且object的使用很多時候很友善，但是MATLAB目前對于OO的實作效率很低，在效率關鍵的代碼中應該慎用objects。
7. 如果需要設計類，應該盡可能采用普通的property，而避免靈活但是效率很低的dependent property。如非确實必要，避免重載subsref和subsasgn函數，因為這會全面接管對于object的接口調用，往往會帶來非常巨大的開銷（成千上萬倍的減慢），甚至使得本來幾乎不是問題的代碼成為性能瓶頸。