首先,提供一个解决方案:
原来scipy.optimize.least_squares
支持设置jac_sparsity
参数以利用雅可比矩阵的结构。
least_squares
函数与minimize
略有不同,因此代价函数需要重写为返回残差:
def residuals(x, a, b):
return np.sum(a * x.reshape(a.shape), axis=1) - b
雅可比矩阵具有块对角稀疏结构,如下所示:
jacs = sp.sparse.block_diag([np.ones((1, 40), dtype=bool)]*500)
并调用优化程序:
res = sp.optimize.least_squares(residuals, np.zeros(500*40),
jac_sparsity=jacs, args=(a, b))
x = res.x.reshape(500, 40)
但它真的更快吗?
%timeit opt1_loopy_min(a, b) # 1 loop, best of 3: 2.43 s per loop
%timeit opt2_loopy_min_start(a, b) # 1 loop, best of 3: 2.55 s per loop
%timeit opt3_loopy_lsq(a, b) # 1 loop, best of 3: 13.7 s per loop
%timeit opt4_dense_lsq(a, b) # ValueError: array is too big
%timeit opt5_jacs_lsq(a, b) # 1 loop, best of 3: 1.04 s per loop
结论:
- 如果不排序,原始解决方案(
opt1
)和重复使用起始点(opt2
)之间没有明显的区别。
- 循环运行
least_squares
(opt3
)比循环运行 minimize
(opt1
、opt2
)要慢得多。
- 问题太大,不能简单地使用
least_squares
运行,因为雅可比矩阵不适合放在内存中。
- 利用
least_squares
中雅可比稀疏结构(opt5
)似乎是最快的方法。
这是时间测试环境:
import numpy as np
import scipy as sp
def cost(x, a, b):
return np.sum((np.sum(a * x.reshape(a.shape), axis=1) - b)**2)
def residuals(x, a, b):
return np.sum(a * x.reshape(a.shape), axis=1) - b
a = np.random.randn(500, 40)
b = np.arange(500)
def opt1_loopy_min(a, b):
x = []
x0 = np.zeros(a.shape[1])
for i in range(a.shape[0]):
res = sp.optimize.minimize(cost, x0, args=(a[None, i], b[None, i]))
x.append(res.x)
return np.stack(x)
def opt2_loopy_min_start(a, b):
x = []
x0 = np.zeros(a.shape[1])
for i in range(a.shape[0]):
res = sp.optimize.minimize(cost, x0, args=(a[None, i], b[None, i]))
x.append(res.x)
x0 = res.x
return np.stack(x)
def opt3_loopy_lsq(a, b):
x = []
x0 = np.zeros(a.shape[1])
for i in range(a.shape[0]):
res = sp.optimize.least_squares(residuals, x0, args=(a[None, i], b[None, i]))
x.append(res.x)
return x
def opt4_dense_lsq(a, b):
res = sp.optimize.least_squares(residuals, np.zeros(a.size), args=(a, b))
return res.x.reshape(a.shape)
def opt5_jacs_lsq(a, b):
jacs = sp.sparse.block_diag([np.ones((1, a.shape[1]), dtype=bool)]*a.shape[0])
res = sp.optimize.least_squares(residuals, np.zeros(a.size), jac_sparsity=jacs, args=(a, b))
return res.x.reshape(a.shape)