標籤：

Matlab求一個矩陣中所有元素的平方和，兩種循環寫法為什麼性能相差很大？

01-05

機器學習中經常要求矩陣所有元素的平方和

測試代碼

function test()
A=rand(7777,7777);
fprintf("Loop1 :%f
",timeit(@()bench_loop1(A)));
fprintf("Loop2 :%f
",timeit(@()bench_loop2(A)));
fprintf("Direct:%f
",timeit(@()bench_sum(A)));
end

function s=bench_loop1(A)
v=0;
[~,N]=size(A);
for i=1:N
v=v+A(:,i).^2;
end
s=sum(v);
end

function s=bench_loop2(A)
s=0;
[~,N]=size(A);
for i=1:N
s=s+sum(A(:,i).^2);
end
end

function s=bench_sum(A)
s=sum(sum(A.^2));
end

結果

&>&> test
Loop1 :0.056253
Loop2 :0.138789
Direct:0.214043

分析及疑惑

看了標題一定會有人問我為什麼不向量化，所以我也用了一個最向量化（也是最直接最常見）的寫法來對比，事實表明在這裡向量化反而是最慢的，原因應該是A.^2產生了一個巨大的中間矩陣，在分配內存上花了很多時間。所以Matlab也並不是什麼時候都是向量化比循環好

但兩種循環寫法性能差一倍多就完全搞不懂了

顯然不是上面所說的中間變數占內存的問題，而且更快的第一種寫法實際上使用了兩倍於第二種的內存（雖然都很小），空間局部性也不如第二種好

調用sum()的overhead也不可能是原因，即使有overhead也不至於大到相當於幾千次加法的時間

如果是解釋器發現了循環可以並行化的話，不至於智障到第一種寫法能發現，第二種寫法就發現不了吧……

而且由於Matlab程序無法反彙編/反位元組碼，也不能通過這種途徑得知到底發生了什麼

順帶用Python的numpy寫了一模一樣的代碼（代碼略），用CPython運行的結果如下

$ python test.py
Loop1:
0.08826357725066399
Loop2:
0.10564035799665879
Direct:
0.2684035420699761

第一種寫法稍快，但差距遠沒有Matlab中的大，另外第二種寫法還顯著地比Matlab里快，什麼鬼？（兩邊的矩陣庫都是mkl）

試驗了一下多種寫法的速度，從慢到快排列。

A為7777階隨機方陣，時間為100次運行平均用時，所用Matlab版本為R2013a。

  • sum(sum(A.^2))：0.273 s
  • sum(A(:).^2)：0.272 s，跟前一種其實沒有區別

  • dot(A(:), A(:))：0.251 s，跟前一種也幾乎沒有區別
  • 題主的loop2：0.160 s
  • 題主的loop1：0.100 s
  • A(:)" * A(:)：0.031 s，明顯比其它寫法快

我猜第二種比第一種慢的原因大概是每次需要一個連續的 A(:,i) 大小的臨時空間用於存儲中間計算結果，之後再將這個結果輸入給 sum 函數；而第一種可能每個元素平方後加到 v 上的，不需要連續的那麼大的臨時空間。

你把第一種改成：

t = A(:,i).^2;
v = v + t;

速度會慢很多，而對第二種做類似修改幾乎沒有影響。

除此之外的差距可能就是 sum 的 overhead 所致。

總之看不到編譯結果只能靠猜了，所以太深究也沒可靠的支持。可能你這個版本知道的到下個版本就不一樣了，況且最近更新了 JIT，還有很多不完善的地方。

更新：

@立黨，本來想放在評論裡邊的，結果發現有點長了，就更新一下：

首先，你這個測試是放在 command window 中執行的吧，這樣會很慢的。建議你放在函數中執行：

function test
tic;
for n=1:1:7000
A=zeros(7000,1);
A=0;
end
toc;
tic;
A=zeros(7000,7000);
toc;
時間已過 0.022421 秒。
時間已過 0.000221 秒。

對比你的測試結果：

Elapsed time is 0.033199 seconds.
Elapsed time is 0.106985 seconds.

可以看到兩種都更快了，但是第二種加速非常明顯（我們電腦配置不同應該也不會產生500倍的差距吧，何況第一種並沒有快多少）以至於直接扭轉了對比結果。

其次，會不會每次都銷毀這點還不好說，JIT 稍微處理得好一些的話應該就不會在每次循環中都做銷毀操作（比方說內存足夠時完全可以在函數結束時再銷毀或者一定次數後銷毀），所以「每次都銷毀」這點應該也只是你的猜測，最好說明一下，否則可能會誤導別人。

另外，你的第一個例子是想簡化模型來對比，不過簡化的有點過了，第一個例子中 b=sum(A) 並沒有開闢新的空間複製A，只是一個求和（頂多在求和過程中會聲明一些臨時空間用於存儲多線程求和的結果，這些臨時空間要比A的小得多），而 A = zeros(7000,1) 只是聲明空間，沒有求和，所以這兩個運算對比體現不出來你想表達的意義

題主的例子中之所以我猜測需要複製 A 是因為他用了索引運算，索引運算和直接引用是非常不同的：

function test
A = rand(4000);
tic
a = sum(A);
toc
tic
b = sum(A(:,:));
toc
tic
C = A(:,:);
toc

時間已過 0.006529 秒。
時間已過 0.056680 秒。
時間已過 0.045250 秒。

可以看到第二種方法用了索引（即使是索引了整個 A ）之後時間增加非常多；第三組是一個對照，可以看到單索引 A 需要的時間已經超出了 a=sum(A) 需要的全部時間，也就是說 sum(A) 的時候不太可能是先去複製 A

通常情況下只有索引結果需要用臨時變數存儲（但似乎題主的第一個例子中處理得更好一些），而如果直接使用整個變數則通常只是引用，即使傳入函數也一樣，除非你在函數內部試圖修改輸入變數MATLAB才會複製一份輸入：

function test
a = rand(5000);
tic, b = fun1(a); toc
tic, c = fun2(a); toc
function b = fun1(a)
b = a + 1;
function b = fun2(a)
a(1) = 0;
b = a + 1;
時間已過 0.038851 秒。
時間已過 0.122280 秒。

不過有種情況例外，就是如果你的輸入和輸出再函數內部以及調用的時候都是一樣的，這樣反而會更快一些。這裡就不多說了

原因前面的人都說了，bench_loop1不需要用一個臨時變數存儲A(:,i).^2的結果，而bench_loop2需要。

@Falccm的答案中指出，A(:,:)這樣的索引會複製A中所有數據。

而A(:)是不複製的，即數據是共享的。這一點我在另一個問題中提到過。

此外，MATLAB優化最好的運算（之一）——矩陣乘法——是值得信賴的。所以，

s = A(:)." * A(:);

bench_loop1比bench_loop2快一倍的原因，是因為在N次循環中，每次sum都需要開闢N-by-1內存空間，運算完s=s+sum(XXX)之後再銷毀。開闢和銷毀N-by-1的內存空間也需要O(N)的時間。測試如下：

tic;
A=zeros(7000,1);
toc;

tic;
b=sum(A);
toc;

Elapsed time is 0.005106 seconds.
Elapsed time is 0.004475 seconds.

用 Julia 跑的結果：

function test()
A=rand(7777,7777);

@time bench_loop1(A)
@time bench_loop2(A)
@time bench_sum(A)
@time bench_vec(A)
end

function bench_loop1(A)
v=0
~,N=size(A)
for i in 1:N
v=v+A[:,i].^2
end
s=sum(v)
return s
end

function bench_loop2(A)
s=0
~,N=size(A)
for i in 1:N
s=s+sum(A[:,i].^2)
end
return s
end

function bench_sum(A)
s=sum(sum(A.^2))
return s
end

function bench_vec(A)
return A[:]"*A[:]
end

test()

在 Win10 x64, julia 0.4.6, i5-4200H 2.8GHz, 12 GB 上是

2.136738 seconds (92.30 k allocations: 1.354 GB, 17.70% gc time)

1.626556 seconds (76.75 k allocations: 924.642 MB, 5.60% gc time)

1.435639 seconds (4 allocations: 461.439 MB, 0.09% gc time)

0.836528 seconds (6 allocations: 922.878 MB, 2.27% gc time

在 OS X 10.11.1, julia 0.4.6, i5-4258U 2.4GHz, 16 GB 上是

2.418516 seconds (92.30 k allocations: 1.354 GB, 14.08% gc time)

1.976998 seconds (76.75 k allocations: 924.642 MB, 2.65% gc time)

2.168147 seconds (4 allocations: 461.439 MB, 0.02% gc time)

0.688714 seconds (6 allocations: 922.878 MB, 0.89% gc time)

同一台 Windows 機器上跑 MATLAB R2015a 的結果（最後是 A(:)"*A(:)）：

Loop1: 0.151757

Loop2: 0.235973

Direct: 0.392878

Vectorized: 0.048278

最後一個快得很顯著。

我同意立黨的說法，覺得本質上還是在循環體里進行函數調用會涉及到參數內存空間的複製與銷毀，而且由於是動態類型語言，在調用的時候還會涉及類型檢查之類的額外開銷，所以對於性能苛刻的應用還是盡量寫簡單的函數體。

至於經cpython優化的python，基本成了強類型語言了，其循環是在c層次實現的，不會有類型檢查之類的開銷，至於調用過程是參數空間複製還是指針層面實現我拿不太准，如果是指針層面實現的那速度快也是應該的。

你問兩種循環方法差距為什麼那麼大我不知道

但是我知道你的向量化方法有問題:

A=rand(7777);

tic

v=0;

[~,N]=size(A);

for i=1:N

v=v+A(:,i).^2;

end

s=sum(v);

toc

tic

b=sum(sum(A.^2));

toc

tic

c=A(:)"*A(:);

toc

Elapsed time is 0.036093 seconds.

Elapsed time is 0.132542 seconds.

Elapsed time is 0.021972 seconds.

明顯還是向量化快嘛

有人說我MATLAB這部分測的不對，那我先刪了吧。以後有空仔細測一下再補上。

@Colliot 提到Julia。在向量化方面Julia和MATLAB剛好相反，由於強大的即時編譯器，for循環是非常高效的，向量化反而會因為動態內存分配浪費時間。因此Julia裡面最高效的實現應該是

function sum(A::Matrix{Float64})
v = 0.0
for i in eachindex(A)
v += A[i]^2
end
return v
end

寫在一個函數里，不要用全局REPL模式。這個代碼的動態內存開銷是0。Julia裡面0 allocation才是最快的。

測試如下：

function mysum(A::Matrix{Float64})
v = 0.0
for i in eachindex(A)
v += A[i]^2
end
v
end

mysum2(A::Matrix{Float64}) = A[:]"*A[:]

function test()
A = rand(7777,7777)
v = 0.0
@time mysum(A)
@time mysum2(A)
end

test()

0.331807 seconds
0.413919 seconds (6 allocations: 922.878 MB, 11.04% gc time)

順帶分享三個Julia技巧，官方Documentation上都可以找到。

  1. 函數最後如果要返回值的話，return關鍵字是可以不寫的，當然為了清楚寫上也挺好的。
  2. 單行的函數可以寫成類似數學函數定義f(x)=a^2的單行形式，不需要寫到function關鍵字裡面。
  3. eachindex是最高效的遍歷迭代器，不要再用什麼length(A), size(A)這種土辦法啦。

歡迎交流

mma：

In[1]:= A = RandomReal[1, {7777, 7777}];
Total[A^2, 2] // AbsoluteTiming
#.# @Flatten@A // AbsoluteTiming

Out[2]= {0.298461, 2.01564*10^7}
Out[3]= {0.223832, 2.01564*10^7}

第二種比第一種快的原因應該是涉及到matlab的矩陣存儲方式問題導致的定址問題，第一種方式每一行都是在內存中連續一段內存，而第二種方式不同每次一個取的向量A(:,i)中的每個數值在原始內存中不是連續的，需要每次定址；

上Octave結果：

Elapsed time is 0.162482 seconds.

Elapsed time is 0.215184 seconds.

Elapsed time is 0.153517 seconds.

神奇！

補充numpy：（執行10次）

7777 1.225519s

7777 1.598701s

7777 2.113372s

推薦閱讀：

※C 語言指針的指針和二維數組的區別?
※如何生成rest api文檔？
※零基礎學習python需要直接使用linux嗎？
※寫代碼的時候適合聽的音樂有哪些推薦？
※函數式編程有必要學設計模式和演算法嗎？

TAG:Python | 數學 | 編程 | MATLAB | 機器學習 |