仅依赖于NumPy/SciPy的二次规划(QP)求解器？

Question

仅依赖于NumPy/SciPy的二次规划(QP)求解器？

pythonnumpyscipymathematical-optimization

42

我希望学生们在作业中解决二次规划问题时不需要安装额外的软件，如cvxopt等。是否有仅依赖于NumPy/SciPy的Python实现可用？

- flxb

如果您能提供一些关于二次规划的链接，以及一两个示例，将有助于更多人回答这个问题。请更新您的问题，因为我不太确定您所说的QP是什么，虽然我可能知道如何编写您的程序，但我不知道它需要什么。谢谢！ - Ryan Saxe

2

抱歉没有说明清楚。QP是一个特殊的线性代数问题，请参见维基百科（http://en.wikipedia.org/wiki/Quadratic_programming）。 - flxb

9

我觉得很奇怪，一个关于寻找仅依赖于 numpy/scipy 而且不需要额外软件（如 cvxopt）的 Python 实现 QP 求解器的问题，有一个答案推荐使用 cvxopt，而另一个答案（即被采纳的答案）推荐的实际上是绑定到另一种语言的未维护的 Python 库（即非 Python 实现）。 - Mike McKerns

5个回答

11

可能有点晚了，但我发现 CVXOPT - http://cvxopt.org/ - 是常用的免费Python库，用于二次规划。然而，安装并不容易，因为它需要安装其他依赖项。

- IssamLaradji

2

好的，正如你所描述的那样，安装并不容易 :-) 点赞作为我的感谢，但我想我会先尝试其他选项。 - Jim Raynor

2

@JimRaynor 在 OS X 上，我直接用 pip install cvxopt 安装 cvxopt 没有任何问题。 pip 会搞定一切。而且，我已经在几台机器上安装了 cvxopt 。当然，你需要安装编译器，但这也很简单，如果您正在使用 scipy ，则很可能已经安装好了。如果有帮助的话，我使用 Anaconda 作为 Python 分发工具（完全免费），安装 Anaconda 也很简单。无需管理员权限，也没有任何需要配置的内容。只需下载、安装即可使用。 - Amelio Vazquez-Reina

3

这个库是我转用Anaconda的原因之一，因为它很容易管理依赖关系。我无法使用pip安装它。如果你已经安装了Anaconda，则可以使用“conda install -c https://conda.anaconda.org/omnia cvxopt”来安装它，然后就完成了。我使用的是Windows 10和Python 2.7。 - blue_chip

3

请注意，问题明确要求不需要安装“cvxopt”的解决方案。 - divenex

5

我发现了一个好的解决方案并想要分享出来。在NICTA的ELEFANT机器学习工具包中有一个LOQO的Python实现(http://elefant.forge.nicta.com.au，本文发布时依然存在)。请查看 optimization.intpointsolver。这个代码是由Alex Smola编写的，我曾经成功地使用过相同代码的C版本。

- Tom Vacek

2

我不相信这个项目还在活跃。下载链接已经失效，但是这个链接可以使用：https://elefant.forge.nicta.com.au/download/release/0.4/index.html。此外，该项目有一个C++分支，位于http://users.cecs.anu.edu.au/~chteo/BMRM.html，但我认为它也不再活跃。 - Tom Vacek

这个答案中的链接已经失效，而且建议使用的软件不是纯Python（+Numpy/Scipy）。 - divenex

5

这个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,))

# min. 1/2 x^T P x + q^T x with G x <= h
print "QP solution:", solve_qp(P, q, G, h)

您还可以尝试不同的QP求解器（例如Curious提到的CVXOPT），通过更改solver关键字参数，例如solver='cvxopt'或solver='osqp'。

- Tastalian

qpsolvers只是其他二次规划包（如cvxopt）的包装器，它需要编译器进行安装，并且不符合所需的纯Python（+Numpy/Scipy）要求。 - divenex

2

提供了一个纯Python实现的非线性/非凸优化算法，具有高级约束功能，通常只在QP求解器中找到。实际上，提供的约束比大多数QP求解器更加健壮。但是，如果您正在寻找优化算法速度，则以下内容不适用于您。并不慢，但它是纯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 能够处理符号数学，因此定义/输入约束的便利性比使用函数的矩阵和导数要好一些。mystic 依赖于 numpy，如果已安装，则将使用 scipy（但是，scipy 不是必需的）。mystic 利用 sympy 来处理符号约束，但这在一般情况下并不是必需的优化工具。

Optimization terminated successfully.
         Current function value: -2.000000
         Iterations: 3
         Function evaluations: 103

x* = [ 2.  1.]; f(x*) = 2.0

获取 mystic：在这里https://github.com/uqfoundation

- Mike McKerns

建议的解决方案不使用二次规划求解器，而是使用非线性求解器。如果通用的非线性求解器足够好，那么更好的答案是由@ali_m提供的，它只依赖于Numpy/Scipy。 - divenex

@divenex：OP并没有要求QP求解器，他们只是要求一些能够解决仅依赖于numpy/scipy的QP问题的东西...而且，mystic中的求解器基本上只有一个numpy依赖项（请注意没有scipy依赖项！）。从投票数可以看出，@ali_m有一个更直接的解决方案（即使用scipy）——我知道这一点。我的观点是，mystic可以解决QP问题以及非线性问题...所以，值得一提。 - Mike McKerns

我的评论旨在为像我一样来到这个页面寻找QP 求解器的其他用户节省时间，正如问题标题所述。 mystic 不包括这样的 QP 求解器。实际上，您的通用优化器 mystic.fmin_powell 是您的软件包中包含的 Scipy 代码的副本。我建议直接调用 Scipy。 - divenex

@divenex：你错了。虽然fmin_powell代码确实是从scipy.optimize派生出来的，但它已经被修改为接受任意硬性和软性约束（因此，如果需要，它可以在有效地使其成为QP、LP或MIP求解器的空间中执行）。如果创建硬性QP约束，则将约束应用于的任何神秘求解器都只会在QP空间中寻找解决方案。这就是我发帖的全部意图，也是mystic的主要特点之一。如果您对mystic的工作原理有疑问，请随时直接与我联系。 - Mike McKerns

你修改的Powell算法也能优化约束二次或线性函数，并不意味着它“有效地”成为了一个二次规划（QP）求解器。QP求解器的优点在于它利用函数的二次形式，从而实现更快速和更稳健的收敛。缺点是真正的QP求解器，与Powell算法不同，不能解决其他类型的函数，如一般非线性函数。以下是真正的Python QP求解器链接，但没有一个“仅依赖于NumPy/SciPy”，如所请求的https://scaron.info/blog/quadratic-programming-in-python.html - divenex

谢谢提供链接。你可以看一下mystic中约束条件的工作方式，它们本质上是坐标变换，这可能会对你有所帮助。 - Mike McKerns

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- ali_m · Accepted Answer

我对二次规划不是很熟悉，但我认为你可以使用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')

输出：

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

enter image description here