优化中的 Julia 抽象类型？
Posted
技术标签:
【中文标题】优化中的 Julia 抽象类型？【英文标题】：Julia abstract types in optimization? 【发布时间】：2018-09-13 21:25:09 【问题描述】：作为 Julia 的起点，我决定实现一个简单的 Strassen 产品：
@inbounds function strassen_product(A :: ArrayNum, 2, B :: ArrayNum, 2, k = 2) :: ArrayNum, 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...。所以我将参数的类型和返回类型更改为AbstractArrayNum, 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 版。
我的问题是：为什么会这样？我没有得到什么？我的推理是使用视图而不是副本会更快......错了吗？
【问题讨论】：
看起来您的时间是从第一次运行开始的：在这种情况下，主要成本是 JIT 编译。运行函数第一次运行后尝试计时。 您是否只运行一次@time？第一次 JIT 编译时，请注意后续调用速度更快。我建议使用 BenchmarkTools 包中的 @btime 代替这些微基准测试。
计时后续运行，我仍然得到大致相同的计时。请注意，您仍在为输出分配结果：尝试就地完成整个事情。函数参数中的类型说明符也不会影响性能。有关详细信息，请参阅docs.julialang.org/en/v1/manual/performance-tips。
@SimonByrne 他们不是第一次运行，5-10 次运行的“中值”值越大
您使用的是哪个 Julia 版本？对于 Julia 1.0 上的两个版本，我得到的时间几乎相同（@views 版本略有优势）。
【参考方案1】：
是的，查看版本比复制版本花费的时间更长，但分配的内存更少。这就是原因。
使用视图而不是副本并不一定意味着加速（尽管它减少了内存分配）。加速取决于程序中的许多因素。首先，您创建的每个视图都是针对矩阵的一个图块。 Julia 中的矩阵以列为主存储在内存中，这使得 CPU 从内存中检索两列比检索两行更容易，因为列的元素是连续存储的。
矩阵的瓦片不是连续存储在内存中的。创建图块的副本会访问矩阵中的每个必要元素并将它们写入内存中的连续块，而在图块上创建视图仅存储限制。尽管创建副本比在图块上创建视图需要更多的时间和更多的内存访问，但副本会连续存储在内存中，这意味着更容易访问、更容易矢量化和 更容易缓存 CPU 供以后使用使用。
您正在多次访问您创建的切片，并且每次新的访问都在一个完整的递归调用之后进行，每个调用都需要一些时间和内存访问，因此您已经加载到缓存中的切片可能会丢失。因此，对同一个图块的每次新访问都可能需要从内存中完全重新加载。但是完全重新加载 view tile 比完全重新加载 copy tile 需要更多时间。 这就是为什么view 版本比复制版本花费更长的时间，尽管view 版本分配的内存更少并且使用的访问次数更少。
查看文档中的性能提示，Consider using views for slices、Copying data is not always bad。
数组连续存储在内存中，供 CPU 使用 由于缓存，矢量化和更少的内存访问。这些是 与建议以列为主访问数组的原因相同 订单（见上文）。不规则的访问模式和不连续的视图 可以大大减慢数组的计算速度，因为 非顺序内存访问。
将不定期访问的数据复制到一个连续的数组中 对其进行操作可能会导致很大的加速，例如在示例中 以下。这里，一个矩阵和一个向量在 800,000 处被访问 它们在相乘之前随机打乱的索引。复制 即使使用普通数组的视图也可以加快乘法速度 复制操作的成本...
不过，尽管如此，差异不应像您的结果所显示的那么大。 此外，4x4 矩阵的乘积甚至不应该花费一毫秒。我相信在你的函数的每次调用中，你也会重新加载你的函数定义，这使得 JIT 编译的版本过时了，而 Julia 已经一次又一次地编译你的函数。您也可能正在创建新函数类型不稳定的情况。我建议你在BenchmarkTools 中使用@btime 来衡量分配和时间。
您还应该使用更大的矩阵进行测试。在 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)
【讨论】：
以上是关于优化中的 Julia 抽象类型？的主要内容，如果未能解决你的问题，请参考以下文章