julia - 计算 julia 中每对矩阵行之间相关性的最有效方法是什么？

Question

假设我有一个名为mat的矩阵:

julia> mat = rand(1:10, 5, 3)
5×3 Matrix{Int64}:
 2  4   3
 5  3  10
 5  7   5
 9  5   7
 4  9   6

我想计算每对 mat 行之间的相关性(例如，cor(mat[1, :], mat[2, :] 等等上), 最后, 实现相关矩阵。我为它写了两个脚本, 我将提供基准测试。但是, 如果我能让它更快, 我会更高兴 (因为我应该在一个大的数据集，比如 2000x20 大小)。

第一种方法

一个非常简单的方法；首先，我用 zeros 创建了一个初始化矩阵，然后尝试在每对行上用计算出的相关性填充它。这不是一个好方法，因为我计算了两次必要的(因为，例如，cor([mat[1, :], mat[3, :]) 等于 cor( [垫[3, :], 垫[1, :])):

using Statistics

function calc_corr(matrix::Matrix)
    n::Int64 = size(matrix, 1)
    corr_mat = zeros(Float64, n, n)

    for (idx1, idx2)=Iterators.product(1:n, 1:n)
        @inbounds corr_mat[idx1, idx2] = cor(
            view(matrix, idx1, :),
            view(matrix, idx2, :)
        )
    end

    return corr_mat
end

第二种方法

计算上三角部分，然后创建对称矩阵，得到完整的相关矩阵:

using LinearAlgebra

function calc_corr2(matrix::Matrix)
    n::Int64 = size(matrix, 1)
    corr_mat = ones(Float64, n, n)

    # find upper triangular indices
    upper_triang_idx = findall(==(1), triu(ones(Int8, n, n), 1))

    for (idx1, idx2)=Tuple.(upper_triang_idx)
        @inbounds corr_mat[idx1, idx2] = cor(
            view(matrix, idx1, :),
            view(matrix, idx2, :)
        )
    end

    corr_mat = Symmetric(corr_mat)
    return corr_mat
end

基准测试

1. 首先在我之前声明的一个小矩阵上(mat):

using BenchmarkTools

@benchmark calc_corr($mat)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.950 μs …   6.210 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.160 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.178 μs ± 289.600 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇▇▂ ▁ ▅█▆▂▁▁▃▆▅▃▁ ▁▁▁▁                                      ▂
  ███████████████████████▆█▇▇█▇▇▆▇▆▇▆▆▄▆▄▄▅▃▁▄▄▄▃▅▃▄▄▄▁▁▁▃▄▁▄ █
  1.95 μs      Histogram: log(frequency) by time      3.62 μs <

 Memory estimate: 256 bytes, allocs estimate: 1.

# ---------------------------------------------------------------------

@benchmark calc_corr2($mat)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.220 μs … 773.080 μs  ┊ GC (min … max): 0.00% … 99.19%
 Time  (median):     1.420 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.698 μs ±   9.921 μs  ┊ GC (mean ± σ):  8.15% ±  1.40%

  █          
  █▇▃▂▃▄▂▂▄▄▃▂▃▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▁
  1.22 μs         Histogram: frequency by time        3.98 μs <

 Memory estimate: 976 bytes, allocs estimate: 7.

测试结果是否相同:

julia> calc_corr(mat) == calc_corr2(mat)
true

1. 在一个大矩阵上:

test_mat = rand(1:10, 2_000, 20);

@benchmark calc_corr($test_mat)
BenchmarkTools.Trial: 8 samples with 1 evaluation.
 Range (min … max):  632.258 ms … 680.094 ms  ┊ GC (min … max): 0.33% … 1.30%
 Time  (median):     646.215 ms               ┊ GC (median):    0.16%
 Time  (mean ± σ):   650.096 ms ±  16.089 ms  ┊ GC (mean ± σ):  0.49% ± 0.60%

  ▁ ▁           ▁  █       ▁                 ▁                ▁  
  █▁█▁▁▁▁▁▁▁▁▁▁▁█▁▁█▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  632 ms           Histogram: frequency by time          680 ms <

 Memory estimate: 30.52 MiB, allocs estimate: 2.

# ---------------------------------------------------------------------

@benchmark calc_corr2($test_mat)
BenchmarkTools.Trial: 14 samples with 1 evaluation.
 Range (min … max):  351.040 ms … 396.431 ms  ┊ GC (min … max): 2.58% … 1.81%
 Time  (median):     357.403 ms               ┊ GC (median):    2.86%
 Time  (mean ± σ):   360.863 ms ±  11.661 ms  ┊ GC (mean ± σ):  2.75% ± 0.80%

    █           
  ▇▁█▇▁▇▇▁▁▁▇▁▁▁▇▇▁▇▁▁▇▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇ ▁
  351 ms           Histogram: frequency by time          396 ms <

 Memory estimate: 99.63 MiB, allocs estimate: 14.

内存不是我现在主要关心的问题，我正在寻找一种方法来使这个过程更快、更优化。如果你创建一个更大的矩阵，比如 10_000x100(所以内存🥵)，速度会很烦人。因此，我正在寻找任何可以帮助我提高此过程速度的建议。

Answer 1

不要重新发明轮子，只要做

cor(test_mat, dims=2) 

这比您的代码快得多。

设置:

test_mat = rand(1:10, 2_000, 20)

现在基准:

julia> @btime calc_corr2($test_mat);
  709.354 ms (16 allocations: 99.63 MiB)

julia> @btime cor(test_mat, dims=2);
  52.679 ms (16 allocations: 30.85 MiB)