numpy - 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 分解进行反转时，我得到了类似的答案。当我使用 dgesvd_ 使用 SVD 进行逆时，我也得到了类似的答案。

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

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

Answer 1

看来numpy.linalg.inv是 scipy.linalg.inv 的精简版根据描述:

This module is a lite version of the linalg.py module in SciPy which contains high-level Python interface to the LAPACK library.

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