优化中的Julia抽象类型？

Question

作为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₁₁等矩阵使用视图，这样就不会发生复制。

所以我在包含索引的4行前面给了@views一个耳光。当然，我得到了一个错误，因为几行之后，递归调用需要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版。

我的问题是:为什么会发生这种情况？我没有得到什么？我的理由是用视图代替拷贝会更快.不对?

Answer 1

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

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

矩阵的块不连续地存储在内存中。创建平铺的副本将访问矩阵中的每个必要元素，并将它们写入内存中的连续块，而在平铺上创建视图只存储限制。虽然创建复制比在块上创建视图花费更多的时间和更多的内存访问，但是副本被连续地存储在内存中，这意味着更容易的访问、更容易的矢量化和更容易的CPU缓存以供以后使用。

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

查看文档中的性能技巧，包括考虑对片使用视图和复制数据并不总是不好的。。

数组被连续地存储在内存中，将自己借给CPU矢量化，而由于缓存而减少了内存访问。这些原因与建议按列-主要顺序访问数组的原因相同(见上文)。不规则的访问模式和非连续的视图可以极大地减缓数组上的计算速度，因为这是因为非顺序的内存访问。 在对不定期访问的数据进行操作之前，将其复制到连续数组中可能会导致很大的加速，如下面的示例所示。在这里，一个矩阵和一个向量是在他们随机调整的指数中的80万，然后再乘以。将视图复制到普通数组中会加快乘法速度，即使要支付复制操作的成本.

然而，这一差异不应该像你的结果所显示的那么大。此外，4x4矩阵的乘积甚至不应该占用毫秒。我相信，在函数的每个调用中，您还会重新加载函数定义，这使得JIT编译版本过时了，而Julia不得不一次又一次地编译您的函数。您还可能正在创建一个新函数不是类型稳定的情况。我建议您在@btime中使用BenchmarkTools来度量分配和时间。

您还应该使用更大的矩阵进行测试。在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

下面是@btime在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)