13 如何写出高性能的Julia

猫叔Rex

发布于 2020-06-30 11:55:21

1.4K0

发布于 2020-06-30 11:55:21

文章被收录于专栏：科学计算

高性能代码

避免全局变量

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

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

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

const VAL = 0

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

global x = rand(10000)

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

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

在函数第一次运行时，由于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

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

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

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

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

code generation

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

@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)

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

再用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)

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

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

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

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

抽象和具体类型

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

但如果函数参数的类型是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

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

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

对于concrete类型：

image

对于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()看一下代码有没有类型转换问题。

运算过程中不存在类型转换的warning。

@benchmark loop_sum()
>>BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     8.868 ms (0.00% GC)
  median time:      8.869 ms (0.00% GC)
  mean time:        8.870 ms (0.00% GC)
  maximum time:     8.925 ms (0.00% GC)
  --------------
  samples:          564
  evals/sample:     1

用向量来完成上面的累加操作

function vec_sum()
    s = 0.0
    k = collect(1:10000)
    for i=1:1000
        s = sum(1 ./(k.^2))
    end
end

@benchmark vec_sum()
>>BenchmarkTools.Trial: 
  memory estimate:  76.48 MiB
  allocs estimate:  3002
  --------------
  minimum time:     19.575 ms (17.75% GC)
  median time:      22.603 ms (21.13% GC)
  mean time:        23.038 ms (22.79% GC)
  maximum time:     59.079 ms (68.95% GC)
  --------------
  samples:          217
  evals/sample:     1

@simd 可以在运算支持被重新recorded时加速运算，加法操作就是可以被加速的

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

@benchmark loop_sum_simd()
>>BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     6.202 ms (0.00% GC)
  median time:      6.236 ms (0.00% GC)
  mean time:        6.240 ms (0.00% GC)
  maximum time:     7.119 ms (0.00% GC)
  --------------
  samples:          802
  evals/sample:     1

说明Julia中向量运算并不会优化速度，这一点在Julia官网也多次说明。

输出预分配

怎么理解呢？就是当我们在操作Array或者其他复杂类型时，我们预先为结果分配一个存储它的Array或其他类型，再进行计算，性能会显著提升

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
   return y
end

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

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

@time loopinc()
>>0.428358 seconds (40.00 M allocations: 1.490 GiB, 10.25% gc time)

@time loopinc_prealloc()
>>0.024745 seconds (6 allocations: 288 bytes)

copy和view

copy就是被复制的变量是原始变量的一份拷贝，存在内存拷贝的操作，而view只是映射关系，不存在内存拷贝

先举个view的例子

A = zeros(3, 3);
@views for row in 1:3
   b = A[row, :]
   b[:] .= row  # A的值也会相应改变
end

分别定义copy和view的函数

fcopy(x) = sum(x[2:end-1]);
@views fview(x) = sum(x[2:end-1]);
x = rand(10^6);

@time fcopy(x)
>>0.031803 seconds (7 allocations: 7.630 MiB, 85.25% gc time)

@time fview(x)
>>0.001218 seconds (6 allocations: 224 bytes)

IO操作优化

在写文件时，对于某个变量，我们可以通过$的方式获取它的值，也可以直接使用变量本身

println(file, "$a $b")

下面的写法更好一些，因为上面这种方式先把表达式转成字符串，再写入文件中。而下面这种方式直接把值写入到文件中

println(file, a, " ", b)

其他性能优化小技巧

避免不必要的Array，比如计算x,y,z的和时，使用x+y+z，不要用sum([x,y,z])
计算复数模值平方时，使用abs2(z)，不要使用abs(z)^2，简而言之，有现成的函数就不要自己写计算过程
用div(x,y)代替trunc(x/y)，用fld(x,y)代替floor(x/y)，用cld(x,y)代替ceil(x/y)，原因也是有现成的函数就不要自己写计算过程

function test_loop(r)
    for i in 1:10000
        for j in length(r)
            r[j] = tanh(r[j])
        end
    end
end

function test_vectorize(r)
    for i in 1:10000
        @. r = tanh(r)
    end
end

r = rand(1000)

@time test_vectorize(r)
>> 0.126987 seconds (4 allocations: 160 bytes)

@time test_loop(r)
>>0.000289 seconds (4 allocations: 160 bytes)

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2019-04-07，如有侵权请联系 cloudcommunity@tencent.com 删除

ide

编程算法

本文分享自傅里叶的猫微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

ide

编程算法

登录后参与评论

0 条评论

热度