如何将np.linalg.solve映射到矩阵数组并保持速度？

Question

如何将np.linalg.solve映射到矩阵数组并保持速度？

3

我有一个线性问题需要解决一定的次数：Ax = B，其中A是一个维度为n的方阵，B是一个维度为n的向量。我需要找到x：

import numpy as np
A = np.random.rand(2,2)
B = np.random.rand(2)
x = np.linalg.solve(A,B)

这很基础。问题在于我想多次解决这个问题。当前的实现方式如下：

import numpy as np
k = 50 # the number of systems to solve
A_list = np.random.rand(k,2,2)
B_list = np.random.rand(k,2)
x = np.array([np.linalg.solve(A,B) for A, B in zip(A_list, B_list)])

但是它非常慢。在这个网站的人们的帮助下，我使用np.newaxis进行智能广播，消除了我的代码中的一个巨大瓶颈。我想知道是否有类似的技巧可以用于这种类型的函数（np.linalg.solve，np.linalg.det等）。

我使用np.vectorize进行测试失败了。

编辑：

输出

>>> import numpy as np
>>> k = 50 # the number of systems to solve
>>> A_list = np.random.rand(k,2,2)
>>> B_list = np.random.rand(k,2)
>>> x = np.array([np.linalg.solve(A,B) for A, B in zip(A_list, B_list)])
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-53-fecc7a7edaf9> in <module>()
----> 1 solution = np.linalg.solve(A_list,B_list)

/usr/lib/python3.3/site-packages/numpy/linalg/linalg.py in solve(a, b)
    309     if one_eq:
    310         b = b[:, newaxis]
--> 311     _assertRank2(a, b)
    312     _assertSquareness(a)
    313     n_eq = a.shape[0]

/usr/lib/python3.3/site-packages/numpy/linalg/linalg.py in _assertRank2(*arrays)
    153         if len(a.shape) != 2:
    154             raise LinAlgError('%d-dimensional array given. Array must be '
--> 155                     'two-dimensional' % len(a.shape))
    156 
    157 def _assertSquareness(*arrays):

LinAlgError: 3-dimensional array given. Array must be two-dimensional

- user1940040

是的，已经确认过了，它确实是一个新的算法实现，适用于1.8。 - alko

啊好的，那么除了升级之外就没有解决方案了。 - user1940040

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- alko · Accepted Answer

Numpy 1.8

从 numpy 版本 1.8 开始，您无需进行任何操作：

np.linalg.solve(A_list, B_list)

示例：

>>> import numpy as np
>>> np.random.seed(11)
>>> k = 10
>>> A_list = np.random.rand(k,2,2)
>>> B_list = np.random.rand(k,2)
>>> solution = np.linalg.solve(A_list,B_list)
>>> all(np.allclose(np.dot(A_list[i, :], solution[i,:]), B_list[i, :])
...         for i in range(A_list.shape[0]))
True

Numpy 1.7及之前版本

在旧版本中，可以尝试使用scipy.linalg.block_diag，但这会引入一些开销，包括内存和速度，并且对于较大的k值，它将输给zip方法：

import scipy.linalg
A = scipy.linalg.block_diag(*A_list)
B = B_list.reshape(-1)
solution = np.linalg.solve(A,B)
solution.reshape(-1,2)

1.8的速度测试

对于k = 2000；seed = 11 ：

>>> timeit('from __main__ import np, A_list, B_list; np.linalg.solve(A_list, B_list)', number = 100)
0.2786309433182055

>>> timeit('from __main__ import np, A_list, B_list; np.array([np.linalg.solve(A,B) for A, B in zip(A_list, B_list)])', number = 100)
8.431871369126554

>>> timeit('from __main__ import np, A, B; np.linalg.solve(A,B)', number = 100)
43.4851636171674712