DennyQi's Log

02 无约束凸优化问题 Unconstrained Convex Optimization

在这一部分我们的目标是求出凸函数的最小值。一般来说,只要我们能解出方程f(x)=0\nabla f(x)=0我们就能求出最小值点。然而很多时候这一方程的封闭解是不存在的,这要求我们用其它的手段来求最小值。

梯度下降(Gradient Descent)

在线性规划的单纯形法中我们注意到每次移动到一个更优值最终能保证我们找到最优解,那在非线性的凸优化中是否存在一个类似的方法?这就是下降方法,我们每次都找到一个更小的值,期待我们最终找到最优值。

f(x0)f(x_0)处,f(x0)-\nabla f(x_0)的方向是函数值下降最快的方向,因此我们能保证每次往负梯度方向移动一小步,函数值是一定会下降的,因此我们期待反复迭代这一过程最终得到最小值。然而步长的确定是一个困难的工作,步长太长会导致反复横跳,步长太短会导致无法收敛到最小值。我们无法对于任意函数给出一个选取最优步长的公式,但我们可以讨论当函数满足一些特殊的性质时,我们能给出步长选取的方法。

当函数满足对于任意x,yx,yf(x)f(y)Lxy\|\nabla f(x)-\nabla f(y)\| \leq L\|x-y\|时,我们称ffLL-smooth函数(这个条件还等价于2f(x)\nabla^2 f(x)的最大特征值的绝对值不超过LL,也等价于f(y)f(x)f(x),yxL2yx2|f(y)-f(x)-\lang \nabla f(x),y-x\rang| \leq \dfrac{L}{2}\|y-x\|^2恒成立)。这个条件保证了函数的一阶导变化不会太快。此时取步长η1L\eta \leq \dfrac{1}{L},我们能够证明每次迭代xk+1=xkηf(xk)x_{k+1}=x_k-\eta\nabla f(x_k),则有f(xk+1)f(xk)t2f(xk)2f(x_{k+1})\leq f(x_k)-\dfrac{t}{2}\|\nabla f(x_k)\|^2成立。这称为下降引理,选取这样的步长我们能保证函数值不断下降,收敛到最小值。

下降引理只给出了单步的下降估计,无法给出整体的下降速率的估计。为了分析梯度下降整体的算法,我们先分析一个连续版本的梯度下降算法:在每点处取梯度方向连成一条曲线,也即曲线x(t)x(t)满足ddtx(t)=f(x(t))\dfrac{d}{dt}x(t)=-\nabla f(x(t)),称为梯度流。对于x(t)x(t)我们可以计算得到f(x(t))f(x)x(0)x222tf(x(t))-f(x^*)\leq \dfrac{\|x(0)-x^*\|_2^2}{2t},其中xx^*是极值点。梯度下降可以看作梯度流下降的离散版本的近似,分析得到f(xk)f(x)x0x222kηf(x_k)-f(x^*)\leq \dfrac{\|x_0-x^*\|^2_2}{2k\eta},也即函数值的收敛可以用迭代次数kk的反比来bound。

一个与迭代次数成反比的估计不是任何时候都够用的bound。例如我们会发现对于二次函数梯度下降的收敛速率达到指数收敛。也就是说,对于满足更苛刻要求的函数,我们可以做出更好的收敛分析。一个重要的性质称为Strong Convexity(强凸性),如果f(x)μ2x2f(x)-\dfrac{\mu}{2}\|x\|^2是凸的,那么称f(x)f(x)是强凸的。直观上可以看到,强凸函数是比二次函数“更凸”的函数。强凸有一些其它等价描述,例如μ\mu-强凸等价于二阶导矩阵的特征值全部μ\geq \mu,也等价于f(y)f(x)+f(y)\geq f(x)+ f(x),yx\lang \nabla f(x),y-x \rang +μ2yx2+\dfrac{\mu}{2}\|y-x\|^2恒成立,根据凸函数的梯度单调性得到它也等价于f(x)f(y),xyμxy2\lang \nabla f(x)-\nabla f(y),x-y\rang\geq \mu\|x-y\|^2恒成立。从中也可以看出强凸一定是严格凸的。对于μ\mu-强凸的LL-smooth函数,梯度下降满足xkx2(1μη)kx0x2\|x_k-x^*\|^2\leq (1-\mu\eta)^k\|x_0-x^*\|^2。可见自变量的距离可以用迭代次数的指数次方来bound。函数值的收敛也是指数的,有f(xk)f(x)f(x_k)-f(x^*) \leq L2(1μη)kx0x2\dfrac{L}{2}(1-\mu\eta)^k\|x_0-x^*\|^2

