如何用numpy来表示这个张量

Question

设x是形状为(a,)的矢量，T是形状为(b, a, a)的张量。如果我想计算(x^T)Tx，我可以用x.dot(w.dot(x).transpose())来做。

例如:

x = np.array([1.,2.,3.,4.,5.])
w = np.array([[[1.,2.,3.,4.,5.],
               [1.,2.,3.,4.,5.],
               [1.,2.,3.,4.,5.]],
              [[1.,2.,3.,4.,5.],
               [1.,2.,3.,4.,5.],
               [1.,2.,3.,4.,5.]]])
x.dot(w.dot(x).transpose())

但是，如果我想将T分解为形状为(b,a,r)、(b,r,a)和r<<a的两个张量P和Q(低秩表达式)，那么T中的每个矩阵(a*a)分解为a*r和r*a，这将减少大量数据。那么我怎么用numpy计算(x^T)PQx呢?

Answer 1

你的例子有问题。

x.shape  (5,)
w.shape  (2,3,5)
x.dot(w.dot(x).transpose())
   ValueError: matrices are not aligned

但是使用你的描述:

`x` `(a,)`,  `T` `(b,a,a)`; `(x^T)Tx`

在考虑复积时，我喜欢使用einsum(爱因斯坦求和)。我想你的x'Tx是:

np.einsum('i,kij,j->k', x, T, x)

T分解为:P (b,a,r), Q (b,r,a);

np.einsum('kir,krj->kij', P,Q) == T

:

np.einsum('i,kir,krj,j->k', x, P, Q, x)

einsum在维数较大时不是最好的，因为k,i,j,r的组合迭代空间可能较大。然而，这是一个思考问题的有用方法。

我认为它可以重写为3 dots:

P1 = np.einsum('i,kir->kr', x, P)
Q1 = np.einsum('krj,j->kr', Q, x)
np.einsum('kr,kr->k', P1, Q1)

---

计算示例:

In [629]: a,b,r = 5,3,2
In [630]: x=np.arange(1.,a+1)
In [632]: P=np.arange(b*a*r).reshape(b,a,r)
In [633]: Q=np.arange(b*a*r).reshape(b,r,a)
In [635]: T=np.einsum('kir,krj->kij',P,Q)
In [636]: P
Out[636]: 
array([[[ 0,  1],
        [ 2,  3],
        [ 4,  5],
        [ 6,  7],
        ...
        [24, 25],
        [26, 27],
        [28, 29]]])
In [637]: Q
Out[637]: 
array([[[ 0,  1,  2,  3,  4],
        [ 5,  6,  7,  8,  9]],
        ...
       [[20, 21, 22, 23, 24],
        [25, 26, 27, 28, 29]]])
In [638]: T
Out[638]: 
array([[[   5,    6,    7,    8,    9],
        [  15,   20,   25,   30,   35],
        [  25,   34,   43,   52,   61],
        [  35,   48,   61,   74,   87],
        ...
        [1105, 1154, 1203, 1252, 1301],
        [1195, 1248, 1301, 1354, 1407],
        [1285, 1342, 1399, 1456, 1513]]])
In [639]: T.shape
Out[639]: (3, 5, 5)
In [640]: R1=np.einsum('i,kij,j->k',x,T,x)
...
In [642]: R1
Out[642]: array([  14125.,  108625.,  293125.])
In [643]: R2=np.einsum('i,kir,krj,j->k',x,P,Q,x)
In [644]: R2
Out[644]: array([  14125.,  108625.,  293125.])
In [645]: P1=np.einsum('i,kir->kr',x,P)
In [646]: Q1=np.einsum('krj,j->kr',Q,x)
In [647]: R3=np.einsum('kr,kr->k',P1,Q1)
In [648]: R3
Out[648]: array([  14125.,  108625.,  293125.])
In [649]: P1
Out[649]: 
array([[  80.,   95.],
       [ 230.,  245.],
       [ 380.,  395.]])
In [650]: Q1
Out[650]: 
array([[  40.,  115.],
       [ 190.,  265.],
       [ 340.,  415.]])

最后一组计算可以用dot

完成

In [656]: np.dot(x,P)
Out[656]: 
array([[  80.,   95.],
       [ 230.,  245.],
       [ 380.,  395.]])
In [657]: np.dot(Q,x)
Out[657]: 
array([[  40.,  115.],
       [ 190.,  265.],
       [ 340.,  415.]])
In [658]: np.dot(np.dot(x,P),np.dot(Q,x).T)
Out[658]: 
array([[  14125.,   40375.,   66625.],
       [  37375.,  108625.,  179875.],
       [  60625.,  176875.,  293125.]])

但是我们只需要最后一个dot的对角线。产品的总和越简单越好:

In [661]: (P1*Q1).sum(axis=1)
Out[661]: array([  14125.,  108625.,  293125.])