我想让学生在作业中解决一个二次元程序,而不需要他们安装额外的软件如cvxopt等。有没有一个只依赖NumPy/SciPy的python实现?
只依赖NumPy/SciPy的二次方程序(QP)求解器?
numpy
/
scipy
and
doesn't
需要额外的软件
like cvxopt
......有一个答案是推荐
cvxopt
,另一个(已接受的答案)是推荐本质上未维护的python与另一种语言的绑定(即非python实现)。
我对二次规划不是很熟悉,但我认为你可以只用
scipy.optimize
的约束最小化算法来解决这种问题。这里有一个例子。
import numpy as np
from scipy import optimize
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
# minimize
# F = x[1]^2 + 4x[2]^2 -32x[2] + 64
# subject to:
# x[1] + x[2] <= 7
# -x[1] + 2x[2] <= 4
# x[1] >= 0
# x[2] >= 0
# x[2] <= 4
# in matrix notation:
# F = (1/2)*x.T*H*x + c*x + c0
# subject to:
# Ax <= b
# where:
# H = [[2, 0],
# [0, 8]]
# c = [0, -32]
# c0 = 64
# A = [[ 1, 1],
# [-1, 2],
# [-1, 0],
# [0, -1],
# [0, 1]]
# b = [7,4,0,0,4]
H = np.array([[2., 0.],
[0., 8.]])
c = np.array([0, -32])
c0 = 64
A = np.array([[ 1., 1.],
[-1., 2.],
[-1., 0.],
[0., -1.],
[0., 1.]])
b = np.array([7., 4., 0., 0., 4.])
x0 = np.random.randn(2)
def loss(x, sign=1.):
return sign * (0.5 * np.dot(x.T, np.dot(H, x))+ np.dot(c, x) + c0)
def jac(x, sign=1.):
return sign * (np.dot(x.T, H) + c)
cons = {'type':'ineq',
'fun':lambda x: b - np.dot(A,x),
'jac':lambda x: -A}
opt = {'disp':False}
def solve():
res_cons = optimize.minimize(loss, x0, jac=jac,constraints=cons,
method='SLSQP', options=opt)
res_uncons = optimize.minimize(loss, x0, jac=jac, method='SLSQP',
options=opt)
print '\nConstrained:'
print res_cons
print '\nUnconstrained:'
print res_uncons
x1, x2 = res_cons['x']
f = res_cons['fun']
x1_unc, x2_unc = res_uncons['x']
f_unc = res_uncons['fun']
# plotting
xgrid = np.mgrid[-2:4:0.1, 1.5:5.5:0.1]
xvec = xgrid.reshape(2, -1).T
F = np.vstack([loss(xi) for xi in xvec]).reshape(xgrid.shape[1:])
ax = plt.axes(projection='3d')
ax.hold(True)
ax.plot_surface(xgrid[0], xgrid[1], F, rstride=1, cstride=1,
cmap=plt.cm.jet, shade=True, alpha=0.9, linewidth=0)
ax.plot3D([x1], [x2], [f], 'og', mec='w', label='Constrained minimum')
ax.plot3D([x1_unc], [x2_unc], [f_unc], 'oy', mec='w',
label='Unconstrained minimum')
ax.legend(fancybox=True, numpoints=1)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('F')
Output:
Constrained:
status: 0
success: True
njev: 4
nfev: 4
fun: 7.9999999999997584
x: array([ 2., 3.])
message: 'Optimization terminated successfully.'
jac: array([ 4., -8., 0.])
nit: 4
Unconstrained:
status: 0
success: True
njev: 3
nfev: 5
fun: 0.0
x: array([ -2.66453526e-15, 4.00000000e+00])
message: 'Optimization terminated successfully.'
jac: array([ -5.32907052e-15, -3.55271368e-15, 0.00000000e+00])
nit: 3
flxb
:
我怀疑这是否非常有效。我认为LOQO的一个实现。二次方程序设计的内部点代码(
citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.39.2191
)会更快。
ali_m
:
你需要你的学生解决的问题有多难?SLSQP在大约1.33msec的时间内解决了我的(公认的相当简单的)例子。它还可以处理任何界限、不等式和平等约束的组合。如果你想使用一个为QP优化的特殊求解器,那么你可能要(A)让你的学生安装额外的依赖,或者(B)自己编写它。
flxb
:
谢谢你的跟进。学生们应该用它来解决一个支持向量机问题,以便与他们应该实现的更有效的算法进行比较。这是一个大约100个变量的凸问题。我可能会实现LOQO,只是觉得我不能成为第一个。
值得在约束定义中加入'jac':(lambda x:-A),以使解算器更加稳健。
ximiki
:
一些有用的Python文档。
scipy.optimize
tutorial/docs (
docs.scipy.org/doc/scipy/reference/tutorial/optimize.html
,
docs.scipy.org/doc/scipy/reference/optimize.html
), scipy-lectures (
scipy-lectures.org/advanced/mathematical_optimization/...
).
IssamLaradji
发布于
2018-11-10
0
人赞同
这可能是一个迟到的答案,但我发现
CVXOPT
。-
http://cvxopt.org/
- 作为常用的免费python库,用于
Quadratic Programming
。然而,它并不容易安装,因为它需要安装其他的依赖项。
嗯,正如你所描述的,它不容易安装 :-)我对你的建议表示感谢,但我想我会先尝试其他方案。
@JimRaynor 我在OS X中直接安装
cvxopt
与
pip install cvxopt
没有问题。就这样了。替换代码2】会照顾到一切。而且我已经在几台机器上安装了
cvxopt
。当然,你需要安装编译器,但这也很简单,如果你正在使用
scipy
,你很可能已经有了它们。如果有帮助的话,我使用Anaconda作为Python发行版(它是完全免费的),安装Anaconda也很简单。你不需要管理员权限,也没有任何需要配置的东西。只要下载它,安装它,就可以使用了。
这个库是我转到Anaconda的原因之一,因为它便于管理依赖关系。我只是无法用pip来安装它。如果你已经有Anaconda,使用
conda install -c https://conda.anaconda.org/omnia cvxopt
就可以了。我在Windows 10和Python 2.7上。
divenex
:
请注意,该问题明确要求求解器能做到
not
需要安装
cvxopt
。
Tom Vacek
发布于
2018-11-10
0
人赞同
我碰到了一个很好的解决方案,想把它拿出来。在NICTA的ELEFANT机器学习工具包中,有一个LOQO的python实现(
http://elefant.forge.nicta.com.au
截至本帖为止)。 请看看optimization.intpointsolver。 这是由Alex Smola编写的,我已经使用了相同代码的C版本,并取得了巨大的成功。
我不相信这个项目是活跃的。 下载链接已被破坏,但这个链接可以使用。
elefant.forge.nicta.com.au/download/release/0.4/index.html
该项目有一个C++分叉,网址是
users.cecs.anu.edu.au/~chteo/BMRM.html
,但我相信它也不活跃。
divenex
:
这个答案中的链接是坏的,建议的软件不是纯Python(+Numpy/Scipy)。
Mike McKerns
发布于
2018-11-10
0
人赞同
替换代码0】为非线性/非凸性优化算法提供了一个纯粹的python实现,具有高级约束功能,通常只在QP求解器中发现。
mystic
实际上提供了比大多数QP求解器更强大的约束。然而,如果你正在寻找优化算法的速度,那么下面的内容就不适合你了。
mystic
并不慢,但它是纯粹的python,而不是python与C的绑定。如果你正在寻找非线性求解器的灵活性和QP约束功能,那么你可能会感兴趣。
Maximize: f = 2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2
Subject to: -2*x[0] + 2*x[1] <= -2
2*x[0] - 4*x[1] <= 0
x[0]**3 -x[1] == 0
where: 0 <= x[0] <= inf
1 <= x[1] <= inf
import numpy as np
import mystic.symbolic as ms
import mystic.solvers as my
import mystic.math as mm
# generate constraints and penalty for a nonlinear system of equations
ieqn = '''
-2*x0 + 2*x1 <= -2
2*x0 - 4*x1 <= 0'''
eqn = '''
x0**3 - x1 == 0'''
cons = ms.generate_constraint(ms.generate_solvers(ms.simplify(eqn,target='x1')))
pens = ms.generate_penalty(ms.generate_conditions(ieqn), k=1e3)
bounds = [(0., None), (1., None)]
# get the objective
def objective(x, sign=1):
x = np.asarray(x)
return sign * (2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2)
# solve
x0 = np.random.rand(2)
sol = my.fmin_powell(objective, x0, constraint=cons, penalty=pens, disp=True,
bounds=bounds, gtol=3, ftol=1e-6, full_output=True,
args=(-1,))
print 'x* = %s; f(x*) = %s' % (sol[0], -sol[1])
需要注意的是,
mystic
可以通用地将LP、QP和高阶平等和不平等约束应用于任何给定的优化器,而不仅仅是一个特殊的QP求解器。其次,
mystic
可以消化符号数学,所以定义/输入约束条件的便利性比处理矩阵和函数的导数要好一些。替换代码0】依赖于
numpy
,如果安装了
scipy
,也会使用它(不过,
scipy
不是必需的)。
mystic
利用
sympy
来处理符号约束,但它在一般情况下也不需要优化。
Output:
Optimization terminated successfully.
Current function value: -2.000000
Iterations: 3
Function evaluations: 103
x* = [ 2. 1.]; f(x*) = 2.0
Get mystic
here: https://github.com/uqfoundation
divenex
:
建议的解决方案没有使用二次编程求解器,而是使用了非线性求解器。如果一个通用的非线性求解器就足够了,一个更好的答案是@ali_m提出的,它只依赖于Numpy/Scipy。
@divenex: OP并没有要求一个QP求解器,他们要求的是解决一个只依赖于
numpy
/
scipy
的QP问题。......好吧,
mystic
中的求解器基本上只有一个
numpy
的依赖性(注意没有
scipy
的依赖性!)。 你可以从投票数中看到,@ali_m有一个更直接的解决方案(即使用
scipy
)--我也知道。 我的观点是,
mystic
可以解决QP问题,也可以解决非线性问题......所以,它是值得一提的。
divenex
:
我的评论是为了给其他像我一样来这个页面寻找QP的用户节省时间。
solver
正如问题标题中所说的那样。
mystic
并不包括这样一个QP求解器。事实上,你的通用优化器
mystic.fmin_powell
是你的软件包中包含的
Scipy
代码的一个副本。我建议人们直接调用
Scipy
。
@divenex:你错了,实际上。虽然
fmin_powell
的代码确实来自于
scipy.optimize
,但它已经被修改为接受任意的硬约束和软约束(并且可以因此它可以在一个空间中执行,如果需要的话,它实际上是一个QP、LP或MIP求解器)。如果你创建了硬的QP约束。
任何
你应用约束条件的神秘求解器只会在QP空间寻找解决方案。这是我帖子的全部内容,也是
mystic
的主要特点之一。如果你对
mystic
的工作方式有疑问,请随时直接与我联系。
divenex
:
你修改过的Powell算法也能优化一个受限的二次或线性函数,但这并不能使它 "有效地 "成为二次编程(QP)求解器。QP求解器的优点是,它利用了函数的二次形式,使收敛速度更快、更稳健。缺点是真正的QP求解器,与Powell算法不同,不能解决其他类型的函数,比如一般的非线性函数。这里有真正的Python QP求解器的链接,但没有一个是按要求 "只依赖于NumPy/SciPy "的
scaron.info/blog/quadratic-programming in-python.html
Tastalian
发布于
2018-11-10
0
人赞同
The
qpsolvers
包似乎也符合这个要求。它只依赖于NumPy,可以通过
pip install qpsolvers
来安装。然后,你可以这样做。
from numpy import array, dot
from qpsolvers import solve_qp
M = array([[1., 2., 0.], [-8., 3., 2.], [0., 1., 1.]])
P = dot(M.T, M) # quick way to build a symmetric matrix
q = dot(array([3., 2., 3.]), M).reshape((3,))
G = array([[1., 2., 1.], [2., 0., 1.], [-1., 2., -1.]])
h = array([3., 2., -2.]).reshape((3,))