为什么这个矩阵复制版本如此缓慢？

Question

为什么这个矩阵复制版本如此缓慢？

4

我注意到julia在矩阵复制时表现出了奇怪的行为。

考虑以下三个函数：

function priv_memcopyBtoA!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  A[1:n,1:n] = B[1:n,1:n]
  return nothing
end

function priv_memcopyBtoA2!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  ii = 1; jj = 1;
  while ii <= n
    jj = 1 #(*)
    while jj <= n
      A[jj,ii] = B[jj,ii]
      jj += 1
    end
    ii += 1
  end
  return nothing
end

function priv_memcopyBtoA3!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  A[1:n,1:n] = view(B, 1:n, 1:n)
  return nothing
end

编辑：1）我测试了代码是否会抛出BoundsError，所以标有jj = 1 #(*)的那行在最初的代码中缺失。测试结果已经来自于修复后的版本，因此它们保持不变。2）我添加了视图变体，感谢@Colin T Bowers解决了两个问题。

看起来这两个函数应该产生大致相同的代码。然而，我得到的结果是：

A = fill!(Matrix{Int32}(2^12,2^12),2); B = Int32.(eye(2^12));

结果

@timev priv_memcopyBtoA!(A,B, 2000)
  0.178327 seconds (10 allocations: 15.259 MiB, 85.52% gc time)
elapsed time (ns): 178326537
gc time (ns):      152511699
bytes allocated:   16000304
pool allocs:       9
malloc() calls:    1
GC pauses:         1

并且

@timev priv_memcopyBtoA2!(A,B, 2000)
 0.015760 seconds (4 allocations: 160 bytes)
elapsed time (ns): 15759742
bytes allocated:   160
pool allocs:       4

并且

@timev priv_memcopyBtoA3!(A,B, 2000)
  0.043771 seconds (7 allocations: 224 bytes)
elapsed time (ns): 43770978
bytes allocated:   224
pool allocs:       7

那是一个非常大的差异。这也很令人惊讶。我预计第一个版本会像memcopy一样难以超越，对于一个大的内存块来说，memcopy很难被超越。

第二个版本有指针算术（getindex），分支条件（<=）和每个赋值都有边界检查的开销。然而，每个赋值只需要约3ns时间。

此外，垃圾收集器消耗的时间对于第一个函数来说是不确定的。如果没有进行垃圾收集，则这种显著差异变得微不足道，但它仍然存在。版本3和版本2之间的差距仍然是2.5倍左右的因素。

那么为什么“memcopy”版本没有“赋值”版本高效呢？

- BlueLemon

2个回答

3

为什么A [1：n，1：n] = view（B，1：n，1：n）或其变体比一组while循环慢？让我们看看A [1：n，1：n] = view（B，1：n，1：n）的操作。 view返回一个包含指向父对象B的指针和计算应该复制的索引信息的迭代器。

A [1：n，1：n] = ... 被解析为调用_setindex！（...）。在此之后，在调用链下进行几次调用后，主要工作由以下代码完成：
.\abstractarray.jl:883;
# In general, we simply re-index the parent indices by the provided ones
function getindex(V::SlowSubArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds r = V.parent[reindex(V, V.indexes, I)...]
    r
end
#.\multidimensional.jl:212;
@inline function next(iter::CartesianRange{I}, state) where I<:CartesianIndex
    state, I(inc(state.I, iter.start.I, iter.stop.I))
end
@inline inc(::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
@inline inc(state::Tuple{Int}, start::Tuple{Int}, stop::Tuple{Int}) = (state[1]+1,)
@inline function inc(state, start, stop)
    if state[1] < stop[1]
        return (state[1]+1,tail(state)...)
    end
    newtail = inc(tail(state), tail(start), tail(stop))
    (start[1], newtail...)
end


getindex接受一个视图V和一个索引I。我们从B获取视图，从A获取索引I。在每一步中，reindex从视图V和索引I计算出用于获取B中元素的索引。它被称为r，并将其返回。最后，r写入A。

每次复制后，inc会将索引I增加到A中的下一个元素，并测试是否完成。请注意，代码来自v0.63，但在master中基本相同。

原则上，该代码可以简化为一组while循环，但它更通用。它适用于B的任意视图和形式为a：b：c的任意切片以及任意数量的矩阵维度。在我们的情况下，大的N是2。

由于函数更加复杂，编译器无法优化它们。也就是说，有一个建议，即编译器应该将它们内联，但它不这样做。这表明所示的函数是非平凡的。

对于一组循环，编译器将最内层循环减少为三个加法（每个指向A和B的指针以及一个用于循环索引的加法）和一条复制指令。

tl;dr A [1：n，1：n] = view（B，1：n，1：n）的内部调用链与多重分派相关联且非平凡，并处理一般情况。这会产生开销。一组while循环已经优化为特殊情况。

请注意，性能取决于编译器。如果我们看一维情况A [1：n] = view（B，1：n），它比while循环更快，因为它矢量化了代码。但是对于高维度N > 2，差异会增大。

- BlueLemon

1

好的！你做得很好。你观察到的关键差异是我们所称的“IndexLinear”和“IndexCartesian”。在这样的代码中，我们需要模拟任意数量的嵌套循环，只使用一个循环。查看SubArray开发者文档获取更多详细信息-希望以后能看到你更多的作品！ - mbauman

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Colin T Bowers · Accepted Answer

首先，您的代码存在一个错误。运行以下代码：

A = [1 2 ; 3 4]
B = [5 6 ; 7 8]
priv_memcopyBtoA2!(A, B, 2)

然后：

julia> A
2×2 Array{Int64,2}:
 5  2
 7  4

在每个内部的 while 循环结束时，你需要将 jj 重新赋值为 1，即：

function priv_memcopyBtoA2!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  ii = 1 
  while ii <= n
    jj = 1
    while jj <= n
      A[jj,ii] = B[jj,ii]
      jj += 1
    end
    ii += 1
  end
  return nothing
end

即使修复了这个bug，你仍然会注意到使用while循环的解决方案更快。这是因为Julia中的数组切片会创建临时数组。所以在这行代码中：

A[1:n,1:n] = B[1:n,1:n]

右侧的操作会创建一个临时的nxn数组，然后将该临时数组分配给左侧。如果想避免临时数组的分配，可以改写为：

A[1:n,1:n] = view(B, 1:n, 1:n)

你会注意到两种方法的时间现在非常接近，虽然while循环仍然略快。总的来说，在Julia中循环速度很快（就像C一样快），显式编写循环通常会得到最优化的编译代码。我仍然期望显式循环比view方法更快。

至于垃圾回收方面的问题，这只是你计时方法的结果。最好使用BenchmarkTools包中的@btime，该包使用各种技巧避免陷入垃圾回收等陷阱。