高效求解多个线性最小二乘系统答案

【问题标题】：Solve a multitude of linear least square system efficiently高效求解多个线性最小二乘系统
【发布时间】：2017-03-23 06:42:22
【问题描述】：

我必须找到 >10^7 方程系统的最佳解决方案，其中每个方程有 2 个变量中的 5 个方程（5 次测量以找到长序列中误差最小的 2 个参数）。下面的代码（通常用来做曲线拟合）做我想要的：

#Create_example_Data
n = 100
T_Arm = np.arange(10*n).reshape(-1, 5, 2)
Erg = np.arange(5*n).reshape(-1, 5)
m = np.zeros(n)
c = np.zeros(n)
#Run
for counter in xrange(n):
     m[counter], c[counter] = np.linalg.lstsq(T_Arm[counter, :, :], 
                                              Erg[counter, :])[0]

不幸的是，它太慢了。有什么办法可以显着加快这段代码的速度吗？我试图对其进行矢量化，但没有成功。使用最后一个解决方案作为初始猜测可能也是一个好主意。使用scipy.optimize.leastsq 也没有加快速度。

【问题讨论】：

什么是Inputlen？是n吗？
n 是方程组的个数，等于 Inputlen，我更正了代码
我认为应该是 xrange(n) 而不是 xrange(len(n)) 因为 n 只是一个整数（在这种情况下为 100）
所以每个counter 值的leastsqr 是独立的？没有办法将其转换为更大的最小平方问题？
你能告诉我们解决所有这些方程需要多少时间，你认为什么时间可以接受？另外：我相信最好的办法是尝试和操作所有这些系统，以便您可以将它们组合在一起并减少要独立解决的系统数量......最后，如果所有其他方法都失败了，您可能会对多线程/处理感兴趣并行解决它们。如果我没记错的话，numpy 应该发布 GIL，所以即使是多线程也应该提供一些好处。

标签： python performance numpy linear-algebra

【解决方案1】：

您可以使用稀疏块矩阵 A，它在其对角线上存储 T_Arm 的 (5, 2) 个条目，并求解 AX = b，其中 b 是由 Erg 的堆叠条目组成的向量。然后用 scipy.sparse.linalg.lsqr(A, b) 求解系统。

要构造 A 和 b，我使用 n=3 进行可视化：

import numpy as np
import scipy
from scipy.sparse import bsr_matrix
n = 3
col = np.hstack(5 * [np.arange(10 * n / 5).reshape(n, 2)]).flatten()
array([ 0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  0.,  1.,  2.,  3.,  2.,
        3.,  2.,  3.,  2.,  3.,  2.,  3.,  4.,  5.,  4.,  5.,  4.,  5.,
        4.,  5.,  4.,  5.])

row = np.tile(np.arange(10 * n / 2), (2, 1)).T.flatten()
array([  0.,   0.,   1.,   1.,   2.,   2.,   3.,   3.,   4.,   4.,   5.,
         5.,   6.,   6.,   7.,   7.,   8.,   8.,   9.,   9.,  10.,  10.,
        11.,  11.,  12.,  12.,  13.,  13.,  14.,  14.])

A = bsr_matrix((T_Arm[:n].flatten(), (row, col)), shape=(5 * n, 2 * n))
A.toarray()
array([[ 0,  1,  0,  0,  0,  0],
       [ 2,  3,  0,  0,  0,  0],
       [ 4,  5,  0,  0,  0,  0],
       [ 6,  7,  0,  0,  0,  0],
       [ 8,  9,  0,  0,  0,  0],
       [ 0,  0, 10, 11,  0,  0],
       [ 0,  0, 12, 13,  0,  0],
       [ 0,  0, 14, 15,  0,  0],
       [ 0,  0, 16, 17,  0,  0],
       [ 0,  0, 18, 19,  0,  0],
       [ 0,  0,  0,  0, 20, 21],
       [ 0,  0,  0,  0, 22, 23],
       [ 0,  0,  0,  0, 24, 25],
       [ 0,  0,  0,  0, 26, 27],
       [ 0,  0,  0,  0, 28, 29]], dtype=int64)

b = Erg[:n].flatten()

然后

scipy.sparse.linalg.lsqr(A, b)[0]
array([  5.00000000e-01,  -1.39548109e-14,   5.00000000e-01,
         8.71088538e-16,   5.00000000e-01,   2.35398726e-15])

编辑：A 的内存并不像看起来那么大：更多关于块稀疏矩阵here。

【讨论】：

真棒答案！
好主意。但是数组 A 包含 (n*5)*(n*2) 个值，它们都是 int64，因此对于 n=1000，需要 80 MB，对于 n=10000，需要 8GB，依此类推。所有这些字节都需要被读取和写入。也许使用 n=100 并分块处理数据比原始解决方案更快。我会试一试的。
@Okapi575 A 不是数组（我只用 A.toarray() 来表示），而是一个稀疏矩阵，所以内存消耗会更低。我尝试了它以查看实现是否有效，但我承认我没有计时。
非常感谢您的好主意。在使用 n=1e5 对其进行计时后，时间消耗比原始代码长 10 倍以上。正如您所说，内存消耗确实较低
@Okapi575 好吧，我不是很乐观。但我将把它留在这里，因为我觉得这个想法很不错，并且使用 numpy 构建 A 很有趣。