为什么 trace 没有保留在 tensordot 中
Why is trace not preserved in tensordot
我正在尝试取三个密度矩阵的张量积并用乘积表示。这些矩阵中的每一个都有迹线 1,理论上,乘积矩阵也应该如此。但是在 numpy 中这样做似乎有一些意想不到的效果。即使将中间数组重塑为二维形式也会给出相同的答案。
In [31]: rho1
Out[31]:
array([[0. , 0. , 0. , 0. , 0. ],
[0. , 0.1, 0. , 0. , 0. ],
[0. , 0. , 0.2, 0. , 0. ],
[0. , 0. , 0. , 0.3, 0. ],
[0. , 0. , 0. , 0. , 0.4]])
In [32]: np.trace(rho1)
Out[32]: 1.0
In [33]: rho2
Out[33]:
array([[0.2, 0. , 0. , 0. , 0. ],
[0. , 0.2, 0. , 0. , 0. ],
[0. , 0. , 0.2, 0. , 0. ],
[0. , 0. , 0. , 0.2, 0. ],
[0. , 0. , 0. , 0. , 0.2]])
In [34]: np.trace(rho2)
Out[34]: 1.0
In [35]: rho3
Out[35]:
array([[0.5, 0. , 0. , 0. , 0. ],
[0. , 0.5, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ]])
In [36]: np.trace(rho3)
Out[36]: 1.0
In [37]: rho = np.tensordot(rho1, np.tensordot(rho2, rho3, axes=0), axes=0)
In [38]: np.trace(rho.reshape(125, 125))
Out[38]: 0.010000000000000002
In [39]: rho = np.tensordot(rho1, np.tensordot(rho2, rho3, axes=0).reshape(25,25), axes=0)
In [40]: np.trace(rho.reshape(125, 125))
Out[40]: 0.010000000000000002
我真的找不到这段代码中的任何错误,所以我觉得我误解了 tensordot 和 reshape 的工作原理。但我并没有真正从文档中得到任何东西。有人可以帮我解决这个问题吗?
np.tensordot 不是张量积,而是张量的点积。张量积类似于外积。
简答:
我猜您正在寻找克罗内克张量积:np.kron()
:
rho = np.kron(rho1, np.kron(rho2, rho3))
而这一次:
np.trace(rho)
Out[22]: 1.0000000000000002
详情:
为什么您的解决方案不起作用?
因为 np.reshape()
没有使用“正确”的维度顺序重塑数组,您最终可能会排列一些维度以获得所需的结果:
rho = np.moveaxis(np.tensordot(rho1, np.moveaxis(np.tensordot(rho1, rho2, axes=0),[0,1,2,3],[0,2,3,1]).reshape((5,5)), axes=0),[0,1,2,3],[0,2,3,1]).reshape((125,125))
会输出正确的结果,但是既然np.kron
存在,就有点大材小用了。
其实你认为 np.tensordot(rho1,rho2,axes=0)
等同于:
np.einsum('ik,jl',rho1,rho2)
但实际上np.tensordot(rho1,rho2,axes=0)
计算:
np.einsum('ij,kl',rho1,rho2)
所以另一种获得正确答案的方法是使用:
np.einsum('ik,jl',rho1,rho2).reshape(5,5)
tensordot 文档确实说 axes=0 给出张量积:
axes = 0 : tensor product :math:`a\otimes b`,
数组的外积。但是你必须注意结果维度是如何排列的。
但是它说:
The shape of the result consists of the non-contracted axes of the
first tensor, followed by the non-contracted axes of the second.
使用 axes=0
,它实际上创建了 rho2[:,:,None]
和 rho3[:,None,:]
,并在正常的 dot
轴上计算乘积和。
In [37]: x=np.tensordot(rho2, rho3, 0)
In [38]: x.shape
Out[38]: (5, 5, 5, 5)
将其与默认 'ijkl' 输出的 einsum
进行比较:
In [40]: y=np.einsum('ij,kl',rho2,rho3)
In [41]: y.shape
Out[41]: (5, 5, 5, 5)
In [42]: np.allclose(x,y)
Out[42]: True
或 dot
等效项:
In [44]: z=np.dot(rho2[:,:,None],rho3[:,None,:])
In [45]: np.allclose(x,z)
Out[45]: True
前两个维度来自rho2
。
如果您真正想要的是 kron
,我们必须首先转置维度,交错 2 个数组中的维度:
In [46]: K=np.kron(rho2,rho3)
In [47]: K.shape
Out[47]: (25, 25)
In [48]: np.trace?
In [49]: np.allclose(x.transpose(0,2,1,3).reshape(25,25),K)
Out[49]: True
正在测试跟踪:
In [50]: K.trace()
Out[50]: 1.0
与 einsum
相同。注意索引的顺序:
In [53]: np.einsum('ij,kl->ikjl',rho2,rho3).reshape(25,25).trace()
Out[53]: 1.0
另一种查看方式是从 4d 数组中提取 2d 对角线值:
In [60]: j = np.arange(5); i=j[:,None]
In [61]: x[i,i,j,j]
Out[61]:
array([[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ]])
In [62]: _.sum()
Out[62]: 1.0
简单地重塑 x
,将这些值放在对角线上:
In [63]: x[i,j,i,j]
Out[63]:
array([[0.1, 0. , 0. , 0. , 0. ],
[0. , 0.1, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ]])
等效于 2 步跟踪:
In [75]: x.trace(0,0,1)
Out[75]:
array([[0.5, 0. , 0. , 0. , 0. ],
[0. , 0.5, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ]])
In [76]: x.trace(0,0,1).trace()
Out[76]: 1.0
我正在尝试取三个密度矩阵的张量积并用乘积表示。这些矩阵中的每一个都有迹线 1,理论上,乘积矩阵也应该如此。但是在 numpy 中这样做似乎有一些意想不到的效果。即使将中间数组重塑为二维形式也会给出相同的答案。
In [31]: rho1
Out[31]:
array([[0. , 0. , 0. , 0. , 0. ],
[0. , 0.1, 0. , 0. , 0. ],
[0. , 0. , 0.2, 0. , 0. ],
[0. , 0. , 0. , 0.3, 0. ],
[0. , 0. , 0. , 0. , 0.4]])
In [32]: np.trace(rho1)
Out[32]: 1.0
In [33]: rho2
Out[33]:
array([[0.2, 0. , 0. , 0. , 0. ],
[0. , 0.2, 0. , 0. , 0. ],
[0. , 0. , 0.2, 0. , 0. ],
[0. , 0. , 0. , 0.2, 0. ],
[0. , 0. , 0. , 0. , 0.2]])
In [34]: np.trace(rho2)
Out[34]: 1.0
In [35]: rho3
Out[35]:
array([[0.5, 0. , 0. , 0. , 0. ],
[0. , 0.5, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ]])
In [36]: np.trace(rho3)
Out[36]: 1.0
In [37]: rho = np.tensordot(rho1, np.tensordot(rho2, rho3, axes=0), axes=0)
In [38]: np.trace(rho.reshape(125, 125))
Out[38]: 0.010000000000000002
In [39]: rho = np.tensordot(rho1, np.tensordot(rho2, rho3, axes=0).reshape(25,25), axes=0)
In [40]: np.trace(rho.reshape(125, 125))
Out[40]: 0.010000000000000002
我真的找不到这段代码中的任何错误,所以我觉得我误解了 tensordot 和 reshape 的工作原理。但我并没有真正从文档中得到任何东西。有人可以帮我解决这个问题吗?
np.tensordot 不是张量积,而是张量的点积。张量积类似于外积。
简答:
我猜您正在寻找克罗内克张量积:np.kron()
:
rho = np.kron(rho1, np.kron(rho2, rho3))
而这一次:
np.trace(rho)
Out[22]: 1.0000000000000002
详情:
为什么您的解决方案不起作用?
因为 np.reshape()
没有使用“正确”的维度顺序重塑数组,您最终可能会排列一些维度以获得所需的结果:
rho = np.moveaxis(np.tensordot(rho1, np.moveaxis(np.tensordot(rho1, rho2, axes=0),[0,1,2,3],[0,2,3,1]).reshape((5,5)), axes=0),[0,1,2,3],[0,2,3,1]).reshape((125,125))
会输出正确的结果,但是既然np.kron
存在,就有点大材小用了。
其实你认为 np.tensordot(rho1,rho2,axes=0)
等同于:
np.einsum('ik,jl',rho1,rho2)
但实际上np.tensordot(rho1,rho2,axes=0)
计算:
np.einsum('ij,kl',rho1,rho2)
所以另一种获得正确答案的方法是使用:
np.einsum('ik,jl',rho1,rho2).reshape(5,5)
tensordot 文档确实说 axes=0 给出张量积:
axes = 0 : tensor product :math:`a\otimes b`,
数组的外积。但是你必须注意结果维度是如何排列的。
但是它说:
The shape of the result consists of the non-contracted axes of the
first tensor, followed by the non-contracted axes of the second.
使用 axes=0
,它实际上创建了 rho2[:,:,None]
和 rho3[:,None,:]
,并在正常的 dot
轴上计算乘积和。
In [37]: x=np.tensordot(rho2, rho3, 0)
In [38]: x.shape
Out[38]: (5, 5, 5, 5)
将其与默认 'ijkl' 输出的 einsum
进行比较:
In [40]: y=np.einsum('ij,kl',rho2,rho3)
In [41]: y.shape
Out[41]: (5, 5, 5, 5)
In [42]: np.allclose(x,y)
Out[42]: True
或 dot
等效项:
In [44]: z=np.dot(rho2[:,:,None],rho3[:,None,:])
In [45]: np.allclose(x,z)
Out[45]: True
前两个维度来自rho2
。
如果您真正想要的是 kron
,我们必须首先转置维度,交错 2 个数组中的维度:
In [46]: K=np.kron(rho2,rho3)
In [47]: K.shape
Out[47]: (25, 25)
In [48]: np.trace?
In [49]: np.allclose(x.transpose(0,2,1,3).reshape(25,25),K)
Out[49]: True
正在测试跟踪:
In [50]: K.trace()
Out[50]: 1.0
与 einsum
相同。注意索引的顺序:
In [53]: np.einsum('ij,kl->ikjl',rho2,rho3).reshape(25,25).trace()
Out[53]: 1.0
另一种查看方式是从 4d 数组中提取 2d 对角线值:
In [60]: j = np.arange(5); i=j[:,None]
In [61]: x[i,i,j,j]
Out[61]:
array([[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ],
[0.1, 0.1, 0. , 0. , 0. ]])
In [62]: _.sum()
Out[62]: 1.0
简单地重塑 x
,将这些值放在对角线上:
In [63]: x[i,j,i,j]
Out[63]:
array([[0.1, 0. , 0. , 0. , 0. ],
[0. , 0.1, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ]])
等效于 2 步跟踪:
In [75]: x.trace(0,0,1)
Out[75]:
array([[0.5, 0. , 0. , 0. , 0. ],
[0. , 0.5, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. ]])
In [76]: x.trace(0,0,1).trace()
Out[76]: 1.0