Julia 代码性能优化

避免全局变量

全局变量的值、类型，都可能变化。这使得编译器很难优化使用全局变量的代码。应尽量使用局部变量，或者把变量当做参数传递给函数。

对性能至关重要的代码，应放入函数中。

声明全局变量为常量可以显著提高性能：

const DEFAULT_VAL = 0

使用非常量的全局变量时，最好在使用时指明其类型，这样也能帮助编译器优化：

global x
y = f(x::Int + 1)

写作功能是很好的编程风格。它带来了更多可重用的代码，并阐明了正在做的步骤，以及它们的输入和输出是什么。

使用 @time 来衡量性能并且留心内存分配

衡量计算性能最有用的工具是 @time 宏. 下面的例子展示了良好的使用方式

julia> function f(n)
           s = 0
           for i = 1:n
               s += i/2
           end
           s
        end
f (generic function with 1 method)
 
julia> @time f(1)
elapsed time: 0.008217942 seconds (93784 bytes allocated)
0.5
 
julia> @time f(10^6)
elapsed time: 0.063418472 seconds (32002136 bytes allocated)
2.5000025e11

在第一次调用时 (@time f(1)), f 会被编译. (如果你在这次会话中还 没有使用过 @time, 计时函数也会被编译.) 这时的结果没有那么重要. 在 第二次调用时, 函数打印了执行所耗费的时间, 同时请注意, 在这次执行过程中 分配了一大块的内存. 相对于函数形式的 tic 和 toc, 这是 @time 宏的一大优势.

出乎意料的大块内存分配往往意味着程序的某个部分存在问题, 通常是关于类型 稳定性. 因此, 除了关注内存分配本身的问题, 很可能 Julia 为你的函数生成 的代码存在很大的性能问题. 这时候要认真对待这些问题并遵循下面的一些个建 议.

另外, 作为一个引子, 上面的问题可以优化为无内存分配 (除了向 REPL 返回结 果), 计算速度提升 30 倍

julia> @time f_improved(10^6)
elapsed time: 0.00253829 seconds (112 bytes allocated)
2.5000025e11

你可以从下面的章节学到如何识别 f 存在的问题并解决.

在有些情况下, 你的函数可能需要为本身的操作分配内存, 这样会使得问题变得 复杂. 在这种情况下, 可以考虑使用下面的工具之一来甄别问题, 或者将函数拆分, 一部分处理内存分 配, 另一部分处理算法 (参见 预分配内存).

工具 —

Julia 提供了一些工具包来鉴别性能问题所在 :

• Profiling 可以用来衡量代码的性能, 同时鉴别出瓶颈所在. 对于复杂的项目, 可以使用 ProfileView <https://github.com/timholy/ProfileView.jl> 扩展包来直观的展示分析 结果.
• 出乎意料的大块内存分配, – @time, @allocated, 或者 -profiler - 意味着你的代码可能存在问题. 如果你看不出内存分配的问题, -那么类型系统可能存在问题. 也可以使用 --track-allocation=user 来 -启动 Julia, 然后查看 *.mem 文件来找出内存分配是在哪里出现的.
• TypeCheck <https://github.com/astrieanna/TypeCheck.jl`_ 可以帮助找 出部分类型系统相关的问题. 另一个更费力但是更全面的工具是 code_typed. 特别留意类型为 Any 的变量, 或者 Union 类型. 这些问题可以使用下面的建议解决.
• Lint 扩展包可以指出程序一 些问题.

避免使用抽象类型参数的容器

当使用参数化类型，包括数组，最好避免可能的抽象类型参数化。

参照下面的代码：

a = Real[]    # typeof(a) = Array{Real,1}
if (f = rand()) < .8
    push!(a, f)
end

因为a是一个抽象类型Real的数组，所以它必须能够保存任何实际值。由于Real对象可以是任意大小和结构，所以必须将a表示为指向单独分配的Real对象的指针数组。因为f将永远是一个float64类型，我们应该这样做：

a = Float64[] # typeof(a) = Array{Float64,1}

这将创建一个64位浮点值的连续块，可以有效地操纵这些浮点值。

也可以参照这里：参数化类型.

类型声明

在 Julia 中，编译器能推断出所有的函数参数与局部变量的类型，因此声名变量类型不能提高性能。然而在有些具体实例中，声明类型还是非常有用的。

type Foo
    field
end

function foo(a::Array{Any,1})
    x = a[1]::Int32
    b = x+1
    ...
end

function with_keyword(x; name::Int = 1)
    ...
end

把函数拆开

把一个函数拆为多个，有助于编译器调用最匹配的代码，甚至将它内联。

举个应该把“复合函数”写成多个小定义的例子：

function norm(A)
    if isa(A, Vector)
        return sqrt(real(dot(A,A)))
    elseif isa(A, Matrix)
        return max(svd(A)[2])
    else
        error("norm: invalid argument")
    end
end

如下重写会更精确、高效：

