不知道这个标题有没有意义。通常单位矩阵是一个二维矩阵,如
In [1]: import numpy as np
In [2]: np.identity(2)
Out[2]:
array([[ 1., 0.],
[ 0., 1.]])
并且没有第三维。
Numpy 可以给我全零的 3D 矩阵
In [3]: np.zeros((2,2,3))
Out[3]:
array([[[ 0., 0., 0.],
[ 0., 0., 0.]],
[[ 0., 0., 0.],
[ 0., 0., 0.]]])
但我想要一个“3D 单位矩阵”,因为前 2 个维度上的所有对角线元素都是 1。例如,对于形状 (2,2,3) 应该是
array([[[ 1., 1., 1.],
[ 0., 0., 0.]],
[[ 0., 0., 0.],
[ 1., 1., 1.]]])
有什么优雅的方式来生成这个吗?
从 2d 单位矩阵开始,您可以通过以下两个选项来制作 “3d 单位矩阵”:
import numpy as np
i = np.identity(2)
选项 1:沿第三维堆叠 2d 单位矩阵
np.dstack([i]*3)
#array([[[ 1., 1., 1.],
# [ 0., 0., 0.]],
# [[ 0., 0., 0.],
# [ 1., 1., 1.]]])
选项 2:重复值然后重塑
np.repeat(i, 3, axis=1).reshape((2,2,3))
#array([[[ 1., 1., 1.],
# [ 0., 0., 0.]],
# [[ 0., 0., 0.],
# [ 1., 1., 1.]]])
选项 3:创建一个零数组,并使用 高级索引:
将 1 分配给对角线元素(第一维和第二维)的位置shape = (2,2,3)
identity_3d = np.zeros(shape)
idx = np.arange(shape[0])
identity_3d[idx, idx, :] = 1
identity_3d
#array([[[ 1., 1., 1.],
# [ 0., 0., 0.]],
# [[ 0., 0., 0.],
# [ 1., 1., 1.]]])
时间:
%%timeit
shape = (100,100,300)
i = np.identity(shape[0])
np.repeat(i, shape[2], axis=1).reshape(shape)
# 10 loops, best of 3: 10.1 ms per loop
%%timeit
shape = (100,100,300)
i = np.identity(shape[0])
np.dstack([i] * shape[2])
# 10 loops, best of 3: 47.2 ms per loop
%%timeit
shape = (100,100,300)
identity_3d = np.zeros(shape)
idx = np.arange(shape[0])
identity_3d[idx, idx, :] = 1
# 100 loops, best of 3: 6.31 ms per loop
一种方法是初始化一个
2D
单位矩阵,然后将其广播到3D
。因此,以 n
作为沿前两个轴的长度,以 r
作为最后一个轴的长度,我们可以这样做 -
np.broadcast_to(np.identity(n)[...,None], (n,n,r))
运行示例让事情变得更清楚 -
In [154]: i = np.identity(3); i # Create an identity matrix
Out[154]:
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])
# Extend it to 3D. This helps us broadcast to reqd. shape later on
In [152]: i[...,None]
Out[152]:
array([[[ 1.],
[ 0.],
[ 0.]],
[[ 0.],
[ 1.],
[ 0.]],
[[ 0.],
[ 0.],
[ 1.]]])
# Broadcast to (n,n,r) shape for the 3D identity matrix
In [153]: np.broadcast_to(i[...,None], (3,3,3))
Out[153]:
array([[[ 1., 1., 1.],
[ 0., 0., 0.],
[ 0., 0., 0.]],
[[ 0., 0., 0.],
[ 1., 1., 1.],
[ 0., 0., 0.]],
[[ 0., 0., 0.],
[ 0., 0., 0.],
[ 1., 1., 1.]]])
这种方法可以提高性能,因为它只是生成单位矩阵的视图。因此,在这种形式下,输出将是一个只读数组。如果您需要一个拥有自己的内存空间的可写数组,只需在此处附加一个
.copy()
即可。
断言性能,这是一个时序测试,用于创建形状为
3D
的单位矩阵:(100, 100, 300)
-
In [140]: n,r = 100,300
In [141]: %timeit np.broadcast_to(np.identity(n)[...,None], (n,n,r))
100000 loops, best of 3: 9.29 µs per loop
与@Psidom类似,使用高级
np.einsum
%%timeit
shape = (100,100,300)
identity_3d = np.zeros(shape)
np.einsum('iij->ij', identity_3d)[:] = 1
1000 loops, best of 3: 251 µs per loop
%%timeit
shape = (100,100,300)
identity_3d = np.zeros(shape)
idx = np.arange(shape[0])
identity_3d[idx, idx, :] = 1
1000 loops, best of 3: 320 µs per loop
%timeit np.broadcast_to(np.identity(100)[...,None], (100,100,300)).copy()
100 loops, best of 3: 12.1 ms per loop
我假设您想要写入单位矩阵的
copy()
,因为 anyarray.dot(identity)
更容易通过 np.broadcast_to(anyarray[..., None], a.shape + (300,))
计算。如果你真的只想要一个裸露的identity
矩阵,那么@Divakar的解决方案比其他任何解决方案都快得多,
%timeit np.broadcast_to(np.identity(100)[...,None], (100,100,300))
The slowest run took 5.16 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 20.8 µs per loop
但是
broadcast_to(anyarray[..., None], . . . )
可能比这更快。
我喜欢@Daniel F 的回答,但多看它给了我一个想法。由于最初的问题要求一种优雅的方式来做到这一点,我相信有一个非常优雅的 einsum 解决方案:
np.einsum("ij,ik->ijk", np.eye(2), np.ones((2, 3)))
这当然效率不高。