短短几年,由 MIT CSAIL 实验室开发的编程语言Julia已然成为编程界的新宠,尤其在科学计算领域炙手可热。很大部分是因为这门语言结合了 C 语言的速度、Ruby 的灵活、Python 的通用性,以及其他各种语言的优势于一身。那么你知道为什么Julia的速度能做到那么快吗?这并不是因为更好的编译器,而是一种更新的设计理念,Julia在开发之初就将这种理念纳入其中,而这也是关注“人生苦短”的Python所欠缺的。
\n为什么要选择Julia?因为它比其他脚本语言更快,它在具备Python、MATLAB、R语言开发速度的同时,又能生成与C语言和Fortran一样快的代码。
\n但Julia新手对这种说法可能会有点怀疑。
\n- \n
- \n
为什么其他脚本语言不也提升一下速度?Julia可以做到的,为什么其他脚本语言做不到?
\n\n - \n
你能提供基准测试来证明它的速度吗?
\n\n - \n
这似乎有违“天底下没有免费的午餐”的道理。它真的有那么完美吗?
\n\n
很多人认为Julia运行速度很快,因为它是即时编译(JIT)型的(也就是说,每条语句都使用编译的函数来运行,这些函数要么在使用之前进行即时编译,要么在之前已经编译过并放在缓存中)。这就引出了一个问题:Julia是否提供了比Python或R语言(MATLAB默认使用JIT)更好的JIT实现?因为人们在这些JIT编译器上所做的工作比Julia要多得多,所以我们凭什么认为Julia这么快就会超过这些编译器?但其实这完全是对Julia的误解。
\n我想以一种非常直观的方式说明,Julia的速度之所以快,是因为它的设计决策。Julia的的核心设计决策是通过多重分派实现专门化的类型稳定性,编译器因此可以很容易地生成高效的代码,同时还能够保持代码的简洁,让它“看起来就像一门脚本语言”。
\n但是,在本文的示例中,我们将看到Julia并不总是像其他脚本语言那样,我们必须接受“午餐不全是免费”的事实。
\n要看出它们之间的区别,我们只需要看看基本的数学运算。
\nJulia中的数学运算
\n一般来说,Julia中的数学运算与其他脚本语言中的数学运算看起来是一样的。它们的数字都是“真正的数字”,比如Float64就是64位浮点数或者类似于C语言中的“double”。Vector{Float64}与C语言double数组的内存布局是一样的,都可以很容易地与C语言进行互操作(实际上,在某种意义上,“Julia是构建在C语言之上的一个层”),从而带来更高的性能。
\n使用Julia进行一些数学运算:
\na = 2+2\nb = a/3\nc = a÷3 #\\div tab completion, means integer division\nd = 4*5\nprintln([a;b;c;d])\n\n[4.0, 1.33333, 1.0, 20.0]\n
\n
我在这里使用了Julia的unicode制表符补全功能。Julia允许使用unicode字符,这些字符可以通过制表符实现Latex风格的语句。同样,如果一个数字后面跟着一个变量,那么不需要使用*运算符就可以进行乘法运算。例如,下面的Julia的代码是合法的:
\nα = 0.5\n∇f(u) = α*u; ∇f(2)\nsin(2π)\n\n-2.4492935982947064e-16\n
\n
类型稳定性和代码内省
\n类型稳定性是指一个方法只能输出一种可能的类型。例如:*(::Float64,::Float64) 输出的类型是Float64。不管你给它提供什么参数,它都会返回一个Float64。这里使用了多重分派:“*”操作符根据它看到的类型调用不同的方法。例如,当它看到浮点数时,就会返回浮点数。Julia提供了代码自省宏,可以看到代码被编译成什么东西。因此,Julia不只是一门普通的脚本语言,还是一门可以让你处理汇编的脚本语言!和其他很多语言一样,Julia被编译成LLVM (LLVM是一种可移植的汇编格式)。
\n@code_llvm 2*5\n\n; Function *\n; Location: int.jl:54\ndefine i64 @\u0026quot;julia_*_33751\u0026quot;(i64, i64) {\ntop:\n %2 = mul i64 %1, %0\n ret i64 %2\n}\n
\n
这段代码的意思是:执行一个浮点数乘法操作,然后返回结果。我们也可以看一下汇编代码。
\n@code_native 2*5\n\n\t.text\n; Function * {\n; Location: int.jl:54\n\timulq\t%rsi, %rdi\n\tmovq\t%rdi, %rax\n\tretq\n\tnopl\t(%rax,%rax)\n;}\n
\n
“*”函数被编译成与C语言或Fortran中完全相同的操作,这意味着它可以达到相同的性能(尽管它是在Julia中定义的)。因此,Julia不仅可以“接近”C语言,而且实际上可以得到相同的C语言代码。那么在什么情况下会发生这种情况?
\nJulia的有趣之处在于,上面的这个问题其实问得不对,正确的问题应该是:在什么情况下代码不能被编译成像C语言或Fortran那样?这里的关键是类型稳定性。如果一个函数是类型稳定的,那么编译器就会知道函数在任意时刻的类型,就可以巧妙地将其优化为与C语言或Fortran相同的汇编代码。如果它不是类型稳定的,Julia必须进行昂贵的“装箱”,以确保在操作之前知道函数的类型是什么。
\n这是Julia与其他脚本语言之间最关键的不同点。
\n好的方面是Julia的函数(类型稳定)基本上就是C语言或Fortran的函数,因此“^”(乘方)运算速度很快。那么,类型稳定的^(::Int64,::Int64)会输出什么?
\n2^5\n\n32\n
\n
2^-5\n\n0.03125\n
\n
这里我们会得到一个错误。为了确保编译器可以为“^”返回一个Int64,它必须抛出一个错误。但在MATLAB、Python或R语言中这么做是不会抛出错误的,因为这些语言没有所谓的类型稳定性。
\n如果没有类型安全性会怎样?让我们看一下代码:
\n@code_native ^(2,5)\n\n\t.text\n; Function ^ {\n; Location: intfuncs.jl:220\n\tpushq\t%rax\n\tmovabsq\t$power_by_squaring, %rax\n\tcallq\t*%rax\n\tpopq\t%rcx\n\tretq\n\tnop\n;}\n
\n
现在,我们来定义自己的整数乘方运算。与其他脚本语言一样,我们让它变得更“安全”:
\nfunction expo(x,y)\n if y\u0026gt;0\n return x^y\n else\n x = convert(Float64,x)\n return x^y\n end\nend\n\nexpo (generic function with 1 method)\n
\n
现在运行一下看看行不行:
\nprintln(expo(2,5))\nexpo(2,-5)\n\n32\n0.03125\n
\n
再来看看汇编代码。
\n@code_native expo(2,5)\n\n\t.text\n; Function expo {\n; Location: In[8]:2\n\tpushq\t%rbx\n\tmovq\t%rdi, %rbx\n; Function \u0026gt;; {\n; Location: operators.jl:286\n; Function \u0026lt;; {\n; Location: int.jl:49\n\ttestq\t%rdx, %rdx\n;}}\n\tjle\tL36\n; Location: In[8]:3\n; Function ^; {\n; Location: intfuncs.jl:220\n\tmovabsq\t$power_by_squaring, %rax\n\tmovq\t%rsi, %rdi\n\tmovq\t%rdx, %rsi\n\tcallq\t*%rax\n;}\n\tmovq\t%rax, (%rbx)\n\tmovb\t$2, %dl\n\txorl\t%eax, %eax\n\tpopq\t%rbx\n\tretq\n; Location: In[8]:5\n; Function convert; {\n; Location: number.jl:7\n; Function Type; {\n; Location: float.jl:60\nL36:\n\tvcvtsi2sdq\t%rsi, %xmm0, %xmm0\n;}}\n; Location: In[8]:6\n; Function ^; {\n; Location: math.jl:780\n; Function Type; {\n; Location: float.jl:60\n\tvcvtsi2sdq\t%rdx, %xmm1, %xmm1\n\tmovabsq\t$__pow, %rax\n;}\n\tcallq\t*%rax\n;}\n\tvmovsd\t%xmm0, (%rbx)\n\tmovb\t$1, %dl\n\txorl\t%eax, %eax\n; Location: In[8]:3\n\tpopq\t%rbx\n\tretq\n\tnopw\t%cs:(%rax,%rax)\n;}\n
\n
这是一个非常直观的演示,说明了Julia通过使用类型推断获得了比其他脚本语言更高的性能。
\n核心思想:多重分派+类型稳定性=\u0026gt;速度+可读性
\n类型稳定性是Julia区别于其他脚本语言的一个关键特性。事实上,Julia的核心思想是这样的:
\n多重分派允许一种语言将函数调用分派给类型稳定的函数。
\n这就是Julia的核心思想,现在让我们花点时间深入了解一下。如果函数内部具有类型稳定性(也就是说,函数内的任意函数调用也是类型稳定的),那么编译器就会知道每一步的变量类型,它就可以在编译函数时进行充分的优化,这样得到的代码基本上与C语言或Fortran相同。多重分派在这里可以起到作用,它意味着“*”可以是一个类型稳定的函数:对于不同的输入,它有不同的含义。但是,如果编译器在调用“*”之前能够知道a和b的类型,那么它就知道应该使用哪个“*”方法,这样它就知道c=a*b的输出类型是什么。这样它就可以将类型信息一路传下去,从而实现全面的优化。
\n我们从中可以学到一些东西。首先,为了实现这种级别的优化,必须具有类型稳定性。大多数语言为了让用户可以更轻松地编码,都没有在标准库中提供这种特性。其次,需要通过多重分派来专门化类型函数,让脚本语言语法“看上去更显式”一些。最后,需要一个健壮的类型系统。为了构建非类型稳定的乘方运算,我们需要使用转换函数。因此,要在保持脚本语言的语法和易用性的同时实现这种原始性能必须将语言设计成具有多重分派类型稳定性的语言,并提供一个健壮的类型系统。
\nJulia基准测试
\nJulia官网提供的基准测试只是针对编程语言组件的执行速度,并没有说是在测试最快的实现,所以这里存在一个很大的误解。R语言程序员一边看着使用R语言实现的Fibonacci函数,一边说:“这是一段很糟糕的代码,不应该在R语言中使用递归,因为递归很慢”。但实际上,Fibonacci函数是用来测试递归的,而不是用来测试语言的执行速度的。
\nJulia使用了类型稳定函数的多重分派机制,因此,即使是早期版本的Julia也可以优化得像C语言或Fortran那样。非常明显,几乎在所有情况下,Julia都非常接近C语言。当然,也有与C语言不一样的地方,我们可以来看看这些细节。首先是在计算Fibonacci数列时C语言比Julia快2.11倍,这是因为这是针对递归的测试,而Julia并没有完全为递归进行过优化。Julia其实也可以加入这种优化(尾递归优化),只是出于某些原因他们才没有这么做,最主要是因为:可以使用尾递归的地方也可以使用循环,而循环是一种更加健壮的优化,所以他们建议使用循环来代替脆弱的尾递归。
\nJulia表现不太好的地方还有rand_mat_stat和parse_int测试。这主要是因为边界检查导致的。在大多数脚本语言中,如果你试图访问超出数组边界的元素就会出错,Julia默认情况下也会这么做。
\nfunction test1()\n a = zeros(3)\n for i=1:4\n a[i] = i\n end\nend\ntest1()\n\nBoundsError: attempt to access 3-element Array{Float64,1} at index [4]\n\nStacktrace:\n [1] setindex! at ./array.jl:769 [inlined]\n [2] test1() at ./In[11]:4\n [3] top-level scope at In[11]:7\n
\n
不过,你可以使用@inbounds宏来禁用这个功能:
\nfunction test2()\n a = zeros(3)\n @inbounds for i=1:4\n a[i] = i\n end\nend\ntest2()\n
\n
这样你就获得了与C语言或Fortran一样的不安全行为和执行速度。这是Julia的另一个有趣的特性:默认情况下是一个安全的脚本语言特性,在必要的时候禁用这个功能,以便获得性能提升。
\n严格类型
\n除了类型稳定性,你还需要严格类型。在Python中,你可以将任何东西放入数组中。而在Julia中,你只能将类型T放入Vector{T}中。Julia提供了各种非严格的类型,例如Any。如果有必要,可以创建Vector{Any},例如:
\na = Vector{Any}(undef,3)\na[1] = 1.0\na[2] = \u0026quot;hi!\u0026quot;\na[3] = :Symbolic\na\n\n3-element Array{Any,1}:\n 1.0 \n \u0026quot;hi!\u0026quot; \n :Symbolic\n
\n
Union是另一个不那么极端的抽象类型,例如:
\na = Vector{Union{Float64,Int}}(undef,3)\na[1] = 1.0\na[2] = 3\na[3] = 1/4\na\n\n3-element Array{Union{Float64, Int64},1}:\n 1.0 \n 3 \n 0.25\n
\n
这个Union只接受浮点数和整数。不过,它仍然是一个抽象类型。接受抽象类型作为参数的函数无法知道元素的类型(在这个例子中,元素要么是浮点数,要么是整数),这个时候,多重分派优化在这里起不到作用,所以Julia此时的性能就不如其他脚本语言。
\n所以我们可以得出一个性能原则:尽可能使用严格类型。使用严格类型还有其他好处:严格类型的Vector{Float64}实际上与C语言或Fortran是字节兼容的,所以不经过转换就可以直接用在C语言或Fortran程序中。
\n不免费的午餐
\n很明显,Julia为了在保持脚本语言特征的同时实现性能目标,做出了非常明智的设计决策。但是,它也为此付出了一些代价。接下来,我将展示Julia的一些奇特的东西及其相应的工具。
\n性能是可选的
\n之前已经说明了Julia提供了多种方法来提升性能(比如@inbounds),但我们不一定要使用它们。你也可以编写类型不稳定的函数,虽然与MATLAB、R语言、Python一样慢,但你绝对可以这么做。在对性能要求没有那么高的地方,可以将其作为一个可选项。
\n检查类型稳定性
\n由于类型稳定性非常重要,Julia为我们提供了一些工具,用来检查一个函数是不是类型稳定的,其中最重要的是@code_warntype宏。让我们用它来检查一个类型稳定的函数:
\n@code_warntype 2^5\n\nBody::Int64\n│220 1 ─ %1 = invoke Base.power_by_squaring(_2::Int64, _3::Int64)::Int64\n│ └── return %1\n
\n
请注意,它将函数中所有变量都显示为严格类型。那么expo会是怎样的?
\n@code_warntype expo(2,5)\n\nBody::Union{Float64, Int64}\n│╻╷ \u0026gt;2 1 ─ %1 = (Base.slt_int)(0, y)::Bool\n│ └── goto #3 if not %1\n│ 3 2 ─ %3 = π (x, Int64)\n│╻ ^ │ %4 = invoke Base.power_by_squaring(%3::Int64, _3::Int64)::Int64\n│ └── return %4\n│ 5 3 ─ %6 = π (x, Int64)\n││╻ Type │ %7 = (Base.sitofp)(Float64, %6)::Float64\n│ 6 │ %8 = π (%7, Float64)\n│╻ ^ │ %9 = (Base.sitofp)(Float64, y)::Float64\n││ │ %10 = $(Expr(:foreigncall, \u0026quot;llvm.pow.f64\u0026quot;, Float64, svec(Float64, Float64), :(:llvmcall), 2, :(%8), :(%9), :(%9), :(%8)))::Float64\n│ └── return %10\n
\n
请注意,可能的返回值是%4和%10,它们是不同的类型,因此返回类型被推断为Union{Float64,Int64}。为了准确地追踪这种不稳定性发生的位置,我们可以使用Traceur.jl:
\nusing Traceur\n@trace expo(2,5)\n\n┌ Warning: x is assigned as Int64\n└ @ In[8]:2\n┌ Warning: x is assigned as Float64\n└ @ In[8]:5\n┌ Warning: expo returns Union{Float64, Int64}\n└ @ In[8]:2\n\n32\n
\n
在第2行,x被分配了一个Int,而在第5行又被分配了一个Float64,因此它被推断为Union{Float64, Int64}。第5行是我们放置显式转换调用的地方,这样我们就确定了问题所在的位置。
\n处理必要的类型不稳定性
\n首先,我已经证明了某些在Julia会出错的函数在其他脚本语言中却可以“读懂你的想法”。在很多情况下,你会发现你可以从一开始就使用不同的类型,以此来实现类型稳定性(为什么不直接使用2.0^-5?)。但是,在某些情况下,你找不到合适的类型。这个问题可以通过转换来解决,但这样会失去类型稳定性。你必须重新考虑你的设计,并巧妙地使用多重分派。
\n假设我们有一个Vector{Union{Float64,Int}}类型的a,并且可能遇到必须使用a的情况,需要在a的每个元素上执行大量操作。在这种情况下,知道给定元素的类型将带来性能的大幅提升,但由于类型位于Vector{Union{Float64,Int}}中,因此无法在下面这样的函数中识别出类型:
\nfunction foo(array)\n for i in eachindex(array)\n val = array[i]\n # do algorithm X on val\n end\nend\n\nfoo (generic function with 1 method)\n
\n
不过,我们可以通过多重分派来解决这个问题。我们可以在元素上使用分派:
\nfunction inner_foo(val)\n # Do algorithm X on val\nend\n\ninner_foo (generic function with 1 method)\n
\n
然后将foo定义为:
\nfunction foo2(array::Array)\n for i in eachindex(array)\n inner_foo(array[i])\n end\nend\n\nfoo2 (generic function with 1 method)\n
\n
因为需要为分派检查类型,所以inner_foo函数是严格类型化的。因此,如果inner_foo是类型稳定的,那么就可以通过专门化inner_foo来提高性能。这就导致了一个通用的设计原则:在处理奇怪或非严格的类型时,可以使用一个外部函数来处理逻辑类型,同时使用一个内部函数来处理计算任务,实现最佳的性能,同时仍然具备脚本语言的通用能力。
\nREPL的全局作用域性能很糟糕
\nJulia全局作用域的性能很糟糕。官方的性能指南建议不要使用全局作用域。然而,新手可能会意识不到REPL其实就是全局作用域。为什么?首先,Julia是有嵌套作用域的。例如,如果函数内部有函数,那么内部函数就可以访问外部函数的所有变量。
\nfunction test(x)\n y = x+2\n function test2()\n y+3\n end\n test2()\nend\n\ntest (generic function with 1 method)\n
\n
在test2中,y是已知的,因为它是在test中定义的。如果y是类型稳定的,那么所有这些工作就可以带来性能的提升,因为test2可以假设y是一个整数。现在让我们来看一下在全局作用域里会发生什么:
\na = 3\nfunction badidea()\n a + 2\nend\na = 3.0\n\n3.0\n
\n
因为没有使用分派来专门化badidea,并且可以随时更改a的类型,因此badidea在编译时无法进行优化,因为在编译期间a的类型是未知的。但是,Julia允许我们声明常量:
\nconst a_cons = 3\nfunction badidea()\n a_cons + 2\nend\n\nbadidea (generic function with 1 method)\n
\n
请注意,函数将使用常量的值来进行专门化,因此它们在设置后应该保持不变。
\n在进行基准测试时会出现这种情况。新手会像下面这样对Julia进行基准测试:
\na = 3.0\n@time for i = 1:4\n global a\n a += i\nend\n\n0.000006 seconds (4 allocations: 64 bytes)\n
\n
但是,如果我们将它放在一个函数中,就可以实现优化。
\nfunction timetest()\n a = 3.0\n @time for i = 1:4\n a += i\n end\nend\ntimetest() # First time compiles\ntimetest()\n\n0.000001 seconds\n0.000000 seconds\n
\n
这个问题非常容易解决:不要在REPL的全局作用域内进行基准测试或计算执行时间。始终将代码放在函数中,或将它们声明为const。
\n结 论
\n速度是Julia的设计目标。类型稳定性和多重分派对Julia编译的专门化起到了关键的作用。而要达到如此精细的类型处理水平,以便尽可能有效地实现类型稳定性,并在不完全可能的情况下实现性能优化,需要一个健壮的类型系统。
\n英文原文:
\nhttps://ucidatascienceinitiative.github.io/IntroToJulia/Html/WhyJulia
\n