高效解决多个线性最小二乘系统

Question

高效解决多个线性最小二乘系统

3

我需要找到最佳解决方案，用于求解超过10^7个包含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也无法加快速度。

- Okapi575

Inputlen是什么？它是n吗？ - TuanDT

我认为应该使用xrange(n)而不是xrange(len(n))，因为n只是一个整数（在这种情况下为100）。 - TuanDT

1

那么每个counter值的leastsqr是独立的吗？没有将其转换为更大的leastsqr问题的方法吗？ - hpaulj

1

你能告诉我们解决所有这些方程需要多长时间，你觉得什么时间是可以接受的吗？另外：我认为最好的方法是尝试操纵所有这些系统，将它们归为一组并减少独立求解的系统数量... 最后，如果一切都失败了，你可能会对多线程/处理感兴趣，以并行地解决它们。如果我没记错的话，NumPy应该会释放GIL，所以即使是多线程也会带来一些好处。 - Bakuriu

@Bakuriu：目前我在这个循环中每1e6次评估需要大约一分钟的时间，这是不切实际的。我梦想着能够通过使用numba或矢量化numpy来实现100倍的速度提升。多进程或使用GPU也可以是一个选择，但首先需要优化代码...... - Okapi575

显示剩余5条评论

1个回答

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

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

- P. Camilleri · Accepted Answer

你可以使用稀疏块矩阵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在内存中并不像看起来那么大：有关块稀疏矩阵的更多信息在此处。