Julia教程13 如果写出高性能的Julia代码

最新推荐文章于 2025-03-26 12:57:44 发布

原创最新推荐文章于 2025-03-26 12:57:44 发布 · 1.7k 阅读

10 ·

CC 4.0 BY-SA版权

文章标签：

#Julia #Julia教程 #性能优化

Julia 专栏收录该内容

14 篇文章

订阅专栏

文章目录

高性能代码

在网易云课堂上直接搜索：Julia教程，就可以找到，教程的全名是：Julia教程从入门到进阶

这是国内第一个免费的完整的Julia视频教程，非常适合Julia的入门。有兴趣的朋友可以去学习一下。

教程链接：
Julia教程

欢迎关注微信公众号：Quant_Times

在这里插入图片描述

文章目录

高性能代码

高性能代码

避免全局变量

全局变量的值和类型随时都会发生变化。这使编译器难以优化使用全局变量的代码。变量应该是局部的，或者尽可能作为参数传递给函数。

任何注重性能或者需要测试性能的代码都应该被放置在函数之中。

把全局变量声明为常量可以巨大的提升性能。

const VAL = 0

如果必须要声明全局变量，可以在使用它的地方标注他们的类型来优化效率。

global x = rand(10000)

function lp1()
    s = 0.0
    for i in x
        s += i
    end
end

使用@time来估算代码运行时间

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-3pky5osj-1592662804747)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-1.png)]

在函数第一次运行时，由于jit的原因，需要预热，第二次再运行的结果是真正的代码运行时间。

以下的@time结果都是不包含jit的时间。

在图中，我们不仅可以知道代码运行时间，还可以看出代码占用的内存大小。从内存大小上，我们也可以大概看出程序是不是存在问题。

比如函数lp1()申请了733KB的内存，仅仅是计算64位浮点数加法，说明肯定是有多余的空间可以节省的。

下面我们用传递函数参数和指定变量类型的方式运行

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

function lp3(x)
    s = 0.0
    for i in x
        s += i
    end
end

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-W6LFYfDu-1592662804761)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-2.png)]

可以看出，运行时间大为减少，但还是占用了160Byte的内存，这是由于在全局作用域中运行@time导致的。

如果我们把测试代码放置到函数之中，就不存在这个问题。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MrG6scLK-1592662804769)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-3.png)]

对于更加正式的性能测试，可以使用BenchmarkTools.jl包，这个包会多次评估函数的性能以降低噪声。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CGPvp1Z7-1592662804773)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-4.png)]

可以看出，指定变量类型的方式是最快的。

code generation

在REPL中，我们输入@code来查看用于code generation的宏

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-OEcYYwi9-1592662804775)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-5.png)]

@code_lowered: Julia底层的运行过程；
@code_typed: 程序运行时type的变化；
@code_llvm： llvm编译器的运行过程；
@code_native：生成程序运行的机器语言；
@code_warntypes：查看程序运行时是否有类型上的warning。

首先我们定义一个简单的函数

cal(a, b, c) = a + b * c
cal(1, 2, 3)
>>7
cal(1.0, 2.0, 3.0)
>>7.0
cal(1, 2, 3.0)
>>7.0

用@code_lowerd查看底层运行过程

@code_lowered cal(1, 2, 3)
@code_lowered cal(1.0, 2.0, 3.0)
@code_lowered cal(1, 2, 3.0)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vMfVUCM2-1592662804777)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-6.png)]

可以看出，三个函数运行过程一样，必须是一样的，只是类型不同，但乘和加的过程还是相同的。

再用code_typed查看运行时类型变化

@code_typed cal(1, 2, 3)
@code_typed cal(1.0, 2.0, 3.0)
@code_typed cal(1, 2, 3.0)

可以看出，前两个运行过程一致，只是输入变量参数不同，计算的函数也相应变化；而第三个函数运行时，多了把两个Int型的变量转为Float变量的过程。

再用@code_llvm查看llvm编译器的运行过程

@code_llvm cal(1, 2, 3)
@code_llvm cal(1.0, 2.0, 3.0)
@code_llvm cal(1, 2, 3.0)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-gMn49poM-1592662804782)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-9.png)]

可以看出，llvm编译器在对前两个函数运行时，过程基本一致，第三个函数又多了很多步骤。

同样的，用@code_native查看机器语言的运行过程。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-nhiSphZy-1592662804783)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-10.png)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-m605tGT3-1592662804785)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-11.png)]

跟@code_llvm的结果类似，前两个函数过程基本一致，第三个函数多出很多步骤。

