標籤：

為什麼在 Mathematica 中使用循環是低效的？

01-08

經常看到這句話，但是不明白個中緣由，在 Wolfram 官方文檔里好像也沒找到專門的分析。

因為在Mathematica中使用循環確實是低效的。。。。。。

深層次的原因涉及到Mathematica的底層實現所以我不太懂，但是至少從下面幾個例子可以看出Mathematica里確實有很多比循環更好的方法

  1. 求和

    首先舉一個最簡單的求和例子，求sum_{i=1}^{n}{sin{i}}的值。為了測試運行時間取n=10^6

    一個剛接觸Mathematica的同學多半會這樣寫

    sum = 0;
    For[i = 1, i &<= 10^6, i++, sum += Sin[N@i]]; (*其中N@i的作用是把整數i轉化為浮點數，類似於C里的double*) sum

    為了便於計時用Module封裝一下，運行時間是2.13秒，如下圖

    然後一個有一定Mathematica經驗的同學多半會知道同樣作為循環的Do速度比For快，於是他可能會這麼寫

    sum = 0;
    Do[sum += Sin[N@i], {i, 1, 10^6}];
    sum

    如下圖，用時1.37秒，比For快了不少

    當然了知道Do速度比For快的同學不太可能不知道Sum函數，所以上面其實是我口胡的，他應該會這麼寫

    Sum[Sin[N@i], {i, 1, 10^6}]

    如下圖，同樣的結果，只用了不到0.06秒

    如果這位同學還知道Listable屬性並且電腦內存不算太小的話，他也可能會這麼寫

    Tr@Sin[N@Range[10^6]]

    如下圖，只用了不到0.02秒，速度超過For循環的100倍

    當然了這只是一個最簡單的例子，而且如果數據量更大的話最後一種方法就不能用了。但是這也足以說明在求和時用循環是低效的，無論是內置的Sum函數還是向量化運算，在效率上都遠遠高於循環

    （這部分模仿了不同程序員如何編寫階乘函數這篇文章，強烈推薦對Mathematica有興趣的同學去看看）

  2. 迭代

    接下來舉一個迭代的例子，x_{n+1}=a x_{n}(1-x_{n})（即Logistic map），取a=3.5,x_{0}=0.5，為了測試運行時間同樣取n=10^6

    還是先用For循環的做法

    x = 0.5;
    For[i = 1, i &<= 10^6, i++, x = 3.5 x (1 - x); ]; x

    如下圖，運行時間2.06秒

    （Do循環和For類似，篇幅所限這裡就不寫了，有興趣的同學可以自行嘗試）

    然後看看內置的Nest函數

    Nest[3.5 # (1 - #) , 0.5, 10^6]

    如下圖，用時0.02秒，又是將近兩個數量級的效率差異

    當然了Nest的使用涉及到純函數，對於Mathematica初學者來說可能有些難以理解，而且一些比較複雜的迭代不太容易寫成Nest的形式，但是在迭代時Nest（還包括Fold）的效率確實要好於循環
  3. 遍歷列表

    4. 依然舉一個簡單的例子：求一個列表中偶數的個數。為測試生成10^6個1到10之間的隨機整數

    list = RandomInteger[{1, 10}, 10^6];
    (*生成10^6個隨機整數*)

    如果用For循環的話代碼是這樣的

    num = 0；
    For[i = 1, i &<= 10^6, i++, If[EvenQ@list[[i]], num++] ]; num

    如下圖，用時1.73秒

    保留上面的思路，單純的將For循環改為Scan （相當於沒有返回結果的Map），代碼如下

    num = 0；
    Scan[If[EvenQ@#, num++] , list];
    num

    如下圖，用時0.91 秒

    （Do循環用時1.00秒左右，篇幅所限就不傳圖了）

    摒棄循環的思路，用其他內置函數寫

    Count[list, _?EvenQ] // AbsoluteTiming
    (*直接用Count數出list中偶數的個數*)
    Count[EvenQ /@ list, True] // AbsoluteTiming
    (*用Map對list中的每個數判斷是否偶數，然後用Count數出結果中True的個數*)
    Select[list, EvenQ] // Length // AbsoluteTiming
    (*選取list中的所有偶數，然後求結果列表長度*)
    Count[EvenQ@list, True] // AbsoluteTiming
    (*利用EvenQ的Listable屬性直接判斷list的每個數是否偶數，然後數出結果中True的個數*)
    Sum[Boole[EvenQ@i], {i, list}]
    (*對list中的每個元素判斷是否偶數，將結果相加*)

    結果如下圖

    這個遍歷的例子舉得不算特別恰當，但也能說明一些問題了：Mathematica中內置了許多神奇的函數，其中大部分只要使用得當效率都比循環高（而且不是一點半點）。就算非要用循環，也要記得（任何能用Do代替For的時候）Do比For快，（遍歷列表時）Scan比Do快

  4. 用向量(矩陣)運算代替循環

    這個例子來自如何用 Python 科學計算中的矩陣替代循環？ - Kaiser 的回答，我只是把代碼從Python翻譯成了Mathematica而已。選這個例子是因為它有比較明確的物理意義，而且效率對比非常明顯

    代碼如下

    AbsoluteTiming[
    n = 100;
    u = unew = SparseArray[{{1, _} -&> 1}, {n, n}] // N // Normal;
    For[k = 1, k &<= 3000, k++, For[i = 2, i &< n, i++, For[j = 2, j &< n, j++, unew[[i, j]] = 0.25 (u[[i + 1, j]] + u[[i - 1, j]] + u[[i, j + 1]] + u[[i, j - 1]]) ] ]; u = unew; ]; u1 = u; ] (*用三重循環，迭代3000次*) ArrayPlot[u1, DataReversed -&> True, ColorFunction -&> "TemperatureMap"]
    (*用ArrayPlot繪圖*)

    AbsoluteTiming[
    n = 100;
    u = SparseArray[{{1, _} -&> 1}, {n, n}] // N // Normal;
    Do[
    u[[2 ;; -2, 2 ;; -2]] =
    0.25 (u[[3 ;; -1, 2 ;; -2]] + u[[1 ;; -3, 2 ;; -2]] +
    u[[2 ;; -2, 3 ;; -1]] + u[[2 ;; -2, 1 ;; -3]]),
    {k, 1, 3000}];
    u2 = u;
    ]
    (*用矩陣運算，迭代3000次*)
    ArrayPlot[u2, DataReversed -&> True, ColorFunction -&> "TemperatureMap"]
    (*用ArrayPlot繪圖*)

    運行結果For循環用時136秒，矩陣運算用時不足0.5秒，且兩者答案完全一樣。在演算法完全相同的情況下兩種寫法有著超過200倍的效率差距

    （圖片太長了這裡就不直接顯示了，鏈接放在下面）

    http://pic4.zhimg.com/65d66161f4c674b1149651c32f69935f_b.png

===========================我是結尾的分隔線===============================

這個答案其實從一開始就跑題了，還寫了這麼長的目的就在於希望讓大家切實地感受到循環的低效並安利一下Mathematica中其它高效的方法。正如wolray的答案中說的，既然選擇了使用Mathematica就應該多利用些MMA獨有的美妙函數，畢竟如果只是用循環的話C和Fortran之類的語言效率比MMA不知高到哪裡去了。。。。。。

當然我也不是讓大家就不用循環了，畢竟很多時候循環的直觀性和易讀性帶來的便利遠遠比那點效率重要。只是希望大家在循環之前能稍稍想一下，自己的目的是不是一定要用循環？可不可以用內置函數代替循環？就像上面的幾個例子，將循環換成內置函數程序的簡潔性和效率都大幅提高，長此以往相信你一定會愛上MMA的~

  • 題外話——關於用編譯提速循環

    在MMA中如果一定要使用循環又對效率有一定要求的話，可以選擇使用編譯，效率能有極大的提高。比如上面的第4個例子使用Complie編譯過後的Do循環用時只有1.86秒，速度提升了將近100倍。如果電腦中有C編譯器的話還可以在Compile中加入CompilationTarget -&> "C"選項，速度還能有所提升。編譯過後的代碼如下：

    In[10]:= cf = Compile[{{n, _Integer}, {times, _Integer}},
    Module[{u},
    u = ConstantArray[0., {n, n}];
    u[[1]] = ConstantArray[1., n];
    Do[
    Do[u[[i, j]] =
    0.25 (u[[i + 1, j]] + u[[i - 1, j]] + u[[i, j + 1]] +
    u[[i, j - 1]]),
    {i, 2, n - 1}, {j, 2, n - 1}
    ], {k, 1, times}];
    u
    ]
    ];
    u3 = cf[100, 3000]; // AbsoluteTiming
    ArrayPlot[u3, DataReversed -&> True, ColorFunction -&> "TemperatureMap"]

    Out[11]= {1.86055, Null}

    前3個例子也都可以通過編譯提速很多，這裡就不放代碼了，有興趣的同學可以自己動手試一試，如果遇到問題歡迎在評論中與我交流。

    需要注意的是編譯有很多注意事項，這裡推薦一篇寫的很好的教程，編譯中常見的問題裡面都有很好的講解：怎樣編譯（Compile）/編譯的通用規則/學會這6條，你也會編譯

    但是一般來講編譯很麻煩，而且再怎麼編譯效率也很難趕上直接用C，所以個人並不特別建議MMA初學者學習編譯。

大概是和Mathematica語言創造者的想法有關，以下只是我的猜測：

一般編程語言的流程式控制制包括：順序、條件分支、循環語句。以前還有goto語句，後來漸漸不用了，但是並不影響語法的表達能力。

而Mathematica語言的創造者可能認為，天下循環大抵如此，都逃不出那幾種範式，並做了歸納總結（即Sum、Count、Select、Nest等等），然後保留For、Do、While只是為了兼容。

你是否想過，在你寫for(i=0; i&<10; i++)這些話時，一個簡單的想法為何要寫得這麼長，是否可以抽象出一種模式？如果你這麼想了，其實就是那些內置函數，Mathematica認為它們是本原的，所以（做了底層優化使其變得）快。

後註：就我實際使用看，90%以上的情境，需要寫循環代碼的地方，邏輯都不會很複雜，編譯器會優化得很好的。

Nov.11更

關於這個問題目前有了許多新認識，頂樓的無影冬瓜列舉了許多MMA里循環對比內建函數的例子。不過最近偶然注意到，對於內建函數，無論是什麼高級語言，例如MATLAB、R，向量化的操作跑起來同樣會比單寫循環快上不少，所以這並不是MMA才有的特點。

大家可以看下這個問題 為什麼Mathematica里矩陣索引這麼慢，有沒有什麼辦法迴避或者改進？

在MMA里，Do循環和Table函數的寫法幾乎一模一樣，區別無非在於前者按過程式風格來更新變數，後者是直接生成矩陣（階數大了可能對內存更敏感），執行速度上Do循環只是稍慢。

相比起循環，更值得注意的是MMA里的矩陣元索引，感覺這才是它的真正短板，尤其數據量大了之後，無論是用循環直接索引，還是用其他函數繞過索引，跑起來的速度都會被MATLAB遠遠甩到後面。

經常能看見MMA里循環低效的說法，現在看來，我想主要還是針對那些習慣了C或者MALTAB這種凡事都要寫循環的語言的同學。對於複雜的數據操作，速度的瓶頸恐怕主要在於矩陣元索引上。

————————————

Mathematica所有操作矩陣相關的函數，例如Apply, Map, Table, Thread等，都經過了非常好的底層優化，換句話說，這些內建函數本身就是封裝好的循環，精心打扮之後向你投懷入抱，何不欣然領了這份情，安心享用之。

喜歡循環的話，用python足夠你進行大部分科學計算和數據處理了。你在mma里寫循環，一來需要額外的編譯，二來也是對庫函數的浪費。mma到底是一款費用高昂的商業軟體，底層的效率交給他們就好，作為主人，不玩點更刺激更上層的遊戲，怎麼對得起人家的鮮美胴體。

我在用mma之前用過一年多的matlab，自以為對矩陣操作瞭然於胸，世間萬物都是活色生香的matrix，只要抓住了行的循環和列的循環，就抓住了整個夏天（什麼鬼…

直到後來接觸mma，被她函數式向量化的操作哲學狠狠的扇到牆角，這才猛然醒悟到：

愛一個人，就應該懷著讚美，享用她最性感的地方。而不是抱怨她為什麼不像隔壁王老二的媳婦那樣，胸大無腦。

所有的循環，可以用Listable,map等列表操作代替，且速度比循環快很多倍。

mma的函數式編程精髓也在此。

原因並不清楚。但是根據Mathematica自己的教程的說法，其數據結構的核心是列表（list）。因此我想，其本意是期望大家使用list這種數據結構的。

很多人可能都是像我一樣，有C語言或者其它面向過程、面向對象的編程經驗，碰到問題首先想到循環。殊不知，在函數式編程之中，有極大的不同。循環，是能不用就不用的。

用一個例子來說明，既用來回答這個問題，也記錄一下我的學習過程。

題目：求1 + (1 + 2) + (1 + 2 + 3) +…… + (1 + 2 +…… + n)

以下定義n=5000；

如果用類似於C語言的模式來寫的話，代碼是這樣的：

sum = 0;

For[i = 1, i &<= n, i++,

For[j = 1, j &<= i, j++,

sum = sum + j

]

] // AbsoluteTiming

sum

簡單易懂，對吧？但是在MMA裡面慢到離譜。用時16.8秒。

用Do來替換，是這樣的，快了一些。用時7.7秒。

sum = 0;

Do[

sum = sum + j, {i, 1, n}, {j, 1, i}

] // AbsoluteTiming

sum

然而這都都不是函數式編程啊！函數式編程，是這樣的：

Total[Total /@ Table[j, {i, 1, n}, {j, 1, i}]] // AbsoluteTiming

用時2.99秒

Total[Flatten[Range[Range[n]]]] // AbsoluteTiming

用時3.03秒

Total[Total[Range[#]] /@ Range[n]] // AbsoluteTiming

用時0.18秒

以上幾個，可能不好理解的主要是這裡：Range[Range[n]]

注意：Range[ Range[ 3 ] ]

將產生這樣一個列表：{ {1}， {1，2}，{1，2，3}}

這樣你就能看懂了。

最後，超級大殺器：

Total[Accumulate[Range[n]]] // AbsoluteTiming

用時0秒

結論：儘可能地用內置函數，很快。

很長一段時間內，我的想法是和 @軒瀟 類似的……不，應該說還沒他/她想得那麼細，只是模模糊糊地覺得這大概是Mathematica的創造者們不鼓勵使用循環才造就了這個局面，直到前些日子吧里有人問了個差不多的問題才好好思考了一下，得出的結論是，傳統的循環語句在Mathematica里的自動編譯更難。更多內容可參考這裡：為什麼MMA用FOR循環的效率特別低？_mathematica吧

沒用過該軟體，僅從寫程序的角度來回答這個問題，若有偏頗或者錯誤請提醒或者摺疊。

程序都有空間與時間之間的矛盾，也就是說某種數據結構及演算法在面對不同類型問題時不可能總是最省內存又算得最快的，但在面對某一類問題時是最好的。

同一問題，該軟體有多種寫法，不同寫法在內部由不同演算法和數據結構實現。你用某種寫法，就是在暗示軟體可以優化的方向，比如並行化，規約等等。

for循環對於許多問題而言比較接近我們的第一感覺，佔用內存空間較少，在速度方面都比較慢。如果你會用for實現規約，循環展開，不一定比其他演算法慢。

剛剛開始學習mma，做題過程中，如果遇到在百萬級別的數據中挑選某些符合條件的數，先把數據Table出來再Select確實比用For循環之類快些，但是感覺速度差別並不明顯。

如果數據到了千萬級別甚至更多，直接把數據Table出來分分鐘炸內存死機重啟。但是循環結果幾乎不佔多少內存。

我也很想知道所謂數組操作的正確姿勢啊。

推薦閱讀：

※如何實現任意精度計算？
※學習 Mathematica 有什麼推薦書籍？
※Mathematica 直接計算二重極限（double limit）的格式是什麼？
※你對即將發布的 Mathematica 11 有什麼期待?
※如何寫出易讀的 Mathematica 代碼？

TAG:編程 | WolframMathematica |