首发于策略算法工程师之路

策略算法工程师之路-凸二次优化

洪九(李戈)

技术改变生活，让平凡人活的更从容!

1.凸优化

凸优化问题（OPT，convex optimization problem）指定义在凸集中的凸函数的最优化的问题。尽管凸优化的条件比较苛刻，但仍然在机器学习领域有十分广泛的应用。这主要得益于凸优化的优良性质:

1).凸优化问题的局部最优解就是全局最优解。

2).很多非凸问题可以等价或近似转化为凸优化问题（比如拉格朗日对偶）。

3).凸优化问题的研究较为成熟。求解方法也较为成熟。

首先来了解几个概念性的问题。

什么是凸集?

定义: C 是凸集，如果对任意的 x,y\in C 和任意的 \theta \in R 满足 0\leq \theta \leq 1 时， \theta x +(1-\theta) y \in C 恒成立。

直观来讲，任取一个集合中的两点连成一条线段，如果这条线段完全落在该集合中，那么这个集合就是凸集。

什么是凸函数?

定义:定义在 R^{n}->R 上的函数 f 是凸函数，如果它的定义域 D(f) 是一个凸集且对任意的 x,y\in D 和 0 \leq \theta \leq 1 , f(\theta x + (1-\theta y)) \leq \theta x f(x) + (1-\theta x)f(y) 恒成立。直观来讲就是形状看上去“凹”下去的函数，注意可不是看上去“凸”的。

什么是非凸函数?

凸优化的形式化定义:

其中， x 为决策变量， f(x) 和 g_{i}(x) 均为凸函数， h_{j}(x) 为仿射(线性)函数。根据 f(x) 和 g(x) 的不同形式，又可以将凸优化细分为以下几种不同的类型。

1).线性规划

其中目标函数 f(x) 和不等式约束 g(x) 都为仿射(线性)函数。

2).二次规划

其中目标函数 f(x)为凸二次型， g(x) 不等式约束为仿射(线性)函数。

3).二次约束的二次规划

其中目标函数 f(x) 和不等式约束 g(x) 都为凸二次型。本章我们主要关注的是“二次规划”。

说起运筹建模，在实践中工程师通常有意将问题建模为“凸优化”问题，其主要原因是“别的我也不会!”，有些时候与自己的智商妥协也不是一件很丢人的事情而是常态。

上图中是阿里的中间件专家沈询(花名)，代表作TDDL(分布式数据库中间件)。当年转正答辩技术环节选的就是TDDL，并且选了其中最难的编译原理部分讲解。因为笃定评委们都不懂，所以就可以随便讲了。即便讲的不对大家也听不出来且还能凸显自己的技术深度。

言归正传，曾经在一次分享中有人问沈询“解决分布式事务最好的方案是什么?”。当时沈询给出的回答是“没有事务”。不愧曾是我的偶像(现在也是)，我喜欢这个答案。

简单思考过关于“科学思维”和“工程思维”在面对“难点”时处理思路的不同。在面对“难点”时，工程同学首先应该想到的是如何“逃避”，解决不了的问题就想办法绕过去，因为我们的首要目标是解决业务问题。而科学家在面对“难点”时，一定是想办法攻克，因为科学家的首要目标为“探索未知世界”。

2.梯度下降法

还是要先简单说下凸优化中常见的“梯度下降法”！

大多数时候我们仅仅知道目标函数是“凸”的，但并不知道最优值(最大值或最小值)在哪。

如前文所说，由于凸优化具有“局部最优解等于全局最优解”的优良性质，因此只要任意选定一个“初始解”，然后沿着“使得目标函数 f(x) 不断减小” 的方向移动，当我们走到“局部最高点”的时候也就走到了“全局最高点”。

将上述讨论用数学语言表达。我站在当前位置 x_{k} ，选择了方向 d_{k} ，然后沿着选择的方向行走了 a_{k} 距离，最终到了新的位置 x_{k+1} 。总结成如下公式:

我们希望新的位置的函数值 f(x_{k+1}) 小于(或大于) f(x_{k}) ，只有这样才能一步一步逼近我们的终极目标 min f(x) 。图中好像画反了(成山顶了)!

所以这里至关重要的是方向 d_{k} 的选择。选择错了越努力越痛苦!

注: \bigtriangledown f(x_{k}) 是 f(x_{k}) 的梯度向量，如下:

3.无约束凸二次优化

3.1 简单凸二次型

min \frac{x^{2}+y^{2}}{2}

针对 f(x,y)=\frac{x^{2}+y^{2}}{2} ,显然最小值0在 (0,0) 处取得。出于讲解需要这里分别采用无约束最优性条件和梯度下降法求解。

1).无约束最优性条件

所谓无约束最优性条件，是指最优化问题的最优解所要满足的必要条件或充分条件。对于无约束最优性条件通常包含“一阶条件”和“二阶条件”。

一阶条件：如果x¯是局部极小点，则在该点梯度为0。

二阶条件：如果x¯是局部极小点，则在该点Hessian矩阵半正定；如果梯度为0

