为了说明我关于广义ufuncs的观点，以及消除egg和ham上的循环的能力，请考虑以下代码：

import numpy as np
A = np.random.randn(4,4,10,10)
AI = np.linalg.inv(A)
#check that generalized ufuncs work as expected
I = np.einsum('xyij,xyjk->xyik', A, AI)
print np.allclose(I, np.eye(10)[np.newaxis, np.newaxis, :, :])

字符串

@yevgeniy你是对的，用scipy有效地解决多个独立的线性系统A x = B有点棘手（假设A数组在每次迭代中都发生变化）。
例如，这里有一个用于求解1000个形式为Ax = b的系统的基准，其中A是10x10矩阵，b是10元素向量。令人惊讶的是，将所有这些放入一个块对角矩阵并调用scipy.linalg.solve一次的方法确实在密集和稀疏矩阵中都较慢。

import numpy as np
from scipy.linalg import block_diag, solve
from scipy.sparse import block_diag as sp_block_diag
from scipy.sparse.linalg import spsolve

N = 10
M = 1000 # number of coordinates 
Ai = np.random.randn(N, N) # we can compute the inverse here,
# but let's assume that Ai are different matrices in the for loop loop
bi = np.random.randn(N)

%timeit [solve(Ai, bi) for el in range(M)]
# 10 loops, best of 3: 32.1 ms per loop

Afull = sp_block_diag([Ai]*M, format='csr')
bfull = np.tile(bi, M)

%timeit Afull = sp_block_diag([Ai]*M, format='csr')
%timeit spsolve(Afull, bfull)

# 1 loops, best of 3: 303 ms per loop
# 100 loops, best of 3: 5.55 ms per loop

Afull = block_diag(*[Ai]*M) 

%timeit Afull = block_diag(*[Ai]*M)
%timeit solve(Afull, bfull)

# 100 loops, best of 3: 14.1 ms per loop
# 1 loops, best of 3: 23.6 s per loop

字符串
具有稀疏阵列的线性系统的解更快，但是创建该块对角阵列的时间实际上非常慢。至于密集阵列，在这种情况下它们只是速度较慢（并且需要大量RAM）。
也许我遗漏了一些关于如何有效地使用稀疏数组的内容，但是如果你保持for循环，你可以做两件事来优化。
从纯python看scipy.linalg.solve的源代码：删除不必要的测试并分解循环外的所有重复操作。例如，假设你的数组不是对称的，我们可以做

from scipy.linalg import get_lapack_funcs

gesv, = get_lapack_funcs(('gesv',), (Ai, bi))

def solve_opt(A, b, gesv=gesv):
    # not sure if copying A and B is necessary, but just in case (faster if arrays are not copied)
    lu, piv, x, info = gesv(A.copy(), b.copy(), overwrite_a=False, overwrite_b=False)
    if info == 0:
        return x
    if info > 0:
        raise LinAlgError("singular matrix")
    raise ValueError('illegal value in %d-th argument of internal gesv|posv' % -info)

%timeit [solve(Ai, bi) for el in range(M)]
%timeit [solve_opt(Ai, bi) for el in range(M)]

# 10 loops, best of 3: 30.1 ms per loop
# 100 loops, best of 3: 3.77 ms per loop

型
这导致6.5倍的加速。
如果你需要更好的性能，你必须在Cython中移植这个for循环，并直接在C中接口gesv BLAS函数，如here所讨论的，或者更好地使用Scipy 0.16中的Cython API for BLAS/LAPACK。

编辑：正如@Eelco Hoogendoorn提到的，如果你的A矩阵是固定的，有一个更简单，更有效的方法。

除了直接使用blas函数的@rth部分，您还可以使用Numba。

@numba.njit()
def numba_solve(A, b):
    return np.linalg.solve(A, b)
#first execution to trigger compilation, otherwise compiling is included in the execution time
numba_solve(Ai[0], bi[0])

字符串
这里我们使用numba库，它通过使用即时（JIT）编译来编译所有可能的部分，从而加速Python函数。与solve_opt by rth类似，它确实大大加快了numpy.linalg.solve()周围的测试，并且是我测试过的最快的方法（见最后的表格）。
另外：如果你的矩阵是奇异的（或者可以是奇异的），并且你想得到最小二乘解scipy.linalg.lstsq是相当快的。Numba或直接使用Blas函数并没有快多少。
以下是建议的想法和其他一些想法的结果：
| 解算器| solver |
| --| ------------ |
| 2.4 numba_solve(A, b)个| numba_solve(A, b) |
| 3.6 solve_opt(A, b,gesv=gesv)个| solve_opt(A, b,gesv=gesv) |
| 4.8 np.linalg.solve(A, b)个| np.linalg.solve(A, b) |
| 13.5稀疏块诊断求解器（通过rth）| sparse blockdiag solver (by rth) |
| 21.3 scipy.linalg.lstsq(A, b,lapack_driver="gelsy")[0]个| scipy.linalg.lstsq(A, b,lapack_driver="gelsy")[0] |
| 22.1 numba_lstsq(A, b)| numba_lstsq(A, b) |
| 23.3 lstsq_opt(A, b,gelss=gelss)| lstsq_opt(A, b,gelss=gelss) |
| 24.2 np.linalg.lstsq(A, b, rcond=-1)[0]个| np.linalg.lstsq(A, b, rcond=-1)[0] |
| 24.2 scipy.linalg.solve(A, b)| scipy.linalg.solve(A, b) |
| 26.1 np.linalg.lstsq(A, b, rcond=None)[0]| np.linalg.lstsq(A, b, rcond=None)[0] |
| 41.3 scipy.linalg.lstsq(A, b,lapack_driver="gelss")[0]| scipy.linalg.lstsq(A, b,lapack_driver="gelss")[0] |
| 42.6 np.linalg.pinv(A) @ b| np.linalg.pinv(A) @ b |
| 43.7 scipy.linalg.lstsq(A, b,lapack_driver="gelsd")[0]| scipy.linalg.lstsq(A, b,lapack_driver="gelsd")[0] |
| 61.2 scipy.optimize.root(lambda x: A @ x - b, np.zeros(10), jac=lambda x: A, method='lm', tol=1e-1)| scipy.optimize.root(lambda x: A @ x - b, np.zeros(10), jac=lambda x: A, method='lm', tol=1e-1) |
| 64.6 scipy.linalg.pinv(A) @ b| scipy.linalg.pinv(A) @ b |
| scipy.optimize.root(lambda x: A @ x - b, np.zeros(10), jac=lambda x: A, method='hybr', tol=1e-1)| scipy.optimize.root(lambda x: A @ x - b, np.zeros(10), jac=lambda x: A, method='hybr', tol=1e-1) |
| 块诊断求解器（通过rth）| blockdiag solver (by rth) |
完整源代码：https://gist.github.com/te2be/4a1c4fb4e436ecaa61d74e09788e429f

**编辑：**在循环中使用Numba甚至可以更快。在这里，你必须使用Numba类型的List，因为在这种情况下，Numba在python列表上的循环效率有问题。

@numba.njit(cache=True)
def loop_numba_solve(Ai, bi):
    for A, b in zip(Ai, bi):
        x = np.linalg.solve(A, b)
    return x
nb_Ai = numba.typed.List()
nb_bi = numba.typed.List()
for A, b in zip(Ai, bi):
    nb_Ai.append(A)
    nb_bi.append(b)
loop_numba_solve(nb_Ai, nb_bi)

型
导致：
| 解算器| solver |
| --| ------------ |
| 1.4 loop_numba_solve(nb_Ai, nb_bi)| loop_numba_solve(nb_Ai, nb_bi) |