数据科学和机器学习中的优化理论与算法（下）

数据科学和机器学习中的优化理论与算法（下）

数据科学和机器学习当前越来越热，其中涉及的优化知识颇多。很多人在做机器学习或者数据科学时，对其中和优化相关的数学基础，包括随机梯度下降、ADMM、KKT 条件，拉格朗日乘数法、对偶问题等，不太了解，一知半解地用，用着用着就出错了。

本文希望从基础知识的角度，尽可能全地对最简单的优化理论和算法做一个小结。内容涵盖以下几个方面：优化简介、无约束优化、线搜索方法、信赖域方法、共轭梯度方法、拟牛顿方法、最小二乘问题、非线性方程、约束优化理论、非线性约束优化算法、二次规划、罚方法…

通过本文档的学习，相信你会掌握数据科学和机器学习中用到的优化基础知识，以后再遇到优化问题，就不会再困惑了。

建议收藏，方便时时查阅。

共轭梯度方法

线性共轭梯度方法

共轭梯度（Conjugate Gradient，简称 CG）方法，最开始在十九世纪五十年代被 Hestenes 和 Stiefel 提出用来迭代地带正定系数的线性方程组。假设 $A$ 对称正定，我们要求解，
$$Ax=b$$
这也可以被等价写为极小化下面的一个凸问题：
$$\min \varphi (x)\equiv \frac{1}{2}{x}^{T}Ax-b^{T}x$$
$\varphi$的梯度也就等于线性系统的残差，
$$\nabla \varphi (x)=Ax-b\equiv r(x)$$
在迭代到$x=x_k$，当前步的残差也可以写为：$r_k=A{x}_k-b$。

共轭方向法

给一个定义，如果对于一个非零向量集 { p 0 , p 1 , ⋯   , p t } \left\{p_{0}, p_{1}, \cdots, p_{t}\right\} {p0​,p1​,⋯,pt​}，对于正定矩阵 $A$ 满足，
$$p_i^{T}A{p}_j=0,\text{ for all }i\ne j$$
我们就说这个向量集是共轭的，其实也就是在 $A$ 意义下的一个正交。

共轭方向之间是相互独立的，任何一个向量可以写成它们之间的一个线性组合 $x={\sum }_{i=1}^n{\alpha }^i{p}_i$，代入到极小化函数 $\varphi$ 中，可得，
$$\varphi (x)=\sum _{i=1}^n\frac{{\left({\alpha }^i\right)}^2}{2}{p}_i^{T}A{p}_i-{\alpha }^i{p}_i^{T}b$$
要极小化 $\varphi$，显然有，
$${\alpha }^i=\frac{p_i^{T}b}{p_i^{T}A{p}_i}$$
也就是说，只要知道了共轭方向，就可以通过他们的线性组合，将最优解表达出来。细细地思考，将这些共轭方向做线性组合，无非就是从原点出发，一次沿一个共轭方向去走一定的步长。那么，不从原点出发，而从某个 $x_0$ 出发，每次也依次沿着共轭方向走，是否也能走到一个最优值点呢？答案是肯定的。

考虑下述的共轭方向方法，给定初值点 $x_0\in {\Re }^n$ 和一个共轭方向集合 { p 0 , p 1 , ⋯   , p n − 1 } \left\{p_{0}, p_{1}, \cdots, p_{n-1}\right\} {p0​,p1​,⋯,pn−1​}，让我们按如下方式进行迭代
$$x_{k+1}=x_k+{\alpha }_k{p}_k$$
这里的 ${\alpha }_k$ 其实是二次函数 $\varphi$ 沿着 $x_k+\alpha p_k$ 方向的一维极小值，显示地给出表达，就是
$${\alpha }_k=-\frac{r_k^T{p}_k}{p_k^{T}A{p}_k}$$
有相关定理表明，以上迭代过程收敛到最优解 $x^{\ast }$，并且具有n步终止性。以上方法，就是传说中的共轭方向法。

那么，问题来了，如何给出一组共轭方向呢？一个容易想到的是，直接取矩阵 $A$ 的 n 个彼此正交的特征向量作为共轭方向。问题在于，实现构造共轭方向，计算量太大，从而可以考虑在迭代中产生。事实上，迭代地去产生和利用共轭方向，就是共轭梯度方法。

共轭梯度法

共轭梯度法在迭代中生成共轭方向，只需要用到前一点出的搜索方向，在内存和计算的开销上，也尤为经济。共轭方向的产生，是将当前步的共轭方向，写成当前步的残差和前一步的共轭方向（搜索方向）的线性组合，
$$p_k=-r_k+{\beta }_k{p}_{k-1}$$
由共轭梯度的定义 $p_k$ 必须满足 $p_{k-1}^{T}A{p}_k=0$，易得，
$${\beta }_k=\frac{r_k^{T}A{p}_{k-1}}{p_{k-1}^{T}A{p}_{k-1}}$$
选择 $p_0$ 是负梯度方向，并且依次执行沿各个共轭梯度方向一维的极小化过程，就是所说的共轭梯度方法。算法流程如下图。

CG 算法具有 n 步终止性，它也是也给 Krylov 子空间方法。为了写成更加经济的程序，我们可以对上述算法的某些变量表达进行化简。

