如何在Julia中有效地计算矩阵的每对列的乘积?

Question

我想计算一个矩阵的每一对列的乘积:

┌─     ─┐         ┌─┐ ┌──┐ ┌──┐
│ 1 3 4 │ ────►   │3│ │4 │ │12│
│ 2 1 5 │         │2│ │10│ │5 │
└─     ─┘         └─┘ └──┘ └──┘

我试着:

julia> using Combinatorics
julia> prods = collect(
Iterators.map(
idx->arr[:, idx[1]] .* arr[:, idx[2]],
combinations(axes(arr, 2), 2)
)
)
3-element Vector{Vector{Int64}}:
[3, 2]
[4, 10]
[12, 5]

由于编译器每次都要编译匿名函数，我复制了6次并收集结果，这对我来说效率不高。最好的方法是什么?

Answer 1

试试这个base julia

julia> col = size(a, 2);
julia> [a[:,i] .* a[:,j] for i in 1:col for j in (i+1):col]
3-element Vector{Vector{Int64}}:
[3, 2]
[4, 10]
[12, 5]

如果你想要一个矩阵，试试reduce

julia> reduce(hcat, [a[:,i] .* a[:,j] for i in 1:col for j in (i+1):col])
2×3 Matrix{Int64}:
3   4  12
2  10   5

julia> a = [[1, 2] [3, 1] [4, 5]];

Answer 2

一个经典的循环将比组合运算快5倍。

function prod_col(arr)
m, n = size(arr)
out = similar(arr, m, n*(n-1)÷2)
c = 1
for i = 1:n, j = i+1:n
for k = 1:m
out[k,c] = arr[k,i] * arr[k,j]
end
c += 1
end
out
end

用10 × 10矩阵测试得到:

arr = rand(0:9,10,10)
@benchmark prod_col($arr)
@benchmark prod_cols_pairs($arr)
BenchmarkTools.Trial: 10000 samples with 195 evaluations.
Range (min … max):  526.667 ns … 48.783 μs  ┊ GC (min … max):  0.00% … 97.62%
Time  (median):       1.969 μs              ┊ GC (median):     0.00%
Time  (mean ± σ):     2.328 μs ±  3.950 μs  ┊ GC (mean ± σ):  21.19% ± 11.86%
▅▆█▄▂▂                                                       ▁
██████▄▃▁▁▁▁▁▁▁▁▄▅▄▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▄▆█ █
527 ns        Histogram: log(frequency) by time      33.4 μs <
Memory estimate: 3.69 KiB, allocs estimate: 1.
BenchmarkTools.Trial: 10000 samples with 3 evaluations.
Range (min … max):   7.967 μs …  1.640 ms  ┊ GC (min … max):  0.00% … 98.25%
Time  (median):      8.767 μs              ┊ GC (median):     0.00%
Time  (mean ± σ):   10.902 μs ± 39.061 μs  ┊ GC (mean ± σ):  13.31% ±  3.74%
▅█▆▄▃▄▄▂                                                    ▂
██████████████▇▆▅▆▆▅▄▄▄▄▄▄▁▃▄▄▁▁▁▄▁▁▃▃▄▃▃▁▁▃▁▁▁▁▄▃▁▁▁▁▁▁▃▃▃ █
7.97 μs      Histogram: log(frequency) by time      33.7 μs <
Memory estimate: 22.88 KiB, allocs estimate: 191.

Answer 3

我认为以下是一个更好的方法:

julia> using Combinatorics
julia> function prod_cols_pairs(mat::Matrix{T})::Matrix{T} where T
idxs = combinations(axes(mat, 2), 2)
arr[:, getindex.(idxs, 1)] .* arr[:, getindex.(idxs, 2)]
end;
julia> prod_cols_pairs(arr)
2×3 Matrix{Int64}:
3   4  12
2  10   5
# Another example
julia> arr = [
1 3 4 2
2 1 5 2
]
2×4 Matrix{Int64}:
1  3  4  2
2  1  5  2
julia> prod_cols_pairs(arr)
2×6 Matrix{Int64}:
3   4  2  12  6   8
2  10  4   5  2  10

1. I avoidcollecting.
2. 我不再定义一个匿名函数。
3. 我将拷贝数从6个减少到2个。我认为这是一个很好的改进，因为如果我的输入有6列，这将导致副本数量从38减少到2(想象一下如果我有50列的矩阵…)。如果我用@views，它可以降到零。然而，我不知道在大型Matrixes的情况下，如果我合并它，性能会如何。

如果我说错了，请纠正我。