首页 > 解决方案 > 将 scipy.linalg.eigh 等矩阵函数应用于高维数组

问题描述

我是 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 次,因此需要提高效率。任何帮助将不胜感激?

标签: pythonarraysnumpyscipy

解决方案


您可以使用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]]))

推荐阅读