时间戳

首页 > 解决方案 > 优化中的 Julia 抽象类型？

arrays - 优化中的 Julia 抽象类型？

问题描述

作为 Julia 的起点，我决定实现一个简单的 Strassen 产品：

@inbounds function strassen_product(A :: Array{Num, 2}, B :: Array{Num, 2}, k = 2) :: Array{Num, 2} where {Num <: Number}

    A_height, A_width = size(A)
    B_height, B_width = size(B)

    @assert A_height == A_width == B_height == B_width  "Matrices are noth both square or of equal size."
    @assert isinteger(log2(A_height))                   "Size of matrices is not a power of 2."

    if A_height ≤ k
        return A * B
    end

    middle = A_height ÷ 2

    A₁₁, A₁₂ =  A[1:middle, 1:middle],     A[1:middle, middle+1:end]
    A₂₁, A₂₂ =  A[middle+1:end, 1:middle], A[middle+1:end, middle+1:end]

    B₁₁, B₁₂ = B[1:middle, 1:middle],     B[1:middle, middle+1:end]
    B₂₁, B₂₂ = B[middle+1:end, 1:middle], B[middle+1:end, middle+1:end]

    P₁ = strassen_product(A₁₁ + A₂₂, B₁₁ + B₂₂)
    P₂ = strassen_product(A₂₁ + A₂₂, B₁₁    )
    P₃ = strassen_product(A₁₁,       B₁₂ - B₂₂)
    P₄ = strassen_product(A₂₂,       B₂₁ - B₁₁)
    P₅ = strassen_product(A₁₁ + A₁₂, B₂₂      )
    P₆ = strassen_product(A₂₁ - A₁₁, B₁₁ + B₁₂)
    P₇ = strassen_product(A₁₂ - A₂₂, B₂₁ + B₂₂)

    C₁₁ = P₁ + P₄ - P₅ + P₇
    C₁₂ = P₃ + P₅
    C₂₁ = P₂ + P₄
    C₂₂ = P₁ + P₃ - P₂ + P₆

    return [ C₁₁ C₁₂ ;
             C₂₁ C₂₂ ]

end

一切都很好。实际上，我喜欢这样的不安全优化的整个想法@inbounds，实际上确实会在很大程度上影响性能，最多不会影响几毫秒。

现在，为了进一步优化，因为我没有 for 循环，将使用视图来处理这些A₁₁等矩阵，这样就不会发生复制。

所以我@views在包含索引的4行前面拍了拍。当然我得到了一个错误，因为几行之后递归调用不需要Array{...}参数SubArray{...}。所以我将参数的类型和返回类型更改为AbstractArray{Num, 2}. 这次它起作用了，因为它AbstractArray是数组类型的基本类型，但是......性能急剧下降，实际上慢了 10 倍，分配更多。

我的测试用例是这样的：

A = rand(1:10, 4, 4)
B = rand(1:10, 4, 4)

@time C = strassen_product(A, B)

使用@views+时AbstractArray：

0.457157 seconds (1.96 M allocations: 98.910 MiB, 5.56% gc time)

使用无视图版本时：

0.049756 seconds (126.92 k allocations: 5.603 MiB)

差异是巨大的，应该更快的版本慢了 9-10 倍，分配量增加了大约 15 倍，空间几乎是其他版本的 20 倍。

编辑：这不是这两种情况的第一次运行，而是~10 次测试运行的最“中值”值。不是第一次运行，当然也不是最小或最大峰值。

编辑：我使用的是 1.0 版。

我的问题是：为什么会发生这种情况？我没有得到什么？我的推理是使用视图而不是副本会更快......错了吗？

标签： arraysoptimizationjulia

解决方案

是的，查看版本比复制版本花费更长的时间，但分配的内存更少。这就是为什么。

使用视图而不是副本并不一定意味着加速（尽管它减少了内存分配。）加速取决于程序中的许多事情。首先，您创建的每个视图都是针对矩阵的一个图块。Julia 中的矩阵以列为主存储在内存中，这使得从内存中检索两列比检索两行更容易，因为列的元素是连续存储的。

矩阵的瓦片不是连续存储在内存中的。创建图块的副本会访问矩阵中的每个必要元素并将它们写入内存中的连续块，而在图块上创建视图仅存储限制。尽管创建副本比在图块上创建视图需要更多的时间和更多的内存访问，但副本连续存储在内存中，这意味着更容易访问、更容易矢量化和更容易缓存供 CPU以供以后使用。

您正在多次访问您创建的切片，并且每次新访问都在完全递归调用之后进行，每个调用都需要一些时间和内存访问，因此您已经加载到缓存中的切片可能会丢失。因此，对同一个图块的每次新访问都可能需要从内存中完全重新加载。但是，视图图块的完全重新加载比复制图块的完全重新加载需要更多时间。这就是为什么view版本比复制版本花费的时间更长，尽管view版本分配的内存更少并且使用的访问次数更少。

看看文档中的性能提示，都考虑使用切片的视图，复制数据并不总是坏的。

