Julia 笔记系列:


唠唠闲话

本篇介绍 Julia 广播,向量化与代码性能,对应课程第四讲,主要内容如下:

课程讲义为 Pluto 文件,Pluto 的安装及介绍参看Julia 初学者指南(一) | 安装、配置及编译器


广播

广播与 for 循环:

  1. Matlab 的广播本质是执行 for 循环,用起来快是因为使用 C 语言代码加速。
  2. 在 MATLAB/Python 这类语言中,for 循环执行效率低,广播运算允许我们将 for 循环的执行从一门缓慢的语言转移到一门高效的语言上执行。
  3. 对于 Julia 来说, 因为存在高效的 for 循环, 所以不用广播代码也很快, 使用广播只是可以让代码变得更简洁易读。

语法

  1. XY 尺寸一致时, f.(X, Y) 等价于 map(f, X, Y)

    1
    2
    g(x, y) = x + y
    @show g.([1, 2], [3, 4])

    深度截图_选择区域_20211104222059

  2. XY 的维数一致,且尺寸不同的维度有一个的长度为 1,此时长度为 1 的那个维度会被复制到尺寸一致。

    1
    2
    3
    4
    5
    @show x1 = reshape(1:6,2,3)
    @show y1 = reshape(1:2,2,1)
    @show g.(x1,y1)
    @show y2 = repeat(y1, 1, 3)
    @show g.(x1,y2)

    深度截图_选择区域_20211104222906

  3. 高维及多变量的规则类似,维数相等时,广播要求各维度上,尺寸要么相同,要么取 1,这一来展开模式不会有歧义。

    1
    2
    3
    4
    5
    g(x,y,z) = z + y + z
    @show x1 = reshape(1:6,2,3)
    @show y1 = reshape(1:2,2,1)
    @show z1 = reshape(1:3,1,3)
    g.(x1,y1,z1)

    深度截图_选择区域_20211104223411

  4. XY 的维数不一致时,维数小的矩阵补尺寸为 1 的维度后,再回到前边规则。

    1
    2
    3
    x4 = reshape(1:6,3,2) # 3x2
    y4 = [7, 8, 9] # 3 ---> 3x1
    @show g.(x4, y4) # 3x2

    深度截图_选择区域_20211104223748

注:计算矩阵元素之和时,用 for 循环逐列操作可以调用 @simd 并行加速;如果逐行操作使用 @simd 只没有加速效果。

调整尺寸

这两个函数常用于修改数据尺寸:

  • reshape: 在不改变内存顺序的情况下调整尺寸
  • permutedims: 交换维度 (同时会改变内存顺序)
  1. Julia 中的矩阵按列存储,查看方式

    1
    2
    3
    x = [1 2 3
    4 5 6]
    @show x[:]

    深度截图_选择区域_20211104224037

  2. 换言之,二维数组存储顺序为

    a1,1a2,1an,1a1,2a2,2an,2\begin{align*} &a_{1,1}\rightarrow a_{2,1} \rightarrow \cdots \rightarrow a_{n,1}\\ &a_{1,2}\rightarrow a_{2,2} \rightarrow \cdots \rightarrow a_{n,2}\\ &\qquad\cdots\rightarrow\cdots \end{align*}

  3. reshape 不改变内存顺序
    深度截图_选择区域_20211104224319

  4. reshape 查看高维数据的存储顺序,靠左维度先增大。
    深度截图_选择区域_20211104224828

  5. permutedims 交换维度,改变内存

    1
    2
    3
    x = [1 2 3
    4 5 6]
    permutedims(x,(2,1))

    深度截图_选择区域_20211104225600


代码性能

向量化编程

向量化编程就是尽可能避免显式 for 循环的代码,背后逻辑是为了尽可能将 for 循环从低效的 Python/MATLAB 端转移到高效的 C/Fortran 端, 从而尽可能少地触发这些动态语言的性能瓶颈。

  1. 我们将采用与广播类似规则的函数称为向量化函数,比如

    1
    2
    f(x::Real, y::Real) = x * y # 标量函数
    f(X::AbstractArray, Y::AbstractArray) = X .* Y # 向量化函数
  2. 向量化代码的问题

    • 向量化代码的实现需要底层 C/Fortran 代码进行支撑。如果你所关心的问题恰好没有人在 C/Fortran 下给出高效实现的话, 那么就需要你自己来做了。
    • 向量化代码的中间结果是数组而非标量, 会带来一定的额外内存开销。
    • 向量化代码可能会阻碍一些本可以进行的性能优化,或带来了一些不必要的代码。
    • 相比于标量代码来说,向量化代码既不容易阅读,也不容易写对。

Fused Dot

  1. A .* B .+ C 这种运算非常普遍,它背后有两种可能的实现方式

    1
    2
    3
    4
    5
    6
    ## 方法一,分两次运行
    tmp = A .* B
    tmp .+ C
    ## 方法二,使用一次广播
    f(a, b, c) = a * b + c
    f.(A, B, C)
  2. 第二种方式可以避免中间矩阵 tmp 的不必要内存开销,因此 Julia 提供了一个内置的 fused dots 机制:当整个运算都是点运算时,Julia 会试图构造类似 f 的标量形式的函数,然后对函数整体进行广播。同时,Julia 也提供了一个 @. 宏用来辅助代码书写。

  3. 举个例子,方法 1 运行时间长,方法 2,3 等价

    1
    2
    3
    4
    5
    6
    7
    8
    A,B,C = rand(100,100),rand(100,100),rand(100,100)
    muladd_v1(A, B, C) = A .* B + C
    # 以下两种写法是完全等价的
    muladd_v2(A, B, C) = A .* B .+ C
    muladd_v3(A, B, C) = @. A * B + C
    @btime muladd_v1($A, $B, $C)
    @btime muladd_v2($A, $B, $C)
    @btime muladd_v3($A, $B, $C);

    深度截图_选择区域_20211104231536

可以看到,第一种方法分配了更多内存,且计算用时更长。

延伸阅读
Dot Syntax for Vectorizing Functions
More Dots: Syntactic Loop Fusion in Julia

view

不必要的内存开销一直都是性能优化致力于解决的问题。

  1. 使用切片获取数据时,创建了新对象

    1
    2
    3
    4
    x = [1,2,3,4]
    y = x[1:2]
    y[1] = 0
    @show x y

    20211104232221

  2. 使用宏 @view 可以引用原对象

    1
    2
    3
    4
    x = [1,2,3,4]
    y = @view x[1:2]
    y[1] = 0
    @show x y

    深度截图_选择区域_20211104232308

  3. 当我们不需要修改数据时,使用 @veiw 可以减少内存开销,继而加快运行速度。

Ps:Python 中切片也是创建新对象,不过没有类似 Julia 的处理工具。