碾壓Python!為什麼Julia速度這麼快?

weixin_33806914發表於2019-03-09

短短几年,由 MIT CSAIL 實驗室開發的程式語言Julia已然成為程式設計界的新寵,尤其在科學計算領域炙手可熱。很大部分是因為這門語言結合了 C 語言的速度、Ruby 的靈活、Python 的通用性,以及其他各種語言的優勢於一身。那麼你知道為什麼Julia的速度能做到那麼快嗎?這並不是因為更好的編譯器,而是一種更新的設計理念,Julia在開發之初就將這種理念納入其中,而這也是關注“人生苦短”的Python所欠缺的。

為什麼要選擇Julia?因為它比其他指令碼語言更快,它在具備Python、MATLAB、R語言開發速度的同時,又能生成與C語言和Fortran一樣快的程式碼。

但Julia新手對這種說法可能會有點懷疑。

  1. 為什麼其他指令碼語言不也提升一下速度?Julia可以做到的,為什麼其他指令碼語言做不到?

  2. 你能提供基準測試來證明它的速度嗎?

  3. 這似乎有違“天底下沒有免費的午餐”的道理。它真的有那麼完美嗎?

很多人認為Julia執行速度很快,因為它是即時編譯(JIT)型的(也就是說,每條語句都使用編譯的函式來執行,這些函式要麼在使用之前進行即時編譯,要麼在之前已經編譯過並放在快取中)。這就引出了一個問題:Julia是否提供了比Python或R語言(MATLAB預設使用JIT)更好的JIT實現?因為人們在這些JIT編譯器上所做的工作比Julia要多得多,所以我們憑什麼認為Julia這麼快就會超過這些編譯器?但其實這完全是對Julia的誤解。

我想以一種非常直觀的方式說明,Julia的速度之所以快,是因為它的設計決策。Julia的的核心設計決策是通過多重分派實現專門化的型別穩定性,編譯器因此可以很容易地生成高效的程式碼,同時還能夠保持程式碼的簡潔,讓它“看起來就像一門指令碼語言”。

但是,在本文的示例中,我們將看到Julia並不總是像其他指令碼語言那樣,我們必須接受“午餐不全是免費”的事實。

要看出它們之間的區別,我們只需要看看基本的數學運算。

Julia中的數學運算

一般來說,Julia中的數學運算與其他指令碼語言中的數學運算看起來是一樣的。它們的數字都是“真正的數字”,比如Float64就是64位浮點數或者類似於C語言中的“double”。Vector{Float64}與C語言double陣列的記憶體佈局是一樣的,都可以很容易地與C語言進行互操作(實際上,在某種意義上,“Julia是構建在C語言之上的一個層”),從而帶來更高的效能。

使用Julia進行一些數學運算:

a = 2+2b = a/3c = a÷3 #\\div tab completion, means integer divisiond = 4*5println([a;b;c;d])[4.0, 1.33333, 1.0, 20.0]

我在這裡使用了Julia的unicode製表符補全功能。Julia允許使用unicode字元,這些字元可以通過製表符實現Latex風格的語句。同樣,如果一個數字後面跟著一個變數,那麼不需要使用*運算子就可以進行乘法運算。例如,下面的Julia的程式碼是合法的:

α = 0.5∇f(u) = α*u; ∇f(2)sin(2π)-2.4492935982947064e-16

型別穩定性和程式碼內省

型別穩定性是指一個方法只能輸出一種可能的型別。例如:*(::Float64,::Float64) 輸出的型別是Float64。不管你給它提供什麼引數,它都會返回一個Float64。這裡使用了多重分派:“*”操作符根據它看到的型別呼叫不同的方法。例如,當它看到浮點數時,就會返回浮點數。Julia提供了程式碼自省巨集,可以看到程式碼被編譯成什麼東西。因此,Julia不只是一門普通的指令碼語言,還是一門可以讓你處理彙編的指令碼語言!和其他很多語言一樣,Julia被編譯成LLVM (LLVM是一種可移植的彙編格式)。

