scipy未优化，返回“由于精度损失，未必达到所需误差”

Question

我有一段代码，目的是想要最小化一个对数似然函数。

#!/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还大，而28是用a=0.01, alpha=0.5, beta=0.6时得到的。可能这个问题可以通过选择更好的初始值来避免，但如果要自动选择初始值，我该怎么做呢？

---

Nelder-Mead、TNC和SLSQP可以直接替换使用，其他方法就不行了。

Answer 1

我知道我来得有点晚，但我通常会连续做三次优化。首先，我使用Nelder-Mead方法来接近目标。如果不先接近目标，我会遇到很多溢出错误。接着，我把上一步得到的结果（res.x）复制到下一步优化的起始参数中。我发现Powell方法最可靠，通常效果也不错。但是，我还会再用Nelder-Mead方法进行一次最小化，以避免陷入局部最小值。
通常，在使用Powell方法后，效果提升不大。

Answer 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

Answer 3

遇到同样的警告，我通过重写对数似然函数来解决这个问题，把 log(params) 和 log(data) 作为参数传进去，而不是直接用 params 和 data。

这样，我就尽量避免在似然函数或雅可比矩阵中使用 np.log()。

Answer 4

还有一个解决办法（对我有效）就是把你的函数（和梯度）缩放到更接近0的值。举个例子，我遇到的问题是需要评估6万点的对数似然值。这意味着我的对数似然值是一个非常大的数字。从概念上讲，对数似然值就像一座非常尖锐的山。

一开始，梯度值很大（就像在爬这座尖锐的山），然后变得适中，但从来没有小于BGFS算法中的默认gtol参数（这是所有梯度必须低于的阈值，才能结束计算）。而且，这个时候我基本上已经找到了正确的值（因为我使用的是生成的数据，所以我知道真实的值）。

发生的事情是我的梯度大约是6万乘以平均个体梯度值，即使平均个体梯度值很小，比如小于1e-8，6万乘以1e-8依然大于gtol。所以即使我已经找到了答案，我也从来没有满足这个阈值。

从概念上讲，由于这座非常尖锐的山，算法在做小步前进，但却总是“越过”了真正的最小值，始终没有达到平均个体梯度 << 1e-8，这意味着我的梯度从来没有低于gtol。

有两个解决方案：

1) 将你的对数似然值和梯度缩放一个因子，比如1/n，其中n是样本数量。

2) 缩放你的gtol：例如"gtol": 1e-7 * n

Answer 5

我照着你的例子试了一下。看起来如果你一直使用BFGS这个求解器，经过几次迭代后，mu + alpha * r会出现一些负数，这就是你收到运行时警告的原因。

我能想到的最简单的解决办法就是换成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