对于二次函数f(x)=12xQxf(x)=\dfrac{1}{2}x^\top Qx,其中QQ是半正定的,根据强凸和LL-smooth的等价条件,我们立马得到f(x)f(x)是一个λmin(Q)\lambda_{\min}(Q)-强凸以及λmax(Q)\lambda_{\max}(Q)-smooth的函数。于是可以取η=1λmax(Q)\eta = \dfrac{1}{\lambda_{\max}(Q)},根据上面的分析得到梯度下降的收敛速率正比于(1λmin(Q)λmax(Q))k\left(1-\dfrac{\lambda_{\min}(Q)}{\lambda_{\max}(Q)}\right)^k。事实上,对于二次函数,经过更精确的分析我们发现可以取η=2λmin(Q)+λmax(Q)\eta = \dfrac{2}{\lambda_{\min}(Q)+\lambda_{\max}(Q)},这样得到收敛速率正比于(λmax(Q)λmin(Q)λmax(Q)+λmin(Q))2k\left(\dfrac{\lambda_{\max}(Q)-\lambda_{\min}(Q)}{\lambda_{\max}(Q)+\lambda_{\min}(Q)}\right)^{2k},把λmax(Q)λmin(Q)\dfrac{\lambda_{\max}(Q)}{\lambda_{\min}(Q)}记为κ\kappa,称为Condition Number,可见二次函数的收敛速率是关于(κ1κ+1)2\left(\dfrac{\kappa-1}{\kappa+1}\right)^2的。对于一般的二阶可导函数f(x)f(x),我们可以在局部泰勒展开到二阶来逼近,这样就可以把它视为二次函数,对它用梯度下降时收敛速率取决于(κ(2f(x))1κ(2f(x))+1)2\left(\dfrac{\kappa(\nabla^2 f(x))-1}{\kappa(\nabla^2 f(x))+1}\right)^2

线搜索(Line Search)

梯度下降的步长选取依赖于函数的光滑程度。而对于非LL-smooth的函数,或是难以确定LLLL-smooth函数,我们就无法使用上面的分析。此时如何选取梯度下降的步长呢?

一个直接的想法是,由于在梯度方向的整条直线上函数依然凸函数, 那么直接下降到这个一元函数的最小值点。这就称为精确线搜索(Exact Line Search),它每次搜索给定“线”上的最小值,然后做迭代。(注意到在这样的迭代过程中,每次到达新的一点后沿原方向的梯度分类就势必是0,不然函数值在该方向上会更小,因此我们每次前进的方向都是互相垂直的。)

如何找到线上的最小值点?这本质上就在问如何找到一个一元凸函数的极小值,一元凸函数的导函数是单调函数,因此我们可以直接对导函数二分。一个更好的方法是用牛顿迭代法求根。牛顿法求根,每次在一点处用切线对函数局部近似,然后跳到切线的零点,不断迭代。在大多数函数上,牛顿迭代法的表现是比较优秀的。

如果一个函数满足μ\mu-强凸和LL-smooth,应用精确线搜索做梯度下降,分析得到收敛速率f(xk)f(x)(1μL)k(f(x0)f(x))f(x_k)-f(x^*) \leq(1-\dfrac{\mu}{L})^k(f(x_0)-f(x^*))。从结果上来看,我们得到了一个和一般梯度一样的收敛速率,因为我们在梯度下降时就会选择η=1L\eta=\dfrac{1}{L}。但好处在于我们不需要提前算出μ\mu或是LL,只需要应用线搜索每一次在梯度方向上用牛顿法求极值,我们就会自动以这样的速率收敛。

