高效求解大量线性最小二乘系统
Solve a multitude of linear least square system efficiently
我必须找到 >10^7 方程系统的最佳解决方案,每个方程有 5 个方程,每个方程有 2 个变量(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 是由 [=12 的堆叠条目组成的向量=].然后用 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。
我必须找到 >10^7 方程系统的最佳解决方案,每个方程有 5 个方程,每个方程有 2 个变量(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 是由 [=12 的堆叠条目组成的向量=].然后用 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。