首先，我们得不加证明地引用残差和共轭梯度的的正交性质。
$$r_k^T{r}_i=0,\forall i\ne k$$
$$r_k^T{p}_i=0,\forall i\ne k$$
那么，根据
$$r_k^T{p}_k=r_k^T\left(-r_k+{\beta }_k{p}_{k-1}\right)=-r_k^T{r}_k$$
以及
$$r_{k+1}=A{x}_{k+1}-b=A\left(x_k+{\alpha }_k{p}_k\right)-b=r_k+{\alpha }_{k}A{p}_k$$
我们可以得到，
$${\alpha }_k=-\frac{r_k^T{p}_k}{p_k^{T}A{p}_k}=\frac{r_k^T{r}_k}{p_k^{T}A{p}_k}$$
和
$${\beta }_{k+1}=\frac{r_{k+1}^{T}A{p}_k}{p_k^{T}A{p}_k}=\frac{r_{k+1}^T\left(r_{k+1}-r_k\right)}{p_k^T\left(r_{k+1}-r_k\right)}=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$$
那么，我们就可以得到一个 CG 算法的实用形式，编程参考这个来写，就是极好了。如图。

定理表明，如果 $A$ 有 $r$ 个不同的特征值，那么 CG 算法至多迭代 $r$ 步终止。线性 CG 方法的收敛速度为，
$${\Vert x_k-x^{\ast }\Vert }_A\le 2{\left(\frac{\sqrt{\kappa (A)}-1}{\sqrt{\kappa (A)}+1}\right)}^k{\Vert x_0-x^{\ast }\Vert }_A$$
这里的 $\kappa (A)=\Vert A{\Vert }_2{\Vert A^{-1}\Vert }_2=\frac{{\lambda }_n}{{\lambda }_1}$。

从收敛率上可以看得出，当 $A$ 的最大最小特征值差距比较大的时候，收敛就会很慢，怎么办？可以用预优方法，也就有了 PCG。所谓的预条件（Preconditioning）CG 方法，就是对变量先做一个变换，$\hat{x}=Cx$，极小化问题变为
$$\hat{\varphi }(\hat{x})=\frac{1}{2}{\hat{x}}^T\left(C^{-T}A{C}^{-1}\right)\hat{x}-{\left(C^{-T}b\right)}^T\hat{x}$$
等价于求解
$$\left(C^{-T}A{C}^{-1}\right)\hat{x}=C^{-T}b$$
我们希望能找到一个合适的 $C$，使得 $C^{-T}A{C}^{-1}$ 的条件数比较小，这就是预条件 CG 方法。

非线性共轭梯度方法

非线性共轭梯度方法用来求解非线性问题：
$$\min f(x)$$
非线性 CG 基本上就是仿这线性 CG 来的。因为我们没有了矩阵 $A$，所以原有 CG 的 ${\alpha }_k$ 的表达就没法用了，我们用原有思想，使用最小化 $\varphi$ 沿着 $p_k$ 对应的步长作为 ${\alpha }_k$。另外，残差 $r$，在非线性 CG 中也没有了，我们用 $f$ 的梯度来替代，于是，就有了如图所示的非线性 CG 方法，根据 $\beta$ 的选择，此时也叫 FR 方法。

在相同的框架下，我们修改 ${\beta }_{k+1}$，就可以得到不同的非线性 CG 方法。
PR 方法：
$${\beta }_{k+1}^{PR}=\frac{\nabla f_{k+1}^T\left(\nabla f_{k+1}-\nabla f_k\right)}{{\Vert \nabla f_k\Vert }^2}$$
DY 方法（戴-袁）：
$${\beta }_{k+1}^{DY}=\frac{\nabla f_{k+1}^T\nabla f_{k+1}}{{\left(\nabla f_{k+1}-\nabla f_k\right)}^T{p}_k}$$
HS 方法：
$${\beta }_{k+1}^{HS}=\frac{\nabla f_{k+1}^T\left(\nabla f_{k+1}-\nabla f_k\right)}{{\left(\nabla f_{k+1}-\nabla f_k\right)}^T{p}_k}$$
观察可发现，这四种方法只不过是两个分子和两个分母不同的组合而已。

当然，这当中还有一些别的有趣的问题，包括重启动、二次终止性、全局收敛性和信赖域子问题的截断 CG 方法等等，由于时间关系，就不再多说。

拟牛顿方法

拟牛顿和割线方程

拟牛顿法是一种以牛顿法为基础设计的，求解非线性方程组或连续的最优化问题函数的零点或极大、极小值的算法。当牛顿法中所要求计算的雅可比矩阵或 Hessian 矩阵难以甚至无法计算时，拟牛顿法便可派上用场。

考虑模型问题
$$m_k(p)=f_k+\nabla f_k^{T}p+\frac{1}{2}{p}^T{B}_{k}p$$
最小化之，可得拟牛顿步长：
$$p_k=-B_k^{-1}\nabla f_k$$
那么，新的迭代为：
$$x_{k+1}=x_k+{\alpha }_k{p}_k$$
问题是，这里的${\alpha }_k$ 和$B_k$如何选取？一般来说，对于 Hessian，以下这种近似关系是要满足的。
$${\nabla }^2{f}_{k+1}\left(x_{k+1}-x_k\right)\approx \nabla f_{k+1}-\nabla f_k$$
因为$B_k$是 Hessian 的逼近，所以需要满足"割线方程"：
$$B_{k+1}{s}_k=y_k$$
这里 $s_k=x_{k+1}-x_k={\alpha }_k{p}_k,y_k=\nabla f_{k+1}-\nabla f_k$。事实上，只要满足曲率条件 $s_k^T{y}_k>0$，那么割线方程总是会被满足的，而曲率条件总是被Wolfe或者强Wolfe条件所保证。

割线方程的解是不唯一的，因为未知数太多了。所以要确定 $B_k$，我们还需要一些额外的条件，不同的附加条件，就对应了不同的方法。

DFP 方法

