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, Jul 28 '17 at 22:42
Also never ever use matrix inverse unless you need to refer to it. For matrix algebra always use solve. — percusse, Jul 28 '17 at 22:49
@percusse I'm not sure what you mean by ____x algorithms. Could you elaborate? — The Dude, Jul 29 '17 at 21:15

score 1 · Answer 1 · answered Oct 30 '18 at 16:54

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 Answers1