解决scipy中的多个独立优化问题

Question

我需要最小化一个成本函数，对于大量（1000个）不同的输入。很明显，这可以通过循环 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 中的元素是分组独立的。 因此，尽管块对角矩阵就足够了，但它计算完整的海森矩阵。 这很慢，会导致计算机内存溢出。

是否有任何方法可以通知 minimize 或另一个优化函数有关问题结构，以便它可以在单个函数调用中解决多个独立优化问题？ （类似于 Matlab 的 fsolve 支持的某些选项。）

- MB-F

---

“minimize”有许多可选参数。你觉得有哪些是有意义的？ - hpaulj

没有一个看起来能满足我的需求，但我不能排除我在寻找错误的东西 :) - MB-F

Answer 1

首先，提供一个解决方案：

原来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

结论：

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)

- MB-F

Answer 2

我猜循环最小化会造成很大的开销。

猜错了。最小化函数所需的时间远远超过任何循环开销。对于这个问题，没有向量化魔法可言。

可以通过使用更好的最小化起点来节省一些时间。首先，将参数排序，使连续循环具有相似的参数。然后使用上一个最小化的结束点作为下一个最小化的起点：

a = np.sort(np.random.randn(500, 40), axis=0)   # sorted parameters
b = np.arange(500)   # no need for np.array here, np.arange is already an ndarray

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

这在我的测试中可以节省30-40％的执行时间。

另一个较小的优化是预先分配适当大小的ndarray作为结果x值，而不是使用列表和append方法。循环之前：x = np.zeros((500, 40))；循环内部：x[i, :] = res.x。

- user6655984

---

改变起点是个好主意，但你现在是对a中的每一列进行排序，这会破坏行的结构。你能否提出一个不会破坏行结构但仍能提高性能的排序标准？（顺便说一句，我关心的不是循环开销，而是在minimize内部重复执行的检查和初始化开销。） - MB-F

正确的，使用向量参数排列它们使行之间靠近是一个单独的问题（值得另一个问题或搜索）。我这样做只是为了说明，因为数字无论如何都是随机的。我仍然认为大部分时间将花费在minimize内的迭代上。调整minimize的参数（例如，要求更少的精度）可能会产生大量的时间节省。 - user6655984