Matrix Inversion in CBLAS/LAPACK vs Python

Question

The matrix I am trying to invert is:

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

The true inverse is:

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

Using Python's numpy.linalg.inv, I get the correct answer. One of my routines for matrix inverse uses dgetri_, it is:

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;
  }

Using this routine, I get:

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

Where e1 and e2 and small numbers on the order of machine precision 1e-16. Now perhaps dgetri_ is not the best to use. However, when I invert using QR decomposition via zgeqrf_ and zungqr_, I get a similar answer. When I use dgesvd_ for inverse using SVD, I get a similar answer as well.

Python seems to use a routine called _umath_linalg.inv. So I have a few questions:

What does that routine do?
What CBLAS/LAPACK routine can I use to invert this matrix and get a result like CBLAS/LAPACK (such that the e1 and e2 get replaced by proper zeros)?

You need to use ?GESV or ?GESVX for the inverses with full precision not with decomposition based solvers. LU is faster but as you see loses precision. Scipy linalg.solve has recently switched to ___x routines for condition number checks. — percusse
Also never ever use matrix inverse unless you need to refer to it. For matrix algebra always use solve. — percusse
@percusse I'm not sure what you mean by ____x algorithms. Could you elaborate? — The Dude

Pedro H. N. Vieira Pedro H. N. Vieira · Accepted Answer · 2018-10-30T16:54:28

It seems that numpy.linalg.inv is a lite version of the scipy.linalg.inv as per the description:

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

Looking at scipy.linalg.inv, it does a call to getrf, then getri.

Matrix Inversion in CBLAS/LAPACK vs Python

1 Answers