了解了code generation对于我们后续的讲解有所帮助。

抽象和具体类型

当我们定义一个函数时，如果函数参数的类型是固定的，比如是一个Int64的Array[1,2,3,4],那他们在内存中会连续存放；

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zuxbaBoc-1592662804787)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-12.png)]

但如果函数参数的类型是Any，那么内存中连续存放只是他们的“指针”,会指向其实际的位置。这样一来，数据存取就慢下来了。

看下面的例子。

import Base: +, *
struct IntNum
    n::Int64
end

struct AnyNum
    n::Number
end

+(a::IntNum, b::IntNum) = IntNum(a.n + b.n)
*(a::IntNum, b::IntNum) = IntNum(a.n * b.n)

+(a::AnyNum, b::AnyNum) = AnyNum(a.n + b.n)
*(a::AnyNum, b::AnyNum) = AnyNum(a.n * b.n)

可以用@code_native +(IntNum(1), IntNum(2))和@code_native +(AnyNum(1), AnyNum(2))来对比两个函数的机器语言的长度，可以明显看出AnyNum的操作要比IntNum多很多操作。

function calc_Int_time(a,b,c,n)
    for i in 1:n
        cal(IntNum(a), IntNum(b), IntNum(c))
    end
end

function calc_Any_time(a,b,c,n)
    for i in 1:n
        cal(AnyNum(a), AnyNum(b), AnyNum(c))
    end
end

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CxtHWEjl-1592662804792)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-14.png)]

从程序运行时所占的内存大小也可以看出IntNum要比AnyNum少很多。
所以，计算concrete类型会比计算abstract类型要节省时间。

我们可以使用@code_warntype来查看运行的函数中是否有abstract类型

