使用 Python 进行优化 (scipy.optimize)

Question

我正在尝试使用 Python 的 scipy.optimize 最大化以下功能。但是，经过多次尝试，它似乎不起作用。函数和我的代码粘贴在下面。感谢您的帮助！

问题

Maximize [sum (x_i / y_i)**gamma]**(1/gamma)
subject to the constraint sum x_i = 1; x_i is in the interval (0,1). 

x是选择变量的向量； y 是参数向量； gamma 是一个参数。 xs 的总和必须为 1。并且每个x 必须在区间 (0,1) 内。

代码

def objective_function(x, y):
    sum_contributions = 0
    gamma = 0.2

    for count in xrange(len(x)):
        sum_contributions += (x[count] / y[count]) ** gamma
    value = math.pow(sum_contributions, 1 / gamma)
    return -value

cons = ({'type': 'eq', 'fun': lambda x: np.array([sum(x) - 1])})

y = [0.5, 0.3, 0.2]
initial_x = [0.2, 0.3, 0.5]

opt = minimize(objective_function, initial_x, args=(y,), method='SLSQP',  
constraints=cons,bounds=[(0, 1)] * len(x))

Answer 1

有时，无论出于何种原因，数值优化器都不起作用。我们可以对问题进行稍微不同的参数化，它就会起作用。 （并且可能工作得更快）

例如，对于(0,1) 的边界，我们可以有一个变换函数，使得(-inf, +inf) 中的值在被变换后将在(0,1) 中结束

我们可以用等式约束做一个类似的技巧。例如，我们可以将维度从 3 减少到 2，因为 x 中的最后一个元素必须是 1-sum(x)。

如果还是不行，我们可以切换到不需要衍生信息的优化器，比如Nelder Mead。

还有Lagrange multiplier。

In [111]:

def trans_x(x):
    x1 = x**2/(1+x**2)
    z  = np.hstack((x1, 1-sum(x1)))
    return z

def F(x, y, gamma = 0.2):
    z = trans_x(x)
    return -(((z/y)**gamma).sum())**(1./gamma)
In [112]:

opt = minimize(F, np.array([0., 1.]), args=(np.array(y),),
               method='Nelder-Mead')
opt
Out[112]:
  status: 0
    nfev: 96
 success: True
     fun: -265.27701747828007
       x: array([ 0.6463264,  0.7094782])
 message: 'Optimization terminated successfully.'
     nit: 52

结果是：

In [113]:

trans_x(opt.x)
Out[113]:
array([ 0.29465097,  0.33482303,  0.37052601])

我们可以通过以下方式对其进行可视化：

In [114]:

x1 = np.linspace(0,1)
y1 = np.linspace(0,1)
X,Y = np.meshgrid(x1,y1)
Z = np.array([F(item, y) for item 
              in np.vstack((X.ravel(), Y.ravel())).T]).reshape((len(x1), -1), order='F')
Z = np.fliplr(Z)
Z = np.flipud(Z)
plt.contourf(X, Y, Z, 50)
plt.colorbar()

Answer 2

即使是棘手的这个问题也有点过时了，我想添加一个替代解决方案，这可能对将来遇到这个问题的其他人有用。

我们的问题是可以通过分析解决的。您可以从写下（等式约束）优化问题的拉格朗日开始：

L = \sum_i (x_i/y_i)^\gamma - \lambda (\sum x_i - 1)

通过将此拉格朗日的一阶导数设置为零来找到最优解：

0 = \partial L / \partial x_i = \gamma x_i^{\gamma-1}/\y_i - \lambda
=> x_i \propto y_i^{\gamma/(\gamma - 1)}

利用这种洞察力，可以通过以下方式简单有效地解决优化问题：

In [4]:
def analytical(y, gamma=0.2):
    x = y**(gamma/(gamma-1.0))
    x /= np.sum(x)
    return x
xanalytical = analytical(y)
xanalytical, objective_function(xanalytical, y)
Out [4]:
(array([ 0.29466774,  0.33480719,  0.37052507]), -265.27701765929692)

CT Zhu 的解决方案很优雅，但它可能违反第三坐标的正性约束。对于gamma = 0.2，这在实践中似乎不是问题，但对于不同的 gamma，您很容易遇到问题：

In [5]:
y = [0.2, 0.1, 0.8]
opt = minimize(F, np.array([0., 1.]), args=(np.array(y), 2.0),
               method='Nelder-Mead')
trans_x(opt.x), opt.fun
Out [5]:
(array([ 1.,  1., -1.]), -11.249999999999998)

对于与您的问题具有相同概率单纯形约束但没有解析解的其他优化问题，可能值得研究投影梯度方法或类似方法。这些方法利用了将任意点投影到该集合上的快速算法这一事实，请参阅https://en.wikipedia.org/wiki/Simplex#Projection_onto_the_standard_simplex。

（要查看完整的代码和更好的方程式渲染，请查看 Jupyter notebook http://nbviewer.jupyter.org/github/andim/pysnippets/blob/master/optimization-simplex-constraints.ipynb）