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

Posted

技术标签:

【中文标题】仅依赖于 NumPy/SciPy 的二次规划 (QP) 求解器?【英文标题】:Quadratic Program (QP) Solver that only depends on NumPy/SciPy? 【发布时间】:2013-06-05 06:37:35 【问题描述】:

我希望学生在作业中解决二次程序,而不必安装 cvxopt 等额外软件。是否有仅依赖于 NumPy/SciPy 的 python 实现?

【问题讨论】:

如果您可以提供一些关于二次程序的含义的链接,并且可能提供一两个示例,它将允许更多人回答这个问题。请更新您的问题,因为我不太确定您所说的 QP 是什么意思,而且我可能知道如何编写您的程序,尽管我不知道它需要什么。谢谢! 抱歉没有澄清。 QP 是一个特殊的线性代数问题,参见 Wikipedia (en.wikipedia.org/wiki/Quadratic_programming)。 我觉得奇怪的是,要求 python 实现的 QP 求解器 依赖于 numpy/scipy不需要需要额外的软件像cvxopt...有一个推荐cvxopt的答案和另一个(公认的答案)推荐基本上未维护的python绑定到另一种语言(即非-python 实现)。 【参考方案1】:

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

【讨论】:

我怀疑这是否非常有效。我认为 LOQO 的实现:二次规划的内点代码 (citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.39.2191) 会更快。 您需要学生解决的问题有多难? SLSQP 在大约 1.33 毫秒内解决了我的(诚然相当简单的)示例。它还可以处理边界、不等式和等式约束的任意组合。如果您决定使用针对 QP 优化的特定求解器,那么您可能必须 (A) 让您的学生安装额外的依赖项,或者 (B) 自己编写。 感谢您的跟进。学生应该使用它来解决支持向量机问题,并将其与他们应该实现的更有效的算法进行比较。这是一个包含大约 100 个变量的凸问题。我可能会实施 LOQO,只是觉得我不能成为第一个。 值得将 'jac':(lambda x:-A) 添加到约束定义中,以使求解器更加健壮。 我试图从头开始实现一些基本的机器学习算法。 SVM 在待办事项列表上,但我没有信心将其拉出。阅读您的答案后,我设法编写了自己的 svm (github.com/Sacry/mla_sani/blob/master/mla_sani/supervised/…),它按预期工作。非常感谢您的回答,非常感谢。【参考方案2】:

这可能是一个迟到的答案,但我发现CVXOPT - http://cvxopt.org/ - 作为Quadratic Programming 的常用免费python 库。但是,它并不容易安装,因为它需要安装其他依赖项。

【讨论】:

好吧,正如你所描述的,安装起来并不容易:-) 感谢您的建议,但我想我会先尝试其他选项。 @JimRaynor 我在 OS X 中直接用pip install cvxopt 安装cvxopt 没有问题。就是这样。 pip 负责一切。我已经在几台机器上安装了cvxopt。当然,您需要安装编译器,但这也很简单,如果您使用的是scipy,您很可能已经安装了它们。如果有帮助,我使用 Anaconda 作为 Python 发行版(完全免费),安装 Anaconda 也很简单。您不需要管理员权限,也不需要配置任何东西。只需下载并安装它,就可以使用了。 这个库是我切换到 Anaconda 以便于管理依赖项的原因之一。我只是无法用 pip 安装它。如果您已经拥有 Anaconda,请使用 conda install -c https://conda.anaconda.org/omnia cvxopt 并完成。我使用的是 Windows 10 和 Python 2.7。 请注意,该问题明确要求求解器需要安装cvxopt【参考方案3】:

我遇到了一个很好的解决方案,并想把它拿出来。在 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 users.cecs.anu.edu.au/~chteo/BMRM.html 有该项目的 C++ 分支,但我不相信它也处于活动状态。 此答案中的链接已损坏,建议的软件不是纯 Python (+Numpy/Scipy)【参考方案4】:

mystic 提供非线性/非凸优化算法的纯 Python 实现,具有通常仅在 QP 求解器中才有的高级约束功能。 mystic 实际上提供了比大多数 QP 求解器更强大的约束。但是,如果您正在寻找优化算法速度,那么以下内容不适合您。 mystic 并不慢,但它是纯 Python,而不是与 C 的 Python 绑定。如果您正在寻找非线性求解器中的灵活性和 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

【讨论】:

建议的解决方案不使用二次规划求解器,而是使用非线性求解器。如果一个通用的非线性求解器就足够了,一个更好的答案是 @ali_m 它只依赖于 Numpy/Scipy @divenex:OP 没有要求 QP 求解器,他们要求解决仅依赖于 numpy/scipy... 的 QP 问题。在mystic 中基本上只有numpy 依赖(注意没有scipy 依赖!)。您可以通过投票数看到@ali_m 有一个更直接的解决方案(即使用scipy)——我知道这一点。我的观点是,mystic 可以解决 QP 问题以及非线性问题......所以,值得一提。 我的评论旨在为像我一样访问此页面的其他用户节省时间,以寻找 QP 求解器,如问题标题中所述。 mystic 不包括这样的 QP 求解器。实际上,您的通用优化器mystic.fmin_powell 是您的包中包含的Scipy 代码的副本。我建议直接拨打Scipy @divenex:实际上,你错了。虽然fmin_powell 代码确实源自scipy.optimize,但它已被修改为接受任意硬约束和软约束(因此,如果需要,它可以在有效地使其成为 QP、LP 或 MIP 求解器的空间中执行)。如果您创建硬 QP 约束,任何您应用约束的神秘求解器将仅在 QP 空间中寻找解决方案。这是我帖子的全部要点,也是mystic 的主要功能之一。如果您对mystic 的工作方式有任何疑问,请随时直接与我联系。【参考方案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)

您还可以通过更改 solver 关键字参数来尝试不同的 QP 求解器(例如 Curious 提到的 CVXOPT),例如 solver='cvxopt'solver='osqp'

【讨论】:

qpsolvers 只是其他二次编程包(如cvxopt)的包装器,它需要编译器进行安装,而不是按要求使用纯 Python(+Numpy/Scypy)

以上是关于仅依赖于 NumPy/SciPy 的二次规划 (QP) 求解器?的主要内容,如果未能解决你的问题,请参考以下文章

最优化学习 常见凸优化问题的类型

Python都有哪些可以做带约束的二次线性规划的包

是否有“增强的”numpy/scipy dot 方法?

Haskell 中的二次规划

Jumpserver 深度集成OpenLdap的二次开发

Jumpserver 深度集成OpenLdap的二次开发