Approximation of the multiplicative matrix inverse (linalg.inv ())

I am trying to compute the multiplicative inverse of a large matrix (~ >40,000x40,000). This can be done with e.g. numpy.linalg.inv or scipy.linalg.inv. Unfortunately, the calculation fails on the HPC to which I have access.

numpy.linalg.inv(I-A) or scipy.linalg.inv(I-A) is equivalent to (I-A)^-1 which can be approximated via I + A + AA + AAA + AAAA ... or I + A + A^2 + A^3 + A^4 ....

As I cannot run linalg.inv I wrote a function to approximate it. However, running my function with the full matrix proved to be extremely slow. I'm therefore wondering whether my code is simply inefficient/flawed.

import numpy as np
import scipy.sparse
from scipy import linalg

# Transactions
T = scipy.sparse.csr_matrix(
np.array([
[8, 5],
[4, 2]
])
)

# Total output
x = np.array(
[16,12]
)

# Technical coefficients
A = scipy.sparse.csr_matrix(T / x)

def getL(
A, # Technical coefficient matrix
log = False,
iterations = 25 # Iterations
):
"""
Approximate (I - A)^-1.

This function is an alternative to the execution of np.linalg.inv(I - A).

L = (I-A)^-1
L = (I-A)^-1 ≈ I + A + AA + AAA + AAAA ... ≈ I + A + A^2 + A^3 + A^4 ...

Parameters
----------
A : Matrix (e.g. np.array or scipy.sparse.csr_matrix)
Technical coefficient matrix.
log : Boolean, optional
Print production layer sum to file. The default is False.
iterations: Integer, optional
Number of iterations.

Returns
-------
L : Matrix (e.g. np.array or scipy.sparse.csr_matrix)
Leontief inverse. Approximation of (I - A)^-1.

"""

# Zeroth production layer
I = scipy.sparse.identity(
A.shape[0],
format = 'csc'
)

L = I.copy()

# First production layer
layer = A.copy()

# ... add the ensuing production layer until
# L_i.sum() is less than 0.001% of L.sum()
i = 0
while i <= iterations:

L += layer
layer = layer @ A

if log: print(f"... layer {i} ...")

i += 1

# Log
if log:

print(f"getL coverage [%]: {1-(layer.sum()/L.sum())}")

return L

# Calculate L

L_getL = getL(A, iterations = 50, log = True)

L_SciPy = scipy.linalg.inv(
np.identity(A.shape[0])
- A.todense()
)

# Compare the approaches

print(L_getL.todense())
print(L_SciPy)
• What is known about $A$? You should realize that for an arbitrary $A$ this process may not converge at all.
– vnp
Jun 17 at 20:17

Guess number one: you've started with a sparse matrix (good), but then you're failing to call into the sparse matrix linalg methods (almost certainly ungood).

Also avoid a dense T / x; you can stay sparse by multiplying with a sparse diagonal matrix from x.

Try:

import numpy as np
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg

T = scipy.sparse.csc_matrix(  # Transactions
np.array([
[8, 5],
[4, 2],
])
)
x = np.array([16, 12])           # Total output
A = T * scipy.sparse.diags(1/x)  # Technical coefficients

L_dense = scipy.linalg.inv(
np.eye(*A.shape) - A.todense()
)
L_sparse = scipy.sparse.linalg.inv(
scipy.sparse.eye(*A.shape, format='csc') - A
)

print(L_dense)
print(L_sparse.todense())
[[2.66666667 1.33333333]
[0.8        1.6       ]]
[[2.66666667 1.33333333]
[0.8        1.6       ]]

Note the use of CSC format instead of your CSR; SciPy asks for this during the inverse for efficiency reasons.