假设我有一个名为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([mat[3, :], mat[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
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
- 在大矩阵上:
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(因此内存🥵)，速度将令人讨厌。因此，我正在寻找任何可以帮助我实现这个过程的更高速度的建议。
不要做无谓的重复，只管去做
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)