convex optimization

1. 凸优化问题

在求解问题时，我们常常会遇到求函数 $f: \mathbb{R}^n \rightarrow \mathbb{R}$ 的最小值（或最大值）问题，当该函数属于凸优化问题时，我们能够在多项式时间内求解出问题的答案，因此我们常常拿凸优化问题做文章。下面我们介绍凸优化问题的相关概念。

2. 相关概念

凸集：一个集合 $C$ ，对于任意 $x \in C$ , $y \in C$ ，如果有 $0 \leq \theta \leq 1$ :

$\theta x + (1-\theta)y \in C$

那么称集合C为凸集(Convex Set)。其物理意义是，集合内任意两点连线上的点仍在集合内。对于 $\lambda_1 x_1 + ... + \lambda_m x_m = x$ 称为凸组合（convex combination），其中 $\sum_i^m \lambda_i = 1$ 并且 $\lambda_i \geq 0$ 。

多面体（polyhedra）。或者称为polyhedral set，是指对于矩阵 $A \in \mathbb R^{m \times n}$ , 向量 $b \in \mathbb R ^m$ 。有 $\{x \in \mathbb R^n, Ax \leq b \}$ 。 polyhedra是由多个hyper plane围成的（多个hyper space的交集）。

需要注意，hyper plane和hyper space都是凸集。

当polyhedra被额外限定 $x\geq 0$ 时，其方向有较为特殊的性质：对于 $P=\{x\in \mathbb R:Ax \leq b, x \geq 0\}$ ，其方向 $d$ 满足： $\{Ad \leq 0, d \geq 0, d \neq 0\}$

射线(Ray)，用 $x_0 \in \mathbb R^n$ 表示点， $d \in \mathbb R^n$ 向量表示方向。则用 $x_0$ 和 $d$ 表示的ray，称为ray，记成 $\{x \mid x = x_0 + \lambda d, \lambda \geq 0\}$ 。需要注意，包含ray的图集是unbounded convex set。

一类特殊的unbounded convex set是锥(cone)。对于任意的 $x \in C$ ，如果有 $\lambda \geq 0$ ，则 $\lambda x\in C$ ，则称 $C$ 是cone。需要注意cone未必都是凸集。需要注意，每个cone都包含它的原点。

极点(extreme points)， $C$ 为凸集，如果不存在 $x_1$ , $x_2$ ( $x_1 \neq x_0, x_2 \neq x_0$ ),使得 $x_0 = \lambda x_1 + (1-\lambda)x_2$ ，对于 $\lambda \in (0,1)$ ，则称 $x_0$ 为极点。

对于凸集，其极点都在边界点（boundary points）上。

方向，凸集 $C$ ，如果对于所有的 $x_0 \in C$ ，如果有从 $x_0$ 出发的射线，对于任意 $\lambda \geq 0$ 有，

$\{x:x=x_0 +\lambda d, \lambda \geq 0 \} \subset C$

极方向(extreme direction)，如果一个方向 $d$ ，不能用其他两个不相同的方向 $d_1$ , $d_2$ 表示，则称该方向为极方向。用数学符号表示，对于 $C \subset \mathbb R^n$ 。不存在 $d_1$ , $d_2$ ，使得对于标量 $\lambda_1>0$ ， $\lambda_2>0$ ：

$d = \lambda_1 d_1 + \lambda_2 d_2$

需要注意，一个方向是极方向的充要条件是 $d$ 是多面体 $D$ 的一个极点。

卡拉西奥多里定理(Caratheodory Characterization Theorem)，对于正象限的多面体（线性规划中的可行域），有有限个、且至少又一个极点。或者对于 $n$ 维空间，用最多 $n+1$ 个点就能表示集合内的任意点。

对于非空unbounded多面体 $P$ ，如果 $P$ 有极点 $x_1...x_k$ ，极方向 $d_1...d_l$ ，如果 $x \in P$ ，则存在常数 $\lambda_1...\lambda_k$ 和 $\mu_1...\mu_l$ ，使得：