@code_llvm 2*5; Function *; Location: int.jl:54define i64 @\u0026quot;julia_*_33751\u0026quot;(i64, i64) {top:  %2 = mul i64 %1, %0  ret i64 %2}

這段程式碼的意思是:執行一個浮點數乘法操作,然後返回結果。我們也可以看一下彙編程式碼。

@code_native 2*5\t.text; Function * {; Location: int.jl:54\timulq\t%rsi, %rdi\tmovq\t%rdi, %rax\tretq\tnopl\t(%rax,%rax);}

“*”函式被編譯成與C語言或Fortran中完全相同的操作,這意味著它可以達到相同的效能(儘管它是在Julia中定義的)。因此,Julia不僅可以“接近”C語言,而且實際上可以得到相同的C語言程式碼。那麼在什麼情況下會發生這種情況?

Julia的有趣之處在於,上面的這個問題其實問得不對,正確的問題應該是:在什麼情況下程式碼不能被編譯成像C語言或Fortran那樣?這裡的關鍵是型別穩定性。如果一個函式是型別穩定的,那麼編譯器就會知道函式在任意時刻的型別,就可以巧妙地將其優化為與C語言或Fortran相同的彙編程式碼。如果它不是型別穩定的,Julia必須進行昂貴的“裝箱”,以確保在操作之前知道函式的型別是什麼。

這是Julia與其他指令碼語言之間最關鍵的不同點。

好的方面是Julia的函式(型別穩定)基本上就是C語言或Fortran的函式,因此“^”(乘方)運算速度很快。那麼,型別穩定的^(::Int64,::Int64)會輸出什麼?

2^532
2^-50.03125

這裡我們會得到一個錯誤。為了確保編譯器可以為“^”返回一個Int64,它必須丟擲一個錯誤。但在MATLAB、Python或R語言中這麼做是不會丟擲錯誤的,因為這些語言沒有所謂的型別穩定性。

如果沒有型別安全性會怎樣?讓我們看一下程式碼:

@code_native ^(2,5)\t.text; Function ^ {; Location: intfuncs.jl:220\tpushq\t%rax\tmovabsq\t$power_by_squaring, %rax\tcallq\t*%rax\tpopq\t%rcx\tretq\tnop;}

現在,我們來定義自己的整數乘方運算。與其他指令碼語言一樣,我們讓它變得更“安全”:

function expo(x,y)    if y\u0026gt;0        return x^y    else        x = convert(Float64,x)        return x^y    endendexpo (generic function with 1 method)

現在執行一下看看行不行:

println(expo(2,5))expo(2,-5)320.03125

再來看看彙編程式碼。

@code_native expo(2,5)\t.text; Function expo {; Location: In[8]:2\tpushq\t%rbx\tmovq\t%rdi, %rbx; Function \u0026gt;; {; Location: operators.jl:286; Function \u0026lt;; {; Location: int.jl:49\ttestq\t%rdx, %rdx;}}\tjle\tL36; Location: In[8]:3; Function ^; {; Location: intfuncs.jl:220\tmovabsq\t$power_by_squaring, %rax\tmovq\t%rsi, %rdi\tmovq\t%rdx, %rsi\tcallq\t*%rax;}\tmovq\t%rax, (%rbx)\tmovb\t$2, %dl\txorl\t%eax, %eax\tpopq\t%rbx\tretq; Location: In[8]:5; Function convert; {; Location: number.jl:7; Function Type; {; Location: float.jl:60L36:\tvcvtsi2sdq\t%rsi, %xmm0, %xmm0;}}; Location: In[8]:6; Function ^; {; Location: math.jl:780; Function Type; {; Location: float.jl:60\tvcvtsi2sdq\t%rdx, %xmm1, %xmm1\tmovabsq\t$__pow, %rax;}\tcallq\t*%rax;}\tvmovsd\t%xmm0, (%rbx)\tmovb\t$1, %dl\txorl\t%eax, %eax; Location: In[8]:3\tpopq\t%rbx\tretq\tnopw\t%cs:(%rax,%rax);}

