将 scipy.linalg.eigh 等矩阵函数应用于高维数组

Question

我是 numpy 的新手，但作为一名工程师，我使用 python 已经有一段时间了。 我正在编写一个程序，该程序当前将应力张量存储为另一个 NxM 数组中的 3x3 numpy 数组，该数组表示通过时间和墙壁厚度的值，所以总的来说它是一个 NxMx3x3 numpy 数组。我想有效地计算这个更大数组中每个 3x3 数组的特征值和向量。到目前为止，我已经尝试使用“fromiter”，但这似乎不起作用，因为函数返回 2 个数组。我也试过 apply_along_axis 也不起作用，因为它说内部 3x3 不是方阵？我可以通过列表理解来做到这一点，但这似乎并不适合使用列表。

使用列表推导计算特征值的示例

import numpy as np
from scipy import linalg
a=np.random.random((2,2,3,3))
f=linalg.eigvalsh
ans=np.asarray([f(x) for x in a.reshape((4,3,3))])
ans.shape=(2,2,3)

我认为这样的东西会起作用，但我已经玩过它并且无法让它起作用：

np.apply_along_axis(f,0,a)

顺便说一句，2x2 位可能高达 5000x100，并且此代码重复了 ~50x50x200 次，因此需要提高效率。任何帮助将不胜感激？

Answer 1

您可以使用numpy.linalg.eigh。它接受像您的示例 a 这样的数组。

这是一个例子。首先，创建一个 3x3 对称数组：

In [96]: a = np.random.random((2, 2, 3, 3))

In [97]: a = a + np.transpose(a, axes=(0, 1, 3, 2))

In [98]: a[0, 0]
Out[98]: 
array([[0.61145048, 0.85209618, 0.03909677],
       [0.85209618, 1.79309413, 1.61209077],
       [0.03909677, 1.61209077, 1.55432465]])

计算所有 3x3 数组的特征值和特征向量：

In [99]: evals, evecs = np.linalg.eigh(a)

In [100]: evals.shape
Out[100]: (2, 2, 3)

In [101]: evecs.shape
Out[101]: (2, 2, 3, 3)

看看a[0, 0]的结果：

In [102]: evals[0, 0]
Out[102]: array([-0.31729364,  0.83148477,  3.44467813])

In [103]: evecs[0, 0]
Out[103]: 
array([[-0.55911658,  0.79634401,  0.23070516],
       [ 0.63392772,  0.23128064,  0.73800062],
       [-0.53434473, -0.55887877,  0.63413738]])

验证是否与分别计算a[0, 0]的特征值和特征向量相同：

In [104]: np.linalg.eigh(a[0, 0])
Out[104]: 
(array([-0.31729364,  0.83148477,  3.44467813]),
 array([[-0.55911658,  0.79634401,  0.23070516],
        [ 0.63392772,  0.23128064,  0.73800062],
        [-0.53434473, -0.55887877,  0.63413738]]))