无约束问题的最小化
最近在看"Machine Learning A Probability Perspective"的逻辑回归，因为算法比较多，一起看书的伙伴们就决定把"Convex Optimization"讲优化算法的内容看一下。
在阅读以下内容之前，要记住问题的出发点，即给定无约束的凸问题该用什么方法求得最优解（这里假设问题是求最小值）？各种方法之间有什么性质。
也许最先想到的方法就是对变量求梯度，让它等于(0)的变量值就是最优解。没错，这样的方法很简单，但是当求导过程过于复杂的时候该怎么办呢？也许你还知道梯度下降，牛顿法等等。如果你对这些方法没又听说过，或者还想知道更多，那接着阅读下去吧，当然最好还是建议你看一下原书。
强凸及其性质
函数强凸是说在凸集( m{S})中，对于(forall mathbf{x}in m S)均有下式成立：
其中(m>0)，将上式带入到函数的泰勒展开式
得到
当固定(mathbf{x})时，等式右边为(mathbf{y})的二次项。为了求与(mathbf{y})无关仅与(mathbf{x})有关的(f(mathbf{y}))下界，在等式右边对(mathbf y)求梯度，得到：
带入得到
因为(eqref{eq:2})对所有的(mathbf{y})均成立，令(mathbf y=mathbf x^*)，可以得到
因此可以得知当(mathbf{x})的梯度范数越小时，其与靠近最优点(mathbf{x^*})，(eqref{eq:2})式子同样可以得到( abla f(mathbf x))的下界
可以看到，当某一点的梯度较小时，作为其下界(2m(f(mathbf x)-p^*))也应该很小，因此(mathbf{x})越接近最优点，而这与最优点梯度为0的常识相符合。因此也可以用(Vert abla f(mathbf x)Vert_2^2leqepsilon)当作收敛条件。
将(mathbf x^*)带入到(eqref{eq:1})中，再利用柯西施瓦茨不等式
带入可得
因为(f(mathbf{x}) > f(mathbf{x^*}))，故(-Vert abla f(mathbf x) Vert_2Vertmathbf x^*-mathbf x Vert_2+frac{m}{2}Vertmathbf x^* - mathbf xVert_2^2 leq0)，整理得到
即通过(Vert abla f(mathbf x)Vert_2)可以得到(mathbf x)和(mathbf x^*)的差距上界。
一般而言( m S)被当做是下子集(sublevel sets)，是有界的，故存在(M>0)使得(forall mathbf xin m S)满足
使用相同的方法，可以得到：
定义(kappa = frac{M}{m})为( m S)的condition number，condition number会在之后的收敛性分析中扮演重要的角色。( abla^2 f(mathbf x))是实对称矩阵，可以看出( extit condition\, number)可以看作( abla^2f(mathbf x))在(forall mathbf x in m S)中最大特征值和最小特征值的比值。
根据泰勒展开(alpha geq f(mathbf{y})approx f(mathbf{x^*})+frac{1}{2}(mathbf{y}-mathbf{x^*})^T abla^2f(mathbf{z})(mathbf{y}-mathbf{x^*}))，可以看出(mathbf{y})大约在一个以(mathbf{x^*})为中心的椭圆内。我们可以以( abla^2 f(mathbf{z}))的特征值描述下子集(S)的形状。下子集最长半轴(r_{max} = 1/sqrt{lambda_{min}( abla^2 f(mathbf{z}))})，最短半轴(r_{min} = 1/sqrt{lambda_{max}( abla^2 f(mathbf{z}))})。根据cond的定义，可见其为(r_{max}^2/r^2_{min})。于是可以集合解释为，设凸集合(C)的沿某一方向上的宽度为
集合(C)在不同方向上的最大和最小宽度如下：
而(condition\, number)为(W_{max}^2/W_{min}^2)，从而可以看出(condition\, number)越小说明集合(C)越圆，反之则其各项异性比较大。
对于(alpha)-sublevel子集(C_{alpha}={mathbf yvert f(mathbf y)leqalpha})，根据(eqref{eq:2})和(eqref{eq:3})式可得
再利用(alpha geq f(mathbf y))，则(alpha geq p^*+frac{m}{2}Vertmathbf x^*-mathbf y Vert_2^2)。当(mathbf y)趋近于(mathbf x^*)时，(p^*+frac{M}{2}Vertmathbf x^*-mathbf y Vert_2^2)趋近于(p^*)，此时(alpha geq p^*+frac{M}{2}Vertmathbf x^*-mathbf y Vert_2^2)可以求得部分(mathbf y)的范围。上述两个不等式可以得到：
可见(alpha)下子集在(B_{inner})和(B_{outer})之间，则根据几何解释( extbf{cond}(C_{alpha})leq frac{M}{m})。这样，(condition \, number)和下子集的形状联系了起来。
当(alpha)趋近于(mathbf x^*)时，根据泰勒展开，其梯度可以认为是0，得到下子集(C_alpha)
可以看到此时的下子集(C_alpha)接近为( abla^2 f(mathbf x^*))的椭球，此时的( extbf{cond}(C_{alpha})=kappa{ abla^2 f(mathbf x^*)})。
下降方法（Descent methods）求解框架
下降方法是一种迭代算法，给定上一部的(mathbf{x}^{k-1})其通过寻找新的(mathbf{x}^{k})来使得(f(mathbf{x}^{k})<f(mathbf{x}^{k-1}))。(mathbf{x^t})与(mathbf {x}^{k-1})的一般关系为(mathbf{x}^k=mathbf{x}^{k-1}+t^{k-1}Deltamathbf{x}^{k-1})。可见下降方法的主要内容为寻找(Deltamathbf{x}^{k-1})和(t^{k-1})，以保证(f(mathbf x^{k})<f(mathbf{x}^{k-1}))。
对于凸集合( m S)和凸函数(f)，当(forall mathbf{y} in m S)，( abla f(mathbf x^{k-1})(mathbf y-mathbf x^{k-1})geq 0)时，则(f(mathbf y)geq f(mathbf x^k))，故要寻找的(Delta mathbf x^{k-1})必然与( abla f(mathbf x^{k-1}))成负角度。
Descent methods的一般流程如下：
- 给定初始点(mathbf{x})
- 重复下列过程
- 决定搜索方向(Delta mathbf{x})
- 决定步长(t) (Line search)
- 对(mathbf{x})进行更新
- 检查是否收敛，如果收敛则停止迭代
有两个Line Search方法，一个是exact line search另外一个是不精确的backtracting line search。
-
exact line search 即是说，在给定搜索方向(Delta mathbf x)后，选择使得(f(mathbf x^{k-1})+tDelta mathbf{x}^{t-1})最小的步长(t)，在上式求解较为简单时该方法可以应用。
-
backtracking line search 相对于exact line search来说应用更为普遍。其算法流程如下
- 选择(0<alpha<0.5, 0<eta<1)，以及搜索方向(Delta mathbf{x})（对于(alpha)为什么要小于0.5请看下一节）
- 当(f(mathbf{x}+tDelta mathbf{x})>f(mathbf x)+alpha t abla f(mathbf x)^T Delta mathbf x)时，令(t:=eta t)
当(t)足够小的时候，(f(mathbf x+tDelta mathbf x)approx f(mathbf x)+t abla f(mathbf x)^T Deltamathbf x <f(mathbf x)+alpha t abla f(mathbf x)^T Deltamathbf x)，因此上述的过程肯定是收敛的。
上图是书中关于(alpha)和(t)的解释：在确定( abla f(mathbf x)^T Delta mathbf x)后，(f(mathbf x)+t abla f(mathbf x)^T Deltamathbf x)即为关于(t)的函数，(alpha abla f(mathbf x)^T Delta mathbf x)即为斜率，(alpha)越小，则直线越平缓。在确定(alpha)和(eta)后，从图中可以看出在(t_0)点，(f(mathbf x)+t abla f(mathbf x)^T Deltamathbf x)与(f(mathbf x+tDelta mathbf x))相交，因此可能的(t)取值为1或者((eta \, t_0, \, t_0])，从而(t geq min{1,\, eta\, t_0})。
梯度下降方法
对于下降方法中的搜索方向(Delta mathbf x)，一个自然的选择为负梯度方向，即(- abla f(mathbf x))（将(f(mathbf x+tDelta mathbf x))一次展开得到(f(mathbf x)+t abla f(mathbf x)^TDelta mathbf x)在(Delta mathbf x)取(- abla f(mathbf x))时可以保证(f(mathbf x+tDelta mathbf x)leq f(mathbf x))）。下面对梯度下降法的收敛速度进行分析。
-
exact line search搜索的收敛速度
将(Delta mathbf x=- abla f(mathbf x))带入到(eqref{eq:3})式子中，得到
[f(mathbf x-t abla f(mathbf x))leq f(mathbf x)-tVert abla f(mathbf x)Vert_2^2+frac{M}{2}t^2Vert abla f(mathbf x)Vert_2^2 ]通过对(t)求最小值，得到
[f(mathbf x-t abla f(mathbf x))leq f(mathbf x)-frac{1}{2M}Vert abla f(mathbf x)Vert_2^2 ]通过(eqref{eq:2})式可以得到，(Vert abla f(mathbf x)Vert_2^2geq2m(f(mathbf x) - p^*))，在上式的两边同时减去(p^*)得到
[f(mathbf x-t abla f(mathbf x)) - p^* leq (1-frac m M)left(f(mathbf x)-p^* ight) ]通过上式可以看到其为线性收敛，当(frac m M)越接近1时，收敛速度越快。而通过(cond)的几何解释，(frac m M)与(alpha)下子集的形状有关，如果其各方向上差异性较小则越接近1，这也是为什么在优化前先对数据进行归一化处理的原因。
-
backtracting line search
首先对上一节中(alpha)的取值进行说明。在(0leq t leq frac 1 M)时，下式成立
[-t + frac M 2 t^2 leq -frac 1 2 t ]将(Delta mathbf x=- abla f(mathbf x))带入得到
[egin{align*} fleft(mathbf x-t abla f(mathbf x) ight)& leq f(mathbf x)-tVert abla f(mathbf x)Vert_2^2+frac M 2 Vert abla f(mathbf x)Vert_2^2 \ & leq f(mathbf x) - frac t 2 Vert abla f(mathbf x)Vert_2^2 \ & leq f(mathbf x) - alpha t Vert abla f(mathbf x)Vert_2^2 end{align*} ]因此搜索步长(t)取值为1或(tgeq frac eta M)。采取跟exact line search一样的分析步骤，在上式两边同时减去(p^*)得到
[egin{align*} fleft(mathbf x-t abla f(mathbf x) ight)-p^* & leq f(mathbf x)-p^*-min{alpha, \, frac {alphaeta}{M}} abla f(mathbf x)Vert_2^2 \ &leq (1-min{2malpha, \, frac {2malphaeta}{M}})(f(mathbf x)-p^*) end{align*} ] -
例子
下面结合书上的一个例子应用exact line search来进行阐述，设一个二次的目标函数为
[f(mathbf x)=frac 1 2left(x_1^2+gamma x_2^2 ight) ]式中(gamma>0)，定义域为(mathbf{R}^2)，目标函数(f)的Hessian矩阵为常量，其Hessian矩阵的特征根为(1)和(gamma)，因此condition number为(frac {max{1, \, gamma}}{min{1\, gamma}}=max{gamma, frac 1 gamma})
，选定初始点(mathbf x_0=(gamma, 1))，其负梯度方向为((gamma \, gamma))，通过计算可以得到下降幅度(t=frac 2 {(1+gamma)})，得到(mathbf x_1=(gammafrac {gamma-1}{gamma+1}, -frac {gamma-1}{gamma+1})，通过归纳可以得到(mathbf x_k=(gammaleft(frac{gamma-1}{gamma+1} ight)^k, left(-frac {gamma-1}{gamma+1} ight)^k)，从而
[f(mathbf x_k)=frac {gammaleft(gamma+1 ight)} 2left(frac{gamma-1}{gamma+1} ight)^{2k}=left(frac{gamma-1}{gamma+1} ight)^{2k}f(mathbf x_0) ]
Steepest descent method
不同于[gradient descent](# 梯度下降方法)中使用负梯度作为下降方向，Steepest descent中使用(Delta mathbf x=underset{v}{min}{- abla f(mathbf x)^T v })作为下降方向。定义normolized steepest descent direction (Delta mathbf x_{nsd} = argminleft{ abla f(mathbf x)^Tv vert leftVert v ightVert=1 ight})，通过对(Delta mathbf x_{nsd})进行变换并记住对偶范数的定义(left Vert z ight Vert_* = underset{v}{sup}left {z^Tv vert leftVert v ightVertleq 1 ight })得到(Delta mathbf x_{sd}=Delta mathbf x_{nsd}left Vert abla f(mathbf x) ight Vert_*)，从而( abla f(mathbf x)^TDelta mathbf x_{sd}= abla f(mathbf x)^TDelta mathbf x_{nsd}left Vert abla f(mathbf x) ightVert_*=-left Vert abla f(mathbf x) ight Vert_*^2)
当使用二范数时，Steepest descent的方向为负梯度方向，与gradient descent的方向一致。
-
Steepest descent for quadratic norm
quadratic norm定义为(left Vert z ight Vert_P=z^TPz)，其中(P in mathbf S_{++}^n)，通过解下述不等式
[egin{align*} & min_{v} abla f(mathbf x)^Tv\ s.t qquad &v^TPvleq 1 end{align*} ]得到
[egin{align*} Delta mathbf x_{nsd}&=-left( abla f(mathbf x)^T P^{-1} abla f(mathbf x) ight)^{-frac 1 2} P^{-1} abla f(mathbf x)\ Delta mathbf x_{sd}&=-P^{-1} abla f(mathbf x) end{align*} ]根据(Delta mathbf x_{nsd})的定义，可以看出其几何意义如下，图中阴影的椭圆表示(v^TPvleq1)
除了上述的集合解释外，还可以通过坐标变换得到(Delta mathbf x_{nsd})。
令(ar{mathbf x}=P^{frac 1 2}mathbf x)，(ar{f}(ar {mathbf x})=f(P^{-frac 1 2}mathbf {ar{x}}))，
通过这样的变换就将quadratic norm 变换成了二范数了。(z^TPz=leftVert ar{z} ightVert_2^2)。此时(Delta mathbf{ar{x}}_{sd}=- abla f(mathbf{ar{x}})=-P^{-frac 1 2} abla f(mathbf x))，(Delta mathbf x_{sd}=P^{-frac 1 2}Delta mathbf{ar{x}}_{sd}=-P^{-1} abla f(mathbf x))。
根据上述变换(ar{f}(mathbf x)=f(P^{-frac 1 2}mathbf {ar{x}}))，对(ar{mathbf {x}})进行求Hessian矩阵，得到(P^{-frac{1}{2}} abla^2 f(mathbf x)P^{-frac{1}{2}})，当(P= abla^2f(mathbf x))时，变换后的函数有更好的condition number。也就是说，如果({mathbf xvert mathbf x^{T}Pmathbf xleq1})如果与sublevel set形状相似的话，descent method在变换后的变量下会取得很好的效果。
-
Steepest descent for (l_1) norm
在(l_1) norm下，(Delta mathbf x_{nsd}=minleft{ abla fleft(mathbf x ight)^Tvvert left Vert v ight Vert_1leq1 ight})，可以得到(Delta mathbf x_{nsd}=left(0,dots,-signleft(frac {partial f(mathbf x)}{partial x_i} ight)e_i,dots,0 ight))，其中(i)为其绝对值最大的分量。这样，每次下降的方向只涉及一个分量，该算法有时被称为coordinate descent algorithm。
可以根据下图对Steepest descent在(l_1) norm的情况下进行理解
-
Converge Analysis
由于任何范数都可以被二范数定界，故假设(gamma,ar{gamma}inleft(0,1 ight])，使得
[egin{align*} Vert mathbf xVert &geq gamma Vert mathbf xVert_2\ Vert mathbf xVert_* &geq gamma Vert mathbf xVert_2\ end{align*} ]同样根据强凸，得到
[egin{align*} fleft(mathbf x+tDeltamathbf x_{sd} ight)&leq fleft(mathbf x ight)+t abla f(mathbf x)^TDeltamathbf x_{sd}+frac{M}{2}t^2VertDelta mathbf x_{sd}Vert_2^2\ &leq f(mathbf x)-tleft Vert abla fleft(mathbf x ight) ight Vert_*^2+frac{M}{2gamma^2}t^2left Vert abla fleft(mathbf x ight) ight Vert_*^2 end{align*} ]在(t=frac {gamma^2} M)时右式取得最小值，带入得到
[f(mathbf x + tDelta mathbf x)leq f(mathbf x)-frac{gamma^2}{2M}Vert abla fleft(mathbf x ight)Vert_*^2 ]发现(t=frac{gamma^2}{M})时满足backtracking line search的条件。因此采用backtracking line search的时候，步长(t=minleft{1, frac{eta gamma^2}{M} ight })，因此
[egin{align*} f(mathbf x+tDelta mathbf x_{sd})&leq f(mathbf x)-alpha ar{gamma}^2minleft{1, frac{eta gamma^2}{M} ight}leftVert abla f(mathbf x) ightVert_2^2\ & leq f(mathbf x)-2malpha ar{gamma}^2minleft{1, frac{eta gamma^2}{M} ight}left(f(mathbf x-p^*) ight) end{align*} ]由此可见收敛性。
Newton's method
-
Newton's step
定义(Delta mathbf x_{nt}=- abla^2fleft(mathbf x ight)^{-1} abla fleft(mathbf x ight))为newton's step，正定的Hessian矩阵意味着( abla fleft(mathbf x ight)Delta mathbf x_{st}leq0)，因此(f(mathbf x + tDelta mathbf x_{sd})leq f(mathbf x))。
对函数(f)的二阶泰勒展开求最小值，可以发现(v=Delta mathbf x_{nt})时取得。
[min_{v} f(mathbf x)+ abla f(mathbf x)^Tv +frac{1}{2}v^T abla ^2 f(mathbf x) v ]
可以理解为在(mathbf x点对函数)(f)进行二次逼近时，逼近函数在(mathbf x+Delta mathbf x_{st})点取到最小值。
根据[Steepest descent](# Steepest descent method)，如果quadratic norm中的(P= abla^2 f(mathbf x))，可以得到(Delta mathbf x_{nt}=Delta mathbf x_{nsd})
对函数(f)得到函数进行一阶泰勒展开
[ abla f(mathbf x+v)approx abla f(mathbf x)+ abla ^2 f(mathbf x)v ]根据极值条件，在最优点其导函数值为(0)，此时(v=Delta mathbf x_{st})。
牛顿方法可以理解为对一阶近似的导函数求极值，得到极值点，并在该点再次进行上述过程。
Newton's step具有仿射不变性，对原变量(mathbf x)左乘非奇异阵(T)进行变换，令(mathbf y=T mathbf x,quad ar{f}(mathbf x)=f(Tmathbf x))，在(ar{f}(mathbf x))中Newton's step为
[egin{align*} Deltaar{mathbf x}&=-T^{-1} abla^2 f(mathbf y)^{-1} abla f(mathbf y)\ &=T^{-1}Delta mathbf x end{align*} ]经过变换的函数(ar{f})在(mathbf x)的Newton's step为(f)在(y)点的newton's step再左乘(T^{-1})。
定义(lambda(mathbf x)=left( abla fleft(mathbf x ight)^T abla ^2 fleft(mathbf x ight)^{-1} abla fleft(mathbf x ight) ight)^{frac 1 2})为newton decrement。
可以看到(lambda(mathbf x)^2)为newton step在( abla^2 f(mathbf x))的quadratic norm，再结合Steepest descent中的(Delta mathbf x_{nsd})，(lambda(mathbf x))与其(Delta mathbf x_{nsd})中的单位化部分相一致。
-
Newton's method
Newton's method 有damped和pure之分，其区别为步长的不同，pure newton method的步长(t)固定为(1)，而damped newton method的步长$t为利用Line search确定，其算法流程如下：
- 确定初始点(mathbf x^{(0)})和收敛条件(epsilon)
- 计算(lambda(mathbf x))并与(epsilon)进行比较，如果未满足收敛条件则仅需进行
- 计算(Delta mathbf x_{st})，并根据backtracking line search or exact line search 计算步长(t)
- 对(mathbf x)进行更新，(mathbf x^{(new)}=mathbf x^{(old)}+tDelta mathbf x_{st})
Newton's method的收敛性分析比较繁琐，这里不再进行分析。后续会补上quasi-Newton methods。