這是一個非常直觀的演示,說明了Julia通過使用型別推斷獲得了比其他指令碼語言更高的效能。

核心思想:多重分派+型別穩定性=\u0026gt;速度+可讀性

型別穩定性是Julia區別於其他指令碼語言的一個關鍵特性。事實上,Julia的核心思想是這樣的:

多重分派允許一種語言將函式呼叫分派給型別穩定的函式。

這就是Julia的核心思想,現在讓我們花點時間深入瞭解一下。如果函式內部具有型別穩定性(也就是說,函式內的任意函式呼叫也是型別穩定的),那麼編譯器就會知道每一步的變數型別,它就可以在編譯函式時進行充分的優化,這樣得到的程式碼基本上與C語言或Fortran相同。多重分派在這裡可以起到作用,它意味著“*”可以是一個型別穩定的函式:對於不同的輸入,它有不同的含義。但是,如果編譯器在呼叫“*”之前能夠知道a和b的型別,那麼它就知道應該使用哪個“*”方法,這樣它就知道c=a*b的輸出型別是什麼。這樣它就可以將型別資訊一路傳下去,從而實現全面的優化。

我們從中可以學到一些東西。首先,為了實現這種級別的優化,必須具有型別穩定性。大多數語言為了讓使用者可以更輕鬆地編碼,都沒有在標準庫中提供這種特性。其次,需要通過多重分派來專門化型別函式,讓指令碼語言語法“看上去更顯式”一些。最後,需要一個健壯的型別系統。為了構建非型別穩定的乘方運算,我們需要使用轉換函式。因此,要在保持指令碼語言的語法和易用性的同時實現這種原始效能必須將語言設計成具有多重分派型別穩定性的語言,並提供一個健壯的型別系統。

Julia基準測試

Julia官網提供的基準測試只是針對程式語言元件的執行速度,並沒有說是在測試最快的實現,所以這裡存在一個很大的誤解。R語言程式設計師一邊看著使用R語言實現的Fibonacci函式,一邊說:“這是一段很糟糕的程式碼,不應該在R語言中使用遞迴,因為遞迴很慢”。但實際上,Fibonacci函式是用來測試遞迴的,而不是用來測試語言的執行速度的。

Julia使用了型別穩定函式的多重分派機制,因此,即使是早期版本的Julia也可以優化得像C語言或Fortran那樣。非常明顯,幾乎在所有情況下,Julia都非常接近C語言。當然,也有與C語言不一樣的地方,我們可以來看看這些細節。首先是在計算Fibonacci數列時C語言比Julia快2.11倍,這是因為這是針對遞迴的測試,而Julia並沒有完全為遞迴進行過優化。Julia其實也可以加入這種優化(尾遞迴優化),只是出於某些原因他們才沒有這麼做,最主要是因為:可以使用尾遞迴的地方也可以使用迴圈,而迴圈是一種更加健壯的優化,所以他們建議使用迴圈來代替脆弱的尾遞迴。

Julia表現不太好的地方還有rand_mat_stat和parse_int測試。這主要是因為邊界檢查導致的。在大多數指令碼語言中,如果你試圖訪問超出陣列邊界的元素就會出錯,Julia預設情況下也會這麼做。

function test1()    a = zeros(3)    for i=1:4        a[i] = i    endendtest1()BoundsError: attempt to access 3-element Array{Float64,1} at index [4]Stacktrace: [1] setindex! at ./array.jl:769 [inlined] [2] test1() at ./In[11]:4 [3] top-level scope at In[11]:7

不過,你可以使用@inbounds巨集來禁用這個功能:

function test2()    a = zeros(3)    @inbounds for i=1:4        a[i] = i    endendtest2()

這樣你就獲得了與C語言或Fortran一樣的不安全行為和執行速度。這是Julia的另一個有趣的特性:預設情況下是一個安全的指令碼語言特性,在必要的時候禁用這個功能,以便獲得效能提升。

