如何实现八度内任意阶张量的张量积？

Question

我无法访问 matlab，所以我正在尝试使用 octave 进行一些操作。您将如何有效地实现以下公式中描述的张量积？

我对任意阶张量 a 和 b 的处理方法如下

% Tensor product
function out = tp(a,b)
  if isvector(a)
    da = prod(size(a));
  else
    da = size(a);
  endif
  if isvector(b)
    db = prod(size(b));
  else
    db = size(b);
  endif
  out = reshape(a(:)*(b(:)'),[da,db]);
endfunction

if 语句只是为了捕捉 a 或 b 是向量的情况。我不知道这是否是一种有效的方法，因为我通常不编程而且我是八度音阶的新手。你的方法是什么？

我会使用Frobenius范数，见下图，看看是否与显式计算有区别。

下面是一些用于检查实现的显式计算。它工作正常，但我想问一下，是否有更好的方法来为任意阶张量执行此操作。谢谢！

% Frobenius norm
function out = nf(a)
  out = sqrt(a(:)'*a(:));
endfunction

% Tests for (m,n)
% (1,1)
disp("(1,1)")
a = rand(4,1);
b = rand(7,1);
c1 = tp(a,b);
c2 = zeros(4,7);
for i1=1:4 for i2=1:7
  c2(i1,i2) = a(i1)*b(i2);
endfor endfor
size(c1)
size(c2)
nf(c1-c2)

%(1,2)
disp("(1,2)")
a = rand(4,1);
b = rand(7,3);
c1 = tp(a,b);
c2 = zeros(4,7,3);
for i1=1:4 for i2=1:7 for i3=1:3
  c2(i1,i2,i3) = a(i1)*b(i2,i3);
endfor endfor endfor
size(c1)
size(c2)
nf(c1-c2)

%(2,1)
disp("(2,1)")
a = rand(4,2);
b = rand(1,3);
c1 = tp(a,b);
c2 = zeros(4,2,3);
for i1=1:4 for i2=1:2 for i3=1:3
  c2(i1,i2,i3) = a(i1,i2)*b(i3);
endfor endfor endfor
size(c1)
size(c2)
nf(c1-c2)

%(3,2)
disp("(3,2)")
a = rand(4,2,5);
b = rand(7,3);
c1 = tp(a,b);
c2 = zeros(4,2,5,7,3);
for i1=1:4 for i2=1:2 for i3=1:5 for i4=1:7 for i5=1:3
  c2(i1,i2,i3,i4,i5) = a(i1,i2,i3)*b(i4,i5);
endfor endfor endfor endfor endfor
size(c1)
size(c2)
nf(c1-c2)

编辑

我看了一下tensorlab包，看到下面Metahominid的回答，太棒了。出于好奇，我想检查我的实现、Andras Deak 的实现（见下面他的回答）和 tensorlab 包之间的时间性能。

% See Andras Deak answer
function c=tensorprod(a, b)
   b_inj = reshape(b, [ones(1,ndims(a)), size(b)]);
   c = a.*b_inj;
end

% Tests
a = rand(10,11,12);
b = rand(9,8,7);
tic; c1=outprod(a,b); t1=toc % tensorlab, see Metahominid's answer
tic; c2=tp(a,b); t2=toc % my approach
tic; c3=tensorprod(a,b); t3=toc % Andras Deak's approach
disp("Check size")
size(c1)
size(c2)
size(c3)
disp("Check Frobenius norm")
frob(c1) % from tensorlab
nf(c2)
disp("Check equality of elements")
nf(c1-c2)
nf(c1-c3)
disp("Compare time performance relative to tp(a,b)")
t1/t2
t3/t2

outprod 的 tensorlab 实现的计算时间 t1（对应于我的 tp）与 t2 的计算时间之比是考虑到的维度大约为 2-4（至少在我的计算机上是这样）。情况确实如此，因为在我的实现中，我不检查输入中的任何错误，也不捕获任何未定义的情况。 t3/t2 将 Andras Deak 的方法与我的方法进行比较，观察到几乎相同的结果。请不要误会我的意思，我不是想吹牛，而是想对可能对此感兴趣的人做一些最后的评论。结论：如果您需要一些用于小型张量的东西，我的简单实现可能对您有用，如果您需要更多东西，您绝对应该看看 tensorlab（请参阅下面的 Metahominid 的回答）。感谢您的回答和参考！

Answer 1

有个东西叫Tensorlab, which for all I can tell is usable with Octave.。您只需获取 link 即可下载它。

编辑：它同时具备这两个功能，而且速度会快得多。

[H,Heff] = hankelize(linspace(0,1,1000),'order',3);
tic; disp(frob(H)); toc; % Using the dense tensor H
tic; disp(frob(Heff)); toc; % Using the efficient representation of H
3.2181e+03
Elapsed time is 0.026401 seconds.
3.2181e+03
Elapsed time is 0.000832 seconds.

outprod(T1,T2)

是您想要的其他命令。

Answer 2

你得到的是广义克罗内克积。该定义非常适合在 Octave 中使用 array broadcasting 来实现。为此，您只需要将与另一个数组的维度一样多的前导单例维度注入到您的一个数组中。这已经足够了，因为每个数组都有无限数量的隐式尾随维度。

function c=tensorprod(a, b)
   b_inj = reshape(b, [ones(1,ndims(a)), size(b)]);
   c = a.*b_inj;
end

如果 a 的大小为 (i,j,k)，b 的大小为 (m,n)，b_inj 的大小为 (1,1,1,m,n)，并且 a已经与大小 (i,j,k,1,1) 隐式兼容。因此，将这两个数组按元素相乘即可得到所需的结果。

证明它应该按照您希望的方式工作：

octave:29> a = rand(2,3);
octave:30> b = rand(4,5);
octave:31> c = tensorprod(a,b);
octave:32> size(c)
ans =

   2   3   4   5

octave:33> c(1,3,2,3) == a(1,3)*b(2,3) % indices chosen by fair dice roll
ans =  1

如果你想以不同的方式处理向量（即你希望两个向量的张量积是一个二维矩阵），你需要自己处理这种特殊情况，因为 Octave 处理向量的方式是 row/column 矩阵。不过，这只是一个微不足道的并发症。

如何实现八度内任意阶张量的张量积？

How to implement tensor product for arbitrary order tensors in octave?

octave

tensor