Why the matrix inversion function in numpy and scipy returns different results with big quadratic matrices?

Question

Lets say I define a big quadratic matrix (e.g. 150x150). One time it is a numpy array (matrix A), one time it is a scipy sparse array (matrix B).

import numpy as np
import scipy as sp

from scipy.sparse.linalg import spsolve

size = 150
A = np.zeros((size, size))
length = 1000
# Set some random numbers at random places
A[np.random.randint(0, size, (2, length)).tolist()] = \
    np.random.randint(0, size, (length, ))
B = sp.sparse.csc_matrix(A)

Now I calculate the inverse of both matrices. For matrix B I use both methods for calculating the inverse (sp.sparse.linalg.inv and spsolve).

epsilon = 10.**-8 # Is needed to prevent singularity of each matrix

inv_A = np.linalg.pinv(A+np.eye(size)*epsilon)
inv_B = sp.sparse.linalg.inv(B+sp.sparse.identity(size)*epsilon)
inv_B2 = spsolve(B+sp.sparse.identity(size)*epsilon, sp.sparse.identity(size))

To check if both inverses of A and B are equal, I will sum up the differences squared.

# Is not equal zero, question: Why?
# Sometimes very small (~+-10**-27), sometimes very big (~+-10**5)
print("np.sum((inv_A - inv_B)**2): {}".format(np.sum((inv_A - inv_B)**2)))
# Should be zero
print("np.sum((inv_B - inv_B2)**2): {}".format(np.sum((inv_B - inv_B2)**2)))

The problem is this: If I use small matrices, e.g. 10x10, the error between numpy as scipy inverse functions is very small (approx. ~+-10**-32). But I need the sparse version for big matrices (e.g. 500x500).

Am I doing here something wrong, or is there any possibility to calc the correct inverse of a sparse matrix in python?

Do you mean the possible relative error for my problem or the relative error between the two matrices? — PiMathCLanguage
I mean the error you computed divided by the squared Euclidean length of one of the two inverses you compare. — Paul Panzer
np.allclose is a convenient tool for checking float equality. — hpaulj
Actually, one likely problem is that your matrices are near-singular by construction. Inverting an almost singular matrix makes it very difficult for the algorithm, I suppose. Typically, in this regime the small unavoidable errors of machine arithmetic accumulate and blow up more readily than with "normal" operands. Why don't you try a more well-behaved example? — Paul Panzer

Paul Panzer Paul Panzer · Accepted Answer · 2017-02-07T13:09:15

The answer to your headline question is: Because of your unfortunate choice of example matrices. Let me elaborate.

Machine precision is limited, therefore floating point arithmetic will rarely be 100% accurate. Just try

>>> np.linspace(0, 0.9, 10)[1:] == np.linspace(0.1, 1, 10)[:-1]
array([ True,  True,  True,  True,  True, False,  True,  True,  True], dtype=bool)

Normally, that's no problem because errors are too small to notice.

However, for many computations, there are some inputs that are difficult to handle and may overstretch a numerical algorithm. This certainly applies to matrix inversion, and you were unlucky enough to pick such difficult inputs.

You can actually check whether a matrix might be 'ill-conditioned' by looking at its singular values see for example here. Here are the matrix condition numbers for several matrices generated with your script (size=200; a well-behaved matrix has a value much closer to 1)

971899214237.0
5.0134186641e+12
36848.0807109
958492416768.0
1.66615247737e+16
1.42435766189e+12
1954.62614384
2.35259324603e+12
5.58292606978e+12

Switching to well-behaved matrices your results should drastically improve.

Why the matrix inversion function in numpy and scipy returns different results with big quadratic matrices?

1 Answers