$x = \sum_{i=1}^{k} \lambda_i x_i + \sum_{j=1}^{l} \mu_j d_j$

凸函数，如果对于函数 $f: \mathbb{R}^n \rightarrow \mathbb{R}$ ，其定义域为图集，并且对于 $0\leq \theta \leq 1$ ：

$f(\theta x + (1-\theta) y) \leq \theta f(x) + (1-\theta) f(y)$

并且对于 $0\leq \theta \leq 1$ , $x_1 \neq x_2$ ，有：

$f(\theta x + (1-\theta) y) < \theta f(x) + (1-\theta) f(y)$

一些性质：1.如果 $f(x)$ 是凹函数(concave function)，则 $-f(x)$ 是凸函数(convex function)。2. $f_1(x)，f_2(x)$ 均为凸函数，则对于 $\forall \lambda_1,\lambda_2>0$ ，有 $\lambda_1 f_1(x) + \lambda_2 f_2(x)$ 仍是凸函数; $Max(f_1(x),f_2(x))$ 仍是凸函数。并且上述两个定义均可以推广到有限个凸函数的组合。

对于非空集合 $S$ ,以及函数 $f:S\rightarrow$ ： $x^* \in R^n, d\in R^n, d \neq 0$ ，对任意的 $\lambda>0$ 有 $x^* +\lambda d \in S$ ，则定义方向导数(dirctional derivative): $f^{'}(x^{*};d) = \lim_{\lambda>0+}\frac{f(x^{*}+\lambda d) - f(x^{*})}{\lambda}$

对于非空集合 $S$ ,以及函数 $f:S\rightarrow$ ， $f$ 的上镜图(epigraph) $epi f$ 是 $R^{n+1}$ 的子集: $\{(x,y), x\in S, y\in R, y \geq f(x)\}$ ，下镜图为 $<$ 。

性质： $f$ 是凸函数的充要条件是：上镜图为凸集。

次梯度(subgradients): 凸集 $S$ ，凸函数 $f$ ，凸函数在集合 $S$ 内的一个点 $x^{*}$ 处的次梯度: $\xi$ :

$f(x) \geq f(x^{*}) + \xi^{T}(x-x^{*})$

注意到， $\Delta_x^2 f(x) \geq 0$

3. Linear Programming

对于LP问题，其标准形式表达如下:

$min\ C^{T}X$ $s.t. AX =b$ $X \geq 0$

定理1，对于存在direction的polyhedra，其存在a finite optimal solution的条件: $C^{T}d_{j} \geq 0$ ，其中 $d_j$ 是extreme direction。

Algorithm1，（原始）单纯形法:

找到一组BFS
判断其是否最优
如果是,停止迭代；否则，调整BFS,继续迭代。

由 $AX=b$ ，可以得到:

$(B\ N） \begin{bmatrix} X_B \\ X_N \end{bmatrix} = b$

即 $BX_B + NX_N = b$ , 由于 $B$ 可逆，故 $X_B = B^{-1}b - B^{-1}NX_N$ 。
$z = CX = (C_B C_N) \begin{bmatrix} X_B \\ X_N \end{bmatrix} = C_BX_B +C_NX_N = C_B(B^{-1}b - B^{-1}NX_N) + C_NX_N = C_{B}B^{-1}b + (C_{N} - C_{B}B^{-1}N)X_N$

其中 $C_{N} - C_{B}B^{-1}N$ 作为进基检验数；出基检验数为 $\frac{B^{-1}b}{B^{-1}p_i}$ , 其中 $p_i$ 为进基列；最小的出基检验数对应原本基的行，出基；

4. Non-linear Programming

NLP, 非线形规划，这里的非线性既可以是目标函数（比如垄断市场下企业利润的函数：价格乘以产量，其中价格又是产量的函数，故利润是产量的二次函数），也可以是约束(无约束，线性约束，二次约束，凸约束，非凸约束等)；

可行方向(Feasible direction)，对于点 $x \in \Omega$ , 如果存在 $a^*$ ，有 $\forall a \in(0,a^*)$ , $x+ad \in \Omega$ , 则称d是可行方向。由定义可知，在集合 $\Omega$ 内的点，任意方向都是可行方向；而对于boundary上的点，并不是所有的方向都是可行方向。

一阶导数必要条件（是局部或全局最优解的必要条件），对于任意的可行方向 $d$ ：

1. 对于有约束问题(约束 $\Omega \in \mathbb{R}^n$ )： $\Delta f(x^*)d \geq 0$
2. 对于无约束问题(约束 $\mathbb{R}^n$ )： $\Delta f(x^*) = 0$

二阶导数必要条件（是局部或全局最优解的必要条件）：

1. 对于有约束问题(约束 $\Omega \in \mathbb{R}^n$ )：1） $\Delta f(x^*)d \geq 0$ 或者 2) if $\Delta f(x^*)d = 0, then, \ d^T \Delta^2f(x^*)d \geq 0$ ，其中 $\Delta^2f(x^*)$ 叫做海森矩阵 $H$ 。

2. 对于无约束问题(约束 $\mathbb{R}^n$ )：1） $\Delta f(x^*) = 0$ 或者 2) $\forall d: \ d^T \Delta^2f(x^*)d \geq 0$ (即海森矩阵是半正定的)。

