Blog Content

cvxopt二次规划

Python 统计学-科学计算   2015-12-19 11:29:10

---

二次规划公示参考官网： http://cvxopt.org/userguide/coneprog.html?highlight=quadratic

二次规划问题的目标函数为二次函数，而约束条件则均为线性约束条件。 

上式中P为一对称矩阵。如果P是一个对称半正定矩阵，则上式中的目标函数是一个凸函数。目标函数为凸函数的 二次规划问题，如果其可行域非空的话，则它的任何一个局域最优解都是全局最优解。该类问题也是二次规划中应用较为广泛的 一类。二次规划问题可以通过内点法求解。

其中P,q,G,h,A,b为输入矩阵，该问题求解采用QP算法。

cvxopt模块采用的时小于等于形式，因此不等式两边需要同时乘以-1。

关于矩阵的维度，可以参考cvxpot.solvers.qp()的说明 。

min 1/2 * (x**2) + 3*x + 4 * y

subject to : x,y >= 0
x + 3y >= 15
2x + 5y <= 100
3x + 4y <= 80

参考： http://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf 

from cvxopt import matrix, solvers

P = matrix([[1.0, 0.0], [0.0, 0.0]]) # 初始化矩阵 np.diag([1.0,0.0])
q = matrix([3.0, 4.0])  # object: min  x 和y的参数
G = matrix([[-1.0, 0.0, -1.0, 2.0, 3.0], [0.0, -1.0, -3.0, 5.0, 4.0]]) # subject 左边x和y的参数。
h = matrix([0.0, 0.0, -15.0, 100.0, 80.0])   # subject 右边的值。

sol=solvers.qp(P, q, G, h)
print(sol['x'])
print(sol['primal objective'])

minimize  2 * (x**2) + (y**2) + xy + x + y
subject to : x >= 0
y >= 0
x + y = 1

from cvxopt import matrix, solvers

Q = 2*matrix([[2, .5], [.5, 1]])
p = matrix([1.0, 1.0])
G = matrix([[-1.0,0.0],[0.0,-1.0]])
h = matrix([0.0,0.0])
A = matrix([1.0, 1.0], (1,2))
b = matrix(1.0)

sol=solvers.qp(Q, p, G, h, A, b)
print(sol['x'])
print(sol['primal objective'])

参考：http://cs229.stanford.edu/section/cs229-cvxopt.pdf

相关参考：

凸优化中的基本概念
1.1 什么是凸集?

简单来说， 凸集是一个点集， 这个点集有一个性质， 就是在这个集合中任取不同的两个点x和y， 他们之间的线段（包括端点）上的点都属于这个点集，那么就说这个点集是一个凸集。

比如下图中左边的图形是凸集，而右边不是，因为我们可以找到两个点，使它们之间的线段上的点不在集合中

数学上，凸集的定义如下：

给定集合C，∀x,y∈C，0≤θ≤1，如果有

θx+(1−θy)∈C
我们就称集合C是凸集，我们把点θx+(1−θy)称为x和y的凸组合。

1.2 什么是凸函数?

假设有一个函数f:ℝn→ℝ，记其定义域为D(f)，如果D(f)是凸集，且在其中任取两个点x,y，满足以下性质：

f(θx+(1−θy))≤θf(x)+(1−θ)f(y)
那么就称f为凸函数。

注意：定义域是凸集这个要求不是必须的，其出发点只是为了使x，y的凸组合有定义

简单来说，我们在定义域任取两个点x，y， 连接他们得到一条线段，如果这个线段上的点都位于对应函数值上方，我们就说该函数是一个凸函数。

更进一步，如果x≠y且0<θ<1我们称f是严格凸的。如果−f是凸函数，那么f就是凹函数。如果−f是严格凸函数，那么f就是严格凹函数。

1.3 凸函数的等价判别方法

上面我们讲了什么是凸函数，然而这个定义在现实中很难用于判断一个函数是不是凸的，因此介绍几个等价的定义。

1.3.1 一阶近似

假设函数f:ℝn→ℝ是可导函数（也就是说f(x)的梯度∇xf(x)在整个定义域上都存在），则f是凸函数当且仅当 其定义域是凸集，且对于所有的x,y∈(f)有下式成立：

f(y)≥f(x)+∇xf(x)T(y−x)

我们将f(x)+∇xf(x)T(y−x)叫做对f的一阶近似，其物理意义实际上是经过点x的切平面，我们用这个切平面上的点来近似f(y)。这个公式的含义是：如果f是凸函数，那么它的一阶近似值始终位于函数值的下方。

