3

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?

PiMathCLanguage
  • 363
  • 4
  • 15
  • 1
    how large is the relative error? – Paul Panzer Feb 06 '17 at 17:47
  • Do you mean the possible relative error for my problem or the relative error between the two matrices? – PiMathCLanguage Feb 06 '17 at 17:49
  • I mean the error you computed divided by the squared Euclidean length of one of the two inverses you compare. – Paul Panzer Feb 06 '17 at 17:52
  • 2
    `np.allclose` is a convenient tool for checking float equality. – hpaulj Feb 06 '17 at 17:52
  • For 10x10 matrices I get, e.g. 100 times, for np.allclose always True. For 150x150 sometimes True sometimes False. For bigger one, like 200x200, it will return probably more False then True. – PiMathCLanguage Feb 06 '17 at 18:06
  • 1
    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 Feb 06 '17 at 18:21
  • I have uploaded a picture of a matrix, https://pl.vc/1crz9p. In this picture the blue color represent the value -2 and the shades of red are +4, +6, +8. You can see, its a very diagonal matrix. Also when I calculate the inverse with numpy, there are no problems. The matrix is inversible. – PiMathCLanguage Feb 06 '17 at 18:50
  • Are you happy with these explanations? – Paul Panzer Feb 06 '17 at 19:03
  • Your answer make for me sense (and makes me happy). However I cannot really understand, why the numpy inversion is better or more precise than the from scipy (which for me doesn't really make sense). – PiMathCLanguage Feb 06 '17 at 19:39
  • How do you know which is the more precise? Did you multiply them back? Like A A^-1 or A^-1 A? – Paul Panzer Feb 06 '17 at 19:46

1 Answers1

2

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.

Paul Panzer
  • 51,835
  • 3
  • 54
  • 99