Python 仅依赖于 NumPy/SciPy 的二次程序 (QP) 求解器?

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/17009774/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me): StackOverFlow

提示:将鼠标放在中文语句上可以显示对应的英文。显示中英文
时间:2020-08-19 00:12:46  来源:igfitidea点击:

Quadratic Program (QP) Solver that only depends on NumPy/SciPy?

pythonnumpyscipymathematical-optimization

提问by flxb

I would like students to solve a quadratic program in an assignment without them having to install extra software like cvxopt etc. Is there a python implementation available that only depends on NumPy/SciPy?

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

采纳答案by Tom Vacek

I ran across a good solution and wanted to get it out there. There is a python implementation of LOQO in the ELEFANT machine learning toolkit out of NICTA (http://elefant.forge.nicta.com.auas of this posting). Have a look at optimization.intpointsolver. This was coded by Alex Smola, and I've used a C-version of the same code with great success.

我遇到了一个很好的解决方案,想把它拿出来。在 NICTA 的 ELEFANT 机器学习工具包中有一个 LOQO 的 Python 实现(http://elefant.forge.nicta.com.au截至本文)。看看 optimization.intpointsolver。这是由 Alex Smola 编写的,我使用了相同代码的 C 版本并取得了巨大成功。

回答by ali_m

I'm not very familiar with quadratic programming, but I think you can solve this sort of problem just using scipy.optimize's constrained minimization algorithms. Here's an example:

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

enter image description here

在此处输入图片说明

回答by Curious

This might be a late answer, but I found CVXOPT- http://cvxopt.org/- as the commonly used free python library for Quadratic Programming. However, it is not easy to install, as it requires the installation of other dependencies.

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

回答by Mike McKerns

mysticprovides a pure python implementation of nonlinear/non-convex optimization algorithms with advanced constraints functionality that typically is only found in QP solvers. mysticactually provides more robust constraints than most QP solvers. However, if you are looking for optimization algorithmic speed, then the following is not for you. mysticis not slow, but it's pure python as opposed to python bindings to C. If you are looking for flexibility and QP constraints functionality in a nonlinear solver, then you might be interested.

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])

Things to note is that mysticcan generically apply LP, QP, and higher order equality and inequality constraints to any given optimizer, not just a special QP solver. Secondly, mysticcan digest symbolic math, so the ease of defining/entering the constraints is a bit nicer than working with the matrices and derivatives of functions. mysticdepends on numpy, and will use scipyif it is installed (however, scipyis not required). mysticutilizes sympyto handle symbolic constraints, but it's also not required for optimization in general.

需要注意的是,mystic可以将 LP、QP 和高阶等式和不等式约束一般应用于任何给定的优化器,而不仅仅是特殊的 QP 求解器。其次,mystic可以消化符号数学,因此定义/输入约束的难度比使用矩阵和函数的导数要好一些。mystic取决于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 mystichere: https://github.com/uqfoundation

获取mystic此:https://github.com/uqfoundation

回答by Stéphane Caron

The qpsolverspackage also seems to fit the bill. It only depends on NumPy and can be installed by pip install qpsolvers. Then, you can do:

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)

You can also try different QP solvers (such as CVXOPT mentioned by Curious) by changing the solverkeyword argument, for example solver='cvxopt'or solver='osqp'.

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