且Hessian矩阵在x¯微分领域内半正定，则x¯是局部极小点。

令 f(x,y) 的偏导为0:

2).梯度下降法

3.2 线性回归

简单了解下什么是回归任务? 如下一组房屋价格数据:

把房屋的面积(Living area)记为特征 x_{1} ，把居室的数目记为特征 x_{2} ，价格记为 y 。我们假设满足如下关系:

因此可以建立模型:

接下来就是如何确定 \theta _{1} 、 \theta _{2} 和 \theta _{3} 的值。我们希望第i个样本的预测值 h_{\theta}(x) 与真实值 y 的误差 J(\theta) 越小越好，用数学语言表达如下:

再考虑所有样本的总体误差:

根据前面所讲“梯度下降法”我们需要确定每一步的梯度方向:

最后迭代更新参数直到收敛:

参考资料:

1. [机器学习必知必会]从无约束优化到拉格朗日法

https://www.jianshu.com/p/0f368102d0aa

2.回归:线性回归

https://www.jianshu.com/p/3a98f33113ac

3. Machine_learning

https://github.com/uestc-buddy/Machine_learning/blob/master/softmax/softmax.py

4.通过线形回归了解算法流程

https://www.cnblogs.com/xiaoyuanqujing/articles/11656814.html

4.有约束凸二次优化

实际上每个人的精力都是有限的，能行走的路程也不是无穷的。因此如何在有限的生命中成就最好的自己，就是“有约束的最优化问题”。

明知道走不到“山顶”，还要继续走吗? 这是一个值得思考的问题。讲一个关于西西弗的故事。

有约束的凸二次优化问题包含:等式约束凸二次优化、不等式约束凸二次优化和混合约束凸二次优化。

4.1 简单有约束优化

4.1.1 简单等式约束

import pprint
from cvxopt import matrix, solvers
P = matrix([[4.0,1.0],[1.0,2.0]])
q = 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])

result = solvers.qp(P,q,G,h,A,b)
print('x\n',result['x'])

result:
pcost       dcost       gap    pres   dres
 0:  1.8889e+00  7.7778e-01  1e+00  3e-16  2e+00
 1:  1.8769e+00  1.8320e+00  4e-02  2e-16  6e-02
 2:  1.8750e+00  1.8739e+00  1e-03  2e-16  5e-04
 3:  1.8750e+00  1.8750e+00  1e-05  1e-16  5e-06
 4:  1.8750e+00  1.8750e+00  1e-07  1e-16  5e-08
Optimal solution found.
x
[ 2.50e-01]
[ 7.50e-01]

4.1.2 简单不等式约束

from cvxopt import matrix, solvers
P = matrix([[1.0, 0.0], [0.0, 0.0]])
q = matrix([3.0, 4.0])
G = matrix([[-1.0, 0.0, -1.0, 2.0, 3.0], [0.0, -1.0, -3.0, 5.0, 4.0]])
h = matrix([0.0, 0.0, -15.0, 100.0, 80.0])
sol=solvers.qp(P, q, G, h)
print(sol['x'])
print(sol['primal objective'])

    pcost       dcost       gap    pres   dres
 0:  1.0780e+02 -7.6366e+02  9e+02  0e+00  4e+01
 1:  9.3245e+01  9.7637e+00  8e+01  6e-17  3e+00
 2:  6.7311e+01  3.2553e+01  3e+01  6e-17  1e+00
 3:  2.6071e+01  1.5068e+01  1e+01  2e-17  7e-01
 4:  3.7092e+01  2.3152e+01  1e+01  5e-18  4e-01
 5:  2.5352e+01  1.8652e+01  7e+00  7e-17  3e-16
 6:  2.0062e+01  1.9974e+01  9e-02  2e-16  3e-16
 7:  2.0001e+01  2.0000e+01  9e-04  8e-17  5e-16
 8:  2.0000e+01  2.0000e+01  9e-06  1e-16  2e-16
Optimal solution found.
[ 7.13e-07]
[ 5.00e+00]
20.00000617311241

看了上面的代码，是不是觉得很简单。难点不在代码，而是在于将实际优化问题转化为标准形式的过程。对于简单的凸二次优化问题可以借助常用的优化包(比如上文代码使用的cvxopt等)解决。但是在很多实际工程问题中(比如机器学习、定价补贴等)，决策变量和约束条件会非常多。此时通用的求解器基本无法满足需要，还是需要我们在充分了解有约束的凸二次优化求解原理的基础上做面向问题的特殊化处理。

参考资料:

1.Python量化投资教程

https://wizardforcel.gitbooks.io/python-quant-uqer/192.html

4.2 有约束凸二次优化原理

4.2.1 Lagrange乘子法

在面对一个棘手的新问题时，科学家们往往会考虑如何把新问题转化成已知的问题。在本文中无约束的凸优化问题是相对比较成熟的，那有约束的凸优化可以转化成无约束凸优化问题吗? 让我们一起来见证奇迹的时刻。

首先来看等式约束: