性能技巧

在以下部分中,我们将简要介绍一些可以帮助您的 Julia 代码尽可能快地运行的技术。

性能关键代码应在函数内

任何性能关键的代码都应在函数内。由于 Julia 编译器的运作方式,函数内的代码通常比顶级代码运行得快得多。

使用函数不仅对性能很重要:函数更易于重用和测试,并阐明了正在执行的步骤以及它们的输入和输出,编写函数,而不仅仅是脚本 也是 Julia 风格指南的建议。

函数应接受参数,而不是直接对全局变量进行操作,请参阅下一条。

避免无类型全局变量

无类型全局变量的值可能随时更改,这可能导致其类型的更改。这使得编译器难以优化使用全局变量的代码。这也适用于类型值变量,即全局级别的类型别名。只要可能,变量应是局部的,或作为参数传递给函数。

我们发现全局名称通常是常量,将它们声明为常量可以大大提高性能

const DEFAULT_VAL = 0

如果已知全局变量始终具有相同的类型,应注释其类型

可以通过在使用点注释其类型来优化无类型全局变量的使用

global x = rand(1000)

function loop_over_global()
    s = 0.0
    for i in x::Vector{Float64}
        s += i
    end
    return s
end

将参数传递给函数是更好的风格。它可以产生更可重用的代码,并阐明输入和输出是什么。

注意

REPL 中的所有代码都在全局作用域中执行,因此在顶层定义和赋值的变量将是 **全局** 变量。在模块内的顶层作用域定义的变量也是全局变量。

在以下 REPL 会话中

julia> x = 1.0

等效于

julia> global x = 1.0

因此,前面讨论的所有性能问题都适用。

使用 @time 测量性能并注意内存分配

一个用于测量性能的有用工具是 @time 宏。我们在这里重复上面使用全局变量的示例,但这次没有类型注释

julia> x = rand(1000);

julia> function sum_global()
           s = 0.0
           for i in x
               s += i
           end
           return s
       end;

julia> @time sum_global()
  0.011539 seconds (9.08 k allocations: 373.386 KiB, 98.69% compilation time)
523.0007221951678

julia> @time sum_global()
  0.000091 seconds (3.49 k allocations: 70.156 KiB)
523.0007221951678

在第一次调用 (@time sum_global()) 时,函数将被编译。(如果您在此会话中尚未使用 @time,它还会编译计时所需的函数。)您不应该认真对待本次运行的结果。对于第二次运行,请注意,除了报告时间外,它还表明分配了大量的内存。我们这里只是计算 64 位浮点数向量中所有元素的总和,因此不需要分配(堆)内存。

我们应该澄清的是,@time 报告的具体是 *堆* 分配,这通常是可变对象或创建/增长可变大小容器(如 ArrayDict,字符串,或类型不稳定的对象,其类型仅在运行时才知道)所需的。分配(或释放)此类内存块可能需要一个昂贵的系统调用(例如,通过 C 中的 malloc),并且必须跟踪它们以进行垃圾回收。相反,不可变值(如数字(bignum 除外)、元组和不可变 struct)可以更便宜地存储,例如在堆栈或 CPU 寄存器内存中,因此人们通常不会担心“分配”它们的性能成本。

意外的内存分配几乎总是表明代码存在问题,通常是类型稳定性问题或创建了许多小型临时数组。因此,除了分配本身之外,您的函数生成的代码很可能远非最佳。认真对待这些指示,并遵循以下建议。

在本例中,内存分配是由于使用了类型不稳定的全局变量 x,因此如果我们改为将 x 作为参数传递给函数,它将不再分配内存(下面报告的剩余分配是由于在全局作用域中运行 @time 宏),并且在第一次调用后速度明显加快

julia> x = rand(1000);

julia> function sum_arg(x)
           s = 0.0
           for i in x
               s += i
           end
           return s
       end;

julia> @time sum_arg(x)
  0.007551 seconds (3.98 k allocations: 200.548 KiB, 99.77% compilation time)
523.0007221951678

julia> @time sum_arg(x)
  0.000006 seconds (1 allocation: 16 bytes)
523.0007221951678

看到的一次分配来自在全局作用域中运行 @time 宏本身。如果我们改为在函数中运行计时,我们可以看到实际上没有进行分配

julia> time_sum(x) = @time sum_arg(x);

julia> time_sum(x)
  0.000002 seconds
523.0007221951678

在某些情况下,您的函数可能需要在操作过程中分配内存,这会使以上简单的情况复杂化。在这种情况下,请考虑使用以下工具之一来诊断问题,或编写函数的版本,将分配与其算法方面分开(请参阅预分配输出)。

注意

对于更严肃的基准测试,请考虑使用BenchmarkTools.jl 包,它可以多次评估函数以减少噪声。

工具

Julia 及其包生态系统包括可以帮助您诊断问题并提高代码性能的工具

  • 性能分析 允许您测量运行代码的性能,并识别充当瓶颈的行。对于复杂的项目,ProfileView 包可以帮助您可视化性能分析结果。
  • Traceur 包可以帮助您找到代码中常见的性能问题。
  • @time@allocated 或分析器(通过调用垃圾收集例程)报告的意外大内存分配提示您的代码可能存在问题。如果您没有看到分配的其他原因,请怀疑类型问题。您也可以使用 --track-allocation=user 选项启动 Julia 并检查生成的 *.mem 文件以查看有关这些分配发生位置的信息。参见 内存分配分析
  • @code_warntype 生成代码的表示形式,这有助于查找导致类型不确定的表达式。参见下文 @code_warntype

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

在处理参数化类型(包括数组)时,最好尽可能避免使用抽象类型进行参数化。

考虑以下内容

julia> a = Real[]
Real[]

julia> push!(a, 1); push!(a, 2.0); push!(a, π)
3-element Vector{Real}:
 1
 2.0
 π = 3.1415926535897...

因为 a 是抽象类型 Real 的数组,所以它必须能够保存任何 Real 值。由于 Real 对象可以具有任意大小和结构,因此 a 必须表示为指向单独分配的 Real 对象的指针数组。但是,如果我们只允许相同类型的数字(例如 Float64)存储在 a 中,则可以更有效地存储它们

julia> a = Float64[]
Float64[]

julia> push!(a, 1); push!(a, 2.0); push!(a,  π)
3-element Vector{Float64}:
 1.0
 2.0
 3.141592653589793

现在将数字分配到 a 将把它们转换为 Float64,并且 a 将存储为一个连续的 64 位浮点值块,可以有效地进行操作。

如果您无法避免使用带有抽象值类型的容器,有时最好使用 Any 进行参数化以避免运行时类型检查。例如,IdDict{Any, Any} 的性能优于 IdDict{Type, Vector}

另请参见 参数化类型 下的讨论。

类型声明

在许多具有可选类型声明的语言中,添加声明是使代码运行更快的主要方法。在 Julia 中 *并非* 如此。在 Julia 中,编译器通常知道所有函数参数、局部变量和表达式的类型。但是,在某些特定情况下,声明会有所帮助。

避免使用带有抽象类型的字段

可以在不指定字段类型的条件下声明类型

julia> struct MyAmbiguousType
           a
       end

这允许 a 为任何类型。这通常很有用,但它确实有一个缺点:对于 MyAmbiguousType 类型的对象,编译器将无法生成高性能代码。原因是编译器使用对象的类型而不是它们的值来确定如何构建代码。不幸的是,关于 MyAmbiguousType 类型的对象,几乎无法推断出任何信息

julia> b = MyAmbiguousType("Hello")
MyAmbiguousType("Hello")

julia> c = MyAmbiguousType(17)
MyAmbiguousType(17)

julia> typeof(b)
MyAmbiguousType

julia> typeof(c)
MyAmbiguousType

bc 的值具有相同的类型,但它们在内存中的数据底层表示方式却大不相同。即使您仅在 a 字段中存储数字值,但由于 UInt8 的内存表示形式与 Float64 不同,因此 CPU 也需要使用两种不同类型的指令来处理它们。由于所需的信息在类型中不可用,因此此类决策必须在运行时进行。这会降低性能。

您可以通过声明 a 的类型来做得更好。在这里,我们关注的是 a 可能属于几种类型中的一种的情况,在这种情况下,自然解决方案是使用参数。例如

julia> mutable struct MyType{T<:AbstractFloat}
           a::T
       end

这是比以下方法更好的选择

julia> mutable struct MyStillAmbiguousType
           a::AbstractFloat
       end

因为第一个版本从包装器对象的类型指定了 a 的类型。例如

julia> m = MyType(3.2)
MyType{Float64}(3.2)

julia> t = MyStillAmbiguousType(3.2)
MyStillAmbiguousType(3.2)

julia> typeof(m)
MyType{Float64}

julia> typeof(t)
MyStillAmbiguousType

a 字段的类型可以很容易地从 m 的类型确定,但不能从 t 的类型确定。实际上,在 t 中,可以更改 a 字段的类型

julia> typeof(t.a)
Float64

julia> t.a = 4.5f0
4.5f0

julia> typeof(t.a)
Float32

相反,一旦构造 mm.a 的类型就不能更改

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float64

m.a 的类型可以从 m 的类型中得知这一事实(再加上它在函数执行期间无法更改类型的事实)使编译器能够为像 m 这样的对象生成高度优化的代码,但不能为像 t 这样的对象生成优化代码。

当然,所有这些都只在我们使用具体类型构造 m 时才成立。我们可以通过显式地使用抽象类型构造它来打破这一点

julia> m = MyType{AbstractFloat}(3.2)
MyType{AbstractFloat}(3.2)

julia> typeof(m.a)
Float64

julia> m.a = 4.5f0
4.5f0

julia> typeof(m.a)
Float32

在所有实际情况下,此类对象的运行方式与 MyStillAmbiguousType 的对象相同。

比较一个简单函数生成的代码量是非常有启发性的

func(m::MyType) = m.a+1

使用

code_llvm(func, Tuple{MyType{Float64}})
code_llvm(func, Tuple{MyType{AbstractFloat}})

由于长度原因,结果未在此处显示,但您可能希望自己尝试一下。由于在第一种情况下类型是完全指定的,因此编译器无需生成任何代码来在运行时解析类型。这会生成更短、更快的代码。

还应牢记,未完全参数化的类型表现得像抽象类型。例如,即使完全指定的 Array{T,n} 是具体的,但 Array 本身没有给出参数就不是具体的

julia> !isconcretetype(Array), !isabstracttype(Array), isstructtype(Array), !isconcretetype(Array{Int}), isconcretetype(Array{Int,1})
(true, true, true, true, true)

在这种情况下,最好避免使用 a::Array 字段来声明 MyType,而是将字段声明为 a::Array{T,N}a::A,其中 {T,N}AMyType 的参数。

避免使用带有抽象容器的字段

相同的最佳实践也适用于容器类型

julia> struct MySimpleContainer{A<:AbstractVector}
           a::A
       end

julia> struct MyAmbiguousContainer{T}
           a::AbstractVector{T}
       end

julia> struct MyAlsoAmbiguousContainer
           a::Array
       end

例如

julia> c = MySimpleContainer(1:3);

julia> typeof(c)
MySimpleContainer{UnitRange{Int64}}

julia> c = MySimpleContainer([1:3;]);

julia> typeof(c)
MySimpleContainer{Vector{Int64}}

julia> b = MyAmbiguousContainer(1:3);

julia> typeof(b)
MyAmbiguousContainer{Int64}

julia> b = MyAmbiguousContainer([1:3;]);

julia> typeof(b)
MyAmbiguousContainer{Int64}

julia> d = MyAlsoAmbiguousContainer(1:3);

julia> typeof(d), typeof(d.a)
(MyAlsoAmbiguousContainer, Vector{Int64})

julia> d = MyAlsoAmbiguousContainer(1:1.0:3);

julia> typeof(d), typeof(d.a)
(MyAlsoAmbiguousContainer, Vector{Float64})

对于 MySimpleContainer,对象由其类型和参数完全指定,因此编译器可以生成优化函数。在大多数情况下,这可能就足够了。

虽然编译器现在可以完美地完成它的工作,但有一些情况是 *您* 可能希望您的代码根据 a 的 *元素类型* 做不同的事情。通常,实现此目的的最佳方法是将您的特定操作(这里为 foo)包装在单独的函数中

julia> function sumfoo(c::MySimpleContainer)
           s = 0
           for x in c.a
               s += foo(x)
           end
           s
       end
sumfoo (generic function with 1 method)

julia> foo(x::Integer) = x
foo (generic function with 1 method)

julia> foo(x::AbstractFloat) = round(x)
foo (generic function with 2 methods)

这使事情保持简单,同时允许编译器在所有情况下生成优化代码。

但是,在某些情况下,您可能需要为不同的元素类型或 MySimpleContainer 中的 a 字段的 AbstractVector 类型声明外部函数的不同版本。您可以这样做

julia> function myfunc(c::MySimpleContainer{<:AbstractArray{<:Integer}})
           return c.a[1]+1
       end
myfunc (generic function with 1 method)

julia> function myfunc(c::MySimpleContainer{<:AbstractArray{<:AbstractFloat}})
           return c.a[1]+2
       end
myfunc (generic function with 2 methods)

julia> function myfunc(c::MySimpleContainer{Vector{T}}) where T <: Integer
           return c.a[1]+3
       end
myfunc (generic function with 3 methods)
julia> myfunc(MySimpleContainer(1:3))
2

julia> myfunc(MySimpleContainer(1.0:3))
3.0

julia> myfunc(MySimpleContainer([1:3;]))
4

注释从未类型化位置获取的值

使用可能包含任何类型的值的数据结构(Array{Any} 类型的数组)通常很方便。但是,如果您正在使用其中一个结构并且碰巧知道元素的类型,那么与编译器共享此知识会有所帮助

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

在这里,我们碰巧知道 a 的第一个元素将是 Int32。进行这样的注释还有额外的好处,如果该值不是预期类型,它会引发运行时错误,从而有可能更早地捕获某些错误。

如果 a[1] 的类型无法精确确定,则可以通过 x = convert(Int32, a[1])::Int32 来声明 x。使用 convert 函数允许 a[1] 为任何可转换为 Int32 的对象(例如 UInt8),从而通过放松类型要求来提高代码的通用性。请注意,convert 本身在这种情况下需要类型注释才能实现类型稳定性。这是因为编译器无法推断函数的返回值类型,即使是 convert,除非所有函数参数的类型都已知。

如果类型是抽象的或在运行时构造的,则类型注释不会提高(实际上可能会降低)性能。这是因为编译器无法使用注释来专门化后续代码,并且类型检查本身需要时间。例如,在以下代码中

function nr(a, prec)
    ctype = prec == 32 ? Float32 : Float64
    b = Complex{ctype}(a)
    c = (b + 1.0f0)::Complex{ctype}
    abs(c)
end

c 的注释会损害性能。要编写涉及在运行时构造的类型的性能良好的代码,请使用下面讨论的 函数屏障技术,并确保构造的类型出现在内核函数的参数类型中,以便编译器正确地专门化内核操作。例如,在上面的代码片段中,一旦构造 b,它就可以传递给另一个函数 k(内核)。例如,如果函数 kb 声明为 Complex{T} 类型的参数,其中 T 是类型参数,那么在 k 中出现的以下形式的类型注释

c = (b + 1.0f0)::Complex{T}

不会损害性能(但也不会帮助),因为编译器可以在编译 k 时确定 c 的类型。

注意 Julia 避免专门化的时间

作为一个启发式方法,Julia 避免在三种特定情况下自动对参数类型进行 专门化TypeFunctionVararg。如果参数在方法中使用,Julia 将始终进行专门化,但如果参数只是传递给另一个函数,则不会进行专门化。这通常不会对运行时性能产生影响,并且 提高了编译器性能。如果您发现它在您的情况下确实会对运行时性能产生影响,则可以通过在方法声明中添加类型参数来触发专门化。以下是一些示例

这不会进行专门化

function f_type(t)  # or t::Type
    x = ones(t, 10)
    return sum(map(sin, x))
end

但这样做会

function g_type(t::Type{T}) where T
    x = ones(T, 10)
    return sum(map(sin, x))
end

这些不会进行专门化

f_func(f, num) = ntuple(f, div(num, 2))
g_func(g::Function, num) = ntuple(g, div(num, 2))

但这样做会

h_func(h::H, num) where {H} = ntuple(h, div(num, 2))

这不会进行专门化

f_vararg(x::Int...) = tuple(x...)

但这样做会

g_vararg(x::Vararg{Int, N}) where {N} = tuple(x...)

只需要引入一个类型参数来强制进行专门化,即使其他类型不受限制。例如,这也会进行专门化,并且在参数不都具有相同类型时很有用

h_vararg(x::Vararg{Any, N}) where {N} = tuple(x...)

请注意,@code_typed 及其相关函数将始终向您显示专门化代码,即使 Julia 通常不会专门化该方法调用。如果您想查看在更改参数类型时是否会生成专门化,即 Base.specializations(@which f(...)) 是否包含针对所讨论参数的专门化,则需要检查 方法内部

将函数分解为多个定义

将函数编写为许多小的定义允许编译器直接调用最适用的代码,甚至内联它。

以下是一个“复合函数”示例,实际上应该编写为多个定义

using LinearAlgebra

function mynorm(A)
    if isa(A, Vector)
        return sqrt(real(dot(A,A)))
    elseif isa(A, Matrix)
        return maximum(svdvals(A))
    else
        error("mynorm: invalid argument")
    end
end

这可以更简洁、更有效地编写为

mynorm(x::Vector) = sqrt(real(dot(x, x)))
mynorm(A::Matrix) = maximum(svdvals(A))

但是,需要注意的是,编译器在优化掉像 mynorm 示例这样编写的代码中的死分支方面非常有效。

编写“类型稳定”的函数

如果可能,确保函数始终返回相同类型的值会有所帮助。考虑以下定义

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

虽然这看起来很无辜,但问题是 0 是一个整数(类型为 Int),而 x 可能属于任何类型。因此,根据 x 的值,该函数可能会返回两种类型之一的值。这种行为是允许的,在某些情况下可能是理想的。但它可以很容易地通过以下方式修复

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

还有一个 oneunit 函数,以及更通用的 oftype(x, y) 函数,它返回将 y 转换为 x 的类型。

避免改变变量的类型

对于在函数中反复使用的变量,存在类似的“类型稳定性”问题

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

局部变量 x 从整数开始,经过一次循环迭代后变成了浮点数(/ 运算符的结果)。这使得编译器更难以优化循环体。有几种可能的修复方法

  • 使用 x = 1.0 初始化 x
  • 显式声明 x 的类型为 x::Float64 = 1
  • 使用 x = oneunit(Float64) 进行显式转换
  • 使用第一次循环迭代进行初始化,即 x = 1 / rand(),然后循环 for i = 2:10

分离内核函数(又名函数屏障)

许多函数遵循执行一些设置工作,然后运行许多迭代来执行核心计算的模式。如果可能,将这些核心计算放在单独的函数中是个好主意。例如,以下人为编造的函数返回一个随机选择的类型数组

julia> function strange_twos(n)
           a = Vector{rand(Bool) ? Int64 : Float64}(undef, n)
           for i = 1:n
               a[i] = 2
           end
           return a
       end;

julia> strange_twos(3)
3-element Vector{Int64}:
 2
 2
 2

这应该写成

julia> function fill_twos!(a)
           for i = eachindex(a)
               a[i] = 2
           end
       end;

julia> function strange_twos(n)
           a = Vector{rand(Bool) ? Int64 : Float64}(undef, n)
           fill_twos!(a)
           return a
       end;

julia> strange_twos(3)
3-element Vector{Int64}:
 2
 2
 2

Julia 的编译器为函数边界处的参数类型专门化代码,所以在原始实现中,它在循环过程中不知道 a 的类型(因为它是由随机选择的)。因此,第二个版本通常更快,因为内部循环可以作为 fill_twos! 的一部分为 a 的不同类型重新编译。

第二种形式通常也是更好的风格,并且可以导致更多代码复用。

这种模式在 Julia Base 中的多个地方使用。例如,参见 abstractarray.jl 中的 vcathcat,或者 fill! 函数,我们本可以使用它来代替编写自己的 fill_twos!

当处理类型不确定的数据时,例如从可能包含整数、浮点数、字符串或其他内容的输入文件加载的数据,就会出现像 strange_twos 这样的函数。

带值作为参数的类型

假设你想创建一个 N 维数组,每个轴的大小都为 3。这样的数组可以通过以下方式创建

julia> A = fill(5.0, (3, 3))
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

这种方法非常有效:编译器可以确定 AArray{Float64,2},因为它知道填充值的类型(5.0::Float64)和维数((3, 3)::NTuple{2,Int})。这意味着编译器可以为同一个函数中对 A 的任何将来使用生成非常有效的代码。

但现在假设你想编写一个函数,该函数在任意维度上创建一个 3×3×... 数组;你可能会想编写一个函数

julia> function array3(fillval, N)
           fill(fillval, ntuple(d->3, N))
       end
array3 (generic function with 1 method)

julia> array3(5.0, 2)
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

这可以工作,但是(你可以使用 @code_warntype array3(5.0, 2) 为自己验证),问题在于输出类型无法推断:参数 N 是一个类型为 Int,类型推断不(也不可能)提前预测它的值。这意味着使用该函数的输出的代码必须保持谨慎,在每次访问 A 时都要检查类型;这样的代码会非常慢。

