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被合适的零所代替）？

- The Dude

---

你需要使用?GESV或?GESVX来进行完全精度的逆矩阵计算，而不是基于分解的求解器。LU分解速度更快，但会失去精度。Scipy linalg.solve最近已经切换到___x例程以进行条件数检查。 - percusse

永远不要使用矩阵求逆，除非你需要引用它。对于矩阵代数，请始终使用solve函数。 - percusse

@percusse 我不确定你所说的____x算法是什么意思。能否详细解释一下？ - The Dude

我指的是GESVX、POSVX、HESVX等等。 - percusse

Answer 1

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

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

查看scipy.linalg.inv, 它调用了getrf和getri函数。

- Pedro H. N. Vieira