norm(x::Vector) = sqrt(real(dot(x,x)))
norm(A::Matrix) = max(svd(A)[2])

写“类型稳定”的函数

尽量确保函数返回同样类型的数值。考虑下面定义：

pos(x) = x < 0 ? 0 : x

尽管看起来没问题，但是 0 是个整数（ Int 型）， x 可能是任意类型。因此，函数有返回两种类型的可能。这个是可以的，有时也很有用，但是最好如下重写：

pos(x) = x < 0 ? zero(x) : x

Julia 中还有 one 函数，以及更通用的 oftype(x,y) 函数，它将 y 转换为与 x 同样的类型，并返回。这仨函数的第一个参数，可以是一个值，也可以是一个类型。

避免改变变量类型

在一个函数中重复地使用变量，会导致类似于“类型稳定性”的问题：

function foo()
    x = 1
    for i = 1:10
        x = x/bar()
    end
    return x
end

局部变量 x 开始为整数，循环一次后变成了浮点数（ / 运算符的结果）。这使得编译器很难优化循环体。可以修改为如下的任何一种：

• 用 x = 1.0 初始化 x
• 声明 x 的类型： x::Float64 = 1
• 使用显式转换: x = one(T)

分离核心函数

很多函数都先做些初始化设置，然后开始很多次循环迭代去做核心计算。尽可能把这些核心计算放在单独的函数中。例如，下面的函数返回一个随机类型的数组：

function strange_twos(n)
    a = Array(randbool() ? Int64 : Float64, n)
    for i = 1:n
        a[i] = 2
    end
    return a
end

应该写成：

function fill_twos!(a)
    for i=1:length(a)
        a[i] = 2
    end
end
 
function strange_twos(n)
    a = Array(randbool() ? Int64 : Float64, n)
    fill_twos!(a)
    return a
end

Julia 的编译器依靠参数类型来优化代码。第一个实现中，编译器在循环时不知道 a 的类型（因为类型是随机的）。第二个实现中，内层循环使用 fill_twos! 对不同的类型 a 重新编译，因此运行速度更快。

第二种实现的代码更好，也更便于代码复用。

标准库中经常使用这种方法。如 abstractarray.jl 文件中的 hvcat_fill 和 fill! 函数。我们可以用这两个函数来替代这儿的 fill_twos! 函数。

形如 strange_twos 之类的函数经常用于处理未知类型的数据。比如，从文件载入的数据，可能包含整数、浮点数、字符串，或者其他类型。

沿列顺序访问数组

Julia中的多维数组以列优先顺序存储。这意味着数组一次堆叠一列。这可以用vec函数或者语法[:]验证，如下所示（注意数组的顺序是1 3 2 4，不是1 2 3 4）：

julia> x = [1 2; 3 4]
2x2 Array{Int64,2}:
 1  2
 3  4
 
julia> x[:]
4-element Array{Int64,1}:
 1
 3
 2
 4

这种排列数组的约定在许多语言中很常见，如FORTRAN、MATLAB和R（仅举几例）。在数组循环时，记住数组的排序会产生显著的性能影响。要记住的一条经验法则是，使用列主数组时，第一个索引变化最快。本质上，这意味着如果内部循环索引是第一个出现在一个切片表达式中，那么循环速度会更快。

来看下一个特殊的例子。假设我们想编写一个函数，接受一个Vector，并返回一个带有输入向量副本的行或列的Matrix（矩阵）。这里行或列是否填充了这些副本并不重要（也许代码的其余部分可以很容易地进行相应的调整）。我们至少有四个方法做到这一点（除了建议调用内置函数repmat）：

function copy_cols{T}(x::Vector{T})
    n = size(x, 1)
    out = Array(eltype(x), n, n)
    for i=1:n
        out[:, i] = x
    end
    out
end
 
function copy_rows{T}(x::Vector{T})
    n = size(x, 1)
    out = Array(eltype(x), n, n)
    for i=1:n
        out[i, :] = x
    end
    out
end
 
function copy_col_row{T}(x::Vector{T})
    n = size(x, 1)
    out = Array(T, n, n)
    for col=1:n, row=1:n
        out[row, col] = x[row]
    end
    out
end
 
function copy_row_col{T}(x::Vector{T})
    n = size(x, 1)
    out = Array(T, n, n)
    for row=1:n, col=1:n
        out[row, col] = x[col]
    end
    out
end

现在我们每个方法运行1万次来计时：

julia> x = randn(10000);
 
julia> fmt(f) = println(rpad(string(f)*": ", 14, ' '), @elapsed f(x))
 
julia> map(fmt, {copy_cols, copy_rows, copy_col_row, copy_row_col});
copy_cols:    0.331706323
copy_rows:    1.799009911
copy_col_row: 0.415630047
copy_row_col: 1.721531501

