python - NumPy：矩阵相对于自身的雅可比矩阵和对角矩阵的 2d 扩展

Question

矩阵相对于自身的雅可比行列式

我正在实现一个内部自动微分模块，仅使用 NumPy 的本机函数，对于任何类型的矩阵运算，从 2D 数组构造一个 4D 数组，如图片中的那个似乎出现在不同的地方。

我目前的方法很简单：如果给我一个名为 的 k-by-d 矩阵a，我正在做类似的事情

r = np.zeros([*a.shape, *a.shape])
for i, j in np.ndindex(a.shape):
    r[i, j, i, j] = a[i,j]

但是，我意识到这是diag()我输入 2D 数组以构造 4D 数组的位置的概括。NumPy 中是否有任何本机方法支持此功能？我主要担心我拥有的for循环最终会成为瓶颈。

感谢您提前回复！

标签： pythonnumpymatrix

Answer 1

如果所有元素都是 0 或 1（我相信你的图片显示），那么：

r = np.identity(np.multiply.reduce(a.shape)).reshape(a.shape + a.shape)

相反，如果您想要快速等效于您编写的循环，则使用np.einsum()返回广义对角线的可写视图r并将您的数组分配给它：

r = np.zeros(a.shape + a.shape)
np.einsum('ijij->ij', r)[:] = a

例子

a = np.arange(6).reshape((2,3))
>>> a
array([[0, 1, 2],
       [3, 4, 5]])

# Jacobian d/da of a
>>> np.identity(np.multiply.reduce(a.shape)).reshape(a.shape + a.shape)
array([[[[1., 0., 0.],
         [0., 0., 0.]],

        [[0., 1., 0.],
         [0., 0., 0.]],

        [[0., 0., 1.],
         [0., 0., 0.]]],


       [[[0., 0., 0.],
         [1., 0., 0.]],

        [[0., 0., 0.],
         [0., 1., 0.]],

        [[0., 0., 0.],
         [0., 0., 1.]]]])

# set a as the generalized diagonal of a n,m,n,m matrix
r = np.zeros(a.shape + a.shape)
np.einsum('ijij->ij', r)[:] = a
>>> r
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., 0.]],

        [[0., 0., 0.],
         [0., 0., 5.]]]])