我们在割线方程的约束框架下，求解子问题：
$$\begin{array}{ll}{\min }_B& \Vert B-B_k\Vert \\ \text{ s.t. }& B=B^T,B{s}_k=y_k\end{array}$$
得到了 $B_k$ 的迭代：
$$\begin{array}{c}{B}_{k+1}=\left(I-{\gamma }_k{y}_k{s}_k^T\right)B_k\left(I-{\gamma }_k{y}_k{s}_k^T\right)+{\gamma }_k{y}_k{y}_k^T\\ {\gamma }_k=1/y_k^T{s}_k\end{array}$$
显然，在应用拟牛顿方法时，我们没有知道用到 $B_k$，而是用到它的逆 $H_k=B_k^{-1}$，利用关系（Sherman-Morrison-Woodbury公式）：
$${\left(A+a{b}^T\right)}^{-1}=A^{-1}-\frac{A^{-1}a{b}^T{A}^{-1}}{1+b^T{A}^{-1}a}$$
我们可以由 $B_k$ 的迭代得到 $H_k$ 的迭代：
$$H_{k+1}=H_k-\frac{H_k{y}_k{y}_k^T{H}_k}{y_k^T{H}_k{y}_k}+\frac{s_k{s}_k^T}{y_k^T{s}_k}$$
将这个迭代套到拟牛顿方法中，就叫DFP方法。

BFGS 方法

直接求解 $B_k$ 逆 $H_k$ 的子问题：
$$\begin{array}{ll}{\min }_H& \Vert H-H_k\Vert \\ \text{ s.t. }& H=H^T,H{y}_k=s_k\end{array}$$

直接得到 $H_k=B_k^{-1}$ 的迭代，称为 BFGS 方法：
$$\begin{array}{c}{H}_{k+1}=\left(I-{\gamma }_k{s}_k{y}_k^T\right)H_k\left(I-{\gamma }_k{s}_k{y}_k^T\right)+{\gamma }_k{s}_k{s}_k^T\\ {\gamma }_k=1/y_k^T{s}_k\end{array}$$

对称-秩一方法

当然，还有其他的一些 $H_k$ 更新方式，如图。

其中的对称-秩一（SR1）方法和 Broyden 族方法都比较有名，秩一更新依然能保证 $B_k$ 和 $H_k$ 的对称性以及能保持对割线方程的满足，不同的是，SR1 更新，不再保证其正定性。数值效果依然很好，具体的推导就不说了。

步长的选取

步长 ${\alpha }_k$ 如何选取？只要保证 Wolfe 条件，就能保证其收敛性。取精确线搜索步长时，其值为 1。

只列出 BFGS 基本框架如图。

在迭代的过程中，我们不需要显性地去求 $B_k$，我们用到的只是它的逆，所以，我们可以用其逆 $H_k$ 的迭代公式进行搜索迭代。给一个简单的代码示例：

clc
clear

f = @(t)t(1)^2+2*t(2)^2;
x0 = [10,10]';
epsilon = 0.001;
H0 = eye(2);
method = 'BFGS';%DFP
x = quasi_Newton(f,x0,epsilon,H0,method);

function [xk,k] = quasi_Newton(f,x0,epsilon,H0,method)
%使用：quasi_Newton(f,x0,method)

if nargin < 3
    help mfilename;
