CODEWITHSUNDEEP
pythonnumpymatrixscipy
Light Yagmi
Light Yagmi

Reputation: 5235

Generate non-singular sparse matrix in Python

When a sparse matrix is generated by scipy.sparse.rand, it can be singular. In fact,the below code raises an error "RuntimeError: superlu failure (singular matrix?) at line 100 in file scipy/sparse/linalg/dsolve/SuperLU/SRC/dsnode_bmod.c".

dim = 20000
ratio = 0.000133

A = scipy.sparse.rand(dim,dim,ratio)
inv_sparse = scipy.sparse.linalg.inv(A)

Is there a way to generate non-singular sparse matrix?

What I really want to do is comparing performance (process time) of scipy.sparse.linalg.inv with np.linalg.inv. That's why I need generate random sparse matrix which is not singular.

Upvotes: 3

Views: 2274

Answers (2)

Nikolay Frick
Nikolay Frick

Reputation: 2205

As was already suggested in the comments by @hpaulj, you can add an identity matrix to your random matrix to make it invertible.

However, in applications where you need to find approximate solutions to systems of linear equations, just multiply identity matrix by a small factor:

dim = 20000
ratio = 0.000133
A = csc_matrix(scipy.sparse.random(dim, dim, density=ratio), dtype=float)
A = A + scipy.sparse.csr_matrix(np.eye(A.shape[0]) * 1e-8)
A = scipy.sparse.csc_matrix(A)
inv_sparse = scipy.sparse.linalg.inv(A)

Upvotes: 1

francis
francis

Reputation: 9817

The density ratio = 0.000133 of your matrix is very low. It means that about one item out of 7518 is non-null. Hence, the probability of each term to be null is about 7517/7518.

Each row is made of 20000 independent terms. So the probability for a row to be null is (7517/7518)^20000=6.99%. Hence, the probability for a row to be non-null is 1-(7517/7518)^20000=93.0%.

Then, the matrix is made of 20000 rows. The rows can be considered as being independent. Hence, the probability that the matrix does not contain null rows is (1-(7517/7518)^20000)^20000=(93.0%)^20000. This probability is very low.

As the matrix is likely to contain a null row, it is often singular.

Moreover, due to the the limited precision of floating-point numbers, programs often consider ill-conditionned matrices as singular. Indeed, in such cases, the computed inverse would be highly unprecise and meaningless.

Finally, to compare the inverse function, it may be better to use matrices known to be invertible... At least, you could try to increase the density so that the probability of a null row becomes very low.

Upvotes: 2

Related Questions