使用广播 Julia 将向量的向量转换为矩阵答案

【问题标题】：using broadcasting Julia for converting vector of vectors to matrices使用广播 Julia 将向量的向量转换为矩阵
【发布时间】：2021-01-01 15:35:11
【问题描述】：

我是 julia 新手，并且有一个婴儿任务来编写一个将向量向量转换为矩阵的函数。这很容易通过迭代元素来完成。

但是，我了解到广播往往更有效率。但我不确定如何在这里执行，因为 .= 操作无法工作，因为它会将向量读取为 1 x n 数组，因此会尝试在两个不同长度的数组上进行广播。

有没有办法广播？

我的代码在下面

function vecvec_to_matrix(vecvec)
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = zeros(Int64, dim1, dim2)
    for i in 1:dim1
        for j in 1:dim2
            my_array[i,j] = vecvec[i][j]
        end
    end
    return my_array
end

【问题讨论】：

代码中最大的性能瓶颈是内存布局。您将vecvec 向量中的向量存储为my_array 的行，并且在Julia 中按列处理数据效率更高（因此，如果您将它们存储在列中，则使用reduce(hcat, vecvec) 会更快，如注释所示下面）。
广播并不比迭代更有效。无论如何，广播都会在幕后进行迭代。

标签： julia

【解决方案1】：

您的方式已经直观且快速。您可以使用一些@inbounds 来提高性能，仅此而已。 vcat 也很快。我认为在您的情况下不需要广播。你以下是我能想到的各种方法的一些基准

function vecvec_to_matrix(vecvec)
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = zeros(Int64, dim1, dim2)
    for i in 1:dim1
        for j in 1:dim2
            my_array[i,j] = vecvec[i][j]
        end
    end
    return my_array
end

function vecvec_to_matrix2(vecvec::AbstractVector{T}) where T <: AbstractVector
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = Array{eltype(vecvec[1]), 2}(undef, dim1, dim2)
    @inbounds @fastmath for i in 1:dim1, j in 1:dim2
        my_array[i,j] = vecvec[i][j]
    end
    return my_array
end

function vecvec_to_matrix3(vecvec::AbstractVector{T}) where T <: AbstractVector
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = Array{eltype(vecvec[1]), 2}(undef, dim1, dim2)
    Threads.@threads for i in 1:dim1
        for j in 1:dim2
            my_array[i,j] = vecvec[i][j]
        end
    end
    return my_array
end

using Tullio

function using_tullio(vecvec::AbstractVector{T}) where T <: AbstractVector
    dim1 = length(vecvec)
    dim2 = length(vecvec[1])
    my_array = Array{eltype(vecvec[1]), 2}(undef, dim1, dim2)

    @tullio my_array[i, j] = vecvec[i][j]

    my_array
end

function using_vcat(vecvec::AbstractVector{T}) where T <: AbstractVector
    vcat(vecvec...)
end

using BenchmarkTools
vecvec =[rand(Int, 100) for i in 1:100];
@benchmark vecvec_to_matrix(vecvec)
@benchmark vecvec_to_matrix2(vecvec)
@benchmark vecvec_to_matrix3(vecvec)
@benchmark using_tullio(vecvec)
@benchmark using_vcat(vecvec)

有结果

