使用scipy.optimize最小化多变量,可微函数
我试图用scipy.optimize
最小化以下功能:
其梯度是这样的:
(对于那些有兴趣的人来说,这是Bradley-Terry-Luce模型的成对比较的似然函数,与logistic回归非常紧密相关。)
很明显,为所有参数添加常量不会改变函数的值。 因此,我让 theta_1 = 0。下面是实现python中的目标函数和渐变(theta在此处变为x
):
def objective(x):
x = np.insert(x, 0, 0.0)
tiles = np.tile(x, (len(x), 1))
combs = tiles.T - tiles
exps = np.dstack((zeros, combs))
return np.sum(cijs * scipy.misc.logsumexp(exps, axis=2))
def gradient(x):
zeros = np.zeros(cijs.shape)
x = np.insert(x, 0, 0.0)
tiles = np.tile(x, (len(x), 1))
combs = tiles - tiles.T
one = 1.0 / (np.exp(combs) + 1)
two = 1.0 / (np.exp(combs.T) + 1)
mat = (cijs * one) + (cijs.T * two)
grad = np.sum(mat, axis=0)
return grad[1:] # Don't return the first element
以下是cijs
可能的外观示例:
[[ 0 5 1 4 6]
[ 4 0 2 2 0]
[ 6 4 0 9 3]
[ 6 8 3 0 5]
[10 7 11 4 0]]
这是我执行最小化的代码:
x0 = numpy.random.random(nb_items - 1)
# Let's try one algorithm...
xopt1 = scipy.optimize.fmin_bfgs(objective, x0, fprime=gradient, disp=True)
# And another one...
xopt2 = scipy.optimize.fmin_cg(objective, x0, fprime=gradient, disp=True)
但是,它在第一次迭代中总是失败:
Warning: Desired error not necessarily achieved due to precision loss.
Current function value: 73.290610
Iterations: 0
Function evaluations: 38
Gradient evaluations: 27
我无法弄清楚它为什么会失败。 由于此行显示错误:https://github.com/scipy/scipy/blob/master/scipy/optimize/optimize.py#L853
所以这个“沃尔夫线搜索”似乎不成功,但我不知道如何从这里开始......任何帮助表示赞赏!
作为@pv。 作为评论指出,我在计算梯度时犯了一个错误。 首先,我的目标函数梯度的正确(数学)表达式是:
(注意减号)。另外,我的Python实现完全错误,超出了符号错误。 这是我更新的渐变:
def gradient(x):
nb_comparisons = cijs + cijs.T
x = np.insert(x, 0, 0.0)
tiles = np.tile(x, (len(x), 1))
combs = tiles - tiles.T
probs = 1.0 / (np.exp(combs) + 1)
mat = (nb_comparisons * probs) - cijs
grad = np.sum(mat, axis=1)
return grad[1:] # Don't return the first element.
为了调试它,我使用了:
scipy.optimize.check_grad
:显示我的梯度函数产生的结果离很近的(有限差分)梯度很远。 scipy.optimize.approx_fprime
获取值的概念应该看起来像。 看来你可以将它转换成(非线性) 最小平方问题。 通过这种方式,你必须定义每个间隔n
,以建立系数矩阵变量,每个变量采样点的数量。
在这个例子中,我对所有变量使用相同数量的点和相同的间隔:
from scipy.optimize import leastsq
from numpy import exp, linspace, zeros, ones
n = 4
npts = 1000
xs = [linspace(0, 1, npts) for _ in range(n)]
c = ones(n**2)
a = zeros((n*npts, n**2))
def residual(c):
a.fill(0)
for i in range(n):
for j in range(n):
for k in range(npts):
a[i+k*n, i*n+j] = 1/(exp(xs[i][k] - xs[j][k]) + 1)
a[i+k*n, j*n+i] = 1/(exp(xs[j][k] - xs[i][k]) + 1)
return a.dot(c)
popt, pconv = leastsq(residual, x0=c)
print(popt.reshape(n, n))
#[[ -1.24886411 1.07854552 -2.67212118 1.86334625]
# [ -7.43330057 2.0935734 37.85989442 1.37005925]
# [ -3.51761322 -37.49627917 24.90538136 -4.23103535]
# [ 11.93000731 2.52750715 -14.84822686 1.38834225]]
编辑:关于上面构建的系数矩阵的更多细节:
链接地址: http://www.djcxy.com/p/20019.html上一篇: minimizing a multivariate, differentiable function using scipy.optimize