对于concrete类型：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VnAtKBa0-1592662804793)(https://raw.githubusercontent.com/Bounce00/pic/master/Julia course/13-15.png)]

对于abstract类型：

对于有abstract类型的地方，会用红色的标出。

再举一个Julia自带函数的例子。

隐藏的类型转换

在C++中，对每个定义的变量都有其固定的类型，但Julia中由于变量定义时可以缺省参数，经常会注意不到参数类型的转换。

function cums(n)
    r = 0
    for i in 1:n
        r += cos(i)
    end
    r
end

function cums_f(n)
    r = 0.0
    for i in 1:n
        r += cos(i)
    end
    r
end

把两个函数放在一起时，有过编程基础的同学可以很容易看出，第一个函数中r是个Int64类型的，程序运行时需要把Int64转成Float64再进行加法计算。但如果没有把这两个函数放在一起，会不会被很多同学忽略掉？

对比两个函数的类型warning

对比两个函数的运行时间

因此我们在定义变量时，要尽量保持与其后面运算时的类型一致。有几个函数可以帮助我们完成这种定义。

zero(value)
eltype(array)
one(value)
similar(array)

我们在定义变量时，可以使用这4个函数来避免常范的错误。

eg.

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

eg.求一个array的和

x = [1.0, 2.1, 4.5]
sm = zero(x[1])
sm = sum(x)    # 当然可以不用定义sm，直接写这一句，这只是个例子

其他三个的用法基本相同，这里不再举例。

避免拥有抽象类型参数的容器

在定义struct时，我们经常会这样写

mutable struct MyType1
           x::AbstractFloat
       end

但用下面的写法会更好

mutable struct MyType2{T<:AbstractFloat}
           x::T
       end

查看他们的类型

a = MyType1(2.3)
typeof(a)
>>MyType1
b = MyType2(2.3)
typeof(b)
>>MyType2{Float64}

x的类型可以由b的类型决定，但不能由a的类型决定。b.x的类型不会改变，永远都是Float64，而a.x的类型则会改变。

typeof(a.x)

a.x = 2.3f0
typeof(a.x)
>>Float32

b.x = 2.3f0
typeof(b.x)
>>Float64

当某个包含多个变量的Array中，如果我们知道其中某个变量的类型，就明确指定出来

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

如果只能确定x的类型，但a[1]的类型不能确定，可以用convert()来完成

x = convert(Int32, a[1])::Int32

从上面我们讲的这些内容也可以知道优化代码的一个策略：程序越简单越好，让编译器明确知道自己想干什么，而不是让编译器去猜我们的目的

因此可以推断出有可变参数类型的函数肯定不如固定关键字参数类型的函数运行的快。

但这种指定某个变量类型在使用时需要注意一点，就是如果变量类型不是在编译时确定而是在运行时才确定，那会有损性能

function nr(n, prec)
    ctype = prec == 32 ? Float32 : Float64
    b = Complex{ctype}(2.3)
    d = 0
    for i in 1:n
        c = (b + 1.0f0)::Complex{ctype}
        d += c
    end
    d
end

那碰到这种需要在运行时才能确定类型的情况应该怎么处理才能使性能最佳？

function nr_op(n, prec)
    ctype = prec == 32 ? Float32 : Float64
    b = Complex{ctype}(2.3)
    d = opFunc(n, b::Complex{ctype})
end

function opFunc(n, b::Complex{T}) where {T}
    d = 0
    for i in 1:n
        c = (b + 1.0f0)::Complex{T}
        d += c
    end
    d
end

@time nr(10000, 64)
>>0.003382 seconds (20.00 k allocations: 625.188 KiB)

@time nr_op(10000, 64)
>>0.000023 seconds (7 allocations: 240 bytes)

我们将使用这种类型的变量的计算在函数值中完成，因为在函数调用时，会先确定输入参数的类型，这样在计算时就无需再去考虑类型的不确定问题了。如果不使用函数，则需要在每次赋值后再获取变量类型，再将结果转成对应类型。这样花费的时间会多很多。

简言之，要先获取变量类型再计算

参数类型优化

当我们的类型是函数参数时，也要用上面的方式先获取

# 定义一个N维矩阵
function array3(fillval, N)
           fill(fillval, ntuple(d->3, N))
       end

typeof(array3(2.0, 4))
>>Array{Float64,4}

这里的N是矩阵的类型参数，该矩阵是一个N维的。同样存在的问题是，该矩阵变量的类型参数就是N的值，如果我们先获取了N的值后再进行矩阵生成，性能会更好。

function array3(fillval, ::Val{N}) where N
           fill(fillval, ntuple(d->3, Val(N)))
       end

在使用Val()时，一定要注意使用方式，如果使用不当，则比不使用Val()性能更差

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

这样编译器就不知道n是个什么东西，也不知道n的类型，跟我们的初衷南辕北辙。

我们再举个可变参数和关键字参数的例子说明两者的运行时间的差异

function func1(n, x...)
    res = 0
    for i in 1:n
        res += x[2]
    end
    res
end

function func2(n, x, y)
    res = 0
    for i in 1:n
        res += y
    end
    res
end

@time func1(1000000000, 10, 10)
>> 0.408809 seconds (5 allocations: 176 bytes)

@time func2(1000000000, 10, 10)
>> 0.000003 seconds (5 allocations: 176 bytes)

把参数明确后，运行速度快了好几个数量级。

用方法代替函数中的判断条件

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如果写成下面的形式效率会更高一些

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

矩阵优化

在Julia中，多维矩阵是以列优先原则排列，这跟MATLAB中是一样的

x = [1 2; 3 4]

# 把x转换为1维矩阵
x[:]

也就是说，Julia中矩阵的每一列的数据在内存上的地址是连续的，每一行的地址不是连续的，操作连续地址比非连续地址速度要快很多。下面举一个矩阵拷贝的例子。

# 按列拷贝
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

x = randn(10000)
fmt(f) = println(rpad(string(f)*": ", 14, ' '), @elapsed f(x))
map(fmt, Any[copy_cols, copy_rows, copy_col_row, copy_row_col])
>>copy_cols:    0.205115039
 copy_rows:     0.856341406
 copy_col_row:  0.279589601
 copy_row_col:  0.845328085

dot运算

在矩阵操作中，操作符前加上.，看下面的例子

f(x) = 5x.^3 - 2x.^2 + 3x
fdot(x) = @. 5x^3 - 2x^2 + 3x

x = rand(10^6);
@time f(x);
>>0.040295 seconds (20 allocations: 53.407 MiB, 17.06% gc time)

@time fdot(x);
>>0.002195 seconds (8 allocations: 7.630 MiB)

@time f.(x);
0.002047 seconds (8 allocations: 7.630 MiB)

可以看出，点乘操作要快很多。

向量化并不会提高Julia的运行速度

很多用过MATLAB和Python的同学都会觉得向量操作肯定要比循环操作要快很多，但在Julia中并没有这个规则，这一点要由为注意。

function loop_sum()
    s = 0.0
    for i = 1:1000
        s = 0.0
        for k = 1:10000
            s += 1.0/(k*k)
        end
    end
    s
end

先用@code_warntype loop_sum()看一下代码有没有类型转换问题。