特别的，对于无约束问题，二阶导数充分条件（是局部或全局最优解的充分条件）：1) $\Delta f(x^*)d \geq 0$ , 2) $H(x^*)$ 是正定 $^{[1]}$ 的。

我们在这篇文章在讲凸优化（线性规划的约束条件是凸的，目标函数是线性函数，也是凸函数；对于非线形规划，其未必是凸优化问题），我想对于非线形规划一样，凸的意义一样只是值得研究；这里很重要的一点是：如果目标函数是凸的，约束条件也是凸的（或者无约束），那么它具有一些典型的性质。在分析凸优化问题的典型性质之前，我们先讲一下，怎么确定一个问题是凸优化问题（约束是凸的，目标函数是凸函数），约束是否为凸比较容易验证，而目标函数是否为凸函数，其验证方法有多种：
1. 根据定义: $f(\theta x + (1-\theta) y) \leq \theta f(x) + (1-\theta) f(y)$
2. 一阶充要条件： $f(x_2) \geq f(x_1) + \Delta f(x_1)(x_2 - x_1)$
3. 二阶充要条件：Hessian矩阵半正定。如果Hessian矩阵是正定的，则该函数是严格凸函数。

如果确定一个问题属于凸优化问题，那么其具有一下性质：
1. Global optimal solution:任何的局部最优解都是全局最优解（这种全局最优点不一定唯一）。根据我们之前提到的一阶、二阶必要、充分条件，我们可以求得一些局部最优解，可以通过这里的性质来验证其是否是全局最优解；
2. 如果目标函数是严格凸函数，则全局最优点是唯一的。

到此为止，我们似乎已经给出了解决无约束问题的方法：令其梯度等于0，得到平稳点 $x^*$ ；然后使用充分条件（ $H(x^*)>0$ ）进行验证。这样就得到的局部最优解，再证明其是全局最优解或通过比较局部最优解来得到全局最优解。但是很多时候下，我们令梯度等于0，是难以解出 $x^*$ 的。在这种背景下，我们常常采用迭代的方法求解问题，下面对这种方法展开介绍, 下降迭代法:
下降迭代法,通过算法A使得： $x^{(k+1)} = A(x^{(k)})$ ，在该问题中如何设计算法A,即确定下降方向 $d_k$ 和步长 $\lambda$ ；其中对于步长确定方法的不同可以区分出几种不同的方法，其中一种方法是使目标函数值下降最多：

$\min_{\lambda_k} \ f(x^{(k)} + \lambda_k d_k)$

上式的优化问题是以 $\lambda_k$ 为变量的一元函数 $f(x^{(k)} + \lambda_k d_k)$ 的最小值，这一过程一般称为一维搜索。