嚴格型別

除了型別穩定性,你還需要嚴格型別。在Python中,你可以將任何東西放入陣列中。而在Julia中,你只能將型別T放入Vector{T}中。Julia提供了各種非嚴格的型別,例如Any。如果有必要,可以建立Vector{Any},例如:

a = Vector{Any}(undef,3)a[1] = 1.0a[2] = \u0026quot;hi!\u0026quot;a[3] = :Symbolica3-element Array{Any,1}: 1.0         \u0026quot;hi!\u0026quot;      :Symbolic

Union是另一個不那麼極端的抽象型別,例如:

a = Vector{Union{Float64,Int}}(undef,3)a[1] = 1.0a[2] = 3a[3] = 1/4a3-element Array{Union{Float64, Int64},1}: 1.0  3    0.25

這個Union只接受浮點數和整數。不過,它仍然是一個抽象型別。接受抽象型別作為引數的函式無法知道元素的型別(在這個例子中,元素要麼是浮點數,要麼是整數),這個時候,多重分派優化在這裡起不到作用,所以Julia此時的效能就不如其他指令碼語言。

所以我們可以得出一個效能原則:儘可能使用嚴格型別。使用嚴格型別還有其他好處:嚴格型別的Vector{Float64}實際上與C語言或Fortran是位元組相容的,所以不經過轉換就可以直接用在C語言或Fortran程式中。

不免費的午餐

很明顯,Julia為了在保持指令碼語言特徵的同時實現效能目標,做出了非常明智的設計決策。但是,它也為此付出了一些代價。接下來,我將展示Julia的一些奇特的東西及其相應的工具。

效能是可選的

之前已經說明了Julia提供了多種方法來提升效能(比如@inbounds),但我們不一定要使用它們。你也可以編寫型別不穩定的函式,雖然與MATLAB、R語言、Python一樣慢,但你絕對可以這麼做。在對效能要求沒有那麼高的地方,可以將其作為一個可選項。

檢查型別穩定性

由於型別穩定性非常重要,Julia為我們提供了一些工具,用來檢查一個函式是不是型別穩定的,其中最重要的是@code_warntype巨集。讓我們用它來檢查一個型別穩定的函式:

@code_warntype 2^5Body::Int64│220 1 ─ %1 = invoke Base.power_by_squaring(_2::Int64, _3::Int64)::Int64│    └──      return %1

請注意,它將函式中所有變數都顯示為嚴格型別。那麼expo會是怎樣的?

@code_warntype expo(2,5)Body::Union{Float64, Int64}│╻╷ \u0026gt;2 1 ─ %1  = (Base.slt_int)(0, y)::Bool│    └──       goto #3 if not %1│  3 2 ─ %3  = π (x, Int64)│╻  ^  │   %4  = invoke Base.power_by_squaring(%3::Int64, _3::Int64)::Int64│    └──       return %4│  5 3 ─ %6  = π (x, Int64)││╻  Type  │   %7  = (Base.sitofp)(Float64, %6)::Float64│  6 │   %8  = π (%7, Float64)│╻  ^  │   %9  = (Base.sitofp)(Float64, y)::Float64││   │   %10 = $(Expr(:foreigncall, \u0026quot;llvm.pow.f64\u0026quot;, Float64, svec(Float64, Float64), :(:llvmcall), 2, :(%8), :(%9), :(%9), :(%8)))::Float64│    └──       return %10

請注意,可能的返回值是%4和%10,它們是不同的型別,因此返回型別被推斷為Union{Float64,Int64}。為了準確地追蹤這種不穩定性發生的位置,我們可以使用Traceur.jl:

using Traceur@trace expo(2,5)┌ Warning: x is assigned as Int64└ @ In[8]:2┌ Warning: x is assigned as Float64└ @ In[8]:5┌ Warning: expo returns Union{Float64, Int64}└ @ In[8]:232

在第2行,x被分配了一個Int,而在第5行又被分配了一個Float64,因此它被推斷為Union{Float64, Int64}。第5行是我們放置顯式轉換呼叫的地方,這樣我們就確定了問題所在的位置。

