0

In my code, each runtime takes about 40sec, and most of the time (38sec) is spent matrix inversion, i.e., numpy.linalg.inv. Here, the matrix dimension generally is 100-1000.

Thus, I want to reduce the runtime. Are there some effective methods to solve it.

Young Q
  • 31
  • 2
  • It is used to solve some specific physical problems, and the matrix mainly depends on the Hamiltonian of system. – Young Q Mar 10 '23 at 04:34
  • That's for multiple inversions? – hpaulj Mar 10 '23 at 05:21
  • Add more details like code, your system configuration etc else the question might be downvoted due to lack of information. – balu Mar 10 '23 at 07:46
  • If it's a dense matrix, numpy.linalg.inv is_probably_ about as good as you're going to get. (maybe you could get some speedup by being sure that numpy is compiled to call BLAS. See https://stackoverflow.com/questions/9000164/how-to-check-blas-lapack-linkage-in-numpy-and-scipy). In most of numerical systems you don't invert the matrix (which is expensive) you find some operation that's equivalent to inversion which is cheaper. – Elliot Mar 10 '23 at 13:44

2 Answers2

0

Out of curiosity I tried the following :

import  numpy as np
from numpy.linalg import inv
from timeit import timeit


def compute_inverse(a):
    return  inv(a)

a = np.random.rand(1000,1000)

loop = 100

result = timeit('compute_inverse(a)', globals=globals(), number=loop)
print(result / loop)

I got ~0.2 seconds.

My PC config is

RAM - 32GiB

Win10 SSD

Intel Core i7

balu
  • 1,023
  • 12
  • 18
0

If you're trying to find the left hand side of a matrix equation Ax = b, the inversion x=A^-1 b is almost never the way to go about it. In most numerical systems you solve the equation for the left hand side. In python you can use scipy.linalg.solve to do this. (See https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html, referenced in Python: Is there a matlab-like backslash operator?)

For a 1000x1000 random, dense matrix there's about a 2x speedup, but if you have symmetries to exploit, as physical systems often do, you can get that faster.

import numpy as np
from scipy import linalg
from timeit import timeit

mat = np.random.random(1000000).reshape(1000,1000)

rhs = np.random.random(1000)

def lhs_byinversion(mat, rhs):
    matinv = linalg.inv(mat)
    lhs = matinv.dot(rhs)
    return lhs
    
def lhs_bysolving(mat, rhs):
    lhs = linalg.solve(mat, rhs)
    return lhs

lhs1 = lhs_byinversion(mat, rhs)
lhs2 = lhs_bysolving(mat, rhs)

# these are equivalent solutions
print('difference in results')
print(np.sqrt(np.sum((lhs1-lhs2)**2)))
# 1.157200161145036e-12

# speed up depends on the actual matrices being used. If there are symmetries to exploit this can be greater
print('time for 100 inversions')
print(timeit('lhs = lhs_byinversion(mat, rhs)', globals=globals(), number=100))
# 1.7345171420020051
print('time for 100 solutions')
print(timeit('lhs = lhs_bysolving(mat, rhs)', globals=globals(), number=100))
# 0.9372443140018731
Elliot
  • 2,541
  • 17
  • 27