数组连续存储在内存中，有助于 CPU 矢量化，并且由于缓存而减少了内存访问。这些与建议以列优先顺序访问数组的原因相同（见上文）。由于非顺序的内存访问，不规则的访问模式和不连续的视图会大大减慢数组的计算速度。

在对不规则访问的数据进行操作之前将其复制到连续数组中可能会导致很大的加速，例如下面的示例。在这里，一个矩阵和一个向量在相乘之前以 800,000 个随机打乱的索引被访问。将视图复制到普通数组中会加快乘法的速度，即使复制操作的成本也是如此……

但是，这就是说，差异不应像您的结果所暗示的那么大。此外，4x4 矩阵的乘积甚至不应该花费一毫秒。我相信在您的函数的每次调用中，您还会重新加载您的函数定义，这会使 JIT 编译的版本过时，而 Julia 必须一次又一次地编译您的函数。您也可能正在创建新函数类型不稳定的情况。我建议您使用@btimeinBenchmarkTools来衡量分配和时间。

您还应该使用更大的矩阵进行测试。在 4x4 数组上创建 2x2 视图仍会将数据写入内存，其大小与 2x2 副本相当。小数据的时序也容易产生噪音。

@inbounds function strassen_product(A, B, k = 2)

     A_height, A_width = size(A)
#     B_height, B_width = size(B)

#     @assert A_height == A_width == B_height == B_width  "Matrices are noth both square or of equal size."
#     @assert isinteger(log2(A_height))                   "Size of matrices is not a power of 2."

     if A_height ≤ k
         return A * B
     end

    middle = A_height ÷ 2

    A₁₁, A₁₂ =  @view(A[1:middle, 1:middle]), @view(A[1:middle, middle+1:end])
    A₂₁, A₂₂ =  @view(A[middle+1:end, 1:middle]), @view(A[middle+1:end, middle+1:end])
    B₁₁, B₁₂ = @view(B[1:middle, 1:middle]),     @view(B[1:middle, middle+1:end])
    B₂₁, B₂₂ = @view(B[middle+1:end, 1:middle]), @view(B[middle+1:end, middle+1:end])

    P₁ = strassen_product(A₁₁ + A₂₂, B₁₁ + B₂₂)
    P₂ = strassen_product(A₂₁ + A₂₂, B₁₁    )
    P₃ = strassen_product(A₁₁,       B₁₂ - B₂₂)
    P₄ = strassen_product(A₂₂,       B₂₁ - B₁₁)
    P₅ = strassen_product(A₁₁ + A₁₂, B₂₂      )
    P₆ = strassen_product(A₂₁ - A₁₁, B₁₁ + B₁₂)
    P₇ = strassen_product(A₁₂ - A₂₂, B₂₁ + B₂₂)

    C₁₁ = P₁ + P₄ - P₅ + P₇
    C₁₂ = P₃ + P₅
    C₂₁ = P₂ + P₄
    C₂₂ = P₁ + P₃ - P₂ + P₆

    return [ C₁₁ C₁₂ ;
             C₂₁ C₂₂ ]

end
@inbounds function strassen_product2(A, B, k = 2)

    A_height, A_width = size(A)
    #B_height, B_width = size(B)

    #@assert A_height == A_width == B_height == B_width  "Matrices are noth both square or of equal size."
    #@assert isinteger(log2(A_height))                   "Size of matrices is not a power of 2."

    if A_height ≤ k
        return A * B
    end

    middle = A_height ÷ 2

    A₁₁, A₁₂ =  A[1:middle, 1:middle], A[1:middle, middle+1:end]
    A₂₁, A₂₂ =  A[middle+1:end, 1:middle], A[middle+1:end, middle+1:end]

    B₁₁, B₁₂ = B[1:middle, 1:middle],     B[1:middle, middle+1:end]
    B₂₁, B₂₂ = B[middle+1:end, 1:middle], B[middle+1:end, middle+1:end]

    P₁ = strassen_product2(A₁₁ + A₂₂, B₁₁ + B₂₂)
    P₂ = strassen_product2(A₂₁ + A₂₂, B₁₁    )
    P₃ = strassen_product2(A₁₁,       B₁₂ - B₂₂)
    P₄ = strassen_product2(A₂₂,       B₂₁ - B₁₁)
    P₅ = strassen_product2(A₁₁ + A₁₂, B₂₂      )
    P₆ = strassen_product2(A₂₁ - A₁₁, B₁₁ + B₁₂)
    P₇ = strassen_product2(A₁₂ - A₂₂, B₂₁ + B₂₂)

    C₁₁ = P₁ + P₄ - P₅ + P₇
    C₁₂ = P₃ + P₅
    C₂₁ = P₂ + P₄
    C₂₂ = P₁ + P₃ - P₂ + P₆

    return [ C₁₁ C₁₂ ;
             C₂₁ C₂₂ ]

end

这是@btimein的测试Benchmarktools

A = rand(1:10, 256, 256)
B = rand(1:10, 256, 256)

@btime C = strassen_product(A, B); # view version
@btime D = strassen_product2(A, B); #copy version

结果：

438.294 ms (4941454 allocations: 551.53 MiB)
349.894 ms (4529747 allocations: 620.04 MiB)

推荐阅读