一维搜索的方法包括“试探法”（斐波那契法、黄金分割法）、“曲线拟合法”（牛顿法、false position法、cubic、Quadratic法等）。
其中试探法的思想是，每次在区间 $(a_i,b_i)$ 内部取两个点位 $(a_{i+1},b_{i+1})$ ，分别计算 $f(a_{i+1}),f(b_{i+1})$ ,对于下单峰函数，如果 $f(a_{i+1})<f(b_{i+1})$ , 则最小值点必在 $(a_i, b_{i+1})$ 之间，否则在 $(a_{i+1},b_i)$ 内，这样迭代下去。
曲线拟合法的思想是，下面介绍曲线拟合法中的牛顿法：

对于一维方程f(x)(这里的x可以理解上下降迭代法中的 $\lambda$ )，用二阶泰勒展开来表示f(x):

$q(x) = f(x_k) + f'(x_k)(x-x_k) + \frac{1}{2}f''(x_k)(x-x_k)^2$

let $q'(x_{k+1})=0$ , 可以得到:

$x_{k+1} = x_k - \frac{f'(x_k)}{f''(x_k)}$

泰勒二阶展开其实是将 $x_k, x_{k+1}$ 两个相邻变量通过一阶导数、二阶导数联系起来，这里对 $q(x)$ 求导，是为了得到极值点。对上式进行迭代，可以得到一维搜索的最优值，牛顿法是超线性收敛的。由上述表达式可知牛顿法涉及二阶导数矩阵（海森矩阵，矩阵的维度由参数的数目决定，而大部分“复杂的”机器学习模型的参数十分多）的求逆，该运算导致每轮迭代计算消耗巨大。这里引入一种能够以较低的计算代价来近似海森矩阵的逆矩阵的方法：拟牛顿法。

拟牛顿法(Quasi-Newton)的参数更新方程与牛顿法一致，区别在于其采用校正公式(Quasi Update)来更新下一轮需要用的海森矩阵的逆，该校正公式的输入是：1.当前轮次的海森矩阵的逆、2.参数更新的幅度 $s_{n+1} = x_{n+1}-x_{n}$ (该更新是由学习速率、当前轮次的海森矩阵的逆，以及 $x_{n}$ 算出来的)，3.一阶导数的变化幅度 $y_{n+1} = \delta f(x_{n+1}) - \delta f(x_n)$ 。

下面介绍拟牛顿法的推导过程：

首先由二阶泰勒展开可以得到：

$f(x_k + \Delta x) = f(x_k) + g_k \Delta x^T + \frac{1}{2}\Delta x ^T H_k \Delta x$

其对 $x_k$ 的导数为，并类推可得第二个式子：

$f'(x_k + \Delta x) = g_k$ $f'(x_{k-1} + \Delta x) = g_{k-1}$

故相减法，并由中值定理可得下式：

$f'(x_k + \Delta x)-f'(x_{k-1} + \Delta x) = g_k-g_{k-1}\\ H_k(x_k - x_{k-1}) = g_k - g_{k-1}\\ H_k s_k = y_k\\ H_k^{-1} y_k = s_k$

基于上面的基础，我们介绍BFGS和LBFGS，前者全称Broyden–Fletcher–Goldfarb–Shanno，是有四位作者的名字命名的，而LBFGS则是加了前缀:Limited-memory。

BFGS模型给出了三个条件：1，满足上述 $H_k^{-1} y_n = s_n$ 的性质，2，满足对称性，3在满足上述约束的所有结果中，我们取变化幅度最小的那个。因此可以建立一个有约束的优化问题：

$min ||H^{-1}_{k} - H^{-1}_{k-1}||^2\\ s.t. \ H^{-1} y_k = s_k\\ H^{-1} is \ Symmetric$

而该优化问题的解，即为我们需要的海森矩阵的逆的校正公式（或者可以称之为更新方程）：

$H_{k}^{-1} = (I - \rho_k y_k s_k^T) H_{k-1}^{-1} (I - \rho_k s_k y_k^T) + \rho_k s_k s_k^T$

其中 $\rho_n = (y_k^Ts_k)^{-1}$

有了上述一维搜索的基础，我们可以介绍几种常见的用于求解无约束问题的下降迭代法：

最速下降法，下降方向取负梯度方向： $-\Delta f(x^{(k)})$ , $\lambda$ 的取值使用上述的一维搜索法。

有时我们认为最速下降法比较慢，这是因为只有在 $x$ 非常接近局部最优解时，这种“最速”才是有效的，当其远离最优解时，并不是那么有效，基于此，我们提出了几种调整策略，一种策略是调整方向从: $-\Delta f(x_k)$ 调整到 $-\Delta f(x_k) +g$ ；另一种调整是从 $-\Delta f(x_k)$ 到 $-D\Delta f(x_k)$ 。依照这几种调整策略，几种经典的算法是：共轭梯度法（conjugate gradient menthod），变尺度法（Scaling）。

共轭梯度法，如果存在矩阵 $Q$ ，并且 $Q$ 是正定的，使得 $d_1^TQd_2=0$ ，则称向量 $d_1,d_2$ 是共轭的。对于向量组 $d_0...d_{n-1}$ 是Q共轭的，即 $d_i^T Qd_j =0(i \neq j)$ 。满足这种共轭条件的向量组之间是相互独立的，即 $\sum_{i=0}^{n-1} \alpha_id_i=0$ ，仅在 $\alpha_i = 0, \forall i=0..n-1$ 下成立。这里有 $x^* = \sum_{i=0}^{n-1} \alpha_id_i$ , 故 $\alpha_i = \frac{d_i^TQx^*}{d_i^TQd_i} = \frac{d_i^Tb}{d_i^TQd_i}$ (对于 $min \frac{1}{2}x^TQx-b^Tx$ )。基于以上内容构造迭代公式：

$x_{k+1} = x_k + \alpha_k d_k$ $\alpha_k = -\frac{g_k^Td_k}{d_k^TQd_k}\\ g_k = Qx_k -b\\d_{k+1}=-g_{k+1}+\beta_kd_k\\\beta_k=\frac{g_{k+1}^TQd_k}{d_k^TQd_k}$

至此，我们介绍了无约束优化问题的求解方法。而对于非线性规划，很有可能是有约束的，此时，常用的策略就是将其转换为非约束问题或线性规划问题。下面我们介绍有约束优化非线形规划（Constraints NLP）的两种常见形式：等式约束和非等式约束。

等式约束：假设对于该问题共有 $h_1(x) = 0$ .. $h_m(x) = 0$ ，总计 $m$ 个等式约束。对于等式约束约束，我们定义正则点（regular point）: $h(x^*)=0$ (对于非等式约束，其正则点， $h(x^*)=0, g(x^*)\leq=0$ )，

等式约束的一阶必要条件： $x^*$ 是等式约束的正则点，则 $x^*$ 也是局部最优点的必要条件：

$\Delta f(x^*)+\lambda^T\Delta h(x^*)=0$

等式约束的二阶必要条件：对于Lagrange function: $L(x,\lambda)=f(x)+\lambda^Th(x)$ 。 $x^*$ 是等式约束的正则点，则 $x^*$ 也是局部最优点的必要条件：

$L(x^*) = F(x^*) + \lambda^TH(x^*) \geq0 \ (postive \ semi-definite)$

其中， $F(x^*)$ 等是Hessian矩阵。等式约束的二阶充分条件是上式中 $L(x^*)$ 是正定的。

非等式约束（该约束即包括等式约束也包括非等式约束）， $min \ f(x):h(x)=0,g(x) \leq 0$ 。非等式约束的正则点是则满足其等式约束和非等约束约束，且部分非等式约束取等式约束的点。

非等式约束的一阶必要条件（First order necessary conditions or KKT conditions）:

$\Delta f(x^*) + \lambda^T \Delta h(x^*) + u^T \Delta g(x^*) = 0\\u^Tg(x^*) =0\\u\geq0$

非等式约束的二阶必要条件（Second order necessary conditions）:

$x^*$ is a local minimum point and a regular point, then:

$\Delta f(x^*) + \lambda^T\Delta h(x^*) + \mu^T\Delta g(x^*) = 0\\ u^Tg(x^*)=0\\L(x^*) = F(x^*) + \lambda^TH(x^*)+u^TG(x^*)\geq 0$

二阶充分条件，上面的半正定条件改为正定。

上面提供了通过必要条件来求解约束非线形规划问题的思路，类似于之前无约束问题，有时令一阶导等于0，难以解出局部最优解。因此，在此基础上，对于有约束的非线形规划问题，可以使用:可行方向法、惩罚和障碍法

惩罚和障碍法：对于带约束的非线形规划问题：

$min \ f(x) \\ \Omega:\{x:g(x) \leq 0, x\in R^n\}$

思想是移除约束，并把其加到目标函数中，根据方法的不同，可以分为：惩罚(penalty)、障碍(barrier)。惩罚方程：如果满足约束，该方程为0；如果不满足约束该方程为正。障碍方程：如果满足约束，该方程为正；如果该点接近边界(boundary)方程趋向于正无穷。总地来说就是，惩罚法是加大惩罚系数c，是得解从外部惩罚到内部；障碍法是，一直保留在内部，防止出去。

外部惩罚法(Exterior Penalty function method)，当满足约束时，P(x) = 0，对于 $g(x) \leq 0$ 这样的约束，一般惩罚函数可以定义为 $c(\sum_i(max\{0, g_i(x)\})^2)$ ：

$min \ f(x) + cP(x)$

内部障碍法，当满足约束时，B(x)为正，当接近于边界(Boundry)时，B(x)趋向于正无穷。对于 $g(x) \leq 0$ 这样的约束，一般障碍函数可以定义为 $-log(-g(x))$ 。

5. 实践中的数值计算

对于比较简单的优化问题（比如一般形式的线性回归），我们可以通过解析法来求得最优值；而大部分的问题，只能通过一些数值计算的方法来求解。而迭代法是数值计算中最常见的方法。

梯度下降法中的更新策略有多种，下面给出常见的几种：

一般形式的更新：参数 $\theta$ ，以 $\lambda$ 为学习率（步长），以 $d(\theta_{i-1})$ 为更新方向。

$\theta_{i} = \theta_{i-1} - \lambda d(\theta_{i-1})$

Momentum：为了避免更新过程中的大幅度震动，引入动量的概念，每次更新时同时受上一轮梯度、本次梯度的方向的共同影响(可以拆分为走了两步)。

$v_i = \mu v_{i-1} - \lambda d(\theta_{i-1})\\ \theta_{i} = \theta_{i-1} + v_i$

Nesterow Accelerated Gradient(NAG or Nesterow Momentum):Momentum可以理解为走了两本，且两步都是以期初的梯度作为方向，如果对其进行改进，将第二步的梯度改为走完第一步之后的梯度，则可以提高一定的收敛速度：

$(\theta_{i-1})' = \theta_{i-1} - \mu v_{i-1}\\ v_i = \mu v_{i-1} - \lambda d((\theta_{i-1})')\\ \theta_{i} = \theta_{i-1} + v_i$

AdaGrad：除了上述改进迭代方向的更新方程外，还可以改进学习速率 $\lambda$ ，下式中 $\epsilon$ 是非常小的常数，其目的是对分母进行平滑，防止其值为0。

$r = r + (d(\theta_{i-1}))^2\\ \theta_{i} = \theta_{i-1} - \frac{\lambda}{\sqrt{r} + \epsilon}d(\theta_{i-1})$

Adam更新方程则是将Momentum和AdaGrad进行了结合。

注：
[1],一个矩阵是正定的，其充分必要条件是该矩阵左上角各阶主子式都大于0（i=1..到n的行列式大于0）。