python - 解决scipy中的多个独立优化

Question

我需要最小化大量(数千)不同输入的成本函数。显然，这可以通过遍历 scipy.optimize.minimize 来实现。或任何其他最小化例程。这是一个例子:

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)

a = np.random.randn(500, 40)
b = np.array(np.arange(500))

x = []
for i in range(a.shape[0]):
    res = sp.optimize.minimize(cost, np.zeros(40), args=(a[None, i], b[None, i]))
    x.append(res.x)

它为每个 a[i, :] 和 b[i 找到最小化 cost 的 x[i, :] ]，但这很慢。我猜遍历 minimize 会导致相当大的开销。

部分解决方案是同时解决所有 x:

res = sp.optimize.minimize(cost, np.zeros_like(a), args=(a, b))

这比循环还要慢。 minimize 不知道 x 中的元素在分组方面是独立的。因此，考虑到问题结构，尽管 block 对角矩阵就足够了，但它还是计算了完整的 hessian。这很慢并且会溢出我的计算机内存。

有什么方法可以通知minimize 或另一个优化函数关于问题结构，以便它可以在单个函数调用中解决多个独立的优化？ (类似于 certain options supported by Matlab's fsolve 。)

Answer 1

首先是一个解决方案:

结果 scipy.optimize.least_squares支持通过设置 jac_sparsity 参数来利用 jacobian 的结构。

least_squares 函数的工作方式与 minimize 略有不同，因此需要重写成本函数以返回残差:

def residuals(x, a, b):
    return np.sum(a * x.reshape(a.shape), axis=1) - b

雅可比具有 block 对角线稀疏结构，由

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

结论:

1. 原始解决方案(opt1)和不排序的重新使用起点(opt2)没有明显区别。
2. 循环 least_squares (opt3) 比循环 minimize (opt1, opt2).
3. 问题太大，无法天真地使用 least_squares 运行，因为雅可比矩阵不适合内存。
4. 在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)