什么要选择 Julia?因为它比其他脚本语言更快,它在具备 Python、MATLAB、R 语言开发速度的同时,又能生成与 C 语言和 Fortran 一样快的代码。
但 Julia 新手对这种说法可能会有点怀疑。
为什么其他脚本语言不也提升一下速度?Julia 可以做到的,为什么其他脚本语言做不到?
你能提供基准测试来证明它的速度吗?
这似乎有违“天底下没有免费的午餐”的道理。它真的有那么完美吗?
很多人认为 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+2
b = a/3
c = a÷3 #\div tab completion, means integer division
d = 4*5
println([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:54
define i64 @"julia_*_33751"(i64, i64) {
top:
%2 = mul i64 %1, %0
ret i64 %2
}
这段代码的意思是:执行一个浮点数乘法操作,然后返回结果。我们也可以看一下汇编代码。
@code_native 2*5
.text
; Function * {
; Location: int.jl:54
imulq %rsi, %rdi
movq %rdi, %rax
retq
nopl (%rax,%rax)
;}
“*”函数被编译成与 C 语言或 Fortran 中完全相同的操作,这意味着它可以达到相同的性能(尽管它是在 Julia 中定义的)。因此,Julia 不仅可以“接近”C 语言,而且实际上可以得到相同的 C 语言代码。那么在什么情况下会发生这种情况?
Julia 的有趣之处在于,上面的这个问题其实问得不对,正确的问题应该是:在什么情况下代码不能被编译成像 C 语言或 Fortran 那样?这里的关键是类型稳定性。如果一个函数是类型稳定的,那么编译器就会知道函数在任意时刻的类型,就可以巧妙地将其优化为与 C 语言或 Fortran 相同的汇编代码。如果它不是类型稳定的,Julia 必须进行昂贵的“装箱”,以确保在操作之前知道函数的类型是什么。
这是 Julia 与其他脚本语言之间最关键的不同点。
好的方面是 Julia 的函数(类型稳定)基本上就是 C 语言或 Fortran 的函数,因此“^”(乘方)运算速度很快。那么,类型稳定的 ^(::Int64,::Int64) 会输出什么?
2^5
32
2^-5
0.03125
这里我们会得到一个错误。为了确保编译器可以为“^”返回一个 Int64,它必须抛出一个错误。但在 MATLAB、Python 或 R 语言中这么做是不会抛出错误的,因为这些语言没有所谓的类型稳定性。
如果没有类型安全性会怎样?让我们看一下代码:
@code_native ^(2,5)
.text
; Function ^ {
; Location: intfuncs.jl:220
pushq %rax
movabsq $power_by_squaring, %rax
callq *%rax
popq %rcx
retq
nop
;}
现在,我们来定义自己的整数乘方运算。与其他脚本语言一样,我们让它变得更“安全”:
function expo(x,y)
if y>0
return x^y
else
x = convert(Float64,x)
return x^y
end
end
expo (generic function with 1 method)
现在运行一下看看行不行:
println(expo(2,5))
expo(2,-5)
32
0.03125
再来看看汇编代码。
@code_native expo(2,5)
.text
; Function expo {
; Location: In[8]:2
pushq %rbx
movq %rdi, %rbx
; Function >; {
; Location: operators.jl:286
; Function <; {
; Location: int.jl:49
testq %rdx, %rdx
;}}
jle L36
; Location: In[8]:3
; Function ^; {
; Location: intfuncs.jl:220
movabsq $power_by_squaring, %rax
movq %rsi, %rdi
movq %rdx, %rsi
callq *%rax
;}
movq %rax, (%rbx)
movb $2, %dl
xorl %eax, %eax
popq %rbx
retq
; Location: In[8]:5
; Function convert; {
; Location: number.jl:7
; Function Type; {
; Location: float.jl:60
L36:
vcvtsi2sdq %rsi, %xmm0, %xmm0
;}}
; Location: In[8]:6
; Function ^; {
; Location: math.jl:780
; Function Type; {
; Location: float.jl:60
vcvtsi2sdq %rdx, %xmm1, %xmm1
movabsq $__pow, %rax
;}
callq *%rax
;}
vmovsd %xmm0, (%rbx)
movb $1, %dl
xorl %eax, %eax
; Location: In[8]:3
popq %rbx
retq
nopw %cs:(%rax,%rax)
;}
这是一个非常直观的演示,说明了 Julia 通过使用类型推断获得了比其他脚本语言更高的性能。
核心思想:多重分派 + 类型稳定性 =>速度 + 可读性
类型稳定性是 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
end
end
test1()
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
end
end
test2()
这样你就获得了与 C 语言或 Fortran 一样的不安全行为和执行速度。这是 Julia 的另一个有趣的特性:默认情况下是一个安全的脚本语言特性,在必要的时候禁用这个功能,以便获得性能提升。
严格类型
除了类型稳定性,你还需要严格类型。在 Python 中,你可以将任何东西放入数组中。而在 Julia 中,你只能将类型 T 放入 Vector{T}中。Julia 提供了各种非严格的类型,例如 Any。如果有必要,可以创建 Vector{Any},例如:
a = Vector{Any}(undef,3)
a[1] = 1.0
a[2] = "hi!"
a[3] = :Symbolic
a
3-element Array{Any,1}:
1.0
"hi!"
:Symbolic
Union 是另一个不那么极端的抽象类型,例如:
a = Vector{Union{Float64,Int}}(undef,3)
a[1] = 1.0
a[2] = 3
a[3] = 1/4
a
3-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^5
Body::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}
│╻╷ >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, "llvm.pow.f64", 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]:2
32
在第 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
end
end
foo (generic function with 1 method)
不过,我们可以通过多重分派来解决这个问题。我们可以在元素上使用分派:
function inner_foo(val)
# Do algorithm X on val
end
inner_foo (generic function with 1 method)
然后将 foo 定义为:
function foo2(array::Array)
for i in eachindex(array)
inner_foo(array[i])
end
end
foo2 (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()
end
test (generic function with 1 method)
在 test2 中,y 是已知的,因为它是在 test 中定义的。如果 y 是类型稳定的,那么所有这些工作就可以带来性能的提升,因为 test2 可以假设 y 是一个整数。现在让我们来看一下在全局作用域里会发生什么:
a = 3
function badidea()
a + 2
end
a = 3.0
3.0
因为没有使用分派来专门化 badidea,并且可以随时更改 a 的类型,因此 badidea 在编译时无法进行优化,因为在编译期间 a 的类型是未知的。但是,Julia 允许我们声明常量:
const a_cons = 3
function badidea()
a_cons + 2
end
badidea (generic function with 1 method)
请注意,函数将使用常量的值来进行专门化,因此它们在设置后应该保持不变。
在进行基准测试时会出现这种情况。新手会像下面这样对 Julia 进行基准测试:
a = 3.0
@time for i = 1:4
global a
a += i
end
0.000006 seconds (4 allocations: 64 bytes)
但是,如果我们将它放在一个函数中,就可以实现优化。
function timetest()
a = 3.0
@time for i = 1:4
a += i
end
end
timetest() # First time compiles
timetest()
0.000001 seconds
0.000000 seconds
这个问题非常容易解决:不要在 REPL 的全局作用域内进行基准测试或计算执行时间。始终将代码放在函数中,或将它们声明为 const。
结 论
速度是 Julia 的设计目标。类型稳定性和多重分派对 Julia 编译的专门化起到了关键的作用。而要达到如此精细的类型处理水平,以便尽可能有效地实现类型稳定性,并在不完全可能的情况下实现性能优化,需要一个健壮的类型系统。