0
\$\begingroup\$

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)
\$\endgroup\$
1
  • 1
    \$\begingroup\$ What is known about \$A\$? You should realize that for an arbitrary \$A\$ this process may not converge at all. \$\endgroup\$
    – vnp
    Jun 17 at 20:17

1 Answer 1

1
\$\begingroup\$

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.

\$\endgroup\$
4
  • \$\begingroup\$ Fantastic. Thank your for your helpful answer. I'll try it out :) \$\endgroup\$
    – Stücke
    Jun 18 at 20:07
  • \$\begingroup\$ Why's the * in np.eye(*A.shape)? \$\endgroup\$
    – Stücke
    Jun 21 at 9:52
  • 1
    \$\begingroup\$ Tuple unpacking to m, n \$\endgroup\$
    – Reinderien
    Jun 21 at 9:58
  • \$\begingroup\$ Thank you for your swift response! \$\endgroup\$
    – Stücke
    Jun 21 at 10:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.