處理必要的型別不穩定性

首先,我已經證明了某些在Julia會出錯的函式在其他指令碼語言中卻可以“讀懂你的想法”。在很多情況下,你會發現你可以從一開始就使用不同的型別,以此來實現型別穩定性(為什麼不直接使用2.0^-5?)。但是,在某些情況下,你找不到合適的型別。這個問題可以通過轉換來解決,但這樣會失去型別穩定性。你必須重新考慮你的設計,並巧妙地使用多重分派。

假設我們有一個Vector{Union{Float64,Int}}型別的a,並且可能遇到必須使用a的情況,需要在a的每個元素上執行大量操作。在這種情況下,知道給定元素的型別將帶來效能的大幅提升,但由於型別位於Vector{Union{Float64,Int}}中,因此無法在下面這樣的函式中識別出型別:

function foo(array)  for i in eachindex(array)    val = array[i]    # do algorithm X on val  endendfoo (generic function with 1 method)

不過,我們可以通過多重分派來解決這個問題。我們可以在元素上使用分派:

function inner_foo(val)  # Do algorithm X on valendinner_foo (generic function with 1 method)

然後將foo定義為:

function foo2(array::Array)  for i in eachindex(array)    inner_foo(array[i])  endendfoo2 (generic function with 1 method)

因為需要為分派檢查型別,所以inner_foo函式是嚴格型別化的。因此,如果inner_foo是型別穩定的,那麼就可以通過專門化inner_foo來提高效能。這就導致了一個通用的設計原則:在處理奇怪或非嚴格的型別時,可以使用一個外部函式來處理邏輯型別,同時使用一個內部函式來處理計算任務,實現最佳的效能,同時仍然具備指令碼語言的通用能力。

REPL的全域性作用域效能很糟糕

Julia全域性作用域的效能很糟糕。官方的效能指南建議不要使用全域性作用域。然而,新手可能會意識不到REPL其實就是全域性作用域。為什麼?首先,Julia是有巢狀作用域的。例如,如果函式內部有函式,那麼內部函式就可以訪問外部函式的所有變數。

function test(x)    y = x+2    function test2()        y+3    end    test2()endtest (generic function with 1 method)

在test2中,y是已知的,因為它是在test中定義的。如果y是型別穩定的,那麼所有這些工作就可以帶來效能的提升,因為test2可以假設y是一個整數。現在讓我們來看一下在全域性作用域裡會發生什麼:

a = 3function badidea()    a + 2enda = 3.03.0

因為沒有使用分派來專門化badidea,並且可以隨時更改a的型別,因此badidea在編譯時無法進行優化,因為在編譯期間a的型別是未知的。但是,Julia允許我們宣告常量:

const a_cons = 3function badidea()    a_cons + 2endbadidea (generic function with 1 method)

請注意,函式將使用常量的值來進行專門化,因此它們在設定後應該保持不變。

在進行基準測試時會出現這種情況。新手會像下面這樣對Julia進行基準測試:

a = 3.0@time for i = 1:4    global a    a += iend0.000006 seconds (4 allocations: 64 bytes)

但是,如果我們將它放在一個函式中,就可以實現優化。

function timetest()    a = 3.0    @time for i = 1:4        a += i    endendtimetest() # First time compilestimetest()0.000001 seconds0.000000 seconds

這個問題非常容易解決:不要在REPL的全域性作用域內進行基準測試或計算執行時間。始終將程式碼放在函式中,或將它們宣告為const。

結 論

速度是Julia的設計目標。型別穩定性和多重分派對Julia編譯的專門化起到了關鍵的作用。而要達到如此精細的型別處理水平,以便儘可能有效地實現型別穩定性,並在不完全可能的情況下實現效能優化,需要一個健壯的型別系統。

英文原文:

https://ucidatascienceinitiative.github.io/IntroToJulia/Html/WhyJulia

\"image\"

相關文章