julia> @benchmark vecvec_to_matrix(vecvec)
BenchmarkTools.Trial:
  memory estimate:  78.20 KiB
  allocs estimate:  2
  --------------
  minimum time:     12.701 μs (0.00% GC)
  median time:      15.001 μs (0.00% GC)
  mean time:        24.465 μs (10.98% GC)
  maximum time:     3.884 ms (98.30% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark vecvec_to_matrix2(vecvec)
BenchmarkTools.Trial:
  memory estimate:  78.20 KiB
  allocs estimate:  2
  --------------
  minimum time:     8.600 μs (0.00% GC)
  median time:      9.800 μs (0.00% GC)
  mean time:        19.532 μs (12.37% GC)
  maximum time:     3.834 ms (98.82% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark vecvec_to_matrix3(vecvec)
BenchmarkTools.Trial:
  memory estimate:  83.28 KiB
  allocs estimate:  32
  --------------
  minimum time:     8.399 μs (0.00% GC)
  median time:      14.600 μs (0.00% GC)
  mean time:        28.178 μs (11.82% GC)
  maximum time:     8.269 ms (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark using_tullio(vecvec)
BenchmarkTools.Trial:
  memory estimate:  78.20 KiB
  allocs estimate:  2
  --------------
  minimum time:     8.299 μs (0.00% GC)
  median time:      10.101 μs (0.00% GC)
  mean time:        19.476 μs (12.15% GC)
  maximum time:     3.661 ms (98.74% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark using_vcat(vecvec)
BenchmarkTools.Trial: 
  memory estimate:  78.20 KiB
  allocs estimate:  2
  --------------
  minimum time:     5.540 μs (0.00% GC)
  median time:      7.480 μs (0.00% GC)
  mean time:        16.236 μs (15.30% GC)
  maximum time:     876.400 μs (97.85% GC)
  --------------
  samples:          10000
  evals/sample:     5

【讨论】：

vcat 将产生不同形状的结果（它将堆叠向量）。此外，写reduce(vcat, vecvec) 效率更高，但您可能希望transpose(reduce(hcat, x)) 获得正确的形状。
reduce 是不是 reduce 效率更高？
是的，因为没有减少你会飞溅，减少会避免它。仅当您有数千个元素要 splat 时才可能很重要（否则您会将数千个位置参数传递给 vcat）。

【解决方案2】：

如果您的向量很短且大小固定（例如，3 维点列表），那么您应该强烈考虑使用StaticArrays package，然后调用reinterpret。演示：

julia> using StaticArrays

julia> A = rand(3, 8)
3×8 Array{Float64,2}:
 0.153872  0.361708  0.39703   0.405625  0.0881371  0.390133  0.185328  0.585539
 0.467841  0.846298  0.884588  0.798848  0.14218    0.156283  0.232487  0.22629
 0.390566  0.897737  0.569882  0.491681  0.499163   0.377012  0.140902  0.513979

julia> reinterpret(SVector{3,Float64}, A)
1×8 reinterpret(SArray{Tuple{3},Float64,1,3}, ::Array{Float64,2}):
 [0.153872, 0.467841, 0.390566]  [0.361708, 0.846298, 0.897737]  [0.39703, 0.884588, 0.569882]  …  [0.390133, 0.156283, 0.377012]  [0.185328, 0.232487, 0.140902]  [0.585539, 0.22629, 0.513979]

julia> B = vec(copy(ans))
8-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [0.1538721224514592, 0.467840786943454, 0.39056612358281706]
 [0.3617079493961777, 0.8462982350893753, 0.8977366743282564]
 [0.3970299970547111, 0.884587972864584, 0.5698823030478959]
 [0.40562472747685074, 0.7988484677138279, 0.49168126614394647]
 [0.08813706434793178, 0.14218012559727544, 0.499163319341982]
 [0.3901332827772166, 0.15628284837250006, 0.3770117394226711]
 [0.18532803309577517, 0.23248748941275688, 0.14090166962667428]
 [0.5855387782654986, 0.22628968661452897, 0.5139790762185006]

julia> reshape(reinterpret(Float64, B), (3, 8))
3×8 reshape(reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}), 3, 8) with eltype Float64:
 0.153872  0.361708  0.39703   0.405625  0.0881371  0.390133  0.185328  0.585539
 0.467841  0.846298  0.884588  0.798848  0.14218    0.156283  0.232487  0.22629
 0.390566  0.897737  0.569882  0.491681  0.499163   0.377012  0.140902  0.513979

【讨论】：