注意copy_cols比copy_rows快上很多。这是因为copy_cols尊重矩阵的基于列的内存布局，并一次填充一列。此外，copy_col_row比copy_row_col更快，因为它符合我们的规则，第一个元素出现在最内环路。

预分配内存

如果函数返回数组或其他复杂类型，则可能需要分配内存。不幸的是，经常分配和逆向收集垃圾会给系统带来巨大的负荷。

有时您可以通过预先分配输出来绕过在每个函数调用上分配内存的需要。这里有一个小例子，比较

function xinc(x)
    return [x, x+1, x+2]
end
 
function loopinc()
    y = 0
    for i = 1:10^7
        ret = xinc(i)
        y += ret[2]
    end
    y
end

和

function xinc!{T}(ret::AbstractVector{T}, x::T)
    ret[1] = x
    ret[2] = x+1
    ret[3] = x+2
    nothing
end
 
function loopinc_prealloc()
    ret = Array(Int, 3)
    y = 0
    for i = 1:10^7
        xinc!(ret, i)
        y += ret[2]
    end
    y
end

计时结果:

julia> @time loopinc()
elapsed time: 1.955026528 seconds (1279975584 bytes allocated)
50000015000000
 
julia> @time loopinc_prealloc()
elapsed time: 0.078639163 seconds (144 bytes allocated)
50000015000000

预分配还有其他优点，例如允许调用者从一个算法中控制“输出”类型。在上面的例子中，我们可以通过一个子数组而不是一个数组。

在的情况下，预分配可以使你的代码更丑陋，所以可能需要性能测量来判断是否使用。

避免I/O的字符串插值

当向一个文件（或其他I/O设备）写入数据时，形成额外的中间字符串会带来许多负载。比如，这里是错误的方法:

println(file, "$a $b")

应该这样写:

println(file, a, " ", b)

代码的第一个版本形成一个字符串，然后将其写入文件，而第二个版本直接将值写入文件。还要注意，在某些情况下，字符串插值可能比较难读。对比下下面的2个例子：

println(file, "$(f(a))$(f(b))")

与:

println(file, f(a), f(b))

处理有关舍弃的警告

被舍弃的函数，会查表并显示一次警告，而这会影响性能。建议按照警告的提示进行对应的修改。

小技巧

注意些有些小事项，能使内部循环更紧致。

• 避免不必要的数组。例如，不要使用 sum([x,y,z]) ，而应使用 x+y+z
• 对于较小的整数幂，使用 * 更好。如 x*x*x 比 x^3 好
• 针对复数 z ，使用 abs2(z) 代替 abs(z)^2 。一般情况下，对于复数参数，尽量用 abs2 代替 abs
• 对于整数除法，使用 div(x,y) 而不是 trunc(x/y), 使用 fld(x,y) 而不是 floor(x/y), 使用 cld(x,y) 而不是 ceil(x/y).

注释的性能

通过代码的优化可以使得程序变得更好。

• 使用@inbounds消除表达式中的数组边界检查。做这事之前一定要100%确定。如果下标总是出界，你可能会抓狂。
• 在for循环前面编写@simd，它适合于向量化。这个特性是实验性的，可以在Julia的未来版本中改变或消失。

这里有2个例子:

function inner( x, y )
    s = zero(eltype(x))
    for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end
 
function innersimd( x, y )
    s = zero(eltype(x))
    @simd for i=1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end
 
function timeit( n, reps )
    x = rand(Float32,n)
    y = rand(Float32,n)
    s = zero(Float64)
    time = @elapsed for j in 1:reps
        s+=inner(x,y)
    end
    println("GFlop        = ",2.0*n*reps/time*1E-9)
    time = @elapsed for j in 1:reps
        s+=innersimd(x,y)
    end
    println("GFlop (SIMD) = ",2.0*n*reps/time*1E-9)
end
 
timeit(1000,1000)

在一台2.4GHz Intel Core i5 处理器的电脑上, 耗时:

GFlop        = 1.9467069505224963
GFlop (SIMD) = 17.578554163920018

@simd for循环的的范围应该是一维范围。用于累加的变量，如示例中的s，称为缩减变量。通过使用@simd，您定义了循环的几个属性：

• 以任意或重叠的顺序执行迭代是安全的，并特别考虑减少变量。
• 减少变量的浮点运算可以重新排序，没有@simd可能导致不同的结果。
• 一个迭代不会等待另一次迭代以取得进展。

使用@simd只是给编译器许可进行矢量化。是否真的需要这么做取决于编译器。为了真的能带来优化效果，您的循环应该具有以下附加属性：

• 循环必须是嵌套最里面的循环。
• 循环体必须是直线代码。这就是为什么目前数组的访问都需要@inbounds。
• 访问必须有一个步幅模式，不能是“聚集”（随机索引读取）或“散射”（随机索引写入）。
• 步幅应是单位步幅。
• 在一些简单的情况，例如一个循环访问2-3个数组，LLVM自动矢量化可能自动启动，导致@simd没有进一步加速。