对于一般情形,很多时候每次迭代都求出精确最小值的操作还是花费了太多开销了。事实上,我们只希望函数值能在每一次迭代中减小一个充分小的值就好了。困难在于如何刻画“充分”。对于凸函数,我们在xkx_k处选定一个下降方向dkd_k,可以过xkx_kdkd_k方向的切线,那么根据凸函数的一阶条件成立f(xk+ηdk)>f(xk)+ηf(xk),dkf(x_k+\eta d_k)>f(x_k)+\eta \lang\nabla f(x_k),d_k\rang。这时候如果我们向上把切线旋转一个角度(但依然是下降的),那么夹在这两条直线之间的函数值部分就可以认为是充分小的了,因为我们正是用旋转后的切线限制住了函数的下降不至于太小(至少也会下降这条新的直线所割的那么多)。刻画旋转,只需要在斜率上乘一个常数因子α\alpha。我们随意预设一个初始的步长,不妨令η=1\eta=1。不停地迭代ηβη\eta\to \beta \eta,直到f(xk+ηdk)<f(xk)+αηf(xk),dkf(x_k+\eta d_k)<f(x_k)+\alpha \eta \lang\nabla f(x_k),d_k\rang。这个操作就称为“回溯”,因为迭代时步长是指数衰减的,因此我们期待这样的迭代是高效的。这就是回溯线搜索(Backtracking Line Search)。假如函数本身就是LL-smooth的,取α1/2\alpha \leq 1/2,可以证明当步长降到小于1L\dfrac{1}{L}的时候一定已经是下降的了,因此我们的每一次步长都恒大于βL\dfrac{\beta}{L}。进一步的分析也表明,回溯线搜索的收敛速率对于强凸且LL-smooth的函数也是指数级的。

