pytorch 具有雅可比信息的scipy全局优化

vatpfxk5  于 2023-10-20  发布在  其他
关注(0)|答案(2)|浏览(129)

我目前正在研究一个全局优化问题,为了加快计算速度,我想在局部优化中包含雅可比信息。我的函数并不简单,但它应该是可微分的,因此我想使用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})

但是这将不起作用,因为它只告诉局部最小化器包括雅可比矩阵,而不告诉围绕局部最小化器的全局优化器。
有解决此问题的方法吗?或者有人有其他可行的框架吗?

qnzebej0

qnzebej01#

有一种稍微更灵活的方法来指定雅可比矩阵,即提供一个可调用的,而不是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以使其更可靠,但以速度为代价。

xfyts7mz

xfyts7mz2#

最后,我在github上与scipy的一位贡献者就这个问题进行了很好的讨论。
线程可以在这里看到:https://github.com/scipy/scipy/issues/19305
这个主题也引出了尼克·奥德尔的建议。它展示了一个缓存派生结果的巧妙技巧,在某些情况下可能很方便。
但总的来说,结论是,你应该使用的方法取决于你正在使用的特定全局和局部优化器,在某些情况下,你可以编写一些没有浪费计算的代码,而在其他情况下,会有一些浪费的向前或导数计算。

相关问题