1.3.2 二阶近似
假设函数f:ℝn→ℝ二阶可导（即海塞矩阵在定义域上都有定义），则f是凸函数当且仅当 其定义域是凸集且其海塞矩阵半正定，即：

∇2xf(x)⪰0
可能有些同学忘了海塞矩阵长什么样了，这里提一下。假设我们的变量来自n维空间，即x∈ℝn，我们记x=(x1,x2,...,xn)={xi}ni=1，即由n个变量组成的向量。那么海塞矩阵（记为H吧）是一个n×n的方块矩阵，且

Hij=∂2f(x)∂xi∂xj
 也就是说，Hij是f(x)分别对xi和xj进行求导两次得到的。

1.4 凸优化问题

上面已经介绍了凸集和凸函数，是时候到凸优化了吧？ 别急，在介绍凸优化概念之前再啰嗦两句。

1.4.1 水平子集(sublevel sets)

由凸函数的概念出发，我们可以引出水平子集（sublevel set）的概念。假定f(x)是一个凸函数， 给定一个实数α∈ℝ，我们把集合

{x ∈ D(f)|f(x)≤α}
叫做α−水平子集。 也就是说α水平子集是所有满足f(x)≤α的点构成的集合。利用凸函数性质，我们可以证明水平子集也是凸集：

f(θx+(1−θy))≤θf(x)+(1−θ)f(y)≤θα+(1−θ)α=α
水平子集告诉我们，给凸函数添加一个上限，定义域内剩下的点构成的点集还是一个凸集。

1.4.2 仿射函数(affine functions)

数学上，我们把形如

h(x)=Ax+b

的函数叫做仿射函数。其中，An×m，一个向量b∈ℝm。直观上理解，仿射函数将一个n维空间的向量通过线性变换A映射到m维空间，并在其基础上加上向量b，进行了平移。

同理，我们可以证明，点集

{x∈ D(h)|h(x)=0}
是一个凸集，证明略。

1.4.3 凸优化(convex optimization)

那么回到凸优化问题上来， 什么是一个凸优化问题？

一个凸优化问题可以定义为：

minimize f(x)
subject to : x ∈ C  

其中f是一个凸函数，C是一个凸集。根据先前介绍过的水平子集等概念，上面问题又可以等价写为：

minimize f(x)
subject to :
gi(x) <= 0  i = 1,2...m
hi(x) = 0  i = 1,2...p

其中，g(x)是凸函数,h(x)是仿射函数。 也就是说，原约束集C被我们表示为一系列凸集的交集（数学上可以证明，凸集的交集还是凸集）。

1.4.4 局部最优（local optima）和全局最优（global optima）

局部最优：周围小范围 内没有比我小的点。

数学定义：

如果存在R>0，对于所有的z：‖x−z‖2
全局最优：我就是整个定义域中的最小的点。

数学定义：

如果对于定义域内的所有z，有f(x)≤f(z)，则称x是全局最优。

现在回到凸优化问题上， 对于凸优化问题，有一个很重要的结论：

对于凸函数来讲， 局部最优就是全局最优。证明如下：

我们用反证法证明。设x是一个局部最优，但不是全局最优，于是我们假设全局最优是z∗，那么我们有f(x)>f(z∗)
由x的局部最优性质，我们有 ：

存在R>0，对于所有的z：‖x−z‖2我们考虑x和z∗的凸组合：z=θx+(1−θ)z∗，无论z∗在哪里，我们总可以找到一个θ，使得z位于x的邻域内，使得f(x)≤f(z)
另一方面，由凸函数性质，我们有：

f(z)=f(θx+(1−θ)z∗)≤θf(x)+(1−θ)f(z∗)<θf(x)+(1−θ)f(x)=f(x)
由此得f(z)
 1.5 常见凸优化问题

 线性规划
如果f和gi都是仿射函数，则凸优化问题变为了线性规划问题：
公式略，不好打。

二次规划
线性规划中，如果f变为一个凸二次函数，则凸优化问题变为二次规划：
公式略，不好打。

二次约束二次规划
f和gi都是凸二次函数
公式略，不好打。

半定规划
公式略，不好打。

其中，X∈ Sn是一个n维对称方阵，并且我们将它约束为半正定矩阵。C,Ai都是对称矩阵。这和前面的问题有点不太相同，前面是优化一个向量，而这里是优化一个矩阵。


参考：
http://cs229.stanford.edu/section/cs229-cvxopt.pdf

---

上一篇：python数据库orm之sqlalchemy的excute执行方式
下一篇：Sqlalchemy中mysql分组groupby用法 incompatible only_full_group_by解决方法