如何为八度音阶中的任意阶张量实现张量积？

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 作为向量的情况。我不知道这是否是一种有效的方法，因为我通常不编程，而且我是 octave 的新手。你的方法是什么？

我将使用 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）与 tp 的计算时间 t2 的比率是针对 2-4 左右的考虑维度（至少在我的计算机上）。肯定是这种情况，因为在我的实现中，我不检查输入中的任何错误，也不捕获任何未定义的情况。将 Andras Deak 的方法与我的方法进行比较，在 t3/t2 中观察到几乎相同。请不要误会我的意思，我不是想吹牛，而是要对可能对此感兴趣的人发表一些最后的评论。结论：如果你需要一些小张量的东西，我的简单实现可能对你有用，如果你需要更多东西，你一定要看看 tensorlab（见下面 Metahominid 的回答）。感谢您的回答和参考！

Answer 1

您所拥有的是一个通用的克罗内克产品。该定义非常适合在 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 将向量作为行/列矩阵处理的方式.不过，这只是一个微不足道的并发症。

Answer 2

有一个叫做Tensorlab 的东西，据我所知，它可以与Octave. 一起使用。您只需通过获取链接下载它。

编辑：它同时具备这两种功能，而且速度会大大加快。

[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)

是您想要的另一个命令。