我有一个简单的凸问题,我正在尝试加快解决速度。我正在求解 argmin (theta)
其中 theta 和 rt 为 Nx1。
我可以使用cvxpy
轻松解决这个问题
import numpy as np
from scipy.optimize import minimize
import cvxpy
np.random.seed(123)
T = 50
N = 5
R = np.random.uniform(-1, 1, size=(T, N))
cvtheta = cvxpy.Variable(N)
fn = -sum([cvxpy.log(1 + cvtheta.T * rt) for rt in R])
prob = cvxpy.Problem(cvxpy.Minimize(fn))
prob.solve()
prob.status
#'optimal'
prob.value
# -5.658335088091929
cvtheta.value
# matrix([[-0.82105079],
# [-0.35475695],
# [-0.41984643],
# [ 0.66117397],
# [ 0.46065358]])
但是对于较大的 R
这会变得太慢,所以我正在尝试使用 scipy
的 fmin_cg
基于梯度的方法:
goalfun
是一个 scipy.minimize
友好函数,它返回函数值和梯度。
def goalfun(theta, *args):
R = args[0]
N = R.shape[1]
common = (1 + np.sum(theta * R, axis=1))**-1
if np.any( common < 0 ):
return 1e2, 1e2 * np.ones(N)
fun = np.sum(np.log(common))
thetaprime = np.tile(theta, (N, 1)).T
np.fill_diagonal(thetaprime, np.ones(N))
grad = np.sum(np.dot(R, thetaprime) * common[:, None], axis=0)
return fun, grad
确保函数和梯度正确:
goalfun(np.squeeze(np.asarray(cvtheta.value)), R)
# (-5.6583350819293603,
# array([ -9.12423065e-09, -3.36854633e-09, -1.00983679e-08,
# -1.49619901e-08, -1.22987872e-08]))
但是解决这个问题只会产生垃圾,无论方法
、迭代等如何。(唯一产生优化终止成功
的是如果x0
实际上等于最佳theta)
x0 = np.random.rand(R.shape[1])
minimize(fun=goalfun, x0=x0, args=R, jac=True, method='CG')
# fun: 3.3690101669818775
# jac: array([-11.07449021, -14.04017873, -13.38560561, -5.60375334, -2.89210078])
# message: 'Desired error not necessarily achieved due to precision loss.'
# nfev: 25
# nit: 1
# njev: 13
# status: 2
# success: False
# x: array([ 0.00892177, 0.24404118, 0.51627475, 0.21119326, -0.00831957])
即cvxpy
可以轻松处理这个看似无害的问题,但对于非凸求解器来说却是完全病态的。这个问题真的那么令人讨厌吗,还是我错过了什么?有什么替代方案可以加快速度?
最佳答案
我认为问题在于 theta
可能导致 log
参数变为负数。看来您已经发现了这个问题,并且在这种情况下让 goalfun
返回元组 (100,100*ones(N))
,显然,作为一种启发式尝试来建议求解器认为此“解决方案”不是首选。然而,必须施加一个更强的条件,即这个“解决方案”不可行。当然,这可以通过提供适当的约束来完成。 (有趣的是,cvxpy
似乎可以自动处理这个问题。)
这是一个示例运行,无需提供衍生产品。请注意使用可行的初始估计x0
。
np.random.seed(123)
T = 50
N = 5
R = np.random.uniform(-1, 1, size=(T, N))
def goalfun(theta, *args):
R = args[0]
N = R.shape[1]
common = (1 + np.sum(theta * R, axis=1))**-1
return np.sum(np.log(common))
def con_fun(theta, *args):
R = args[0]
return 1+np.sum(theta * R, axis=1)
cons = ({'type': 'ineq', 'fun': lambda x: con_fun(x, R)})
x0 = np.zeros(R.shape[1])
minimize(fun=goalfun, x0=x0, args=R, constraints=cons)
fun: -5.658334806882614 jac: array([ 0.0019, -0.0004, -0.0003, 0.0005, -0.0015, 0. ]) message: 'Optimization terminated successfully.' nfev: 92 nit: 12 njev: 12 status: 0 success: True x: array([-0.8209, -0.3547, -0.4198, 0.6612, 0.4605])
请注意,当我运行此命令时,我收到日志中遇到无效值
警告,表明在搜索中的某个时刻检查了theta
值,这几乎没有影响满足约束条件。然而,结果相当接近cvxpy
。当在 cvxpy.Problem
公式中显式施加约束时,检查 cvxpy
解决方案是否发生变化将会很有趣。
关于python - "pathological"凸函数的快速优化,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/40005320/