5

我想要一种 numpy-sh 向量化特征值计算的方法，这样我就可以给它一个矩阵矩阵，它会返回一个各自特征值的矩阵。

例如，在下面的代码中，B 是由 3x3 矩阵 A 的 4 个副本组成的块 6x6 矩阵。C 是我希望看到的输出，即维度为 (2,2,3) 的数组（因为 A有 3 个特征值）。

这当然是一个非常简化的例子，一般情况下矩阵 A 可以有任意大小（尽管它们仍然是正方形），而矩阵 B 不一定由 A 的副本组成，而是由不同的 A1、A2 等组成（所有大小相同但包含不同元素）。

import numpy as np
A = np.array([[0, 1, 0],
              [0, 2, 0],
              [0, 0, 3]])
B = np.bmat([[A, A], [A,A]])
C = np.array([[np.linalg.eigvals(B[0:3,0:3]),np.linalg.eigvals(B[0:3,3:6])],
              [np.linalg.eigvals(B[3:6,0:3]),np.linalg.eigvals(B[3:6,3:6])]])
4

2 回答 2

4

编辑：如果您使用的是 numpy >= 1.8.0 的版本，则np.linalg.eigvals对您提供的任何数组的最后两个维度进行操作，因此如果您将输入重塑为(n_subarrays, nrows, ncols)数组，则只需调用eigvals一次：

import numpy as np

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

# the input needs to be an array, since matrices can only be 2D.
B = np.repeat(A[np.newaxis,...], 4, 0)

# for arbitrary input arrays you could do something like:
# B = np.vstack(a[np.newaxis,...] for a in input_arrays)
# but for this to work it will be necessary for each element in 
# 'input_arrays' to have the same shape

# eigvals will operate over the last two dimensions of the array and return
# a (4, 3) array of eigenvalues
C = np.linalg.eigvals(B)

# reshape this output so that it matches your original example
C.shape = (2, 2, 3)

如果您的输入数组并非都具有相同的维度，例如input_arrays[0].shape == (2, 2)等，input_arrays[1].shape == (3, 3)那么您只能在具有匹配维度的子集之间对该计算进行矢量化。

如果您使用的是旧版本的 numpy，那么不幸的是，我认为没有任何方法可以对多个输入数组上的特征值计算进行矢量化 - 您只需要在 Python 中循环输入即可。

于 2013-10-19T19:03:51.243 回答
1

你可以做这样的事情

C = np.array([[np.linalg.eigvals(B[i:i+3, j:j+3])
              for i in xrange(0, B.shape[0], 3)]
                for j in xrange(0, B.shape[1], 3)])

也许更好的方法是使用https://stackoverflow.com/a/5078155/1352250block_view中的函数：

B_blocks = block_view(B)
C = np.array([[np.linalg.eigvals(m) for m in v] for v in B_blocks])

更新

正如 ali_m 所指出的，这种方法是一种语法糖，不会减少eigvals大量调用所产生的开销。虽然如果应用它的每个矩阵都很大，那么这个开销应该很小，但对于 OP 感兴趣的 6x6 矩阵，它并不是微不足道的（参见下面的评论；根据 ali_m，可能有一个因子三我上面给出的版本和他发布的使用 Numpy >= 1.8.0 的版本之间的区别）。

于 2013-10-19T20:16:48.920 回答