end
k = 0;
syms t1 t2;
t = [t1,t2]';
fs = f(t);
dfs = gradient(fs);
df = matlabFunction(dfs);
df = @(x) df(x(1),x(2));
df0 = df(x0);
normdf = sqrt(df0'*df0);
H = H0;
xk = x0;
dfk = df0;
while normdf > epsilon
    p = -H*dfk;
    alpha = cal_alpha(H,dfk);
    xk1 = xk + alpha*p;
    dfk1 = df(xk1);
    sk = xk1 - xk;
    yk = dfk1 - dfk;
    eval(['H = ' method '(H,sk,yk);']);
    %H = BFGS(H,sk,yk);
    k = k + 1;
    xk = xk1;
    dfk = dfk1;
    normdf = sqrt(dfk'*dfk);
    xk
end
end
function alpha = cal_alpha(H,dfk)
    %alpha = dfk'*H*dfk/(dfk'*H'*dfk);
    alpha = 1;
end
function H = BFGS(H,sk,yk)
    gammak = 1/(yk'*sk);
    skykT = sk * yk';
    skskT = sk * sk';
    E = eye(2);
    H = (E-gammak*skykT)*H*(E-gammak*skykT) + gammak*skskT;
end
function H = DFP(H,sk,yk)
    Hyk = H*yk;
    ykTsk =  yk'*sk;
    skskT = sk * sk';
    H = H - (Hyk*yk'*H)/(yk'*Hyk) + skskT/ykTsk;
end

其他内容

关于拟牛顿方法，还有狠多有趣的内容，包括但不仅限于其二次终止性、非对称 Rank-1 方法、BFGS 收敛率、BFGS 和 CG 的一个关系等，下面简要提两个。

BB 方法

在牛顿方法中，我们取 $B_k$ 矩阵为单位矩阵，就得到了 BB 方法：
$$x_{k+1}=x_k-D_k^{-1}{g}_k$$
其中，$D_k\equiv {\alpha }_{k}I$，这本质上也是梯度下降方法。步长的选取要让 $D_k$ 尽可能地满足割线方程 $B_k{s}_{k-1}=y_{k-1}$，求解
$$\underset{D=\alpha I}{\min }{\Vert D{s}_{k-1}-y_{k-1}\Vert }_2$$
可得，所谓的BB步长，
$${\alpha }_k^{BB}=\frac{s_{k-1}^T{s}_{k-1}}{s_{k-1}^T{y}_{k-1}}$$
对于凸函数 $f(x)=\frac{1}{2}{x}^{T}Hx+b^{T}x$，BB 步长就是
$${\alpha }_k^{BB}=\frac{g_{k-1}^{T}H{g}_{k-1}}{g_{k-1}^T{g}_{k-1}}=\frac{1}{{\alpha }_{k-1}^{SD}}$$

有限内存 BFGS 方法（L-BFGS）

对于大规模问题，BFGS 方法因其逼近矩阵 $B_k$ 和 $H_k$ 都是稠密的，那么必然消耗比较大。为了解决这个问题，工程上用得比较多的是有限内存（Limited-Memory）的 BFGS 方法，简称 L-BFGS。L-BFGS 的想法比较简单，就是在 $H_k$ 迭代的时候，我们只用到最近几步的曲率信息，而把之前的信息丢弃，然后来执行 BFGS 过程。也就是说，反正 BFGS 算着算着，误差累积，导致 $H_k$ 反正都差了精度，还不如丢弃了，重新从初值开始。这样还能降低 $H_k{g}_k$ 的计算量。L-BFGS 方法用得比较多。

最小二乘问题

问题描述

最小二乘问题，描述的是求解形如
$$\min f(x)=\frac{1}{2}\sum _{j=1}^m{r}_j^2(x)$$
的问题。

如果把 $r(x)$ 写成
$$r(x)={\left(r_1(x),r_2(x),\cdots ,r_m(x)\right)}^T$$
那么最小二乘问题可以简洁地写为：$f(x)=\frac{1}{2}\Vert r(x){\Vert }_2^2$。

导数

$r$ 雅克比矩阵的雅克比矩阵定义为：
$$J(x)=\left[\frac{\partial r_j}{\partial x_i}\right]=\left[\begin{array}{c}\nabla r_1(x{)}^T\\ \nabla r_2(x{)}^T\\ ⋮\\ \nabla r_m(x{)}^T\end{array}\right]$$
如果
$$f=\frac{1}{2}\sum _{j=1}^m{r}_j^2(x)=\frac{1}{2}\Vert r(x){\Vert }^2$$
那么梯度和 Hessian 矩阵可以表达为：
$$\begin{array}{rl}\nabla f(x)& =\sum _{j=1}^m{r}_j(x)\nabla r_j(x)=J(x{)}^{T}r(x)\\ {\nabla }^{2}f(x)& =\sum _{j=1}^m\nabla r_j(x)\nabla r_j(x{)}^T+\sum _{j=1}^m{r}_j(x){\nabla }^2{r}_j(x)\\ & =J(x{)}^{T}J(x)+\sum _{j=1}^m{r}_j(x){\nabla }^2{r}_j(x)\end{array}$$

线性最小二乘问题

对于线性最小二乘问题，
$$f(x)=\frac{1}{2}\Vert Jx-y{\Vert }_2^2$$
$f$ 的一阶导和二阶导为：$\nabla f(x)=J^T(Jx-y),{\nabla }^{2}f(x)=J^{T}J$。极小值点要满足一阶必要性条件 $\nabla f\left(x^{\ast }\right)=0$，那么原问题就等价于求解所谓的正规方程：
$$J^{T}J{x}^{\ast }=J^{T}y$$
$J^{T}J$ 是半正定的厄米特矩阵，如果 $J$ 是列满秩的，那么 $J^{T}J$ 就是正定的，上述方程组的解存在唯一。

除了直接解法，求解最小二乘问题的算法还有 Cholesky 分解，QR 分解、SVD 分解、迭代逼近等等，这里就是不说了。

非线性最小二乘问题

求解非线性最小二乘问题有高斯-牛顿法、Levenberg-Marquardt 方法和截断 CG 方法，以及对应于大残差问题的方法的等等，时间关系，等到下一次再来补充。

非线性方程

这部分想表达的内容有牛顿方法扩展及其优缺点、非精确牛顿方法、Broyden 方法、评价函数等等，下次有空再来补充。

约束优化理论

拉格朗日条件和 KKT

$x$ 对应的积极集是指等式约束指标集和不等式约束指标集中取到等式的那些指标的并：
A ( x ) = E ∪ { i ∈ I ∣ c i ( x ) = 0 } \mathcal{A}(x)=\mathcal{E} \cup\left\{i \in \mathcal{I} | c_{i}(x)=0\right\} A(x)=E∪{i∈I∣ci​(x)=0}
LICQ（linear independence constraint qualication）条件满足是说，积极约束梯度集合
{ ∇ c i ( x ) , i ∈ A ( x ) } \left\{\nabla c_{i}(x), i \in \mathcal{A}(x)\right\} {∇ci​(x),i∈A(x)}
是线性独立的。

定义一般约束优化问题的拉格朗日函数为：
$$\mathcal{L}(x,\lambda )=f(x)-\sum _{i\in \mathcal{E}\cup \mathcal{I}}{\lambda }_i{c}_i(x)$$
假设 $x^{\ast }$ 是问题 $(1)$ 局部解，并且在 $x^{\ast }$ 点LICQ条件成立。那么存在一个拉格朗日乘子向量 ${\lambda }^{\ast },$，它的分量表示为 ${\lambda }_i^{\ast },i\in \mathcal{E}\cup \mathcal{I}$，使以下几个条件在 $\left(x^{\ast },{\lambda }^{\ast }\right)$ 点处成立，
$$\begin{array}{rlrl}{\nabla }_x\mathcal{L}\left(x^{\ast },{\lambda }^{\ast }\right)& =& 0,& \\ c_i\left(x^{\ast }\right)& =& 0,& \forall i\in \mathcal{E}\\ c_i\left(x^{\ast }\right)& \ge & 0,& \forall i\in \mathcal{I}\\ {\lambda }_i^{\ast }& \ge & 0,& \forall i\in \mathcal{I}\\ {\lambda }_i^{\ast }{c}_i\left(x^{\ast }\right)& =& 0,& \forall i\in \mathcal{E}\cup \mathcal{I}\end{array}$$
这几个条件就叫做著名的 Karush-Kuhn-Tucker 条件，简称 KKT 条件，这个方程系统，一般称为 KKT 系统。第一个条件是拉个朗日乘数法的要求，即极小化目标函数，极大化约束函数，极值点处对于各个变量的偏导为零。第二第三个条件分别为等式约束和不等式约束。以上三个条件都比较好记，都是已有的知识。后面两个条件，是对不等式约束而言的，叫做严格互补条件。一言以蔽之，就是对不等式约束指标对应的 ${\lambda }_i$ 必须大等于零，且 ${\lambda }_i$ 和 $c_i$ 之间必须有一个为零。互补性条件是对不等式约束而言的，没有不等式约束，就没有互补性条件。

对偶问题

考虑一个目标函数和约束函数是凹的凸优化问题，
$$\underset{x\in {\mathbb{R}}^n}{\min }f(x)\text{ s.t. }c(x)\ge 0$$
它的对偶问题是，
$$\underset{\lambda \in {\mathbb{R}}^m}{\max }q(\lambda )\text{ s.t. }\lambda \ge 0$$
这里 $q(\lambda )\equiv {\inf }_x\mathcal{L}(x,\lambda )$。对偶问题可以用来推导增广拉格朗日（augmented Lagrangian）方法。

其它

留一些关键词，等忙完手头的活儿，再来补充：局部解、孤立解、光滑性、单不等式约束、两不等式约束、线性化可行方向、切锥、稳定锥、wolfe 对偶…

非线性约束优化算法

积极集方法

求解非线性规划问题的一个挑战在于求解不等式约束。一个重要的逼近解的方法，就是积极集方法。

• 从一个初始的工作集 $\mathcal{W}$ 开始，作为最优积极集 ${\mathcal{A}}^{\ast }$，所谓的最优积极集就是最优值对应的等式约束的指标集。

• 然后求解这样一个问题：工作集中对应的约束取到等式约束。

• 检查是否存在一个拉格朗日乘子使得上述问题的解是原问题的解。如果不是，我们改变工作集，重复以上过程。

其他内容

求解非线性约束优化，除了积极集方法还有内点法，消元法等等，在这里就不讲了。当然，这一节本该还有一些其他内容，比如硬约束、软约束、滤子等等，这里也不说。

二次规划

序列二次规划

我们考虑一个等式约束问题，
$$\begin{array}{rl}& \min f(x)\\ & \text{ s.t. }c(x)=0\end{array}$$
这里 $f:{\mathbb{R}}^n\to \mathbb{R}$ 并且 $c={\left(c_1,\dots ,c_m\right)}^T$，$c_i:{\mathbb{R}}^n\to \mathbb{R}$ 是光滑函数。

我们知道这个问题的拉格朗日函数是：
$$\mathcal{L}(x,\lambda )=f(x)-{\lambda }^{T}c(x)$$
我们使用 $A(x)$ 去表示约束的雅克比矩阵，
$$A(x{)}^T=\left[\nabla c_1(x),\nabla c_2(x),\cdots ,\nabla c_m(x)\right]$$

由 KKT 条件，可得，
$$F(x,\lambda )=\left[\begin{array}{c}\nabla f(x)-A(x{)}^T\lambda \\ c(x)\end{array}\right]=0$$
接下来我们可以使用牛顿方法来求解这个非线性方程。

$F$ 关于 $x$ 和 $\lambda$ 的雅克比矩阵是，
$$F^{\prime }(x,\lambda )=\left[\begin{array}{cc}{\nabla }_{xx}^2\mathcal{L}(x,\lambda )& -A(x{)}^T\\ A(x)& 0\end{array}\right]$$
那么牛顿迭代就是：
$$\left[\begin{array}{c}{x}_{k+1}\\ {\lambda }_{k+1}\end{array}\right]=\left[\begin{array}{l}{x}_k\\ {\lambda }_k\end{array}\right]+\left[\begin{array}{l}{p}_k\\ p_{\lambda }\end{array}\right]$$
这里的 $p_k$ 和 $p_{\lambda }$ 是求解下面的牛顿-KKT 系统得到：
$$\left[\begin{array}{cc}{\nabla }_{xx}^2{\mathcal{L}}_k& -A_k^T\\ A_k& 0\end{array}\right]\left[\begin{array}{l}{p}_k\\ p_{\lambda }\end{array}\right]=\left[\begin{array}{c}-\nabla f_k+A_k^T{\lambda }_k\\ -c_k\end{array}\right]$$
这个方法一般就叫做牛顿-拉格朗日方法，对于 KKT 矩阵是非奇异的情况下还是很好用的。它具有二次收敛性。

我们可以换一个视角来看上面的牛顿迭代。假设我们在迭代过程中定义二次规划问题：
$$\begin{array}{ll}{\min }_p& f_k+\nabla f_k^{T}p+\frac{1}{2}{p}^T{\nabla }_{xx}^2{\mathcal{L}}_{k}p\\ \text{ s.t. }& A_{k}p+c_k=0\end{array}$$
在一些假设满足的条件下（LICQ），这个问题有一个唯一的解 $\left(p_k,l_k\right)$，满足，
$$\begin{array}{rl}{\nabla }_{xx}^2{\mathcal{L}}_k{p}_k+\nabla f_k-A_k^T{l}_k& =0\\ A_k{p}_k+c_k& =0\end{array}$$
仔细观察和对比，我们貌似可以发现，求解这个问题，好像就是在求解前面的牛顿-拉格朗日迭代的方向。将上面求牛顿方向系统的第一式的等式两边同时减去 $A_k^T{\lambda }_k$，就得到
$$\left[\begin{array}{cc}{\nabla }_{xx}^2{\mathcal{L}}_k& -A_k^T\\ A_k& 0\end{array}\right]\left[\begin{array}{c}{p}_k\\ {\lambda }_{k+1}\end{array}\right]=\left[\begin{array}{c}-\nabla f_k\\ -c_k\end{array}\right]$$
对比中，我们可以发现，在迭代中，我们只要求解这个问题，并令 $l_k={\lambda }_{k+1}=p_{\lambda }+{\lambda }_k$，也就解决了牛顿迭代求方向问题，上面的二次迭代问题也就迎刃而解了。这也就是所谓的序列二次规划（SQP）方法，贴一个算法框架，如图。

二次规划的内点法

可行域指的是满足不等式和等式约束的那些点的集合。内点法也叫障碍法，保证迭代的点始终在可行域当中，而跑出不等式约束定义的高高的"围墙障碍"，故叫它内点法。

考虑一个带不等式约束的凸二次规划问题：
$$\begin{array}{ll}{\min }_x& q(x)\equiv \frac{1}{2}{x}^{T}Gx+x^{T}c\\ \text{s.t.}& Ax\ge b\end{array}$$
这里的 $G$ 是对称半正定的，$A$ 是 $m\times n$ 的。这个问题的 KKT 系统可以用松弛的方式写出，即
$$\begin{array}{rl}Gx-A^T\lambda +c& =0\\ Ax-y-b& =0\\ y_i{\lambda }_i& =0,i=1,2,\cdots ,m\\ (y,\lambda )& \ge 0\end{array}$$
我们可以求解这个KKT系统来得到原问题的解。给定一个当前的迭代 $(x,y,\lambda )$，它满足 $(y,\lambda )>0$，我们定义一个互补性度量，
$$\mu =\frac{y^T\lambda }{m}$$
我们想使用原-对偶方法，考虑扰动的 KKT 系统：
$$F(x,y,\lambda ,\sigma ,\mu )=\left[\begin{array}{c}Gx-A^T\lambda +c\\ Ax-y-b\\ \mathcal{Y}\Lambda e-\sigma \mu e\end{array}\right]=0$$
这里
$$\mathcal{Y}=\mathrm{diag}\left(y_1,\cdots ,y_m\right),\Lambda =\mathrm{diag}\left({\lambda }_1,\cdots ,{\lambda }_m\right),e=(1,\cdots ,1{)}^T$$
并且 $\sigma \in [0,1]$。固定 $\mu$ 对上述系统使用牛顿方法，我们可以得到线性系统，
$$\left[\begin{array}{ccc}G& 0& -A^T\\ A& -I& 0\\ 0& \Lambda & \mathcal{Y}\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\\ \Delta \lambda \end{array}\right]=\left[\begin{array}{c}-\left(Gx-A^T\lambda +c\right)\\ -(Ax-y-b)\\ -\Lambda ye+\sigma \mu e\end{array}\right]$$
我们能得到下一步迭代，通过设置
$$\left(x^{+},y^{+},{\lambda }^{+}\right)=(x,y,\lambda )+\alpha (\Delta x,\Delta y,\Delta \lambda )$$
这里 $\alpha$ 的选择使得不等式约束 $\left(y^{+},{\lambda }^{+}\right)>0$ 被满足，且尽可能地满足更多的条件。

其它

等式约束二次规划：对称无限期分解、shur补方法、零空间方法、Krylov方法（GMRES、QMR、LSQR）、CG 方法。

不等式约束：积极集方法，梯度投影方法，内点法。

有界约束优化。

梯度投影方法。

罚方法

$l_1$和二次罚方法

惩罚函数法是求解有约束的最优化问题的一种算法。惩罚函数法的要旨是将一个有约束的最优化问题转化为一系列的无约束问题，这些无约束问题由原问题及罚函数，再加上惩罚因子组成。而且，这些无约束问题的解会收敛于所求问题的解。

一般的约束优化问题是：
min ⁡ x ∈ R n f ( x )  subject to  c i ( x ) = 0 , i ∈ E = { 1 , … , m e } c i ( x ) ≥ 0 , i ∈ I = { m e + 1 , … , m } \begin{array}{cl} {\min _{x \in \mathbb{R}^{n}}} & {f(x)} \\ {\text { subject to }} & {c_{i}(x)=0, \quad i \in \mathcal{E}=\left\{1, \ldots, m_{e}\right\}} \\ {} & {c_{i}(x) \geq 0, \quad i \in \mathcal{I}=\left\{m_{e}+1, \ldots, m\right\}} \end{array} minx∈Rn​ subject to ​f(x)ci​(x)=0,i∈E={1,…,me​}ci​(x)≥0,i∈I={me​+1,…,m}​
怎么构造罚函数？$l_1$ 罚函数定义为：
$${\varphi }_1(x;\mu )=f(x)+\mu \sum _{i\in \mathcal{E}}\left|c_i(x)\right|+\mu \sum _{i\in \mathcal{I}}{\left|c_i(x)\right|}^{-}$$
这里 [ z ] − = max ⁡ { 0 , − z } [z]^{-}=\max \{0,-z\} [z]−=max{0,−z}，$\mu$ 是罚因子。除了 $l_1$ 罚方法，还有二次罚方法。考虑等式约束问题，
$$\underset{x}{\min }f(x)\text{s.t.}{c}_i(x)=0,i\in \mathcal{E}$$
那么它的罚函数就可以写为，
$$Q(x;\mu )\equiv f(x)+\frac{\mu }{2}{c}_i^2(x)$$
一般的二次罚函数方法是，对于一般的约束问题，
$$\underset{x}{\min }f(x)\text{ s.t. }{c}_i(x)=0,i\in \mathcal{E};c_i(x)\ge 0,i\in \mathcal{I}$$
它的二次罚函数定义为，
$$Q(x;\mu )\equiv f(x)+\frac{\mu }{2}{c}_i^2(x)+\frac{\mu }{2}\sum _{i\in \mathcal{I}}{\left({\left[c_i(x)\right]}^{-}\right)}^2$$
这里$[y{]}^{-}:=\max (-y,0)$。

增广拉格朗日乘子法

增广拉格朗日法（Augmented Lagrange Method，简称ALM）是对二次惩罚法的一种改进。为了用无约束的目标函数替代原约束问题，二次惩罚法要求二次惩罚项的系数趋近于无穷（对约束的偏离给予很高的惩罚），但是这种要求会使得替代的目标函数的 Hessian 矩阵趋近于无穷（病态）。这使得替代函数的优化变得很困难，尤其是在最优点的附近。目标函数的行为很诡异，用二阶的泰勒公式逼近的话只有在非常小的区域内有效，因此收敛速度非常的慢。通过牛顿方程来计算下降方向会由于 Hessian 矩阵的病态性而变得很不准确。

为了解决这个问题，在二次惩罚的基础上引入了线性逼近的部分。个人感觉这种方法可以这样感性的理解：线性项和二次惩罚项实际是对约束偏离的一种惩罚，只不过他们有各自的针对性，二次项比较适合惩罚大的偏离，而小的偏离，只有当其系数很大时才起作用。相反，线性项比较适合小的偏离，而大的偏离二次项比它要好。那么这样的互补就能够使得在二次项系数比较小的情况下依然适用。而理论上恰恰可以验证这一点。细心的人就会发现，为什么不直接用拉格朗日乘子法求呢？这就是ALM算法巧妙的地方，因为多出一个二次惩罚项会使得算法的收敛速度很快。而且增广拉格朗日乘子法比拉格朗日乘子法普适性更好，需要的条件更加温和，比如不要求原函数是强凸的，甚至可以是非凸的。而且原函数在趋近于无穷的条件下，拉格朗日乘子法就无能为力了，可以用增广拉格朗日乘子法。究其原因，是因为二次惩罚项具有很好的矫正作用，在原函数非凸的情况下，只要满足一定的条件（二次惩罚项系数足够大），增广拉格朗日函数在最优点处的二阶导是正定的，因此具有严格的局部极小值。

增广拉格朗日方法其实是拉格朗日乘数法和二次罚方法的一个结合，比如，对于等式约束问题，增广拉格朗日函数写为：
$${\mathcal{L}}_A(x,\lambda ;\mu )\equiv f(x)-\sum _{i\in \mathcal{E}}{\lambda }_i{c}_i(x)+\frac{\mu }{2}\sum _{i\in \mathcal{E}}{c}_i^2(x)$$

交替方向乘子法（ADMM）

网上的一些资料根本就没有把ADMM的来龙去脉说清楚，发现只是一个地方简单写了一下流程，别的地方就各种抄，共轭函数，对偶梯度上升什么的，都没讲清楚，给跪了。下面我来讲讲在机器学习中用得很多的ADMM方法到底是何方神圣。

共轭函数

给定函数 $f:{\mathbb{R}}^n\to \mathbb{R}$，那么函数
$$f^{\ast }(y)=\underset{x}{\max }(y^{T}x-f(x))$$
就叫做它的共轭函数。其实一个更直观的理解是：对一个固定的 $y$，将 $y^{T}x$ 看成是一条斜率为 $y$ 的直线，它和 $f(x)$ 关于 $x$ 的距离的最大值，就是 $f^{\ast }(y)$。百度百科有关于这个直观说法的一个解释，看看就明白了。

关于共轭函数有几点重要的说明：

• 不管 $f$ 凸不凸，它的共轭函数总是一个凸函数。

• 如果 $f$ 是闭凸的（闭指定义域是闭的），那么 $f^{\ast \ast }=f$。

• 如果 $f$ 是严格凸的，那么
  $$\nabla f^{\ast }(y)=\underset{z}{\mathrm{argmin}}(f(z)-y^{T}z)$$

• 共轭总是频频出现在对偶规划中，因为极小问题总是容易凑出一个共轭：$-f^{\ast }(y)={\min }_x(f(x)-y^{T}x)$。

关于 $f^{\ast }$ 的凸性，这篇博客给了一个比较直观的图示说明，下面我从数学上，不太严格地做个简单证明。以一维的情况说明。

假设 $f$ 是一个凸函数，下面都不妨考虑函数的最值都不再边界处取到。${\max }_x(yx-f(x))$ 的极值点在 $f_x^{\prime }(x)=y$ 处取到，定义 $g:=(f_x^{\prime }{)}^{-1}$，那么 $x=g(y)$ 可能会是一堆点。则有
$$f^{\ast }(y)=\underset{x}{\max }(yx-f(x))=yg(x)-f(g(y))$$
进而
$$(f^{\ast }(y){)}_y^{\prime }=g(y)+y{g}_y^{\prime }(y)-f_x^{\prime }(g(y))g_y^{\prime }(y)=g(y)$$
那么
$$(f^{\ast }(y){)}_y^{\prime \prime }y=g_y^{\prime }(y)=((f_x^{\prime }{)}^{-1}{)}_y^{\prime }\ge 0$$
故而，$f^{\ast }$ 是凸的。

关于 $f^{\ast \ast }=f$，也是容易证明的。我们假设 $f$ 是闭凸的。${\max }_y(zy-f^{\ast }(y))$ 的值在 $g(y)=z$ 处取到，那么
$$f^{\ast \ast }(z)=\underset{y}{\max }(zy-f^{\ast }(y))=f(g(g^{-1}(z)))=f(z)$$
$z$ 换成 $x$ 就是 $f(x)$。

第三条非常重要，它说明了共轭函数的梯度，其实就是共轭函数取到极大值对应的 $x$ 值，它从 $(f\ast (y){)}_y^{\prime }=g(y)$ 就可以看出来。

对偶梯度上升法

有了上知识的铺垫，我们就可以说清楚对偶上升方法了。以考虑等式约束问题为例（一般约束问题也是类似的流程），假设$f(x)$ 是严格凸的，我们考虑问题：
$$\underset{x}{\min }f(x)\text{ subject to }Ax=b$$
它的拉格朗日对偶问题是：
$$\underset{u}{\max }\underset{x}{\min }(f(x)+u^T(Ax-b))$$
有理论表明，若原问题和对偶问题满足强对偶条件，即对偶函数关于 $u$ 的最大值等价于原优化问题关于 $x$ 的最小。那么原问题和对偶问题对于 $x$ 是同解的。也就是说只要找到使得对偶问题对应最大的 $u$，其对应的 $x$ 就是原优化问题的解，那么我们就解决了原始优化问题。

所以，下面我们来求解这个对偶问题。先把和 $x$ 无关的变量提出 ${\min }_x$，再想办法凑出 $f^{\ast }$，因为我们要用到对偶的性质。
KaTeX parse error: No such environment: split at position 7: \begin{̲s̲p̲l̲i̲t̲}̲ \max _u \min…

那么对偶问题就成了
$$\underset{u}{\max }-f^{\ast }\left(-A^{T}u\right)-b^{T}u$$
这里 $f^{\ast }$ 是 $f$ 的共轭，这里 ${\max }_u$ 后面不加括号，表示它管着下面的所有，下同，不再重述。定义 $g(u)=-f^{\ast }\left(-A^{T}u\right)-b^{T}u$，我们希望能极大化 $g(u)$，一个简单的想法是沿着 $g(u)$ 梯度上升的方向去走。注意到，
$$\partial g(u)=A\partial f^{\ast }\left(-A^{T}u\right)-b$$
因此，利用共轭的性质，
$$\partial g(u)=Ax-b\text{ where }x\in \underset{z}{\mathrm{argmin}}f(z)+u^{T}Az$$
因为 $f$ 是严格凸的，$f^{\ast }$ 是可微的，那么，就有了所谓的对偶梯度上升方法。从一个对偶初值 $u^{(0)}$ 开始，重复以下过程：
$$\begin{array}{rl}& x^{(k)}=\underset{x}{\mathrm{argmin}}f(x)+{\left(u^{(k-1)}\right)}^{T}Ax\\ & u^{(k)}=u^{(k-1)}+t_k\left(A{x}^{(k)}-b\right)\end{array}$$
这里的步长 $t_k$ 使用标准的方式选取的。近端梯度和加速可以应用到这个过程中进行优化。

交替方向乘子法

交替方向乘子法（ADMM）是一种求解具有可分离的凸优化问题的重要方法，由于处理速度快，收敛性能好，ADMM算法在统计学习、机器学习等领域有着广泛应用。ADMM算法一般用于解决如下的凸优化问题:
$$\underset{x,y}{\min }f(x)+g(y)\text{ subject to }Ax+By=c$$
其中的 $f$ 和 $g$ 都是凸函数。

它的增广拉格朗日函数如下：
$$L_p(x,y,\lambda )=f(x)+g(y)+{\lambda }^T(Ax+By-c)+(\rho /2)\Vert Ax+By-c{\Vert }_2^2,\rho >0$$

ADMM算法求解思想和推导同梯度上升法，最后重复迭代以下过程：
$$\begin{array}{rl}{x}^{k+1}& :=\arg \underset{x}{\min }{L}_p(x,y,\lambda )\\ x^{k+1}& :=\arg \underset{y}{\min }{L}_p(x,y,\lambda )\\ {\lambda }^{k+1}& :={\lambda }^k+\rho \left(A{x}^{k+1}+B{y}^{k+1}-c\right)\end{array}$$
上述迭代可以进行简化。

• 第一步简化，通过公式 $\Vert a+b{\Vert }_2^2=\Vert a{\Vert }_2^2+\Vert b{\Vert }_2^2+2a^{T}b$，替换掉拉格朗日函数中的线性项 ${\lambda }^T(Ax+By-c)$ 和二次项 $\rho /2\Vert Ax+By-c{\Vert }_2^2$，可以得到
  $${\lambda }^T(Ax+By-c)+\rho /2\Vert Ax+By-c{\Vert }_2^2=\rho /2\Vert Ax+By-c+\lambda /\rho {\Vert }_2^2-\rho /2\Vert \lambda /\rho {\Vert }_2^2$$
  那么ADMM的过程可以化简如下：
  $$\begin{array}{rl}{x}^{k+1}& :=\arg \underset{x}{\min }(f(x)+\rho /2{\Vert Ax+B{y}^k-c+{\lambda }^k/\rho \Vert }_2^2\\ y^{k+1}& :=\arg \underset{y}{\min }(g(y)+\rho /2{\Vert A{x}^{k+1}+By-c+{\lambda }^k/\rho \Vert }_2^2\\ {\lambda }^{k+1}& :={\lambda }^k+\rho \left(A{x}^{k+1}+B{y}^{k+1}-c\right)\end{array}$$

• 第二步化简，零缩放对偶变量$u=\lambda /\rho$，于是ADMM过程可化简为：
  $$\begin{array}{rl}{x}^{k+1}& :=\arg \min (f(x)+\rho /2{\Vert Ax+B{y}^k-c+u^k\Vert }_2^2\\ y^{k+1}& :=\arg \min (g(y)+\rho /2{\Vert A{x}^{k+1}+By-c+u^k\Vert }_2^2\\ u^{k+1}& :=u^k+\left(A{x}^{k+1}+B{y}^{k+1}-c\right)\end{array}$$

ADMM相当于把一个大的问题分成了两个子问题，缩小了问题的规模，分而治之。

[help me with MathJax]