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

我必须找到大于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并没有加快速度.

解决方法:

您可以使用稀疏块矩阵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上更多.

标签：performance,linear-algebra,python,numpy
来源： https://codeday.me/bug/20191118/2024634.html