加速你的MATLAB開發(3): 使用Profiler找出性能瓶頸

要說我們實習生myc最喜歡又最討厭的事，就是老闆給他定的「提交的程序一定不能比之前慢」這個標準了。

作為一個學習Computer Vision的學生，myc經常和矩陣打交道，之前他照著這篇paper寫了一個簡單的Triangulation的演算法，摘錄部分如下：

function alpha = slowernM1 = rand(3);nM2 = rand(3);npoints1 = rand(3,10000);npoints2 = rand(3,10000);nnfor i = 1:size(points1,2);n n point1 = points1(:,i);n point2 = points2(:,i);n n a1 = inv(M1)*point1;n a2 = inv(M2)*point2;n n A = [a1,-a2];n y = rand(3,1);n n AtA = A*A;n if cond(AtA) > 1/eps(double)n alpha(:,i) = 0;n elsen alpha(:,i) = inv(AtA) * A * y;n end n nendnnend%End of slowern

老闆說，

Why so slow?[1]

作為合格的實習生，myc知道調用MATLAB的神器，Profiler，的時候到了

>>profile on;slow;profile reportn

得到一份運行時間的報告。

首先我們關注的是當中一塊Code Analyzer results， 這一部分是最容易改進並且最有可能提高速度的環節。

這份報告指出，在解這種問題的時候

inv(A)*bn

這種寫法，有可能要比

Abn

更慢，且更不精確。畢竟Ab是MATLAB的招牌嗎

Code Analyzer results還指出，alpha這個變數在循環中的大小是不斷變化的。作為優化速度的常用手段，我們需要在循環前給alpha預留大小。

於是僅僅根據Code Analyzer的建議，myc就將程序改成了如下的樣子：

function alpha = betternM1 = rand(3);nM2 = rand(3);npoints1 = rand(3,10000);npoints2 = rand(3,10000);n% Pre-allocatenalpha = zeros(2,size(points1,2));nnfor i = 1:size(points1,2);n n point1 = points1(:,i);n point2 = points2(:,i);n n a1 = M1 point1;n a2 = M2 point2;n n A = [a1,-a2];n y = rand(3,1);n n AtA = A*A;n if cond(AtA) > 1/eps(double)n alpha(:,i) = 0;n elsen alpha(:,i) = AtA A * y;n end n nendnnend%End of bettern

再次運行Profile

>>profile on;better;profile reportn

用時稍微快了一點(2.63s vs 3.26s)，僅僅是參考給出的建議就可以輕鬆提速20%有沒有

Code Analyzer的建議都沒啦，那我們的優化到此為止了嗎？答案顯然是沒有，接下來我們開始分析報告最上方的表格里靠前的代碼。

這些代碼通常是你程序的瓶頸所在。

myc看到20，23，13，14，和16行耗時好像都挺多的，他決定逐條來進行分析

20行，if判斷，執行了10000次，耗時略久，但是好像暫時沒有什麼好的想法

23行，A * A A * y, 好像和Ay是等價的，這個地方可以優化

13行和14行，都是一個3x3的矩陣左除(mldivide)一個3x1的向量，查看mldivide的文檔,我發現可以將Ab變為AB,其中B為3xN的矩陣，其效果等同於循環做Ab。這裡也可以優化

16行，將兩個向量拼接，似乎也可以優化，但是先放著吧。

改好的程序如下

function alpha = betterThanBetternM1 = rand(3);nM2 = rand(3);npoints1 = rand(3,10000);npoints2 = rand(3,10000);nalpha = zeros(2,size(points1,2));nna1 = M1 points1;na2 = M2 points2;nnfor i = 1:size(points1,2);n n A = [a1(:,i),-a2(:,i)];n y = rand(3,1);n n if cond(A) > 1/eps(double)n alpha(:,i) = 0;n elsen alpha(:,i) = A y;n endn nendnnend%End of betterThanBettern

同樣，運行Profile

>>profile on;betterThanBetter;profile offn

又快了不少有沒有（1.62s vs 2.63s)。顯然我們將以前最大的瓶頸，10000次Ab變成了AB，且只進行一次計算，大幅提高了程序的速度。

這時候再看第一個表格中前幾行的用時，myc是這麼想的

19行憑我的能力估計不能再快了，13，14行都是內置的函數，已經也不能再快了

16行好像仍然是程序的瓶頸，1/eps(double)應該不是很慢，估計就是cond這個函數太慢了，讓我來看看cond的文檔....

看我發現了什麼？有個叫rcond的函數，也可以檢驗矩陣的條件數，但是他沒有cond精確，然而速度會更快一點。我的應用里不需要很精確的條件數，速度比精度更重要。

...

然而rcond好像只能用方形矩陣，看來我為了用rcond,還是需要進行額外的一步運算，不知道這樣會不會比以前更慢？試試

function alpha = bestnM1 = rand(3);nM2 = rand(3);npoints1 = rand(3,10000);npoints2 = rand(3,10000);nalpha = zeros(2,size(points1,2));nna1 = M1 points1;na2 = M2 points2;nnfor i = 1:size(points1,2);n n A = [a1(:,i),-a2(:,i)];n y = rand(3,1);n AtA = A * A;n if rcond(AtA) < eps(double)n alpha(:,i) = 0;n elsen alpha(:,i) = A y;n endn nendnnend%End of bestn

>>profile on;best;profile reportn

結果出來了，真的比之前的還要快！(0.99s vs 1.6s)

看來rcond帶來的速度提升完全抵消了A*A這一步矩陣運算。myc對這個優化結果非常滿意。

myc通過這次練習得出了幾點經驗，用Profiler優化你的程序一般分以下幾步

1. 關注Code Analyzer result給出的提示，基本都是無腦提高程序效率的方式
2. 修復非常明顯可以提升速度的代碼，用最少的修改時間換取最大程度的提速
3. 仔細閱讀文檔，使用更適合的替代函數。
4. 有時候為了提速反而需要增加更多的代碼。
5. 當你覺得無法再優化代碼的時候，就可以停下了，與其花更多的時間鑽研如何將自己的一個程序提高5%的速度，不如把這些時間用在做10%額外的工作上

其實第5點是myc老闆說的……

但是myc覺得，有時為了優化代碼而犧牲了可讀性，也是不行的。如果只是1%左右的速度提升，而別人需要多花1分鐘去讀懂你的代碼，是不是也得不償失呢？

於是他決定將代碼上傳到這裡，讓大家來幫他優化，他自己當然是去旅遊去了……

你能打敗我們實習生的best代碼嗎？

==============更新==============

@Falccm提供了一個更快的方法，完美地解決了循環開始時A和AtA的瓶頸.

function alpha = best nM1 = rand(3); nM2 = rand(3); npoints1 = rand(3,10000); npoints2 = rand(3,10000); ny = rand(size(points1)); nalpha = zeros(2,size(points1,2)); nA = [ reshape(M1 points1, size(M1,1), 1, []) ...n -reshape(M2 points2, size(M2,1), 1, [])]; nEps = eps(double); nAtA = [sum(bsxfun(@times, A, A(:,1,:))); sum(bsxfun(@times, A, A(:,2,:)))]; nfor i = 1:size(points1,2) n if rcond(AtA(:,:,i)) < Eps n alpha(:,i) = 0; n else n alpha(:,i) = A(:,:,i) y(:,i); n end nend n

[1]: 真事，這話是張正友在某年的CVPR對俄羅斯實習生說的

推薦閱讀：