logo标签列表

矩阵/张量三重积是什么？

matlabnumpymatrixmatrix-multiplicationblas
10

我正在研究的算法需要在一些地方计算一种矩阵三重积。该操作需要使用具有相同维度的三个方阵，并生成一个三指标张量。将运算数标记为A、B和C，结果中的(i,j,k)元素为：

X[i,j,k] = \sum_a A[i,a] B[a,j] C[k,a]
在 numpy 中，您可以使用 einsum('ia,aj,ka->ijk', A, B, C) 来计算这个操作。
问题：
  • 这个操作有标准名称吗？
  • 我能用单个 BLAS 调用来计算它吗？
  • 还有其他哪些高度优化的数值 C/Fortran 库可以计算这种类型的表达式？
1
请查看 MATLAB 张量工具箱版本 2.6：http://www.sandia.gov/~tgkolda/TensorToolbox/index-2.6.html - bla
2
你可以在Matlab中使用一个bsxfun和一个矩阵乘法来实现这个功能，我觉得这两个操作都很快速。这样对你来说可以吗？ - Luis Mendo
3个回答
7

介绍和解决方案代码

np.einsum非常强大，但在极少数情况下，如果你可以将矩阵乘法引入计算中，仍然可以超越它。经过几次尝试，似乎可以使用np.dot进行矩阵乘法来提高性能，而不是使用np.einsum('ia,aj,ka->ijk', A, B, C)。

基本思想是将“所有einsum”操作分解成以下np.einsum和np.dot的组合：

  • 对于A:[i,a]和B:[a,j]的求和使用np.einsum得到一个3D数组：[i,j,a]。
  • 然后将此3D数组重塑为2D数组：[i*j,a]，并且第三个数组C[k,a]被转置为[a,k]，以便在这两个数组之间执行矩阵乘法，从而得到[i*j,k]作为矩阵乘积，因为我们在那里失去了索引[a]。
  • 将乘积重塑为一个3D数组：[i,j,k]，作为最终输出。

以下是迄今为止讨论的第一个版本的实现 -

import numpy as np

def tensor_prod_v1(A,B,C):   # First version of proposed method
    # Shape parameters
    m,d = A.shape
    n = B.shape[1]
    p = C.shape[0]
    
    # Calculate \sum_a A[i,a] B[a,j] to get a 3D array with indices as (i,j,a)
    AB = np.einsum('ia,aj->ija', A, B)
    
    # Calculate entire summation losing a-ith index & reshaping to desired shape
    return np.dot(AB.reshape(m*n,d),C.T).reshape(m,n,p)

由于我们要对三个输入数组的a-th索引进行求和，因此可以有三种不同的方法沿着a-th索引进行求和。之前列出的代码是针对(A,B)的。因此，我们还可以有(A,C)和(B,C)，这给我们带来了另外两种变化，如下所示：

def tensor_prod_v2(A,B,C):
    # Shape parameters
    m,d = A.shape
    n = B.shape[1]
    p = C.shape[0]
    
    # Calculate \sum_a A[i,a] C[k,a] to get a 3D array with indices as (i,k,a)
    AC = np.einsum('ia,ja->ija', A, C)
    
    # Calculate entire summation losing a-ith index & reshaping to desired shape
    return np.dot(AC.reshape(m*p,d),B).reshape(m,p,n).transpose(0,2,1)
    
def tensor_prod_v3(A,B,C):
    # Shape parameters
    m,d = A.shape
    n = B.shape[1]
    p = C.shape[0]
    
    # Calculate \sum_a B[a,j] C[k,a] to get a 3D array with indices as (a,j,k)
    BC = np.einsum('ai,ja->aij', B, C)
    
    # Calculate entire summation losing a-ith index & reshaping to desired shape
    return np.dot(A,BC.reshape(d,n*p)).reshape(m,n,p)

根据输入数组的形状，不同的方法会相对于彼此产生不同的加速效果，但我们希望所有方法都比all-einsum方法更好。性能数据列在下一节中。

运行时测试

这可能是最重要的部分，因为我们试图查看所提出的三种方法与问题中最初提出的all-einsum方法相比的加速数值。

数据集＃1（等形状的数组）：

In [494]: L1 = 200
     ...: L2 = 200
     ...: L3 = 200
     ...: al = 200
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [495]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum('ia,aj,ka->ijk', A, B, C)
     ...: 
1 loops, best of 3: 470 ms per loop
1 loops, best of 3: 391 ms per loop
1 loops, best of 3: 446 ms per loop
1 loops, best of 3: 3.59 s per loop

数据集 #2（更大的A）：

