关于单个numpy数组中多个numpy数组的维度初始化的问题
Questions regarding the dimension initialization of multiple numpy arrays within a single numpy array
假设我们有 3 个 Pauli matrices，每个维度为 (2x2)。如下图：
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
现在，如果我将这些每个单独的 (2x2) 矩阵作为另一个 (2x2) 矩阵的条目。说：
A = np.array([[X, 0], [0, Y]])
B = np.array([[X, X], [Y, Y]])
奇怪的是，A 的亮度为 (2x2) - 这正是我想要的 - B 的亮度为 (2, 2, 2, 2) 不管这是什么，如下所示
A = np.array([[X, 0], [0, Y]])
A.shape
Out[131]: (2, 2)
B = np.array([[X, X], [Y, Y]])
B.shape
Out[133]: (2, 2, 2, 2)
另一方面，设 C 为 (1x3) 矩阵，D 为 (1x2) 矩阵，例如
C = np.array([[X, Y, 0]])
D = np.array([[X, Y]])
再次查看初始化矩阵的维度
C = np.array([[X, Y, 0]])
C.shape
Out[136]: (1, 3)
D = np.array([[X, Y]])
D.shape
Out[138]: (1, 2, 2, 2)
所以似乎每当我在这样的数组中初始化数组时，如果有混合数据类型作为条目，即矩阵和整数，如 A 和 C，它给了我合理的形状我想要的，即维度 (2,2)，每个条目的 "hidden" 维度 (2x2)。但是一旦条目只是像 B 和 D 中那样的严格矩阵，它就会给我不敏感的维度，例如 (2, 2, 2, 2)。所以我的问题是：
如何使用严格的 (2, 2) 矩阵作为条目初始化 (n, n) numpy 数组（矩阵），并仍然保留其 (n, n) 维度，即而不是给出我有些奇怪的 numpy 维度 (w, x, y, z)?
我想要这个的原因是因为我正在使用量子力学中的运算符进行计算，使用这些 Pauli 矩阵，例如 X、Y 和 Z，作为量子计算中的量子门。所以如果我有一些状态 rho 这也是一个 (2x2) 矩阵。让
rho = np.array([[1, 0],
[0, 0]])
并设 RHO 为 (2x2) 对角矩阵，其条目为 (2x2) rho 矩阵。
RHO = np.array([[rho, 0],
[0, rho]])
我希望计算类似 np.dot(D, RHO) 的东西，这样它可以给出
np.array([np.dot(X, rho), 0],
[0, np.dot(Y, rho)])
而且我在 python 上检查过，两个 (2x2) 矩阵的点积以 (2x2) 矩阵作为条目，其条目乘法也将全部是点积。
我上面讨论的所有内容的动机是我希望使用这些属性作为向量化我的算法的手段。目前我的算法的一个非常粗略的例子看起来像这样：
for i in (integer_1):
for j in (integer_2):
#do something that involves operations with sums of dot products of many matrices#
并将其矢量化，使其有可能成为
for i in (integer_1):
#do something that involves operations with multiples of matrices with dot product of matrices as its entries#
哪个可能有效，也可能无效！但我很好奇我的这种方法是否会产生加速。
我希望我已经很好地解释了我的问题。
提前致谢。
编辑(1)
I've added latex formatted maths so hopefully you can understand what I'm trying to do.
def compute_channel_operation(rho, operators):
"""
Given a quantum state's density function rho, the effect of the
channel on this state is:
rho -> sum_{i=1}^n E_i * rho * E_i^dagger
Args:
rho (2x2 matrix): Density function
operators (list): List of operators(matrices)
Returns:
number: The result of applying the list of operators
"""
operation = [E@rho@E.conj().T for i, E in enumerate(operators)]
return np.sum(operation, axis=0)
So then I was hoping that instead of using loops, this direct multiplication of tensor method might be quicker as I scale up my simulation, say having to do this 1 million times
这里的 K 应该是任何 (1xn) 维度的列表，即 [I] 或 [I, X] 或 [I, X, Y] 或 [I, X, Y, Z]。我知道这里 X = X^{\dagger} Y 和 Z 也是如此，但我在模拟中会遇到不会出现这种情况的情况。
希望我已经解释清楚了。
(2, 2, 2, 2) 不是 奇怪的维度 ，它只是一个形状为 2x2x2x2
的 4D 张量
您看到 A 和 B 不同形状的原因是因为您使用标量 0 而不是 2x2 零矩阵来设置 A。改成
A = np.array([[X, np.zeros((2, 2))], [np.zeros((2, 2)), Y]])
B = np.array([[X, X], [Y, Y]])
你将得到 2x2x2x2 张量。
或改为
C = np.vstack([
np.hstack([X, np.zeros((2, 2))]),
np.hstack([np.zeros((2, 2)), Y])
])
D = np.vstack([
np.hstack([X, X]),
np.hstack([Y, Y])
])
你会得到 4x4 矩阵。
您还可以使用
从一种形式转换为另一种形式
E = A.transpose(0, 2, 1, 3).reshape(4, 4)
F = B.transpose(0, 2, 1, 3).reshape(4, 4)
np.allclose(C, E) # True
np.allclose(D, F) # True
来回
G = E.reshape(2, 2, 2, 2).transpose(0, 2, 1, 3)
H = F.reshape(2, 2, 2, 2).transpose(0, 2, 1, 3)
np.allclose(A, G) # True
np.allclose(B, H) # True
编辑： 关于你的函数 compute_channel_operation()，如果你不执行列表推导而是向量化操作并传入 3D，你可以大大加快它的速度您所有操作的张量
rho = np.random.rand(2, 2)
operators = [np.random.rand(2, 2) for _ in range(1000)]
operators_tensor = np.asarray(operators) # same as np.random.rand(1000, 2, 2)
def compute_channel_operation(rho, operators):
operation = [E@rho@E.conj().T for i, E in enumerate(operators)]
return np.sum(operation, axis=0)
def compute_channel_operation2(rho, operators):
return np.sum(operators @ rho @ operators.transpose(0, 2, 1).conj(), axis=0)
A = compute_channel_operation(rho, operators)
B = compute_channel_operation(rho, operators_tensor)
C = compute_channel_operation2(rho, operators_tensor)
np.allclose(A, B) # True
np.allclose(A, C) # True
%timeit compute_channel_operation(rho, operators)
# 6.95 ms ± 103 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit compute_channel_operation(rho, operators_tensor)
# 7.53 ms ± 141 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit compute_channel_operation2(rho, operators_tensor)
# 416 µs ± 12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [343]: X = np.array([[0, 1], [1, 0]], dtype=complex)
...: Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
...: Z = np.array([[1, 0], [0, -1]], dtype=complex)
In [344]: X
Out[344]:
array([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]])
In [345]: Y
Out[345]:
array([[ 0.+0.j, -0.-1.j],
[ 0.+1.j, 0.+0.j]])
In [346]: Z
Out[346]:
array([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, -1.+0.j]])
np.array 尝试创建一个多维数值数组，如果失败则创建一个对象 dtype 数组（或引发错误）。
In [347]: np.array([X,0])
Out[347]:
array([array([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]]), 0], dtype=object)
这个数组与列表 [X,0] 基本相同 - 它包含指向对象 X 和数字 0.
的指针或引用
但是给定 2 个（或更多）大小匹配的数组，它会生成更高维的数值数组。例如，具有复杂数据类型的 (2,2,2) 数组。
In [348]: np.array([X,Y])
Out[348]:
array([[[ 0.+0.j, 1.+0.j],
[ 1.+0.j, 0.+0.j]],
[[ 0.+0.j, -0.-1.j],
[ 0.+1.j, 0.+0.j]]])
block，或 concatenate/stack 的某种组合可以构成二维数组。例如一个 (2,4) 复数数组：
In [349]: np.block([X,Y])
Out[349]:
array([[ 0.+0.j, 1.+0.j, 0.+0.j, -0.-1.j],
[ 1.+0.j, 0.+0.j, 0.+1.j, 0.+0.j]])
要创建包含 X 和 Y` 的对象 dtype 数组，我可以这样做：
In [356]: xy = np.empty((2,), object)
In [357]: xy[0]= X
In [358]: xy[1]= Y
In [359]: xy
Out[359]:
array([array([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]]),
array([[ 0.+0.j, -0.-1.j],
[ 0.+1.j, 0.+0.j]])], dtype=object)
此 empty 后接单独赋值是构建对象 dtype 数组的最可靠方法。它绕过 Out[348].
中显示的多维结果
我不知道这些方法是否有助于您实现计算目标。我还没有充分研究你的描述。但请记住，快速 numpy 代码适用于多维数值（包括 complex）数组（例如 Out[348]）。使用对象 dtype 数组的数学是 hit-or-miss，而且几乎总是较慢。
@ 矩阵乘法适用于 X、Y、Out[348]、rho 等，但不适用于 Out[347] 或 RHO。 + 使用对象 dtype 数组，前提是元素本身支持加法。
假设我们有 3 个 Pauli matrices，每个维度为 (2x2)。如下图：
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
现在，如果我将这些每个单独的 (2x2) 矩阵作为另一个 (2x2) 矩阵的条目。说：
A = np.array([[X, 0], [0, Y]])
B = np.array([[X, X], [Y, Y]])
奇怪的是，A 的亮度为 (2x2) - 这正是我想要的 - B 的亮度为 (2, 2, 2, 2) 不管这是什么，如下所示
A = np.array([[X, 0], [0, Y]])
A.shape
Out[131]: (2, 2)
B = np.array([[X, X], [Y, Y]])
B.shape
Out[133]: (2, 2, 2, 2)
另一方面，设 C 为 (1x3) 矩阵，D 为 (1x2) 矩阵，例如
C = np.array([[X, Y, 0]])
D = np.array([[X, Y]])
再次查看初始化矩阵的维度
C = np.array([[X, Y, 0]])
C.shape
Out[136]: (1, 3)
D = np.array([[X, Y]])
D.shape
Out[138]: (1, 2, 2, 2)
所以似乎每当我在这样的数组中初始化数组时，如果有混合数据类型作为条目，即矩阵和整数，如 A 和 C，它给了我合理的形状我想要的，即维度 (2,2)，每个条目的 "hidden" 维度 (2x2)。但是一旦条目只是像 B 和 D 中那样的严格矩阵，它就会给我不敏感的维度，例如 (2, 2, 2, 2)。所以我的问题是：
如何使用严格的 (2, 2) 矩阵作为条目初始化 (n, n) numpy 数组（矩阵），并仍然保留其 (n, n) 维度，即而不是给出我有些奇怪的 numpy 维度 (w, x, y, z)?
我想要这个的原因是因为我正在使用量子力学中的运算符进行计算，使用这些 Pauli 矩阵，例如 X、Y 和 Z，作为量子计算中的量子门。所以如果我有一些状态 rho 这也是一个 (2x2) 矩阵。让
rho = np.array([[1, 0],
[0, 0]])
并设 RHO 为 (2x2) 对角矩阵，其条目为 (2x2) rho 矩阵。
RHO = np.array([[rho, 0],
[0, rho]])
我希望计算类似 np.dot(D, RHO) 的东西，这样它可以给出
np.array([np.dot(X, rho), 0],
[0, np.dot(Y, rho)])
而且我在 python 上检查过，两个 (2x2) 矩阵的点积以 (2x2) 矩阵作为条目，其条目乘法也将全部是点积。
我上面讨论的所有内容的动机是我希望使用这些属性作为向量化我的算法的手段。目前我的算法的一个非常粗略的例子看起来像这样：
for i in (integer_1):
for j in (integer_2):
#do something that involves operations with sums of dot products of many matrices#
并将其矢量化，使其有可能成为
for i in (integer_1):
#do something that involves operations with multiples of matrices with dot product of matrices as its entries#
哪个可能有效，也可能无效！但我很好奇我的这种方法是否会产生加速。 我希望我已经很好地解释了我的问题。 提前致谢。
编辑(1)
I've added latex formatted maths so hopefully you can understand what I'm trying to do.
def compute_channel_operation(rho, operators):
"""
Given a quantum state's density function rho, the effect of the
channel on this state is:
rho -> sum_{i=1}^n E_i * rho * E_i^dagger
Args:
rho (2x2 matrix): Density function
operators (list): List of operators(matrices)
Returns:
number: The result of applying the list of operators
"""
operation = [E@rho@E.conj().T for i, E in enumerate(operators)]
return np.sum(operation, axis=0)
So then I was hoping that instead of using loops, this direct multiplication of tensor method might be quicker as I scale up my simulation, say having to do this 1 million times 这里的 K 应该是任何 (1xn) 维度的列表，即 [I] 或 [I, X] 或 [I, X, Y] 或 [I, X, Y, Z]。我知道这里 X = X^{\dagger} Y 和 Z 也是如此，但我在模拟中会遇到不会出现这种情况的情况。
希望我已经解释清楚了。
(2, 2, 2, 2) 不是 奇怪的维度 ，它只是一个形状为 2x2x2x2
您看到 A 和 B 不同形状的原因是因为您使用标量 0 而不是 2x2 零矩阵来设置 A。改成
A = np.array([[X, np.zeros((2, 2))], [np.zeros((2, 2)), Y]])
B = np.array([[X, X], [Y, Y]])
你将得到 2x2x2x2 张量。
或改为
C = np.vstack([
np.hstack([X, np.zeros((2, 2))]),
np.hstack([np.zeros((2, 2)), Y])
])
D = np.vstack([
np.hstack([X, X]),
np.hstack([Y, Y])
])
你会得到 4x4 矩阵。
您还可以使用
从一种形式转换为另一种形式E = A.transpose(0, 2, 1, 3).reshape(4, 4)
F = B.transpose(0, 2, 1, 3).reshape(4, 4)
np.allclose(C, E) # True
np.allclose(D, F) # True
来回
G = E.reshape(2, 2, 2, 2).transpose(0, 2, 1, 3)
H = F.reshape(2, 2, 2, 2).transpose(0, 2, 1, 3)
np.allclose(A, G) # True
np.allclose(B, H) # True
编辑： 关于你的函数 compute_channel_operation()，如果你不执行列表推导而是向量化操作并传入 3D，你可以大大加快它的速度您所有操作的张量
rho = np.random.rand(2, 2)
operators = [np.random.rand(2, 2) for _ in range(1000)]
operators_tensor = np.asarray(operators) # same as np.random.rand(1000, 2, 2)
def compute_channel_operation(rho, operators):
operation = [E@rho@E.conj().T for i, E in enumerate(operators)]
return np.sum(operation, axis=0)
def compute_channel_operation2(rho, operators):
return np.sum(operators @ rho @ operators.transpose(0, 2, 1).conj(), axis=0)
A = compute_channel_operation(rho, operators)
B = compute_channel_operation(rho, operators_tensor)
C = compute_channel_operation2(rho, operators_tensor)
np.allclose(A, B) # True
np.allclose(A, C) # True
%timeit compute_channel_operation(rho, operators)
# 6.95 ms ± 103 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit compute_channel_operation(rho, operators_tensor)
# 7.53 ms ± 141 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit compute_channel_operation2(rho, operators_tensor)
# 416 µs ± 12 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [343]: X = np.array([[0, 1], [1, 0]], dtype=complex)
...: Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
...: Z = np.array([[1, 0], [0, -1]], dtype=complex)
In [344]: X
Out[344]:
array([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]])
In [345]: Y
Out[345]:
array([[ 0.+0.j, -0.-1.j],
[ 0.+1.j, 0.+0.j]])
In [346]: Z
Out[346]:
array([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, -1.+0.j]])
np.array 尝试创建一个多维数值数组，如果失败则创建一个对象 dtype 数组（或引发错误）。
In [347]: np.array([X,0])
Out[347]:
array([array([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]]), 0], dtype=object)
这个数组与列表 [X,0] 基本相同 - 它包含指向对象 X 和数字 0.
但是给定 2 个（或更多）大小匹配的数组，它会生成更高维的数值数组。例如，具有复杂数据类型的 (2,2,2) 数组。
In [348]: np.array([X,Y])
Out[348]:
array([[[ 0.+0.j, 1.+0.j],
[ 1.+0.j, 0.+0.j]],
[[ 0.+0.j, -0.-1.j],
[ 0.+1.j, 0.+0.j]]])
block，或 concatenate/stack 的某种组合可以构成二维数组。例如一个 (2,4) 复数数组：
In [349]: np.block([X,Y])
Out[349]:
array([[ 0.+0.j, 1.+0.j, 0.+0.j, -0.-1.j],
[ 1.+0.j, 0.+0.j, 0.+1.j, 0.+0.j]])
要创建包含 X 和 Y` 的对象 dtype 数组，我可以这样做：
In [356]: xy = np.empty((2,), object)
In [357]: xy[0]= X
In [358]: xy[1]= Y
In [359]: xy
Out[359]:
array([array([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]]),
array([[ 0.+0.j, -0.-1.j],
[ 0.+1.j, 0.+0.j]])], dtype=object)
此 empty 后接单独赋值是构建对象 dtype 数组的最可靠方法。它绕过 Out[348].
我不知道这些方法是否有助于您实现计算目标。我还没有充分研究你的描述。但请记住，快速 numpy 代码适用于多维数值（包括 complex）数组（例如 Out[348]）。使用对象 dtype 数组的数学是 hit-or-miss，而且几乎总是较慢。
@ 矩阵乘法适用于 X、Y、Out[348]、rho 等，但不适用于 Out[347] 或 RHO。 + 使用对象 dtype 数组，前提是元素本身支持加法。