Julia 中没有内存分配的外积

Question

我正在使用 Julia 1.9.0。

我希望在程序的每次迭代中，张量

T

被向量

A1, A2, ... , AD

的外积覆盖。例如，如果

D

为 5，我将张量

T

更新为

T[i1,i2,i3,i4,i5] .= A1[i1] * A2[i2] * A3[i3] * A4[i4] * A5[i5]

在我之前的问题中，我得到了一个名为

outer

的函数，它从向量

T

中获取

A1, A2, ..., AD

。

outer(v...) = reshape(kron(reverse(v)...),length.(v))

但是，它需要内存分配。

using BenchmarkTools
T = rand(3,4,8,5);
v = [rand(3), rand(4), rand(8), rand(5)];
@btime $T .= outer($v...)
# 2.582 μs (36 allocations: 6.91 KiB)

在我的程序中，我会反复更新

T

。所以内存分配会导致性能下降。我们怎样才能在没有（或更少）内存分配的情况下创建这样一个函数

outer!(T, v)

？

我担心的是，

kron!

不接受多个向量，而

kron

接受。这个事实让我很难开发这个功能

outer!(T, v)

。

a = rand(4); b = rand(3); ab = rand(4*3);
kron!(ab, a, b); # it works!

a = rand(4); b = rand(3); c = rand(4); ABC = rand(4*3*4)
kron!(ABC, a, b, c);
ERROR: MethodError: no method matching kron!(::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64},)

--

我猜

@Tulio

可能是另一种方法，如下所示。

@tullio T[i1,i2,i3,i4.i5] := A1[i1]*A2[i2]*A3[i3]*A4[i4]*A5[i5]

但是我需要实现通用的程序

D

，我认为我不应该使用@Tulio。

Answer 1

将

ABC

的定义调整为大小为 (4,3,4) 的张量，使之成为可能：

julia> a = rand(4); b = rand(3); c = rand(4);

julia> ABC = similar(a, (4,3,4));

julia> for I in CartesianIndices(ABC)
       ABC[I] = prod(getindex.((a,b,c),Tuple(I)))
       end;

for

循环可以包装在一个函数中来模仿

cron

，但在实际代码中它的可读性非常好。