In [497]: L1 = 1000
     ...: L2 = 100
     ...: L3 = 100
     ...: al = 100
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [498]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum('ia,aj,ka->ijk', A, B, C)
     ...: 
1 loops, best of 3: 442 ms per loop
1 loops, best of 3: 355 ms per loop
1 loops, best of 3: 303 ms per loop
1 loops, best of 3: 2.42 s per loop

数据集 #3 (更大的 B) :

In [500]: L1 = 100
     ...: L2 = 1000
     ...: L3 = 100
     ...: al = 100
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [501]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum('ia,aj,ka->ijk', A, B, C)
     ...: 
1 loops, best of 3: 474 ms per loop
1 loops, best of 3: 247 ms per loop
1 loops, best of 3: 439 ms per loop
1 loops, best of 3: 2.26 s per loop

数据集#4（更大的C）:


In [503]: L1 = 100
     ...: L2 = 100
     ...: L3 = 1000
     ...: al = 100
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)

In [504]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum('ia,aj,ka->ijk', A, B, C)
     ...: 
1 loops, best of 3: 250 ms per loop
1 loops, best of 3: 358 ms per loop
1 loops, best of 3: 362 ms per loop
1 loops, best of 3: 2.46 s per loop

数据集 #5（更长的a维度长度）:


In [506]: L1 = 100
     ...: L2 = 100
     ...: L3 = 100
     ...: al = 1000
     ...: 
     ...: A = np.random.rand(L1,al)
     ...: B = np.random.rand(al,L2)
     ...: C = np.random.rand(L3,al)
     ...: 

In [507]: %timeit tensor_prod_v1(A,B,C)
     ...: %timeit tensor_prod_v2(A,B,C)
     ...: %timeit tensor_prod_v3(A,B,C)
     ...: %timeit np.einsum('ia,aj,ka->ijk', A, B, C)
     ...: 
1 loops, best of 3: 373 ms per loop
1 loops, best of 3: 269 ms per loop
1 loops, best of 3: 299 ms per loop
1 loops, best of 3: 2.38 s per loop

结论：我们发现，相比于问题中列出的 all-einsum 方法，采用提出的方法的变体可以实现 8x-10x 的加速。

6

设矩阵大小为nxn。在Matlab中，您可以执行以下操作：

  1. 将A和C组合成一个n^2xn矩阵AC，其中AC的行对应于A和C的所有行的组合。
  2. 用B后置乘以AC。这会给出所需的结果，只是形式不同。
  3. 重新调整大小和排列维度以获得所需的形式的结果。

代码：

AC = reshape(bsxfun(@times, permute(A, [1 3 2]), permute(C, [3 1 2])), n^2, n); % // 1
X = permute(reshape((AC*B).', n, n, n), [2 1 3]);                               %'// 2, 3

使用逐字循环的方法进行检查：

%// Example data:
n = 3;
A = rand(n,n);
B = rand(n,n);
C = rand(n,n);

%// Proposed approach:
AC = reshape(bsxfun(@times, permute(A, [1 3 2]), permute(C, [3 1 2])), n^2, n);
X = permute(reshape((AC*B).', n, n, n), [2 1 3]); %'

%// Loop-based approach:
Xloop = NaN(n,n,n); %// initiallize
for ii = 1:n
    for jj = 1:n
        for kk = 1:n
            Xloop(ii,jj,kk) = sum(A(ii,:).*B(:,jj).'.*C(kk,:)); %'
        end
    end
end

%// Compute maximum relative difference:
max(max(max(abs(X./Xloop-1))))

ans =
    2.2204e-16

最大相对差的数量级为eps，因此结果在数值精度范围内是正确的。
1
总之，使用numpy吧。抱歉，我忍不住要说一句 ;) - Eelco Hoogendoorn
0

我知道这个话题现在有点老了，但它经常被提起。在Matlab中，很难超越Jason Farquhar编写的tprod，可以在这里找到一个MEX文件。

https://www.mathworks.com/matlabcentral/fileexchange/16275-tprod-arbitary-tensor-products-between-n-d-arrays

tprod的工作方式与einsum非常相似，但它仅限于二元操作（2个张量）。这可能并不是真正的限制，因为我怀疑einsum只是执行一系列二元操作。这些操作的顺序很重要，我的理解是einsum只是按照传递数组的顺序执行它们，并且不允许多个中间产品。

tprod也仅限于密集（全）数组。Kolda的Tensor Toolbox（在早期帖子中提到）支持稀疏张量，但其功能比tprod更受限制（它不允许输出中有重复的索引）。我正在努力填补这些空白，但如果Mathworks能够做到这一点，那不是很好吗？

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接