现在,解决此类问题的非常好的方法之一是使用 函数屏障技术。但是,在某些情况下,你可能希望完全消除类型不稳定性。在这种情况下,一种方法是将维数作为参数传递,例如通过 Val{T}()(参见 "值类型"

julia> function array3(fillval, ::Val{N}) where N
           fill(fillval, ntuple(d->3, Val(N)))
       end
array3 (generic function with 1 method)

julia> array3(5.0, Val(2))
3×3 Matrix{Float64}:
 5.0  5.0  5.0
 5.0  5.0  5.0
 5.0  5.0  5.0

Julia 有一个专门化的 ntuple 版本,它接受 Val{::Int} 实例作为第二个参数;通过将 N 作为类型参数传递,你让它的“值”为编译器所知。因此,这个版本的 array3 允许编译器预测返回类型。

然而,利用这些技术可能会出乎意料地微妙。例如,如果你从以下函数调用 array3,它将没有帮助

function call_array3(fillval, n)
    A = array3(fillval, Val(n))
end

在这里,你又创建了同样的问题:编译器无法猜到 n 是什么,所以它不知道 Val(n)类型。试图使用 Val,但使用不正确,很容易在许多情况下导致性能下降。(只有当你有效地将 Val 与函数屏障技巧结合使用,以提高内核函数的效率时,才应该使用上面的代码。)

使用 Val 的正确示例将是

function filter3(A::AbstractArray{T,N}) where {T,N}
    kernel = array3(1, Val(N))
    filter(A, kernel)
end

在这个例子中,N 被作为参数传递,所以它的“值”是编译器已知的。从本质上讲,Val(T) 仅在 T 是硬编码/字面量(Val(3))或在类型域中已指定的情况下才有效。

滥用多重调度的危害(又名,更多关于带值作为参数的类型)

一旦学会了欣赏多重调度,就有一种可以理解的倾向,即过度使用它,试图将其用于所有事情。例如,你可能想象使用它来存储信息,例如

struct Car{Make, Model}
    year::Int
    ...more fields...
end

然后对类似 Car{:Honda,:Accord}(year, args...) 的对象进行调度。

当以下情况之一成立时,这可能会有用

  • 你需要对每个 Car 进行 CPU 密集型处理,并且如果你在编译时知道 MakeModel,并且将使用的不同 MakeModel 的总数不是太大,那么效率将大大提高。
  • 你拥有要处理的相同类型的 Car 的同类列表,以便你可以将它们全部存储在 Array{Car{:Honda,:Accord},N} 中。

当后者成立时,处理这种同类数组的函数可以有效地进行专门化:Julia 提前知道容器中每个元素的类型(容器中的所有对象都具有相同的具体类型),因此 Julia 可以“查找”正确的方法调用当函数被编译时(避免在运行时进行检查),从而为处理整个列表发出有效的代码。

当这些不成立时,你很可能不会得到任何好处;更糟糕的是,由此产生的“类型组合爆炸”将适得其反。如果 items[i+1] 的类型与 item[i] 不同,Julia 必须在运行时查找类型,在方法表中搜索适当的方法,确定(通过类型交集)哪个方法匹配,确定它是否已进行 JIT 编译(如果没有,则进行 JIT 编译),然后进行调用。从本质上讲,你是在要求完整的类型系统和 JIT 编译机制基本上在自己的代码中执行相当于 switch 语句或字典查找的操作。

可以在 邮件列表中 找到比较 (1) 类型调度、(2) 字典查找和 (3)“switch”语句的某些运行时基准。

也许比运行时影响更糟糕的是编译时影响:Julia 将为每个不同的 Car{Make, Model} 编译专门化的函数;如果你有数百或数千种这样的类型,那么每个接受这种对象作为参数的函数(从你可能自己编写的自定义 get_year 函数,到 Julia Base 中的通用 push! 函数)将有数百或数千种变体被编译。这些都增加了已编译代码缓存的大小、方法的内部列表的长度等等。对带值作为参数的过度热情很容易浪费大量的资源。

按照内存顺序访问数组,沿着列

Julia 中的多维数组按列优先顺序存储。这意味着数组是按列依次堆叠的。可以使用 vec 函数或 [:] 语法验证这一点,如下所示(请注意,数组的顺序为 [1 3 2 4],而不是 [1 2 3 4]

julia> x = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> x[:]
4-element Vector{Int64}:
 1
 3
 2
 4

这种数组排序约定在 Fortran、Matlab 和 R 等许多语言中很常见(仅举几例)。列优先排序的替代方案是行优先排序,这是 C 和 Python(numpy)等其他语言采用的约定。记住数组的排序在循环数组时可能会产生重大的性能影响。需要记住的一条经验法则是,对于列优先数组,第一个索引变化最快。从本质上讲,这意味着如果最内层循环索引是切片表达式中出现的第一个元素,那么循环将更快。请记住,使用 : 对数组进行索引是一个隐式循环,它迭代地访问特定维度内的所有元素;例如,提取列可能比提取行更快。

考虑以下人为编造的示例。假设我们要编写一个接受 Vector 并返回一个正方形 Matrix 的函数,该矩阵的行或列填充了输入向量的副本。假设填充行还是列并不重要(也许其余代码可以很容易地相应地进行调整)。我们可能至少可以用四种方法来实现这一点(除了对内置 repeat 的推荐调用之外)

function copy_cols(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[:, i] = x
    end
    return out
end

function copy_rows(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for i = inds
        out[i, :] = x
    end
    return out
end

function copy_col_row(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for col = inds, row = inds
        out[row, col] = x[row]
    end
    return out
end

function copy_row_col(x::Vector{T}) where T
    inds = axes(x, 1)
    out = similar(Array{T}, inds, inds)
    for row = inds, col = inds
        out[row, col] = x[col]
    end
    return out
end

现在我们将使用相同的随机 10000×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_colscopy_rows 快得多。这是预期的,因为 copy_cols 尊重 Matrix 的基于列的内存布局,并一次填充一列。此外,copy_col_rowcopy_row_col 快得多,因为它遵循了我们的经验法则,即切片表达式中出现的第一个元素应该与最内层循环相结合。

预先分配输出

如果你的函数返回一个 Array 或其他一些复杂类型,它可能必须分配内存。不幸的是,分配及其反面,垃圾回收,通常是重大的瓶颈。

有时,你可以在每次函数调用时通过预先分配输出来避免分配内存的需要。作为一个简单的示例,比较

julia> function xinc(x)
           return [x, x+1, x+2]
       end;

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

julia> function xinc!(ret::AbstractVector{T}, x::T) where T
           ret[1] = x
           ret[2] = x+1
           ret[3] = x+2
           nothing
       end;

julia> function loopinc_prealloc()
           ret = Vector{Int}(undef, 3)
           y = 0
           for i = 1:10^7
               xinc!(ret, i)
               y += ret[2]
           end
           return y
       end;

计时结果

julia> @time loopinc()
  0.529894 seconds (40.00 M allocations: 1.490 GiB, 12.14% gc time)
50000015000000

julia> @time loopinc_prealloc()
  0.030850 seconds (6 allocations: 288 bytes)
50000015000000

预分配还有其他优点,例如允许调用者控制算法的“输出”类型。在上面的示例中,如果需要,我们可以传递一个SubArray而不是一个Array

如果将预分配发挥到极致,可能会使代码更难看,因此可能需要进行性能衡量和一些判断。但是,对于“向量化”(逐元素)函数,可以使用方便的语法x .= f.(y)进行就地操作,并使用融合循环和没有临时数组(请参阅用于向量化函数的点语法)。

更多点:融合向量化操作

Julia 有一个特殊的点语法,它将任何标量函数转换为“向量化”函数调用,并将任何运算符转换为“向量化”运算符,具有嵌套“点调用”是融合的特殊属性:它们在语法级别上合并到一个循环中,而不会分配临时数组。如果您使用.=和类似的赋值运算符,结果也可以就地存储在预分配的数组中(见上文)。

在线性代数环境中,这意味着即使定义了vector + vectorvector * scalar等运算,使用vector .+ vectorvector .* scalar也可能更有利,因为生成的循环可以与周围的计算融合。例如,考虑以下两个函数

julia> f(x) = 3x.^2 + 4x + 7x.^3;

julia> fdot(x) = @. 3x^2 + 4x + 7x^3; # equivalent to 3 .* x.^2 .+ 4 .* x .+ 7 .* x.^3

ffdot都计算相同的值。但是,fdot(在@.宏的帮助下定义)在应用于数组时速度明显更快

julia> x = rand(10^6);

julia> @time f(x);
  0.019049 seconds (16 allocations: 45.777 MiB, 18.59% gc time)

julia> @time fdot(x);
  0.002790 seconds (6 allocations: 7.630 MiB)

julia> @time f.(x);
  0.002626 seconds (8 allocations: 7.630 MiB)

也就是说,fdot(x)的速度是f(x)的十倍,并且分配的内存是f(x)的 1/6,因为f(x)中的每个*+运算都会分配一个新的临时数组并在单独的循环中执行。在本例中,f.(x)fdot(x)的速度一样快,但在许多情况下,在表达式中添加一些点比为每个向量化操作定义一个单独的函数更方便。

考虑对切片使用视图

在 Julia 中,类似于array[1:5, :]的数组“切片”表达式会创建该数据的副本(除了在赋值的左侧,其中array[1:5, :] = ...就地分配给array的该部分)。如果您对切片执行了许多操作,这可能会提高性能,因为处理较小的连续副本比索引到原始数组效率更高。另一方面,如果您只对切片执行几个简单的操作,则分配和复制操作的成本可能很高。

另一种方法是创建数组的“视图”,它是一个数组对象(SubArray),实际上就地引用了原始数组的数据,而无需进行复制。(如果您写入视图,它也会修改原始数组的数据。)这可以通过调用view来对单个切片执行,或者更简单地通过将@views放在该表达式之前来对整个表达式或代码块执行。例如

julia> fcopy(x) = sum(x[2:end-1]);

julia> @views fview(x) = sum(x[2:end-1]);

julia> x = rand(10^6);

julia> @time fcopy(x);
  0.003051 seconds (3 allocations: 7.629 MB)

julia> @time fview(x);
  0.001020 seconds (1 allocation: 16 bytes)

请注意,fview版本的函数速度提高了 3 倍,并且内存分配也减少了。

复制数据并不总是坏事

数组在内存中连续存储,这有利于 CPU 向量化和由于缓存造成的内存访问次数减少。这些与建议以列优先顺序访问数组的原因相同(见上文)。不规则的访问模式和非连续视图会导致对数组的计算速度急剧下降,因为内存访问不是连续的。

将不规则访问的数据复制到连续数组中,然后进行重复访问可能会导致速度大幅提升,例如以下示例。这里,一个矩阵在被乘以之前是在随机排序的索引处访问的。即使增加了复制和分配的成本,将数据复制到普通数组中也能加速乘法运算。

julia> using Random

julia> A = randn(3000, 3000);

julia> x = randn(2000);

julia> inds = shuffle(1:3000)[1:2000];

julia> function iterated_neural_network(A, x, depth)
           for _ in 1:depth
               x .= max.(0, A * x)
           end
           argmax(x)
       end

julia> @time iterated_neural_network(view(A, inds, inds), x, 10)
  0.324903 seconds (12 allocations: 157.562 KiB)
1569

julia> @time iterated_neural_network(A[inds, inds], x, 10)
  0.054576 seconds (13 allocations: 30.671 MiB, 13.33% gc time)
1569

只要有足够的内存,将视图复制到数组的成本会被重复矩阵乘法在连续数组上执行带来的速度提升所抵消。

考虑对小型固定大小的向量/矩阵操作使用 StaticArrays.jl

如果您的应用程序涉及许多大小很小的(< 100 个元素)固定大小的数组(即大小在执行之前已知),那么您可能需要考虑使用StaticArrays.jl 包。此包允许您以一种避免不必要的堆分配的方式表示这些数组,并允许编译器针对数组的大小专门化代码,例如通过完全展开向量运算(消除循环)并将元素存储在 CPU 寄存器中。

例如,如果您正在进行 2D 几何计算,您可能需要对包含 2 个分量的向量进行许多计算。通过使用 StaticArrays.jl 中的SVector类型,您可以使用方便的向量表示法和运算,例如在向量vw上使用norm(3v - w),同时允许编译器将代码展开为最小计算,等效于@inbounds hypot(3v[1]-w[1], 3v[2]-w[2])

避免使用字符串插值进行 I/O

当将数据写入文件(或其他 I/O 设备)时,形成额外的中间字符串会带来开销。与其使用

println(file, "$a $b")

不如使用

println(file, a, " ", b)

代码的第一个版本会形成一个字符串,然后将其写入文件,而第二个版本则直接将值写入文件。此外,请注意,在某些情况下,字符串插值可能更难阅读。考虑

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

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

优化并行执行期间的网络 I/O

当并行执行远程函数时

using Distributed

responses = Vector{Any}(undef, nworkers())
@sync begin
    for (idx, pid) in enumerate(workers())
        @async responses[idx] = remotecall_fetch(foo, pid, args...)
    end
end

using Distributed

refs = Vector{Any}(undef, nworkers())
for (idx, pid) in enumerate(workers())
    refs[idx] = @spawnat pid foo(args...)
end
responses = [fetch(r) for r in refs]

速度更快。前者会导致单个网络往返于每个工作器,而后者会导致两次网络调用 - 第一次由@spawnat执行,第二次由fetch(甚至wait)执行。fetch/wait也是串行执行的,导致整体性能下降。

修复弃用警告

一个已弃用的函数会在内部执行查找,以便只打印一次相关的警告。此额外的查找会导致明显的减速,因此应根据警告建议修改所有对已弃用函数的使用。

调整

以下是一些在紧凑的内部循环中可能有所帮助的小要点。

性能注释

有时,您可以通过承诺某些程序属性来启用更好的优化。

  • 使用@inbounds在表达式内消除数组边界检查。在执行此操作之前请务必确认。如果索引超出了范围,您可能会遇到崩溃或静默损坏。
  • 使用@fastmath允许对浮点数进行优化,这些优化对于实数是正确的,但会导致 IEEE 数的差异。在执行此操作时要小心,因为这可能会改变数值结果。这对应于 clang 的-ffast-math选项。
  • for循环之前写入@simd,以承诺迭代是独立的,并且可以重新排序。请注意,在许多情况下,Julia 可以自动向量化代码,而无需@simd宏;它只在无法进行这种转换的情况下才有益,包括允许浮点数重新关联和忽略依赖的内存访问(@simd ivdep)的情况。同样,在断言@simd时要非常小心,因为错误地用依赖迭代注释循环会导致意外结果。特别是,请注意,某些AbstractArray子类型的setindex!本质上依赖于迭代顺序。此功能尚处于实验阶段,可能在 Julia 的未来版本中发生变化或消失。

使用 1:n 索引到 AbstractArray 的通用习惯用法如果 Array 使用非传统的索引是不安全的,并且如果关闭边界检查可能会导致段错误。请改用LinearIndices(x)eachindex(x)(另请参阅具有自定义索引的数组)。

注意

虽然@simd需要直接放置在最内层的for循环之前,但@inbounds@fastmath都可以应用于单个表达式或出现在嵌套代码块中的所有表达式,例如,使用@inbounds begin@inbounds for ...

以下是一个使用@inbounds@simd标记的示例(我们在这里使用@noinline来阻止优化器尝试过于聪明并破坏我们的基准测试)

@noinline function inner(x, y)
    s = zero(eltype(x))
    for i=eachindex(x)
        @inbounds s += x[i]*y[i]
    end
    return s
end

@noinline function innersimd(x, y)
    s = zero(eltype(x))
    @simd for i = eachindex(x)
        @inbounds s += x[i] * y[i]
    end
    return 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/sec        = ", 2n*reps / time*1E-9)
    time = @elapsed for j in 1:reps
        s += innersimd(x, y)
    end
    println("GFlop/sec (SIMD) = ", 2n*reps / time*1E-9)
end

timeit(1000, 1000)

在一台配备 2.4GHz 英特尔酷睿 i5 处理器的计算机上,这将生成

GFlop/sec        = 1.9467069505224963
GFlop/sec (SIMD) = 17.578554163920018

GFlop/sec衡量性能,数字越大越好。)

以下是一个使用所有三种标记类型的示例。该程序首先计算一维数组的有限差分,然后评估结果的 L2 范数

function init!(u::Vector)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds @simd for i in 1:n #by asserting that `u` is a `Vector` we can assume it has 1-based indexing
        u[i] = sin(2pi*dx*i)
    end
end

function deriv!(u::Vector, du)
    n = length(u)
    dx = 1.0 / (n-1)
    @fastmath @inbounds du[1] = (u[2] - u[1]) / dx
    @fastmath @inbounds @simd for i in 2:n-1
        du[i] = (u[i+1] - u[i-1]) / (2*dx)
    end
    @fastmath @inbounds du[n] = (u[n] - u[n-1]) / dx
end

function mynorm(u::Vector)
    n = length(u)
    T = eltype(u)
    s = zero(T)
    @fastmath @inbounds @simd for i in 1:n
        s += u[i]^2
    end
    @fastmath @inbounds return sqrt(s)
end

function main()
    n = 2000
    u = Vector{Float64}(undef, n)
    init!(u)
    du = similar(u)

    deriv!(u, du)
    nu = mynorm(du)

    @time for i in 1:10^6
        deriv!(u, du)
        nu = mynorm(du)
    end

    println(nu)
end

main()

在一台配备 2.7 GHz 英特尔酷睿 i7 处理器的计算机上,这将生成

$ julia wave.jl;
  1.207814709 seconds
4.443986180758249

$ julia --math-mode=ieee wave.jl;
  4.487083643 seconds
4.443986180758249

这里,选项--math-mode=ieee禁用了@fastmath宏,以便我们可以比较结果。

在这种情况下,@fastmath带来的加速效果约为 3.7 倍。这很不寻常,通常情况下,加速效果会更小。(在这个特定的例子中,基准测试的工作集足够小,可以完全放入处理器的 L1 缓存中,因此内存访问延迟不会产生影响,计算时间主要由 CPU 使用率决定。在许多现实世界中的程序中,情况并非如此。)此外,在这种情况下,此优化不会改变结果,一般来说,结果会略有不同。在某些情况下,特别是对于数值不稳定的算法,结果可能会有很大不同。

注释 @fastmath 会重新排列浮点数表达式,例如改变求值的顺序,或假设某些特殊情况(无穷大、NaN)不会发生。在这种情况下(以及在这台特定的计算机上),主要的区别在于 deriv 函数中的表达式 1 / (2*dx) 被提升到循环之外(即在循环之外计算),就好像你写了 idx = 1 / (2*dx) 一样。在循环中,表达式 ... / (2*dx) 就会变成 ... * idx,这在计算上要快得多。当然,编译器实际应用的优化以及由此产生的加速效果在很大程度上取决于硬件。可以使用 Julia 的 code_native 函数来查看生成代码的变化。

请注意,@fastmath 还假设在计算过程中不会出现 NaN,这会导致意外的行为。

julia> f(x) = isnan(x);

julia> f(NaN)
true

julia> f_fast(x) = @fastmath isnan(x);

julia> f_fast(NaN)
false

将次正规数视为零

次正规数(以前称为非正规数)在许多情况下很有用,但在某些硬件上会造成性能损失。调用set_zero_subnormals(true) 允许浮点运算将次正规输入或输出视为零,这可能会在某些硬件上提高性能。调用 set_zero_subnormals(false) 强制执行次正规数的严格 IEEE 行为。

下面是一个次正规数在某些硬件上明显影响性能的例子

function timestep(b::Vector{T}, a::Vector{T}, Δt::T) where T
    @assert length(a)==length(b)
    n = length(b)
    b[1] = 1                            # Boundary condition
    for i=2:n-1
        b[i] = a[i] + (a[i-1] - T(2)*a[i] + a[i+1]) * Δt
    end
    b[n] = 0                            # Boundary condition
end

function heatflow(a::Vector{T}, nstep::Integer) where T
    b = similar(a)
    for t=1:div(nstep,2)                # Assume nstep is even
        timestep(b,a,T(0.1))
        timestep(a,b,T(0.1))
    end
end

heatflow(zeros(Float32,10),2)           # Force compilation
for trial=1:6
    a = zeros(Float32,1000)
    set_zero_subnormals(iseven(trial))  # Odd trials use strict IEEE arithmetic
    @time heatflow(a,1000)
end

这将产生类似于以下内容的输出

  0.002202 seconds (1 allocation: 4.063 KiB)
  0.001502 seconds (1 allocation: 4.063 KiB)
  0.002139 seconds (1 allocation: 4.063 KiB)
  0.001454 seconds (1 allocation: 4.063 KiB)
  0.002115 seconds (1 allocation: 4.063 KiB)
  0.001455 seconds (1 allocation: 4.063 KiB)

请注意,每个偶数迭代都明显更快。

此示例生成许多次正规数,因为 a 中的值变成了指数递减的曲线,随着时间的推移逐渐变平。

将次正规数视为零应该谨慎使用,因为这样做会破坏某些恒等式,例如 x-y == 0 意味着 x == y

julia> x = 3f-38; y = 2f-38;

julia> set_zero_subnormals(true); (x - y, x == y)
(0.0f0, false)

julia> set_zero_subnormals(false); (x - y, x == y)
(1.0000001f-38, false)

在某些应用中,将次正规数清零的替代方法是在其中注入少许噪声。例如,不是将 a 初始化为零,而是将它初始化为

a = rand(Float32,1000) * 1.f-9

@code_warntype

@code_warntype(或其函数变体 code_warntype)有时可以帮助诊断与类型相关的問題。以下是一个例子

julia> @noinline pos(x) = x < 0 ? 0 : x;

julia> function f(x)
           y = pos(x)
           return sin(y*x + 1)
       end;

julia> @code_warntype f(3.2)
MethodInstance for f(::Float64)
  from f(x) @ Main REPL[9]:1
Arguments
  #self#::Core.Const(f)
  x::Float64
Locals
  y::Union{Float64, Int64}
Body::Float64
1 ─      (y = Main.pos(x))
│   %2 = (y * x)::Float64
│   %3 = (%2 + 1)::Float64
│   %4 = Main.sin(%3)::Float64
└──      return %4

解释 @code_warntype 的输出,就像解释其同类 @code_lowered@code_typed@code_llvm@code_native 的输出一样,需要一些练习。你的代码以已被大量处理的形式呈现,以便生成编译后的机器代码。大多数表达式都带有类型注释,用 ::T 表示(其中 T 可以是 Float64,例如)。@code_warntype 最重要的特点是,非具体类型以红色显示;由于本文档是用 Markdown 编写的,没有颜色,因此在本文件中,红色文本用大写字母表示。

在顶部,函数的推断返回值类型显示为 Body::Float64。接下来的几行代表 f 的主体,用 Julia 的 SSA IR 形式表示。编号的方框是标签,代表代码中跳转(通过 goto)的目标。查看主体,你可以看到,首先会调用 pos,并且推断的返回值类型为 Union 类型 Union{Float64, Int64},因为它是非具体类型,所以用大写字母表示。这意味着我们无法根据输入类型确定 pos 的确切返回值类型。但是,无论 yFloat64 还是 Int64y*x 的结果都是 Float64。最终结果是,f(x::Float64) 的输出将不会出现类型不稳定,即使某些中间计算是类型不稳定的。

如何使用这些信息取决于你。显然,最好是将 pos 更改为类型稳定的:如果你这样做,f 中的所有变量都将是具体的,并且其性能将达到最佳。但是,在某些情况下,这种短暂的类型不稳定可能并不重要:例如,如果 pos 从未单独使用,那么 f 的输出类型稳定(对于 Float64 输入)将保护后面的代码免受类型不稳定传播的影响。这在修复类型不稳定很困难或不可能的情况下尤其相关。在这种情况下,上面的技巧(例如,添加类型注释和/或分解函数)是控制类型不稳定“损害”的最佳工具。此外,请注意,即使是 Julia Base 也包含类型不稳定的函数。例如,函数 findfirst 返回数组中找到键的索引,如果未找到则返回 nothing,这是一个明显的类型不稳定。为了便于查找可能重要的类型不稳定,包含 missingnothingUnion 将以黄色突出显示,而不是红色。

以下示例可能有助于你解释标记为包含非叶类型的值的表达式。

  • 函数主体以 Body::Union{T1,T2}) 开头

    • 解释:具有不稳定返回值类型的函数
    • 建议:即使你必须添加注释,也要使返回值类型稳定
  • invoke Main.g(%%x::Int64)::Union{Float64, Int64}

    • 解释:调用类型不稳定的函数 g
    • 建议:修复函数,或者在必要时对返回值进行注释
  • invoke Base.getindex(%%x::Array{Any,1}, 1::Int64)::Any

    • 解释:访问类型不佳的数组的元素
    • 建议:使用类型定义更好的数组,或者在必要时对单个元素访问的类型进行注释
  • Base.getfield(%%x, :(:data))::Array{Float64,N} where N

    • 解释:获取非叶类型的字段。在本例中,x 的类型(例如 ArrayContainer)有一个字段 data::Array{T}。但 Array 也需要维度 N 才能成为具体类型。
    • 建议:使用具体的类型,例如 Array{T,3}Array{T,N},其中 N 现在是 ArrayContainer 的参数。

捕获变量的性能

考虑以下定义内部函数的示例

function abmult(r::Int)
    if r < 0
        r = -r
    end
    f = x -> x * r
    return f
end

函数 abmult 返回一个函数 f,该函数将它的参数乘以 r 的绝对值。分配给 f 的内部函数被称为“闭包”。内部函数也被语言用于 do 块和生成器表达式。

这种代码风格对语言提出了性能挑战。解析器在将其转换为低级指令时,会大幅重组上述代码,将内部函数提取到单独的代码块中。由内部函数及其封闭范围共享的“捕获”变量(例如 r)也会被提取到一个堆分配的“框”中,该框可供内部函数和外部函数访问,因为语言指定内部范围中的 r 必须与外部范围中的 r 相同,即使外部范围(或另一个内部函数)修改了 r 也是如此。

前一段中提到的“解析器”指的是编译的阶段,该阶段在包含 abmult 的模块首次加载时发生,而不是在它首次调用时发生的后续阶段。解析器“不知道” Int 是一个固定类型,也不知道语句 r = -rInt 转换为另一个 Int。类型推断的魔力发生在编译的后续阶段。

因此,解析器不知道 r 具有固定类型(Int),也不知道 r 在创建内部函数后不会改变值(因此不需要框)。因此,解析器会为框发出代码,该框包含一个具有抽象类型(例如 Any)的对象,这需要对 r 的每次出现进行运行时类型调度。可以通过将 @code_warntype 应用于上述函数来验证这一点。框和运行时类型调度都会导致性能下降。

如果捕获的变量用于代码的性能关键部分,那么以下技巧有助于确保其使用是高性能的。首先,如果已知捕获的变量不会改变其类型,那么可以使用类型注释显式声明这一点(在变量上,而不是在右侧)

function abmult2(r0::Int)
    r::Int = r0
    if r < 0
        r = -r
    end
    f = x -> x * r
    return f
end

类型注释部分恢复了由于捕获而导致的性能下降,因为解析器可以将具体类型与框中的对象关联起来。更进一步,如果捕获的变量根本不需要框(因为它在创建闭包后不会被重新分配),那么可以使用 let 块来指示这一点。

function abmult3(r::Int)
    if r < 0
        r = -r
    end
    f = let r = r
            x -> x * r
    end
    return f
end

let 块创建了一个新的变量 r,它的范围仅限于内部函数。第二种技术在存在捕获变量的情况下恢复了完整的语言性能。请注意,这是编译器快速发展的方面,未来的版本可能不需要这种程度的程序员注释来实现性能。在此期间,一些用户贡献的软件包(例如 FastClosures)可以自动插入 let 语句,就像在 abmult3 中一样。

多线程和线性代数

本节适用于多线程 Julia 代码,在每个线程中,它执行线性代数运算。实际上,这些线性代数运算涉及 BLAS/LAPACK 调用,这些调用本身就是多线程的。在这种情况下,必须确保核心不会因两种不同类型的多线程而被过度订阅。

Julia 编译并使用自己的 OpenBLAS 副本进行线性代数运算,其线程数量由环境变量 OPENBLAS_NUM_THREADS 控制。它可以在启动 Julia 时作为命令行选项设置,或者在 Julia 会话期间使用 BLAS.set_num_threads(N) 修改(子模块 BLASusing LinearAlgebra 导出)。可以使用 BLAS.get_num_threads() 访问其当前值。

当用户没有指定任何内容时,Julia 会尝试为 OpenBLAS 线程数量选择一个合理的值(例如,基于平台、Julia 版本等)。但是,通常建议手动检查和设置该值。OpenBLAS 的行为如下

  • 如果 OPENBLAS_NUM_THREADS=1,OpenBLAS 使用调用 Julia 线程,即它“驻留在”运行计算的 Julia 线程中。
  • 如果 OPENBLAS_NUM_THREADS=N>1,OpenBLAS 会创建并管理自己的线程池(总共 N 个)。所有 Julia 线程之间共享一个 OpenBLAS 线程池。

当您以多线程模式启动 Julia,使用 JULIA_NUM_THREADS=X 时,通常建议设置 OPENBLAS_NUM_THREADS=1。鉴于上述行为,将 BLAS 线程数量增加到 N>1 非常容易导致性能下降,尤其是在 N<<X 的情况下。但是,这只是一个经验法则,设置每个线程数量的最佳方法是在您的特定应用程序上进行实验。

替代线性代数后端

作为 OpenBLAS 的替代方案,还存在其他几个后端,它们可以帮助提高线性代数性能。突出示例包括 MKL.jlAppleAccelerate.jl.

这些是外部包,因此我们在此不再详细讨论。请参阅它们各自的文档(尤其是因为它们在多线程方面与 OpenBLAS 的行为不同)。