使用雅可比信息进行 scipy 全局优化

Question

我目前正在研究全局优化问题，为了加快计算速度，我想在局部优化中包含雅可比信息。我的函数并不简单，但它应该是可微分的，因此我想使用 pytorch 进行自动微分。就像这里所做的那样：https://bnikolic.co.uk/blog/pytorch/python/2021/02/28/pytorch-scipyminim.html

import numpy
import scipy, scipy.optimize
import torch


def minim(obs, f, p0):
    """Fit function f to observations "obs" starting at p0"""
    
    def fitfn(pars):
        # NB the require_grad parameter specifying we want to
        # differentiate wrt to the parameters
        pars=torch.tensor(pars,
                          requires_grad=True)
        y=f(pars)
        # Simple least-squares fitting
        res=((obs-y)**2).sum()
        res.backward()
        # Note that gradient is taken from the "pars" variable
        return res.data.cpu().numpy(), pars.grad.data.cpu().numpy()

    res=scipy.optimize.minimize(fitfn,
                                p0,
                                method="BFGS",
                                jac=True) # NB: we will compute the jacobian
    return res

# Points on which observations are made
X=torch.arange(100)

def exfn(p):
    y=torch.sin(p[0]*X)+torch.cos(p[1]*X)
    return y

# Sample observations
yp=exfn(torch.tensor([0.3, 0]) )

# Sample run
minim(yp, exfn, numpy.array([0.34, 0.01]))

因此，上面的工作原理是向最小化器提供关键字 jac=True ，此时它知道该函数不仅会返回该值，还会在求值时返回雅可比矩阵。 问题是我需要使用全局优化器，而框架似乎不支持这一点。至少我看不到任何方法来告诉全局优化器我的函数不仅返回值，还返回雅可比矩阵。

使用双重退火，它看起来像这样：

dual_annealing(f,bounds,args=(x,y0),minimizer_kwargs={'jac':True})

但这不起作用，因为它只告诉本地最小化器包含雅可比，而不是告诉本地最小化器周围的全局优化器。

这个问题有解决办法吗？或者有人有任何可能有效的替代框架吗？

Answer 1

有一种稍微灵活的方式来指定雅可比行列式，即提供可调用的，而不是 True。这使得最小化器可以调用提供函数值的函数，或者调用提供函数雅可比行列式的函数。

以下是如何更改函数以向双退火最小化器的局部搜索步骤提供雅可比行列式。

import numpy
import scipy, scipy.optimize
import torch

def minim_da(obs, f, p0):
    """Fit function f to observations "obs" starting at p0"""
    
    def fitfn(pars):
        # NB the require_grad parameter specifying we want to
        # differentiate wrt to the parameters
        pars=torch.tensor(pars,
                          requires_grad=True)
        y=f(pars)
        # Simple least-squares fitting
        res=((obs-y)**2).sum()
        res.backward()
        # Note that gradient is taken from the "pars" variable
        return res, pars

    def fitfn_val(pars):
        res, pars = fitfn(pars)
        return res.data.cpu().numpy()
    
    def fitfn_grad(pars):
        res, pars = fitfn(pars)
        return pars.grad.data.cpu().numpy()

    bounds = [(0, 1)] * 2
    res = scipy.optimize.dual_annealing(fitfn_val,
                                        bounds,
                                        maxiter=1000,
                                        minimizer_kwargs={'jac':fitfn_grad})
    return res


# Points on which observations are made
X=torch.arange(100)

def exfn(p):
    y=torch.sin(p[0]*X)+torch.cos(p[1]*X)
    return y

# Sample observations
yp=exfn(torch.tensor([0.3, 0]) )

# Sample run
print(minim_da(yp, exfn, numpy.array([0.34, 0.01])))

注意：我为您选择了一个

bounds

参数，因为使用双退火时需要该参数。您可能需要根据您的应用程序更改它。

注意：这并非 100% 可靠。大约 10% 的时间它无法找到最小值。您可以提高

maxiter

使其更加可靠，但代价是速度。