CBLAS/LAPACK 与 Python 中的矩阵反演

Question

我试图反转的矩阵是：

[ 1  0  1]
A = [ 2  0  1]
[-1  1  1] 

真正的反面是：

[-1  1  0]
A^-1 = [-3  2  1]
[ 2 -1  0]

使用Python的numpy.linalg.inv，我得到了正确的答案。我的矩阵逆用法之一dgetri_，它是：

void compute_matrix_inverse_dbl( double* matrix,
int order,
double * inverse )
{
int N, lwork;
int success;
int *pivot;
double* workspace;
//===Allocate Space===//
pivot = malloc(order * order * order * sizeof(*pivot));
workspace = malloc(order * order * sizeof(*workspace));
//===Run Setup===//
N = order;
copy_array_dbl(matrix, order*order, inverse);
lwork = order*order;
//===Factor Matrix===//
dgetrf_(&N,&N,inverse,&N,pivot,&success);
//===Compute Inverse===//
dgetri_(&N, inverse, &N, pivot, workspace, &lwork, &success);
//===Clean Up===//
free(workspace);
free(pivot);
return;
}

使用此例程，我得到：

[-1   1 +-e1 ]
A^-1 = [-3   2   1  ]
[ 2  -1 +-e2 ]

其中 e1 和 e2 以及机器精度为 1e-16 的小数字。现在也许dgetri_不是最好的使用。但是，当我通过zgeqrf_和zungqr_使用QR分解进行反转时，我得到了类似的答案。当我使用 SVD 将dgesvd_用于反转时，我也得到了类似的答案。

Python似乎使用了一个名为_umath_linalg.inv的例程。所以我有几个问题：

• 这个例程有什么作用？
• 我可以使用什么 CBLAS/LAPACK 例程来反转此矩阵并获得类似于 CBLAS/LAPACK 的结果(这样 e1 和 e2 被正确的零替换(？

Answer 1

根据描述，numpy.linalg.inv似乎是scipy.linalg.inv的精简版：

此模块是 SciPy 中 linalg.py 模块的精简版，它 包含到 LAPACK 库的高级 Python 接口。

查看scipy.linalg.inv，它调用getrf，然后getri。