【问题标题】:Julia Array of arrays: (rows -> columns) performanceJulia 数组数组:(行 -> 列)性能
【发布时间】:2020-05-25 05:38:54
【问题描述】:

在这里完成 Julia 新手。

给定一个数组数组,我想组合每个子数组的对应元素。像这样的:

 [2, 7, 9]       [2, 3, 2, 7, 3]
 [3, 5, 4]       [7, 5, 7, 9, 5]
 [2, 7, 7]  ->   [9, 4, 7, 1, 1]
 [7, 9, 1]
 [3, 5, 1]

搜索 stackoverflow 我遇到了几个解决方案,而不是直接循环或列表理解。

julia> a=Vector{Int}[rand(1:10,3) for i=1:5]
5-element Array{Array{Int64,1},1}:
 [2, 7, 9]
 [3, 5, 4]
 [2, 7, 7]
 [7, 9, 1]
 [3, 5, 1]

julia> using BenchmarkTools

julia> @btime a2=mapslices( x -> [x], hcat(a...), dims=2)[:]
  6.174 μs (65 allocations: 3.45 KiB)
3-element Array{Array{Int64,1},1}:
 [2, 3, 2, 7, 3]
 [7, 5, 7, 9, 5]
 [9, 4, 7, 1, 1]

julia> @btime a3=[getindex.(a,i) for i=1:length(a[1])]
  948.087 ns (14 allocations: 768 bytes)
3-element Array{Array{Int64,1},1}:
 [2, 3, 2, 7, 3]
 [7, 5, 7, 9, 5]
 [9, 4, 7, 1, 1]

我的问题是:为什么第二个比第一个快六倍?和hcat有关系吗?

【问题讨论】:

  • 如果您想要第三个更好的,请使用:julia> @btime a3 = [[a[i] for a in $a] for i = 1:3] 193.593 ns (9 allocations: 640 bytes) 3-element Array{Array{Int64,1},1}: [2, 3, 2, 7, 3] [7, 5, 7, 9, 5] [9, 4, 7, 1, 1]。这一切都归结为分配,第一个例子(65 allocations: 3.45 KiB),第二个(14 allocations: 768 bytes),还有我的(9 allocations: 640 bytes)。在 Julia 代码中,最好避免任何不必要的分配以获得最大速度。
  • 不,这只是一个基准测试人工制品。您将a 插入到基准表达式中,而@toylas 没有。如果插值正确,它们具有相似的性能。不过,我更愿意建议:[getindex.(a,i) for i in eachindex(first(a))],在我的笔记本电脑上,它的速度提高了 10%,但更重要的是,它更强大(尽管你仍然应该这样做)。顺便说一句,我认为写[a[i] for a in a] 不是一个好主意,将a 重用作为迭代变量非常令人困惑,事实上我有点惊讶它甚至可以工作。
  • 注意原始时序问题,在非常量全局变量上进行基准测试是一个问题。您应该通过编写$a 而不是a 来插入它们,就像@aboammar 一样,或者将其设为常量。我已经看到这会导致误导性的结果,甚至可以颠倒正在考虑的不同替代方案的顺序,因为由此产生的类型推断中断占主导地位。我没有检查这是否是这样一种情况

标签: arrays julia


【解决方案1】:

正确的基线和基准测试

好的,首先让我们在我的电脑上建立基线。

在我们做任何其他事情之前,我们需要确保我们没有对全局变量进行基准测试。 来自BenchmarkTools readme

如果您要进行基准测试的表达式依赖于外部变量,您应该将$ to "interpolate" 用于avoid the problems of benchmarking with globals 的基准表达式。 本质上,任何插值变量 $x 或表达式 $(...) 在基准测试开始之前都是“预先计算的”......

julia> a=Vector{Int}[rand(1:10,3) for i=1:5];

julia> @btime a2=mapslices( x -> [x], hcat($a...), dims=2)[:];
  6.015 μs (65 allocations: 3.45 KiB)

julia> @btime a3=[getindex.($a,i) for i=1:length($a[1])];
  149.228 ns (6 allocations: 544 bytes)

(如果我没有插值,我会得到和你大致相同的a3999.500 ns (14 allocations: 768 bytes))。

所以a3 不是快 6 倍,而是实际上快了 33 倍。

为什么会有差异?

分配。

与其他操作(所有语言)相比,分配相当慢。 我们可以看到a2 代码比a3 代码分配的更多。

让我们看看分配的位:

a2

  • [x] 为每一列分配一个新的 1 元素数组
  • hcat 分配一个新矩阵,所有内容都连接在一起
  • mapslices 为从矩阵中取出的每个切片分配一个数组
  • mapslice 分配一个数组来保存输出(有趣的是它不做视图,但我检查了)
  • [:] 执行输出的重塑副本。 (替代方案是vec,它返回一个重塑视图)

