不用for循环将矩阵乘以标量数组的Numpy方法

Question

假设我有以下代码：

# define a 3x3 array
A = np.array([[-2, 1, 1], [1,-2,1], [1,1,-2]])
# define a 1D array of scalars to multiply A with
test = np.arange(10)

for i in range(len(test)):
    matrix = scipy.linalg.expm(A*test[i])

我想知道是否有一种不使用for循环的方法来完成这个乘法。我并不是想使用矩阵乘法来计算两个数组的乘积。我将测试数组看作是一组我想一个一个地与A相乘的标量值。必须有某种聪明的NumPy方式可以做到这一点。有什么建议吗？

- sjhaque14

---

从这些数据中制作一个（10,3,3）或（3,3,10）数组很容易。但是，expm接受什么？ - hpaulj

Answer 1

这是一个关于如何在不使用for循环的情况下通过多个标量来缩放矩阵的答案。由于我认为没有办法在不使用for循环的情况下执行矩阵指数，因此本文不考虑矩阵指数，但问题似乎并不涉及指数步骤。

您可以使用3D数组来使用numpy的向量化乘法。

import numpy as np

A = np.array([[[-2, 1, 1], [1,-2,1], [1,1,-2]]])
test = np.arange(10).reshape(-1,1,1)

result = test*A

- abatea

---

是的，expm 只能用于 (m,m) 矩阵。它使用某种版本的 Pade 近似方法，不适合进行“批量”计算。https://en.wikipedia.org/wiki/Matrix_exponential#Computing_the_matrix_exponential - hpaulj

我喜欢使用A * np.arange(10)[:, None, None]的表示法，但结果是相同的。使用None/np.newaxis来扩展数组的维度以进行广播是一种很好的视觉习惯。 - hpaulj

这个方法对于一般问题是有效的，但正如其他人所说，我确实需要一个可以通过expm函数传递的东西。但还是谢谢你的解决方案！ - sjhaque14

Answer 2

scipy.linalg.expm 无法应用于非方阵矩阵：

ValueError: expected a square matrix

我认为最简单的做法是：

list_result = list(map(lambda i: scipy.linalg.expm(A*i), range(10)))

然后如果你只想要一个数组：

np.concatenate(list_result) #for 2D

或者

np.stack(list_result) #for 3D

- ansev

Answer 3

你可以通过堆叠由所需值乘以单位矩阵得到的特殊矩阵来完成技巧。

我将以3x3矩阵为例，但对于任意矩阵，思路都是相同的：

假设我们想要乘以矩阵：

A = np.array([
   [1,0,1],
   [0,1,0],
   [0,0,1] 
])

通过 [2、3 和 5]。不使用循环。

我们可以构建一个特殊的矩阵来连接身份以实现此目的：

[ 2 * np.eye(3), 3 * np.eye(3), 5 * np.eye(3)]

[ 0 0 2 ] [ 0 0 3] [ 0 0 5 ]
[ 0 2 0 ] [ 0 3 0] [ 0 5 0 ]
[ 2 0 0 ] [ 3 0 0] [ 5 0 0 ]

我们可以使用列表推导来实现这个功能：

multiply_values = [2,3,5]

special_matrix = np.concatenate( [ x * eye(3) for x in multiply_values ], axis=1)

# special_matrix is:
# array([[2., 0., 0., 3., 0., 0., 5., 0., 0.],
#        [0., 2., 0., 0., 3., 0., 0., 5., 0.],
#        [0., 0., 2., 0., 0., 3., 0., 0., 5.]])

现在我们可以一步完成矩阵相乘：

result = np.dot( , special_matrix)

# result is
array([[2., 0., 0., 3., 0., 0., 5., 0., 0.],
       [0., 2., 0., 0., 3., 0., 0., 5., 0.],
       [0., 0., 2., 0., 0., 3., 0., 0., 5.]])

可能这个大宽矩阵不是我们需要的。

我们可以通过切片结果来获取部分结果：

result[:,0:3]

# array([[2., 0., 2.],
#       [0., 2., 0.],
#       [0., 0., 2.]])

我们可以使用另一种列表推导式来处理切片：

[ special_matrix[:,x:x+3] for x in [0,3,6] ]

# or in a more general way

[ special_matrix[:,x:x+A.shape[0]] for x in list(range(0,A.shape[0]**2,A.shape[0] ]

- Iñigo González

---

通过广播和逐元素乘法，就不需要使用 dot 和连接的 eyes。 - hpaulj