Julia中大规模科学计算数据的最优数据结构选择

背景 #

从Python转用Julia，看中其性能优势与灵活性，目前正在将重要项目迁移至Julia，但处理大规模数据（数千行/列甚至百万级行）时，始终无法充分发挥Julia的最佳性能。

遇到的具体问题 #

1. StaticArrays.jl的局限性 #

• SMatrix无法支持大尺寸结构：尝试创建5000行1000列的SMatrix时直接崩溃：

julia> rnd = SMatrix{5000, 1000, Float64}(rand(5000, 1000));

[process exited with code 3221225725 (0xc00000fd)]

• SVector向量的运算效率低下：虽然创建SVector的向量速度尚可：

julia> @benchmark [SVector{1000, Float64}(rand(1000)) for _∈1:5000]
BenchmarkTools.Trial: 100 samples with 1 evaluation.
 Range (min … max):  40.097 ms … 76.436 ms  ┊ GC (min … max):  0.00% … 30.38%
 Time  (median):     49.359 ms              ┊ GC (median):    19.95%
 Time  (mean ± σ):   50.296 ms ±  6.551 ms  ┊ GC (mean ± σ):  17.96% ±  8.91%

  ▃            ▃██▆▆▄▁▃  ▃▁
  █▄▆▁▆▆▄▁▄▆▄▇▇████████▆▄██▆▁▄▁▆▄▁▁▄▄▁▁▆▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▄
  40.1 ms         Histogram: frequency by time        72.5 ms <

 Memory estimate: 76.90 MiB, allocs estimate: 5002.

但基于该结构的运算速度极慢，比如实现逐元素相乘的函数：

julia> rnd_sv1 = [SVector{1000, Float64}(rand(1000)) for _∈1:5000];

julia> rnd_sv2 = [SVector{1000, Float64}(rand(1000)) for _∈1:5000];

julia> function dot(sv1, sv2)
           result = Vector{SVector{1000, Float64}}(undef, 5000)
           for idx∈eachindex(sv1)
               result[idx]=SVector{1000, Float64}([sv1[idx][i]*sv2[idx][i] for i∈1:1000])
           end
       end;

julia> @benchmark dot($rnd_sv1, $rnd_sv2)
BenchmarkTools.Trial: 2 samples with 1 evaluation.
 Range (min … max):  3.211 s …   3.232 s  ┊ GC (min … max): 0.27% … 0.00%
 Time  (median):     3.221 s              ┊ GC (median):    0.13%
 Time  (mean ± σ):   3.221 s ± 14.793 ms  ┊ GC (mean ± σ):  0.13% ± 0.19%

  █                                                       █
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  3.21 s         Histogram: frequency by time        3.23 s <

 Memory estimate: 76.90 MiB, allocs estimate: 5002.

2. 常规Array的性能优势 #

常规矩阵的运算速度明显更快，比如矩阵乘法操作：

julia> rnd_mat1 = rand(5000, 1000);

julia> rnd_mat2 = rand(1000, 5000);

julia> @benchmark $rnd_mat1 * $rnd_mat2
BenchmarkTools.Trial: 7 samples with 1 evaluation.
 Range (min … max):  706.033 ms … 857.549 ms  ┊ GC (min … max): 0.06% … 12.98%
 Time  (median):     770.682 ms               ┊ GC (median):    5.71%
 Time  (mean ± σ):   775.272 ms ±  61.905 ms  ┊ GC (mean ± σ):  6.11% ±  6.41%

  █    █ █                  █             █              █    █
  █▁▁▁▁█▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁█ ▁
  706 ms           Histogram: frequency by time          858 ms <

 Memory estimate: 190.73 MiB, allocs estimate: 2.

核心疑问 #

由于无法创建大尺寸SMatrix，且SVector向量无法利用矩阵乘法的性能优势，特此询问：

• 处理百万级行甚至更高维度的大规模科学计算数据时，Julia的最优数据结构是什么？
• 如何写出性能优于Python的Julia代码？

内容的提问来源于stack exchange，提问作者Shayan