a3

  • getindex.(a, i) 为输出的每一列分配一个数组(与mapslice 输入矩阵的内部切片相同)
  • [ ... for ...] 为输出分配一个数组(与 mapslices 输出相同)

    所以我们可以看到a2 中正在进行大量额外分配,而a3 中没有。

    如果我们只有与 hcat 相关的分配。

    既然最初的问题是因为hcat,让我们来看看。

我定义了一个新的基准保存到a4。 它使用eachslice 将视图的(惰性)生成器返回到矩阵的切片中。所以那里的分配可以忽略不计。 为了阻止它变得懒惰,我们collect它。 最终输出是SubArrays 中的Array(而不是Arrays 中的Array),但这很好,它仍然是AbstractArray 的子类型。

julia> @btime a4 = collect(eachslice(hcat($a...), dims=1));
  734.320 ns (13 allocations: 704 bytes)

我们的主要分配是 - hcat - collect 分配输出(与[ ... for ...] 相同)。

所以是的,hcat 有效果,但与大部分差异相差甚远。

喷溅和reduce(hcat, xs)

喷溅是一种代价。它通常很小,直到你开始喷出数百个项目,但因为这是一个微基准,其他一切都非常快,让我们看看它是如何删除它的。

Julia 有一个针对 reduce(hcat, xs) 的优化函数,因为 xs 是一个数组数组。

让我们看看情况如何:

julia> @btime a2_s=mapslices(x -> [x], reduce(hcat, $a), dims=2);
  5.278 μs (59 allocations: 3.17 KiB)

julia> @btime a4_s=collect(eachslice(reduce(hcat, $a), dims=1));
  337.656 ns (8 allocations: 528 bytes)

我们可以看到它有所作为。 但是在a2 的情况下,这并不多,因为hcat 只完成了一次,而x->xmapslices 中的缓慢分配从hcat 复制切片发生了很多次。

我们可以走得更快吗?

不是真的。 a3 是非常理想的代码。 它不分配任何不返回的内容。

想如果我们愿意换成使用StaticArrays 我们可以非常快地得到一些东西。

julia> b = @SVector [@SVector [rand(1:10) for ii in 1:3] for i=1:5];

julia> @btime b3=[getindex.($b,i) for i in 1:length($b[1])];
  36.055 ns (1 allocation: 208 bytes)

静态数组为编译器提供了更多信息。 特别是所有数组的大小,以及它们都不会被改变的承诺。 这意味着它可以: - 展开循环 - 编译时的边界检查 - 将它们分配在堆栈上(而不是堆上) - 可能我忘记了其他一些事情。

这让优化器(在 Julia 和 LLVM 中)变得非常疯狂。 他们被编译成每个输入列(/输出行)基本上 2 个 SSE/AVX 矢量化移动操作,加上少量的固定开销。

julia> @code_native (b->[getindex.(b,i) for i in 1:length(b[1])])(b)
    .section    __TEXT,__text,regular,pure_instructions
; ┌ @ REPL[83]:1 within `#161'
    subq    $136, %rsp
    vmovups (%rdi), %ymm0
    vmovups 32(%rdi), %ymm1
    vmovups 64(%rdi), %ymm2
    vmovups 88(%rdi), %ymm3
    vmovups %ymm3, 88(%rsp)
    vmovups %ymm2, 64(%rsp)
    vmovups %ymm1, 32(%rsp)
    vmovups %ymm0, (%rsp)
    movabsq $5152370032, %rax       ## imm = 0x1331AED70
; │┌ @ generator.jl:32 within `Generator' @ generator.jl:32
    vmovaps (%rax), %xmm0
    vmovups %xmm0, 120(%rsp)
; │└
    movabsq $collect, %rax
    movq    %rsp, %rdi
    vzeroupper
    callq   *%rax
    addq    $136, %rsp
    retq
    nop
; └

【讨论】:

  • 如果用@inbounds 写出循环似乎会快一点,但在 v1.4 上相差 10%。
  • 是的,我试过了,但不确定这是否有意义。此外,鉴于我们仅使用第一个的长度来确定内部长度,并且假设其他的长度相同,这似乎并不是特别安全的东西。它消除了早期代码中的错误的护栏。
  • 我在代码中使用了最短的数组长度,但无论如何性能优势很小。
  • 感谢您如此详细地解释。现在事情对我来说确实有意义。
猜你喜欢
  • 1970-01-01
  • 2016-08-12
  • 1970-01-01
  • 2021-11-12
  • 2018-04-11
  • 1970-01-01
  • 1970-01-01
  • 2021-04-06
  • 1970-01-01
相关资源
最近更新 更多