出于好奇，我尝试了以下方法：

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)

我只有0.2秒。
我的PC配置为
内存-32千兆字节
Win10固态硬盘
英特尔酷睿i7

如果你试图找到一个矩阵方程Ax = b左手，求逆x=A^-1 b几乎永远不是解决问题的方法。在大多数数值系统中，你需要求解方程的左边。在python中，你可以使用scipy.linalg.solve来完成这个任务。（参见https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html，在Python: Is there a matlab-like backslash operator?中引用）
对于一个1000x1000的随机、密集矩阵，大约有2倍的加速，但是如果你有对称性可以利用，就像物理系统经常做的那样，你可以得到更快的速度。

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