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)