LAPACKE C++ Real Matrix Inversion

Question

I am trying to invert a real matrix in C++ LAPACKE. I have the same function for complex matrices and it works. But the real case gives the wrong answer. Here's my function:

void inv(std::vector<std::vector<double>> &ans, std::vector<std::vector<double>> MAT){

    int N = MAT.size();

    int *IPIV = new int[N];

    double * arr = new double[N*N];
    for (int i = 0; i<N; i++){
        for (int j = 0; j<N; j++){
            int idx = i*N + j;
            arr[idx] = MAT[i][j];
        }
    }

    LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);
    LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);

     for (int i = 0; i<N; i++){
        for (int j = 0; j<N; j++){
            int idx = i*N + j;
            ans[i][j] = arr[idx];
        }
    }

    delete[] IPIV;
    delete[] arr;
}

I tried inverting a 24 by 24 matrix of doubles. While the program seems to be almost there, the inverse is not quite there yet and it differs a lot from what python linalg inverse gives me (python is right here because I multiplied the matrix with inverse and result is very close to indentity). In the LAPACKE output, I multiply the matrix by its inverse and I get the diagonals being 1 but the non diagonals go to values up to 0.17, which is huge compared to 0. Is there a way to make the LAPACKE program to give a better result? Thanks!

Is it well-defined matrix? What is the value of matrix determinant? — Severin Pappadeux
2.2864066779666567e+37. It has large values so I think that's why the program is having trouble? — Imran Alkhatib
Yes, that might explain the effect. You could try to multiply input matrix (each term) by (1/25), which will decrease determinant by factor of 25^{-24}=2.8e-34 (if I'm not mistaken), thus making input matrix determinant about 1000. Then compute inverse and multiply back by 1/25. I've put simple python code in the answer — Severin Pappadeux
Taking a look at the source of numpy's linalg.inv(), in particular line 1555 - 1693 -1727 of umath_linalg.c, it solves A.X=Id, using dgesv() instead of calling dgetrf() and dgetri(). Scipy makes use of dgetrf() / dgetri(). Could you try to use scipy's scipy.linalg.inv() in python and see if it also result in a correct output? — francis

Severin Pappadeux Severin Pappadeux · Accepted Answer · 2020-01-24T21:30:57

For matrix with large determinant you could rescale input, compute inverse and then rescale output back. Here is very simple Python example, your scale factor should be 1/25 or so in order to get total multiplier of (1/25)²⁴=2.8e-34, thus making input matrix determinant about 1000.

import numpy as np

scale = 0.5

i = np.array([[1,2],[3,4]])
print(i)
print(np.linalg.det(i))
print("-----------------------------------")
x = np.multiply(i, [scale]) # rescale matrix
print(x)
print(np.linalg.det(x))     # determinant should be less
print("-----------------------------------")
y = np.linalg.inv(x)
print(y)
print(np.linalg.det(y))
print("-----------------------------------")
o = np.multiply(y, [scale]) # rescale matrix
print(o)
print(np.linalg.det(o))
print(np.dot(i, o))

You could play with the scale and verify that whatever value is code returns properly inverted input matrix.

So your code would be

double scale = 1.0/25.0;

for (int i = 0; i<N; i++){
    for (int j = 0; j<N; j++){
        int idx = i*N + j;
        arr[idx] = scale*MAT[i][j];
    }
}

....

for (int i = 0; i<N; i++){
    for (int j = 0; j<N; j++){
        int idx = i*N + j;
        ans[i][j] = scale * arr[idx];
    }
}

LAPACKE C++ Real Matrix Inversion

1 Answers