scipy 未优化并返回“由于精度损失不一定能实现所需的错误”答案

【问题标题】：scipy is not optimizing and returns "Desired error not necessarily achieved due to precision loss"scipy 未优化并返回“由于精度损失不一定能实现所需的错误”
【发布时间】：2014-09-06 04:18:02
【问题描述】：

我有以下代码试图最小化对数似然函数。

#!/usr/bin/python
import math
import random
import numpy as np
from scipy.optimize import minimize

def loglikelihood(params, data):
    (mu, alpha, beta) = params
    tlist = np.array(data)
    r = np.zeros(len(tlist))
    for i in xrange(1,len(tlist)):
        r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
    loglik  = -tlist[-1]*mu
    loglik = loglik+alpha/beta*sum(np.exp(-beta*(tlist[-1]-tlist))-1)
    loglik = loglik+np.sum(np.log(mu+alpha*r))
    return -loglik

atimes = [ 148.98894201,  149.70253172,  151.13717804,  160.35968355,
        160.98322609,  161.21331798,  163.60755544,  163.68994973,
        164.26131871,  228.79436067]
a= 0.01
alpha = 0.5
beta = 0.6
print loglikelihood((a, alpha, beta), atimes)

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'BFGS',args = (atimes,))
print res

它给了我

28.3136498357
./test.py:17: RuntimeWarning: invalid value encountered in log
  loglik = loglik+np.sum(np.log(mu+alpha*r))
   status: 2
  success: False
     njev: 14
     nfev: 72
 hess_inv: array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]])
      fun: 32.131359359964378
        x: array([ 0.01,  0.1 ,  0.1 ])
  message: 'Desired error not necessarily achieved due to precision loss.'
      jac: array([ -2.8051672 ,  13.06962156, -48.97879982])

请注意，它根本没有设法优化参数，并且最小值 32 大于 28，这就是你得到的 a= 0.01, alpha = 0.5, beta = 0.6 。通过选择更好的初始猜测可以避免这个问题，但如果是这样，我怎样才能自动做到这一点？

【问题讨论】：

我认为你想要最大化 LL，而不是最小化它。如果您要最小化平方和，那么您就是在最大化 LL。
@MikeDunlavey 是的。请注意该函数返回处理此问题的-loglik。
@felix 只是一个评论 - 我曾经遇到过与您的症状相同的问题，但原因完全不同。原来我的梯度函数有一个错误，所以当我通过jac 参数将它传递到例程中时，例程无法工作。这些错误很神秘，只有在重新检查我的代码后我才发现了这个错误。也就是说，下面使用Nelder-Mead 的答案确实很有帮助，因为它可以在没有渐变的情况下进行优化，并为我提供了正确的答案，帮助我意识到问题出在我的渐变函数中。
有人能解释一下警告指的是什么precision loss吗？我觉得它很晦涩。
这里好像有答案scicomp.stackexchange.com/questions/2004/…

标签： python math optimization scipy

【解决方案1】：

我复制了您的示例并尝试了一点。看起来如果您坚持使用 BFGS 求解器，经过几次迭代后，mu+ alpha * r 将有一些负数，这就是您获得 RuntimeWarning 的方式。

我能想到的最简单的解决方法是切换到 Nelder Mead 求解器。

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'Nelder-Mead',args = (atimes,))

它会给你这个结果：

28.3136498357
  status: 0
    nfev: 159
 success: True
     fun: 27.982451280648817
       x: array([ 0.01410906,  0.68346023,  0.90837568])
 message: 'Optimization terminated successfully.'
     nit: 92

【讨论】：

谢谢！关于为什么/何时 Nelder-Mead 会是比 BGFS 更好的选择的任何直觉？

【解决方案2】：

注意 log() 函数的负值，解决它们并通过添加惩罚告诉优化器它们是坏的：

#!/usr/bin/python
import math
import random
import numpy as np
from scipy.optimize import minimize

def loglikelihood(params, data):
    (mu, alpha, beta) = params
    tlist = np.array(data)
    r = np.zeros(len(tlist))
    for i in xrange(1,len(tlist)):
        r[i] = math.exp(-beta*(tlist[i]-tlist[i-1]))*(1+r[i-1])
    loglik = -tlist[-1]*mu
    loglik += alpha/beta*sum(np.exp(-beta*(tlist[-1]-tlist))-1)
    argument = mu + alpha * r
    limit = 1e-6
    if np.min(argument) < limit:
        # add a penalty for too small argument of log
        loglik += np.sum(np.minimum(0.0, argument - limit)) / limit
        # keep argument of log above the limit
        argument = np.maximum(argument, limit)
    loglik += np.sum(np.log(argument))
    return -loglik

atimes = [ 148.98894201,  149.70253172,  151.13717804,  160.35968355,
        160.98322609,  161.21331798,  163.60755544,  163.68994973,
        164.26131871,  228.79436067]
a= 0.01
alpha = 0.5
beta = 0.6
print loglikelihood((a, alpha, beta), atimes)

res = minimize(loglikelihood, (0.01, 0.1,0.1), method = 'BFGS',args = (atimes,))
print res

【讨论】：

【解决方案3】：

面对同样的警告，我通过重写对数似然函数以获取 log(params) 和 log(data) 作为参数而不是参数和数据来解决它。

因此，如果可能，我避免在似然函数或雅可比函数中使用np.log()。

【讨论】：

实际上，对于参数值具有完全不同比例的情况（例如我的情况），这听起来是一个非常好的主意。谢谢！

【解决方案4】：

另一个解决方案（对我有用）是将函数（和梯度）缩放到接近 0 的值。例如，当我必须评估 60k 点的对数似然时，我的问题就出现了。这意味着我的对数似然是一个非常大的数字。从概念上讲，对数似然是一个非常非常尖锐的函数。

梯度开始时很大（爬上这座尖峰山），然后变得适度小，但绝不会小于 BGFS 例程中默认的 gtol 参数（这是所有梯度必须低于终止的阈值） .此外，此时我基本上已经得出了正确的值（我使用的是生成的数据，所以我知道真实值）。

发生的事情是我的渐变是大约。 60k * average individual gradient value，即使 average individual gradient value 很小，比如小于 1e-8，60k * 1e-8 > gtol。因此，即使我已经找到了解决方案，我也从未达到过门槛。

从概念上讲，由于这座山峰非常陡峭，算法会迈出小步，但会超过真正的最小值并且从未达到average individual gradient << 1e-8，这意味着我的渐变从未低于gtol。

两种解决方案：

1) 将您的对数似然和梯度按一个因子缩放，例如 1/n，其中 n 是样本数。

2) 扩展您的gtol：例如"gtol": 1e-7 * n

【讨论】：

感谢您分享解决问题的方法。但是，我仍然想知道，为什么不重新缩放就不能收敛？如果我理解正确，您的解决方案相当于扩大 gtol。为什么即使使用真实的对数似然也无法达到所需的精度？

【解决方案5】：

我知道我迟到了，但我连续进行了 3 次优化。首先，我使用 Nelder-Mead 接近。如果没有先接近，我会得到太多溢出错误。然后我将 res.x 复制到下一个优化例程的起始参数。我发现鲍威尔是最可靠的，而且它通常做得很好。但是，然后我再次使用 Nelder-Mead 进行另一个最小化，以避免陷入局部最小值。
通常，使用 Powell 最小化后并没有太大改善。

【讨论】：