牛顿法(Newton's Method)

在精确线搜索时我们应用牛顿迭代法来求单调函数的零点,然后确定梯度下降的步长。现在我们直接用牛顿法来下降。对于函数f(x)f(x),我们直接在xkx_k点处对ff做二阶泰勒展开f(x)f(xk)+f(xk)(xxk)+f(x)\approx f(x_k)+\nabla f(x_k)^\top (x-x_k)+ 12(xxk)f2(xk)(xxk)\dfrac{1}{2}(x-x_k)^\top \nabla f^2(x_k)(x-x_k),右侧是二次函数,我们直接下降到这个二次函数的极小值点xk+1=xk(2f(xk))1f(xk)x_{k+1}=x_k-(\nabla^2 f(x_k))^{-1}\nabla f(x_k),反复迭代。这就是用来求函数极值的牛顿法。(注意我们需要ff的二阶导是可逆的,这对于严格凸函数是始终成立的,因为此时二阶导是正定的。)

在一般情形下,牛顿法是无法保证正确性的。因为泰勒展开只是局部近似,无法保证全局的正确性。因此,我们需要ff满足额外的条件,在这里主要是控制二阶导的变化范围。

如果2f(x)2f(y)Mxy\|\nabla^2 f(x)-\nabla^2f(y)\|\leq M\|x-y\|,就称2f(x)\nabla^2 f(x)MM-Lipschitz的(前者是矩阵范数,A=maxv=1Av=λmax(AA)\|A\|=\max\limits_{\|v\|=1}\|Av\|=\sqrt{\lambda_{\max}(A^\top A)})。对于μ\mu-强凸以及二阶导是MM-Lipschitz的函数,能够证明xk+1xM2μxkx2\|x_{k+1}-x^*\|\leq\dfrac{M}{2\mu}\|x_k-x^*\|^2。注意到等式右侧的平方项,这意味着xkx(M2μ)2k1x0x2k\|x_k-x^*\| \leq \left(\dfrac{M}{2\mu}\right)^{2k-1}\|x_0-x^*\|^{2k},这比梯度下降收敛得快得多,称为二次收敛。然而,这样得收敛只有在x0x<1\|x_0-x^*\|<1时才能成立,因此只对xx^*邻域内才能生效。相反地,梯度下降虽然不能达到二次收敛,但对于全局都是有效的。

尽管牛顿法并不能保证每一步都下降,但我们可以验证对于严格凸函数,dk=(2f(xk))1f(xk)d_k=-(\nabla^2 f(x_k))^{-1}\nabla f(x_k)确实是一个下降方向,因为方向导数f(xk),dk=f(xk)(2f(xk))1f(xk)\lang \nabla f(x_k),d_k\rang=-\nabla f(x_k)^\top(\nabla^2 f(x_k))^{-1}\nabla f(x_k)是一个二次型,而中间的矩阵是正定的,因此方向导数小于0。只不过,牛顿法每次走过的步长太大了。于是一个直观的想法是,我们在原先牛顿法的dkd_k方向上应用回溯线搜索,保证牛顿法下降。这就是阻尼牛顿法(Damped Newton's Method),它给牛顿法的步长一个阻尼来保证下降。于是我们可以把牛顿法分为两个阶段,先采用阻尼牛顿法,等到接近收敛时,采用一般牛顿法直接二次收敛。对于μ\mu-强凸,二阶导是MM-Lipschitz的函数并且二阶导满足存在常数LL使得LI2f(x)LI-\nabla^2 f(x)是正定的,那么存在δ(0,min{μ2M,3(12α)μ2M})\delta\in(0,\min\{\dfrac{\mu^2}{M},3(1-2\alpha)\dfrac{\mu^2}{M}\})γ=αβ2δ2μL2\gamma=\alpha\beta^2\delta^2\dfrac{\mu}{L^2},对于f(xk)>δ\|\nabla f(x_k)\|>\delta时我们用阻尼牛顿法保证f(xk+1)f(xk)γf(x_{k+1})-f(x_k)\leq-\gamma,对于f(xk)δ\|\nabla f(x_k)\|\leq \delta时用一般牛顿法保证二次收敛。

另一个控制二阶导范围的方法是考虑三阶导的范围。对于一元函数ff,如果满足f(x)2f(x)3/2|f'''(x)|\leq 2f''(x)^{3/2}恒成立,就称其为自和谐(Self-concordant)的。它可以等价描述为ddx1f(x)1\left|\dfrac{d}{dx}\dfrac{1}{\sqrt{f''(x)}}\right|\leq 1。对于多元函数,如果沿任意直线都是自和谐的就称多元函数是自和谐的。对于自和谐的严格凸函数,我们可以分析阻尼牛顿法的收敛速率。定义λ(x)=f(x)(2f(x))1f(x)\lambda(x)=\sqrt{\nabla f(x)^\top (\nabla^2f(x))^{-1}\nabla f(x)},那么存在δ(0,1/4]\delta \in(0,1/4],和γ>0\gamma>0(只关于回溯线搜索中的α,β\alpha,\beta),满足:当λ(xk)δ\lambda(x_k)\geq \delta时,阻尼牛顿法保证f(xk+1)f(xk)γf(x_{k+1})-f(x_k)\leq -\gamma;否则,一般牛顿法保证λ(xk+1)22λ(xk)2\lambda(x_{k+1})^2\leq2\lambda(x_k)^2。(另一个好处在于,以上收敛分析是仿射不变的,而关于MM-Lipschitz的分析却不是)

近端梯度下降(Proximal Gradient Descent)

对于不可微的函数,梯度下降法和牛顿法都失效了。现在我们要开始讨论如何求解不可微的凸函数的极小值。

一个有趣的观察是,梯度下降可以“看作”是在使用牛顿法。因为当我们迭代xk+1=xkηf(xk)x_{k+1}=x_k-\eta \nabla f(x_k)时,我们可以把xkηf(xk)x_k-\eta\nabla f(x_k)想象成某个二次函数的极小值点,那么这个二次函数可以构造为f^(x)=f(xk)+f(xk),xxk\hat f(x)=f(x_k)+\lang \nabla f(x_k),x-x_k\rang +12ηxxk2+\dfrac{1}{2\eta}\|x-x_k\|^2。因此梯度下降法可以看作是对f^(x)\hat f(x)做牛顿法。

现在假设f(x)f(x)本身是不可微的,但是可以拆成一个可微函数g(x)g(x)和一个不可微函数h(x)h(x)。(我们要求f,g,hf,g,h都是凸的)既然梯度下降的迭代可以看作求f(xk)+f(xk),xxk+12ηxxk2f(x_k)+\lang \nabla f(x_k),x-x_k\rang+\dfrac{1}{2\eta}\|x-x_k\|^2的最小值,那么现在我们可以用g(xk)+g(xk),xxk+12ηxxk2+h(x)g(x_k)+\lang \nabla g(x_k),x-x_k\rang+\dfrac{1}{2\eta}\|x-x_k\|^2+h(x)来近似,求它的最小值。配方化简后等价于求12ηx(xkηg(xk))2+h(x)\dfrac{1}{2\eta}\|x-(x_k-\eta\nabla g (x_k))\|^2+h(x)的最小值。定义近端梯度算子proxh(y)=argmin12xy2+h(x)\text{prox}_h(y)=\arg\min \dfrac{1}{2}\|x-y\|^2+h(x),可以写出迭代:xk+1=proxηh(xkηg(xk))x_{k+1}=\text{prox}_{\eta h}(x_k-\eta\nabla g(x_k))。这就称为近端梯度下降算法。

近端梯度下降算法有一个更高观点的理解方法:在离散近似梯度流算法中,可以用xk+1=xkηf(xk)x_{k+1}=x_k-\eta\nabla f(x_k)近似,也可以用xk+1=xkηf(xk+1)x_{k+1}=x_k-\eta \nabla f(x_{k+1})近似,后者恰好可以等价地写成xk+1=proxηf(xk)x_{k+1}=\text{prox}_{\eta f}(x_k)(因为(12xxk2+ηf(x))=xxk+ηf(x)\nabla (\dfrac{1}{2}\|x-x_k\|^2+\eta f(x))=x-x_k+\eta f(x))。因此近端梯度下降本质上就是轮换地用两种方法来近似梯度流算法。

我们看到近端梯度下降最关键的就是要能合适的选取h(x)h(x)使得极端梯度算子是可以计算的。对于复杂的h(x)h(x),我们是很难求出极端梯度算子的。这里举一个可以应用近端梯度下降的例子,称为LASSO(Least Absolute Shrinkage and Selection
Operator),在这里f(β)=yXβ22+λβ1f(\beta)=\|y-X\beta\|_2^2+\lambda\|\beta\|_1。我们注意到1-范数是不可微的,因此选取h(β)=λβ1h(\beta)=\lambda\|\beta\|_1,在求proxh(γ)\text{prox}_h(\gamma)我们发现每一维是独立的,对于第ii个坐标分类,可以分类讨论求得proxh(γ)i\text{prox}_h(\gamma)^i的取值当γi>λ\gamma_i>\lambda时为γiλ\gamma_i-\lambda,当γi<λ\gamma_i<-\lambda时为γi+λ\gamma_i+\lambda,其余为0。我们用符号Sλ(γi)\mathcal{S}_\lambda(\gamma_i)来表示这个函数,称为软阈值算子(Soft Thresholding Operator)。应用软阈值算子的近端梯度下降算法求出LASSO的最小值的算法称为ISTA(Iterative Soft-Thresholding Algorithm)。

下面要分析近端梯度下降算法的收敛情况。注意到hh是不可微的,为此我们要引入次梯度(Subgradient)的概念,如果对于每个xx都存在一个向量vxv_x成立f(y)f(x)+vx,yxf(y)\geq f(x)+\lang v_x,y-x \rang,就称vxv_xff的次梯度,记为v\partf(x)v \in \part f(x)。对于凸函数,我们看到次梯度就扮演着一阶条件中梯度的角色。最后我们分析得到对于f(x)=g(x)+h(x)f(x)=g(x)+h(x),如果ggLL-smooth且令η1L\eta \leq \dfrac{1}{L},那么f(xT)f(x)x0x22Tηf(x_T)-f(x^*)\leq\dfrac{\|x_0-x^*\|^2}{2T\eta}。如果ggμ\mu-强凸的,那么有xk+1x2(1μη)xkx2\|x_{k+1}-x^*\|^2 \leq (1-\mu\eta)\|x_k-x^*\|^2。(如果无法求出LL,我们可以用线搜索得到与梯度下降一样的收敛结果)。