优化基本理论与方法(12)约束集上的优化问题
信息来源:网络    时间:2024-08-12 02:22

封面图来自 pixiv.net/artworks/7853

本文中我们介绍两个经典的求解受限优化的方法:梯度投影法和 Frank--Wolfe 方法。

梯度投影法(Gradient Projection method)在 20 世纪 60 年代由 Goldstein、Levitin 和 Polyak 等人发明,读者可以考古 Dimitri P. Bertsekas. (1974). On the Goldstein-Levitin-Polyak gradient projection method. IEEE Conference on Decision and Control. ieeexplore.ieee.org/doc

与无约束的问题相比,此时目标函数的梯度应该被区别对待;因为直接像梯度方法那样进行迭代,得到的结果可能在可行集之外。不过就像次梯度方法那样,我们可以想办法找到梯度的替代品。

在求解 \\mathcal{F}^{1,1}_{L}(\\mathbb{R}^n) 中函数的极小化问题时,在梯度方法中,我们用到了两个重要的性质。第一个性质是沿梯度反向的一步下降使得函数值减小了与梯度范数平方相当的量,

\\begin{align}f\\left(\\bm x-\\frac{1}{L}\
abla f(\\bm x)\\right)\\le f(\\bm x)-\\frac{1}{2L}\\lVert\
abla f(\\bm x)\\rVert^2 \\end{align}

第二个性质是不等式:

\\begin{align}\\langle\
abla f(\\bm x),\\bm x-\\bm x^*\\rangle \\ge \\frac 1L \\lVert\
abla f(\\bm x)\\rVert^2 \\end{align}

我们可以找到继承这两个重要性质的对象。取定 \\gamma > 0 ,定义

\\begin{align}p_Q(\\bar{\\bm x};\\gamma)&=\\operatorname*{argmin}_{\\bm x\\in Q}\\left(f(\\bar{\\bm x})+\\langle \
abla f(\\bar{\\bm x}),\\bm x-\\bar{\\bm x}\\rangle+\\frac{\\gamma}{2}\\lVert\\bm x-\\bar{\\bm x}\\rVert^2\\right)\\\\ g_Q(\\bar{\\bm x};\\gamma)&=\\gamma\\cdot(\\bar{\\bm x}-p_Q(\\bar{\\bm x};\\gamma)) \\end{align}

我们在未来介绍近端梯度方法时还会见到这个形式,它们分别被称作近邻算子和梯度映射。另外,我们可以将这两个定义作如下改写:

\\begin{align}p_Q(\\bar{\\bm x};\\gamma)&=\\operatorname*{argmin}_{\\bm x\\in Q}\\left(f(\\bar{\\bm x})+\\langle \
abla f(\\bar{\\bm x}),\\bm x-\\bar{\\bm x}\\rangle+\\frac{\\gamma}{2}\\lVert\\bm x-\\bar{\\bm x}\\rVert^2\\right)\\\\ &=\\operatorname*{argmin}_{\\bm x\\in Q}\\Bigg(\緻brace{f(\\bar{\\bm x})-\\frac{1}{2\\gamma}\\|\
abla f(\\bar{\\bm x})\\|^2}_{\	ext{与 }\\bm x\	ext{ 无关,可以丢掉}}+\\frac{\\gamma}{2}\\left\\|\\bm x-\\bar{\\bm x}+\\frac{1}{\\gamma}\
abla f(\\bar{\\bm x})\\right\\|^2\\Bigg)\\\\ &=\\operatorname*{argmin}_{\\bm x\\in Q}\\left\\|\\bm x-\\bar{\\bm x}+\\frac{1}{\\gamma}\
abla f(\\bar{\\bm x})\\right\\|^2\\\\ &=\\pi_Q\\left(\\bar{\\bm x}-\\frac{1}{\\gamma}\
abla f(\\bar{\\bm x})\\right) \\end{align}

其中 \\pi_Q 代表在 Q 上的欧几里得投影,定义为

\\begin{align}\\pi_Q(\\bm x_0)=\\operatorname*{argmin}_{\\bm x\\in Q}\\lVert\\bm x-\\bm x_0\\rVert\\end{align}

另外,假设 \\mathrm{dom}\\,f=\\mathbb{R}^n ,那么 p_Q(\\bar{\\bm x};\\gamma) 对任意的 \\bar{\\bm x}\\in\\mathbb{R}^n 都有定义,而不必限定于 \\bar{\\bm x}\\in Q 。另外,根据第 10 篇文章的定理 13,向凸集的投影是唯一的,因此 p_Q(\\bar{\\bm x};\\gamma) 也是唯一的。

有如下定理:

定理 1:令 f\\in\\mathcal{S}_{\\mu,L}^{1,1}(Q)\\gamma\\ge L , 则对 \\forall\\,\\bm x,\\bar{\\bm x}\\in Q ,我们有

\\begin{align}f(\\bm x)\\geq f(p_Q(\\bar{\\bm x};\\gamma))+\\langle g_Q(\\bar{\\bm x};\\gamma),\\bm x-\\bar{\\bm x}\\rangle+\\frac{1}{2\\gamma}\\| g_Q(\\bar{\\bm x};\\gamma)\\|^2+\\frac{\\mu}{2}\\| \\bm x-\\bar{\\bm x}\\|^2 \\end{align}

这个定理的证明见第 13 篇文章中的结论 2。根据该定理,分别令 \\bm x=\\bar{\\bm x}\\bm x=\\bm x^* ,可以得到如下推论:

推论:令 f\\in\\mathcal{S}_{\\mu,L}^{1,1}(Q)\\gamma\\ge L ,且 \\bar{\\bm x}\\in Q ,则

\\begin{align}&f(p_Q(\\bar{\\bm x};\\gamma))\\le f(\\bar{\\bm x})-\\frac{1}{2\\gamma}\\lVert g_Q(\\bar{\\bm x};\\gamma)\\rVert^2\\\\ &\\langle g_Q(\\bar{\\bm x};\\gamma),\\bar{\\bm x}-\\bm x^*\\rangle\\ge \\frac{1}{2\\gamma}\\lVert g_Q(\\bar{\\bm x};\\gamma)\\rVert^2+\\frac{\\mu}{2}\\|\\bar{\\bm x}-\\bm x^*\\|^2+\\frac{\\mu}{2}\\| p_Q(\\bar{\\bm x};\\gamma)-\\bm x^*\\|^2 \\end{align}

这说明 g_Q(\\bar{\\bm x};\\gamma) 有我们希望的类似于梯度的性质,而且对任意的 \\bar{\\bm x} ,一定有 g_Q(\\bar{\\bm x};\\gamma)\\in Q 。自然,我们得到了如下算法:

\\Large \	extbf{梯度投影法}\\\\  \\rule[0pt]{8cm}{0.05em}\\\\ \\begin{align}0.& \	ext{ 选择初始点 }\\bm x_0\\in Q\	ext{ 和参数 }\\gamma >0\\\\ 1.&\	ext{ 第 }k\	ext{ 次迭代 }(k\\ge 0):\\\\ &\\qquad \\bm x_{k+1}=\\bm x_k+\\frac{1}{\\gamma}g_Q(\\bm x_k;\\gamma) \\end{align}

注意到在这个算法中,

\\begin{align}\\bm x_{k+1}=p_Q(\\bm x_k;\\gamma)=\\pi_Q\\left(\\bm x_k-\\frac{1}{\\gamma}\
abla f(\\bm x_k)\\right) \\end{align}

当测试点 \\bm x_k\\in Q 时,由于 \\pi_Q(\\bm x_k)=\\bm x_k ,那么此时 \\begin{align}p_Q(\\bm x_k;\\gamma)=\\bm x_k-\\frac{1}{\\gamma}\
abla f(\\bm x_k)\\end{align}g_Q(\\bm x_k;\\gamma)=\
abla f(\\bm x_k) 。因此,我们可以这样理解梯度投影法:

也就是说,我们通过这种投影操作保证了每一步得到的新测试点都位于 Q 内。

Trendafilov, Nickolay. (2006). The Dynamical System Approach to Multivariate Data Analysis. Journal of Computational and Graphical Statistics. 15. 10.1198/106186006X130828.

该算法的效率分析与无约束条件下的分析十分相似,有如下定理:

定理 2:令 f\\in\\mathcal{S}_{\\mu,L}^{1,1}(Q) ,如果在投影梯度法中, \\begin{align}\\gamma\\ge \\frac{L+\\mu}{2}\\end{align} ,则

\\begin{align}\\|\\bm x_k-\\bm x^*\\|\\le \\left(1-\\frac{\\mu}{\\gamma}\\right)^k\\lVert\\bm x_0-\\bm x^*\\|\\end{align}

证明:设 r_k=\\|\\bm x_k-\\bm x^*\\| ,则

\\begin{align}r_{k+1}^{2}&=\\quad\\left\\|\\pi_{Q}\\left(\\bm x_{k}-\\frac{1}{\\gamma}\
abla f\\left(\\bm x_{k}\\right)\\right)-\\pi_{Q}\\left(\\bm x^{*}-\\frac{1}{\\gamma}\
abla f\\left(\\bm x^{*}\\right)\\right)\\right\\|^{2}\\\\ &\\le\\quad\\left\\|\\bm x_k-\\bm x^*-\\frac1\\gamma(\
abla f(\\bm x_k)-\
abla f(\\bm x^*))\\right\\|^2 \\\\ &=\\quad r_{k}^{2}-\\frac{2}{\\gamma}\\big\\langle\
abla f(\\bm x_{k})-\
abla f(\\bm x^{*}),\\bm x_{k}-\\bm x^{*}\\big\\rangle+\\frac{1}{\\gamma^{2}}\\|\
abla f(\\bm x_{k})-\
abla f(\\bm x^{*})\\|^{2}\\\\ &\\le\\left(1-\\frac{2}{\\gamma}\\cdot\\frac{\\mu L}{\\mu+L}\\right)r_{k}^{2}+\\left(\\frac{1}{\\gamma^{2}}-\\frac{2}{\\gamma}\\cdot\\frac{1}{\\mu+L}\\right)\\|\
abla f(\\bm x_{k})-\
abla f(\\bm x^{*})\\|^{2}\\\\ &\\le\\left(1-\\frac2\\gamma\\cdot\\frac{\\mu L}{\\mu+L}+\\mu^{2}\\left(\\frac1{\\gamma^{2}}-\\frac2\\gamma\\cdot\\frac1{\\mu+L}\\right)\\right)r_{k}^{2}\\\\&=\\left(1-\\frac\\mu\\gamma\\right)^{2}r_{k}^{2}\\end{align}

其中第二行的不等式来自第 10 篇文章的推论 2;第四行的不等式来自于第 4 篇文章的式 (3)

再根据 Lipschitz 连续性,有

\\begin{align}f(\\bm x_k)-f(\\bm x^*)&\\le \\langle \
abla f(\\bm x^*),\\bm x_k-\\bm x^*\\rangle +\\frac{L}{2}\\|\\bm x_k-\\bm x^*\\|^2\\\\ &=\\frac{L}{2}\\|\\bm x_k-\\bm x^*\\|^2\\\\ &\\le \\frac{L}{2}\\left(1-\\frac{\\mu}{\\gamma}\\right)^{2k}\\lVert\\bm x_0-\\bm x^*\\|^2 \\end{align}

可以看到,强凸函数保证了该方法以线性速度收敛;如果选择 \\begin{align}\\gamma=\\frac{\\mu+L}{2}\\end{align} ,那么该方法与普通的梯度方法有相同的收敛速度。

至此我们仅仅分析了强凸函数,其收敛速度为 O(\\log\\frac{1}{\\varepsilon}) ;此外,凸函数时的收敛速度也与梯度方法相同,为 O(\\frac 1\\varepsilon) ;并且,只要 f 满足 Lipschitz 光滑,该方法都是收敛的。我们有如下定理:

定理 3:假设 \
abla fL -Lipschitz 连续的,设 \\{\\bm x_k\\} 是由梯度投影法得到的序列,那么其极限点一定是 f 的稳定点。

证明见 C. T. Kelley. (1999). Iterative Methods for Optimization. SIAM, Philadelphia. 的 Theorem 5.4.6。

最后,值得一提的是,投影梯度法在深度学习中的对抗样本攻击(adversarial attack) 中非常常用。其思路是:在给定干净样本 (\\bar{\\bm x};\\bar{\\bm y})\\bar{\\bm y} 是标签)的条件下,在其上加一个不超过 \\varepsilon 扰动,使得模型得到错误的结果;亦即,在集合 Q=B_2(\\bar{\\bm x},\\varepsilon) 内最大化目标模型的损失函数 \\mathcal{L}(\\bm \	heta;\\bm x;\\bar{\\bm y}) 。作为一个受限的最大化问题,使用投影梯度法是很自然的。读者可以参考论文:Madry, A., Makelov, A., Schmidt, L., Tsipras, D., & Vladu, A. (2017). Towards Deep Learning Models Resistant to Adversarial Attacks. ArXiv, abs/1706.06083。

梯度投影法有一个缺陷,就是“投影”这一操作可能并不容易实现。事实上,投影梯度法最适用于那些“简单”的闭凸集,例如非负象限约束、 n 维盒子约束、单纯形约束、 \\ell_1 球或者 \\ell_2 球约束等,从而让投影的结果(或者 g_Q(\\bar{\\bm x};\\gamma) )可以直接有闭式解。否则,计算投影相当于计算一个最小化问题,这样的开销实在不可接受。

而我们接下来要介绍的算法则完全没有这种烦恼,它直接去除了投影的操作。

Frank--Wolfe (FW) 算法,又称条件梯度(conditional gradient)算法,是解决非线性受限优化问题的一种经典求解方法。首次发表于 Frank, Marguerite & Philip Wolfe. (1956). An algorithm for quadratic programming. Naval Research Logistics。该方法内存消耗低,而且最重要的是无需投影操作,因此同样很受欢迎。

传统的 Frank--Wolfe 算法可以被用来求解问题\\begin{align}\\min_{\\bm x\\in Q}f(\\bm x) \\end{align} ,其中 Q 是一个紧的凸集, f\\in C_{L}^{1,1}(Q) ,即 f 的梯度满足 Lipschitz 条件。所谓紧集(compact set),可以理解为在闭集之上加了一个有界条件(在欧式空间这紧急都等价于有界闭集)。根据前文投影梯度法中 p_Q(\\bm x;\\gamma) 的定义,可以看到投影梯度法是采用了目标函数的局部二次逼近来更新:

\\begin{align}\\bm x_{k+1}&=p_Q(\\bm x_k;\\gamma)\\\\ &=\\operatorname*{argmin}_{\\bm x\\in Q}\\left(f(\\bm x_k)+\\langle \
abla f(\\bm x_k),\\bm x-\\bm x_k\\rangle+\\frac{\\gamma}{2}\\lVert\\bm x-\\bm x_k\\rVert^2\\right) \\end{align}

而如果去掉其二次项只做一次逼近,那么更新就变成了一个线性问题:

\\begin{align*}\\bm{x}_{k+1}&=\\operatorname*{\\arg\\min}\\limits_{\\bm{x}\\in Q}\\big(f(\\bm x_k)+\\langle\
abla f(\\bm x_k),\\bm{x}-\\bm{x}_k\\rangle\\big)\\\\ &=\\operatorname*{\\arg\\min}\\limits_{\\bm{x}\\in Q}\\langle \
abla f(\\bm x_k),\\bm{x}\\rangle \\end{align*}

这一过程通常被称作线性最小化 Oracle(Linear Minimization Oracle, LMO)或者线性最优 Oracle。LMO 也可以看做是求 Q 中与负梯度方向 -\
abla f(\\bm x_k) 内积最大的方向;相比于投影,它的复杂度会降低很多,这就 Frank--Wolfe 算法的核心想法。也就是说,FW 算法免除了投影操作,每次迭代中只需要依赖一个求解限制域上的线性优化问题。对任意给定的集合 Q ,我们很容易能确定计算 LMO 的复杂度,例如如果 X 是一个单纯形或者凸多面体(对应线性约束),那么求解 LMO 就变为了一个线性规划(Linear Programming)问题,用单纯形法或者内点法等可以迅速求解(虽然单纯形法最差需要指数级的运算,但是一般情况下求解器都能很快得到结果)。

近几年来,FW 受到了机器学习和优化领域的广泛关注,主要由于如 下原因:

Frank--Wolfe 算法形式十分简单,我们直接来看其伪代码:

\\begin{align}&\\Large\\mathbf{Frank--Wolfe 算法}\\\\ &\\rule[2pt]{11.5cm}{0.08em}\\\\  &\\begin{array}{cl}\\mathbf{Input}:&\	ext{initial point }\\bm x_0\\in Q,\\\\ &\	ext{upper bound of error }\\delta > 0\\\\ \\end{array}\\\\ &\\rule[2pt]{11cm}{0.05em}\\\\ &\\mathbf{for}\\;k=0,\\dots\\;\\mathbf{do}:\\\\ &\\qquad \\bm s_k=\\operatorname*{argmin}_{\\bm s\\in Q}\\langle \
abla f(\\bm x_k),\\bm s\\rangle\\\\ &\\qquad \\bm d_k=\\bm s_k-\\bm x_k\\\\ &\\qquad g_k=-\\langle\
abla f(\\bm x_k),\\bm d_k\\rangle\\\\ &\\qquad\\mathbf{if}\\; g_k < \\delta:\\\\ &\\qquad\\qquad \\mathbf{return}\\;\\bm x_k\\\\ &\\qquad \\mathbf{end\\;if}\\\\ &\\qquad \	ext{Update step size }\\gamma_k\	ext{, e.g., using line search.}\\\\ & \\qquad \\bm x_{k+1}=\\bm x_k+\\gamma_k\\bm d_{k}\\\\ &\\mathbf{end\\;for}\\\\ \\end{align}

其中 g_k 是一个评价函数,用来与 \\delta 做比较,判断程序何时退出;它的正式名字是 Frank-Wolfe 间隔,我们在“收敛性分析”章节会详细介绍。\\begin{align}\\bm s_k=\\operatorname*{argmin}_{\\bm s\\in Q}\\langle \
abla f(\\bm x_k),\\bm s\\rangle \\end{align} 是该算法的核心,即求解 LMO。不过,求解 LMO 的过程中没有用到步长因素,因此还不能将其结果直接作为 \\bm x_{k+1} 。所以我们要在 \\bm s_k 的方向上按照步长 \\gamma_k 迭代得到 \\bm x_{k+1} ,即 \\bm x_{k+1}=\\bm x_k+\\gamma_k(\\bm s_k-\\bm x_k) 。因此可以将 FW 算法提炼成两个步骤:

  1. 计算 LMO: \\begin{align}\\bm s_k\\in \\operatorname*{argmin}_{\\bm s\\in Q}\\langle \
abla f(\\bm x_k),\\bm s\\rangle \\end{align} ;
  2. \\bm x_{k+1}=\\bm x_k+\\gamma_k(\\bm s_k-\\bm x_k)

步长系数 \\gamma_k 的更新方式有很多种,这里给出常见的两种方式:

第二种看上去十分奇怪,后面我们会解释。

此外,还有一种预先定义步长的方式:

如果限制域是某种特殊的形式,那么我们可以直接得到 \\bm s_k 的闭式解。例如,如果限制域是由某个范数 \\|\\cdot\\| 定义的 Q=\\{\\|\\bm x\\|\\le \	au\\} ,那么 LMO 步骤可以简化成

\\begin{align}\\bm s_k &\\in\\operatorname*{argmin}_{\\|\\bm s\\|\\le \	au}\\langle\
abla f(\\bm x_k),\\bm s\\rangle\\\\ &=-\	au\\operatorname*{argmax}_{\\|\\bm s\\|\\le 1}\\langle \
abla f(\\bm x_k),\\bm s\\rangle\\\\ &=-\	au\\partial \\|\
abla f(\\bm x_k)\\|_* \\end{align}

其中 \\|\\cdot\\|_* 代表对偶范数(dual norm)。换言之,如果我们知道对偶范数的次梯度,那么 FW 算法即可十分简单地进行迭代更新。

何为对偶范数呢?假设 \\|\\cdot\\| 表示 \\mathbb{R}^n 上的范数,其对应的对偶范数有如下定义:

\\begin{align}\\|\\bm x\\|_*=\\sup\\{\\langle\\bm x,\\bm z\\rangle\\,|\\,\\|\\bm x\\|\\le 1\\}\\end{align}

其等价定义为 \\begin{align}\\|\\bm x\\|_*=\\sup_{\\|\\bm x\\|\\le 1}\\langle\\bm x,\\bm z\\rangle=\\sup_{\\|\\bm x\\|=1}\\langle\\bm x,\\bm z\\rangle=\\sup_{\\bm x\
e \\bm 0}\\frac{\\langle\\bm x,\\bm z\\rangle}{\\|\\bm x\\|}\\end{align} 。考察 \\mathbb{R}^n 上的一个线性泛函 f_{\\bm z}(\\bm x)=\\bm z^\	op\\bm x ,其算子范数正是对偶范数 \\begin{align}\\|\\bm z\\|_*=\\sup_{\\bm x\
e \\bm0}\\frac{\\bm z^\	op\\bm x}{\\| \\bm x\\|}\\end{align} 。熟悉泛函分析的同学可以得知,对于 Banach 空间 (\\mathbb{R},\\|\\cdot\\|) ,其对偶空间对 \\|\\cdot\\|_* 诱导的度量也是完备的,即其对偶空间是 (\\mathbb{R},\\|\\cdot\\|_*)

至于范数的次梯度,我们可以证明如下定理

定理 4\\partial \\|\\bm x\\|=\\{\\bm v\\,|\\,\\langle\\bm v,\\bm x\\rangle=\\|\\bm x\\|,\\|\\bm v\\|_*\\le 1\\}

设集合 V(\\bm x)=\\{\\bm v\\,|\\,\\langle\\bm v,\\bm x\\rangle=\\|\\bm x\\|,\\|\\bm v\\|_*\\le 1\\} ,我们仅需证明 V(\\bm x)\\subseteq \\partial \\|\\bm x\\| 以及 V(\\bm x)\\supseteq \\partial \\|\\bm x\\| 即可。首先证明 V(\\bm x)\\subseteq\\partial\\|\\bm x\\| ,对于 \\forall\\,\\bm v\\in V(\\bm x)\\forall\\,\\bm y\\in \\mathbb{R}^n ,有

\\begin{align}\\|\\bm x\\|+\\langle \\bm v,\\bm y-\\bm x\\rangle&=\\|\\bm x\\|+\\langle\\bm v,\\bm y\\rangle-\\langle\\bm v,\\bm x\\rangle\\\\ &=\\|\\bm x\\|+\\langle\\bm v,\\bm y\\rangle-\\|\\bm x\\|\\\\ &=\\langle\\bm v,\\bm y\\rangle\\\\ &\\le \\|\\bm y\\|\\|\\bm v\\|_*\\\\ &\\le \\|\\bm y\\| \\end{align}

因此 \\bm v\\in\\partial\\|\\bm x\\| ;其中第四行的不等号来自 \\begin{align}\\langle\\bm v,\\bm y\\rangle=\\|\\bm y\\|\\left(\\frac{\\bm v^\	op\\bm y}{\\|\\bm y\\|}\\right)\\le \\|\\bm y\\|\\|\\bm v\\|_* \\end{align}

接下来证明 \\partial\\|\\bm x\\|\\subseteq V(\\bm x) 。设 \\bm v\\in\\partial\\|\\bm x\\| ,根据次梯度定义,对 \\forall\\,\\bm y\\in\\mathbb{R}^n ,有

\\begin{align}&\\|\\bm y\\|\\ge\\|\\bm x\\|+\\langle\\bm v,\\bm y-x\\rangle\\\\ \\Leftrightarrow\\;& \\langle\\bm v,\\bm y\\rangle-\\|\\bm y\\|\\le \\langle\\bm v,\\bm x\\rangle-\\|\\bm x\\| \\end{align}

该式对所有 \\bm y 都成立,因此

\\begin{align}\\sup_{\\bm y}\\big(\\langle\\bm v,\\bm y\\rangle-\\|\\bm y\\|\\big)\\le \\langle\\bm v,\\bm x\\rangle-\\|\\bm x\\| \\end{align}

其中左侧正是 \\|\\cdot\\|\\bm v 处的 Fenchel 共轭,熟悉这个概念的同学可以立刻得知其是对偶范数的单位球之示性函数,即

\\sup\\limits_{\\bm y}\\big(\\langle\\bm v,\\bm y\\rangle-\\|\\bm y\\|\\big)=\\left\\{\\begin{array}{cl}0,& \\|\\bm v\\|_*\\le 1\\\\ +\\infty,& \\|\\bm v\\|_*> 1 \\end{array}\\right.

不熟悉的读者也不要慌,这一结论很容易证明(Stephen Boyd, Convex Optimization 的 Example 3.26):当 \\|\\bm v\\|_*\\le 1 时, \\langle\\bm v,\\bm y\\rangle-\\|\\bm y\\|\\le\\|\\bm y\\|\\|\\bm v\\|_*-\\|\\bm y\\|=(\\|\\bm v\\|_*-1)\\|\\bm y\\|\\le 0 ;当 \\|\\bm v\\|_*> 1 时,根据对偶范数定义,我们知道存在 \\bm y\\|\\bm y\\|\\le 1 ,使得 \\langle\\bm v,\\bm y\\rangle> 1 ,因此 \\langle\\bm v,\\bm y\\rangle-\\|\\bm y\\|>0 ,取 \\bm z=t\\bm yt>0 ,可知 \\langle\\bm v,\\bm z\\rangle-\\|\\bm z\\|=t\\big(\\langle\\bm v,\\bm y\\rangle-\\|\\bm y\\|\\big) 可以任意大。

总结一下,我们目前有

\\langle\\bm v,\\bm x\\rangle-\\|\\bm x\\|\\ge\\sup\\limits_{\\bm y}\\big(\\langle\\bm v,\\bm y\\rangle-\\|\\bm y\\|\\big)=\\left\\{\\begin{array}{cl}0,& \\|\\bm v\\|_*\\le 1\\\\ +\\infty,& \\|\\bm v\\|_*> 1 \\end{array}\\right.

如果 \\|\\bm v\\|_*>1 那么 \\langle \\bm v,\\bm x\\rangle-\\|\\bm x\\| 就恒为 +\\infty ,这显然是不可能的。因此,我们有 \\|\\bm v\\|_*\\le 1 ,此时 \\langle\\bm v,\\bm x\\rangle-\\|\\bm x\\|\\ge 0 。那么有

\\begin{align}0 \\le \\langle \\bm v,\\bm x\\rangle-\\|\\bm x\\|  \\le \\|\\bm x\\|\\|\\bm v\\|_*-\\|\\bm x\\|=\\|\\bm x\\|(\\|\\bm v\\|_*-1) \\le 0 \\end{align}

这说明 \\langle\\bm v,\\bm x\\rangle=\\|\\bm x\\| ,所以 \\bm v\\in V(\\bm x) 。定理得证。

定理 5\\ell_p 范数的对偶范数是 \\ell_q 范数,其中 p\\in(1,\\infty) ,满足 \\begin{align}\\frac{1}{p}+\\frac{1}{q}=1\\end{align} 。特别地, \\ell_1 的对偶范数是 \\ell_\\infty ,反过来也成立; \\ell_2 的对偶范数是其自身。

该定理根据霍尔德不等式(H?lder inequality)可以立即得到。特别地,对于 \\ell_1\\ell_\\infty ,我们有

\\begin{align}\\sup\\{\\langle \\bm x,\\bm z\\rangle\\,|\\,\\|\\bm x\\|_{\\infty }\\leq 1\\}=\\sum _{i=1}^{n}\\left| z^{(i)}\\right|=\\|\\bm z\\|_{1}\\end{align}

\\begin{align}\\sup\\{\\langle \\bm x,\\bm z\\rangle\\,|\\,\\|\\bm x\\|_1\\leq 1\\}=\\sup\\left\\{\\left.\\sum _{i=1}^n x^{(i)}z^{(i)}\\right|\\,\\sum_{i=1}^nx^{(i)}=1\\right\\}=\\max_{i=1,\\dots,n}z^{(i)}=\\|\\bm z\\|_\\infty\\end{align}


我们来看几个例子

\\begin{align}&\\min_{\\bm x}f(\\bm x)\\\\ \	ext{s.t.}\\;& \\|\\bm x\\|_1\\le t \\end{align}

有: \\bm s_k\\in-t\\partial \\|\
abla f(\\bm x_k)\\|_\\infty ;由于 \\partial\\|\\bm z\\|_\\infty=\\{\\bm x\\,|\\,\\langle\\bm x,\\bm z\\rangle=\\|\\bm z\\|_\\infty,\\,\\|\\bm x\\|_1\\le 1\\} ,易得 \\bm x=\\bm e_j 或者 \\bm x=-\\bm e_j ,其中 j 是最大的 \\begin{align}\\left|z^{(j)}\\right|\\end{align} 对应的坐标。

FW 更新即为:

\\begin{align}i_k&\\in\\operatorname*{argmax}_{i=1,\\dots p}\\left|(\
abla f(\\bm x_k))^{(i)}\\right|\\\\ \\bm s_{k+1}&=-t\\,\\mathrm{sign}\\left((\
abla f(\\bm x_k))^{(i_k)}\\right)\\bm e_{i_k}\\end{align}

这个形式十分类似于贪心坐标下降(greedy coordinate descent),只不过是步长选择有所不同。相比于在 \\ell_1 球上的投影操作,这个方法更加简单,虽然二者都需要 O(n) 的复杂度。

(投影到 \\ell_1 球上的算法复杂度,请参考 Duchi, J.C., Shalev-Shwartz, S., Singer, Y., & Chandra, T. (2008). Efficient projections onto the l1-ball for learning in high dimensions.International Conference on Machine Learning. icml2008.cs.helsinki.fi;作者提出了两个算法,复杂度分别为 O(n\\log n) (因为要排序)和 O(n) 。)

\\begin{align}&\\min_{\\bm x}f(\\bm x)\\\\ \	ext{s.t.}\\;& \\|\\bm x\\|_p\\le t \\end{align}

\\bm s_k=-t\\partial\\|\
abla f(\\bm x_k)\\|_q ,其中 \\begin{align}\\frac{1}{p}+\\frac{1}{q}=1 \\end{align} 。Frank--Wolfe 更新方向即为:

\\begin{align}\\bm s_{k}^{(i)}&=-t\\mathrm{~sign}\\left((\
abla f(\\bm x_k))^{(i)}\\right)\\left|(\
abla f(\\bm x_k))^{(i)}\\right|^{p/q},\\enspace i=1,\\ldots,n \\end{align}

对于给定的 p ,该操作相比于 \\ell_p 球上的投影操作简单很多,只需要 O(n) 步即可计算完成;除了一些特殊情况(如 p=1, 2,\\infty ), \\ell_p 球上的投影操作通常不能直接求解,需要当成一个优化子问题进行求解。

\\begin{align}&\\min_{\\bm X}f(\\bm X)\\\\ \	ext{s.t.}\\;& \\|\\bm X\\|_{tr}\\le t \\end{align}

其中矩阵 \\bm X 的迹范数表示矩阵的所有奇异值之和,它在机器学习中最为经典的应用莫过于用迹范数最小化来代替秩最小化。其具体定义为:给定任意矩阵 \\boldsymbol{X}\\in\\mathbb{R}^{m\	imes n} ,若 r=\\min\\{m,n\\} ,且 \\sigma_1\\geq\\sigma_2\\geq\\cdots\\geq \\sigma_r 为矩阵 \\boldsymbol{X} 的奇异值,则矩阵 \\boldsymbol{X} 的迹范数为

\\begin{align}\\|\\boldsymbol{X}\\|_{tr}=\\mathrm{tr}(\\bm X^*\\bm X)=\\sum_{i=1}^r \\sigma_i \\end{align}

Given thar singular values are non-negative, the trace norm  \\|\\boldsymbol{\\Theta}\\|_{tr} can be regarded as the \\ell_1 norm on singular values. Thus, using the trace norm  \\|\\boldsymbol{\\Theta}\\|_{tr} tends to produce a sparce solution on singular values, implying that \\boldsymbol{\\Theta} becomes a low-rank matrix. Note that if \\boldsymbol{\\Theta} is squared and diagonal,  \\|\\boldsymbol{\\boldsymbol{\\Theta}}\\|_{tr} is reduced to the \\ell_1 norm of the diagonal entries. Thus, the trace norm can be regarded as an extension of the \\ell_1 norm from vectors to matrices.
—— Masashi Sugiyama. 2015. Introduction to Statistical Machine Learning. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA. Chapter 24.4. p. 278.

可以证明,迹范数的对偶是谱范数(spectral norm),对应 |\\sigma_{max}| ;为了不与 \\ell_2 混淆,我们这里将谱范数记作 \\|\\bm X\\|_{sp} ,则有 \\bm S_k\\in-t\\partial\\|\
abla f(\\bm X_k)\\|_{sp} 。对于 \\partial\\|\\bm X\\|_{sp} ,不失一般性,设 \\sigma_{max}=\\sigma_1 ,并设 \\bm u_1,\\bm v_1 是对应的左右奇异向量\\|\\bm u_1\\|=\\|\\bm v_1\\|=1 ,有

\\begin{align}\\|\\bm X\\|_{sp}=\\sigma_{1}=\\bm u_1^\	op\\bm X\\bm v_1=\\bm X^\	op\\bm u_1\\bm v_1^\	op=\\left\\langle \\bm u_1\\bm v_1^\	op,\\bm X\\right\\rangle \\end{align}

\\bm u_1\\bm v_1^\	op\\in\\partial\\|\\bm X\\|_{sp} 。那么 FW 的更新方向即为 \\bm S_k=-t\\bm u_1\\bm v_1^\	op

相比于需要迹范数上的投影操作,这一操作只需要求最大奇异值对应奇异向量,有简单高效方法求解,例如

反而,投影梯度法求解这个问题,需要对矩阵做全奇异值分解,计算量比 FW 复杂很多。


值得说明的是,我们求解如下约束优化问题时,

\\begin{align}&\\min_{\\bm x}f(\\bm x)\\\\ \	ext{s.t.}\\;& \\|\\bm x\\|\\le t \\end{align} (约束形式)

可能会想到使用拉格朗日对偶来解决:

\\begin{align}&\\min_{\\bm x}f(\\bm x)+\\lambda\\|\\bm x\\| \\end{align} (拉格朗日形式)

如果我们允许 t\\lambda 的值在 [0,\\infty] 变动的话,这两个问题的解通常是相同的(也就是说可以根据某个 t 找到对应的 \\lambda ,反之亦然;尽管如此, \\lambda 也只是作为一个超参数)。在这种情况下,约束形式或拉格朗日形式并无哪个更好,对于给定的统计学习/机器学习问题,只需要选择一个更容易求解、或者效果更好的形式即可。例如,求解拉格朗日形式可以用近端方法(proximal method)。

但必须说明的是,这两种形式并不完全等同;拉格朗日形式只是一种“松弛”的形式。严格来说,约束形式等价于:

\\begin{align}\\min_{\\bm x}f(\\bm x)+ I_{\\{\\bm x|\\,\\|\\bm x\\|\\le t\\}}(\\bm x) \\end{align}

其中 I_{C} 代表集合 C示性函数(characteristic function),  I_{C}(\\bm x)=\\left\\{ \\begin{array}{cl}0,& \\bm x\\in C\\\\ +\\infty,&\\bm x\
ot\\in C \\end{array}\\right.

另外,许多其它的正则项也都有高效的 Frank--Wolfe 更新方式,例如,特殊的多面体(polyhedra) /锥(cone)限制,组合范数正则化,原子范数等。参见:

回忆一下我们所要求解的问题模型\\begin{align}\\min_{\\bm x\\in Q}f(\\bm x) \\end{align},其中 Q 是一个紧的凸集, f\\in C_{L}^{1,1}(Q) ;该问题类的 Orcale 是一阶梯度 Orcale 与 LMO。在正式开始收敛性分析之前,我们先给出两个重要的定义以及一个引理。

稳定点的定义:我们称点 \\bm x^*\\in Q 是该问题一个稳定点(stationary point,也称驻点),如果它满足

\\begin{align}\\langle \
abla f(\\bm x^*),\\bm x-\\bm x^*\\rangle\\ge 0,\\quad\\forall\\,\\bm x\\in Q \\end{align}

(如果你对这个定义感到困惑:驻点不应该是导数等于 0 的点吗?那么请戳:math.stackexchange.com/。)

Frank-Wolfe 间隔的定义:我们称如下定义的 g_k 为 Frank-Wolfe 间隔(Frank-Wolfe gap),

\\begin{align}g_k=\\langle\
abla f(\\bm x_k),\\bm x_k-\\bm s_k\\rangle \\end{align}

注意到根据 \\begin{align}\\bm s_k=\\operatorname*{argmin}_{\\bm s\\in Q}\\langle \
abla f(\\bm x_k),\\bm s\\rangle \\end{align} 的定义,有 \\langle \
abla f(\\bm x_k),\\bm s_k\\rangle\\le \\langle \
abla f(\\bm x_k),\\bm x_k\\rangle ,因此  g_k \\ge 0 。如果 f 是凸函数,那么

\\begin{align}g_k&=\\max_{\\bm s}\\langle\
abla f(\\bm x_k),\\bm x_k-\\bm s\\rangle\\\\ &\\ge \\langle\
abla f(\\bm x_k),\\bm x_k-\\bm x^*\\rangle\\\\ &\\ge f(\\bm x_k)-f(\\bm x^*)\\qquad\\cdots\\cdots(*) \\end{align}

由此可知,Frank--Wolfe 间隔可以被当作一个衡量输出的次优性(suboptimality)指标。有的文献中也称其对偶间隔(duality gap)。

Jaggi, M. (2013). Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization. International Conference on Machine Learning.

下面的定理建立了两个连续迭代点之间目标函数的联系,在类似的证明凸/非凸问题的收敛性中都有重要的作用。

定理 6:设 \\{\\bm x_0,\\bm x_1,\\dots\\} 表示 FW 算法产生的迭代点,对于任意的 \\xi\\in[0,1] ,均有下面的不等式成立

\\begin{align}f(\\bm x_{k+1})\\le f(\\bm x_k)-\\xi g_k+\\frac{1}{2}\\xi^2L\\,\\mathrm{diam}(Q)^2 \\end{align}

其中 \\mathrm{diam}(Q) 代表集合 Q直径(diameter),代表集合内距离最远的两个点的距离(的上界),定义为 \\begin{align}\\mathrm{diam}(Q)=\\sup_{\\bm x,\\bm y\\in Q}\\|\\bm x-\\bm y\\| \\end{align}

证明:根据 f 的 Lipschitz 光滑假设,对于 \\forall \\bm x,\\bm y\\in Q ,有

\\begin{align}f(\\bm y)\\le f(\\bm x)+\\langle \
abla f(\\bm x),\\bm y-\\bm x\\rangle+\\frac{L}{2}\\|\\bm y-\\bm x\\|^2 \\end{align}

代入 \\bm x=\\bm x_k\\bm y=\\bm x_{k+1}=(1-\\gamma)\\bm x_k+\\gamma\\bm s_k ,其中 \\gamma\\in[0,1] 为任意选定的步长,有

\\begin{align}f\\big((1-\\gamma)\\bm x_k+\\gamma\\bm s_k\\big)\\le f(\\bm x_k)+\\gamma\緻brace{\\langle \
abla f(\\bm x_k),\\bm s_k-\\bm x_k\\rangle}_{-g_k}+\緻brace{\\frac{L\\gamma^2}{2}\\|\\bm s_k-\\bm x_k\\|^2}_{\\le \\frac{1}{2}\\gamma^2L\\,\\mathrm{diam}(Q)^2}\\end{align}

这就证明了定理。

进一步,我们调整 \\gamma \\in[0, 1] ,让不等式右边最小化。可以看到这是一个关于 \\gamma 的二次函数,其最优点 \\gamma_k^*

\\begin{align}\\gamma_k^*=\\min\\left\\{\\frac{g_k}{L\\|\\bm x_k-\\bm s_k\\|^2},1\\right\\}\\end{align}

步长选择的策略 2 就来自于此。该步长与 FW 间隔 g_k 成比例,而 g_k 是问题次优性的度量,因此当我们接近最优值时,这个步长自然会变为零,这是一个理想的特性。同时,它的效率也很高。Demyanov 和 Rubinov 在 20 世纪 60 年代首次发表了这一策略。

Alexander Rubinov(1940–2006)(左)和 Vladimir Demyanov(1938–2014)(右),苏联优化领域的两位先驱。图片来自 https://fa.bianp.net/blog/2022/adaptive_fw/

至于线搜索策略,虽然可以得到每次迭代中能够下降最多的最优步长,然而每次迭代中的线搜索可能是一个代价高昂的优化问题;除非线搜索易于求解,否则这种方法并不适用。

我们记 C=L\\mathrm{diam}(Q)^2 ,也有的文献称其为曲率常数(curvature constant)。根据上面定理 6 的证明过程,可以看到其含义是

\\begin{align}C=\\sup\\left\\{\\left.\\frac{2}{\\gamma^2}\\big(f(\\bm y)-f(\\bm x)-\\langle\
abla f(\\bm x),\\bm y-\\bm x\\rangle\\big)\\right|\\,\\gamma\\in[0,1],\\bm x,\\bm s\\in Q,\\bm y=(1-\\gamma)\\bm x+\\gamma\\bm s\\right\\}\\end{align}

其中的 f(\\bm y)-f(\\bm x)-\\langle\
abla f(\\bm x),\\bm y-\\bm x\\rangle 正是布雷格曼散度(Bregman divergence)。

下面我们给出 FW 算法的第一个收敛性结论。注意这里的目标函数只需要满足 Lipschitz 光滑性,而不需要凸性。这一结论最先是由 Simon Lacoste--Julien 证明的。读者可参考:Simon Lacoste--Julien. (2016). Convergence rate of frank--wolfe for non--convex objectives. arXiv:1607.00345.

定理 7:若 f 可微且梯度满足 -Lipschitz 条件,那么如下的最优 Frank Wolfe 上界成立:

\\begin{align}\\min_{0\\le i\\le k}\\bm g_i\\le \\frac{\\max\\{2h_0,C\\}}{\\sqrt{k+1}}\\end{align}

其中 \\begin{align}h_0=f(\\bm x_0)-\\min_{\\bm x\\in Q}f(\\bm x)\\end{align} 表示初始点的全局次优性。

证明:根据定理 6,我们有如下不等式对于 \\forall\\,\\xi\\in[0,1] 成立

\\begin{align}f(\\bm x_{k+1}) \\leq f(\\bm x_k)-\\xi g_k+\\frac{\\xi^2C}2 \\end{align}

对等值右侧求最小化可知 \\begin{align}\\xi^*=\\min\\left\\{\\frac{g_k}{ C},1\\right\\}\\end{align} ;根据 \\xi^* 不同取值情况,我们有如下两种情况:

综合上述两种情况可得

\\begin{align}f(\\bm x_{k+1})\\leq f(\\bm x_k)-\\frac{g_k}2\\min\\left\\{\\frac{g_k}C,1\\right\\}\\end{align}

将上述不等式从 0 累加至 t ,调整不等式左右两边相关项后可得

\\begin{align*}-h_{0}& \\leq f(\\bm x_{k+1})-f(\\bm x_0)\\leq-\\sum_{i=0}^k\\frac{g_i}2\\min\\left\\{\\frac{g_i}C,1\\right\\}\\leq-(k+1)\\frac{g_k^*}2\\min\\left\\{\\frac{g_k^*}C,1\\right\\}\\end{align*}

其中 \\begin{align}g_k^*=\\min_{0\\le i\\le t}g_i\\end{align} 。同样,我们根据 g_k^* 的不同值分两种情况讨论。

因此,在这两种情况下,均有 \\begin{align}g_k^*\\leq\\frac{\\max\\{2h_0,C\\}}{\\sqrt{k+1}}\\end{align} 成立,这就证明了定理。


如果目标函数是凸的,我们可以得到一个更加紧的上界。

定理 8:若 f 为可微凸函数且梯度满足 Lipschitz 条件,那么 FW 算法满足:

\\begin{align}f(\\bm x_k)-f(\\bm x^*)\\leq\\frac{2C}{k+1}\\end{align}

该定理说明这种情况下 FW 算法有 O\\left(\\frac{1}{\\varepsilon}\\right) 的收敛速率。

证明:

e_k=A_k(f(\\bm x_k)-f(\\bm x^*)) ,其中 A_k 表示一个正的常数。那么,对 \\forall\\,\\xi\\in[0,1] ,有如下不等式成立

\\begin{align}e_{k+1}-e_{k}&=A_{k+1}(f(\\bm x_{k+1})-f(\\bm x^*))-A_k(f(\\bm x_k)-f(\\bm x^*))  \\\\ &\\leq A_{k+1}(f(\\bm x_k)-\\xi g_k+\\frac{\\xi^2C}2-f(\\bm x^*)) -A_k(f(\\bm x_k)-f(\\bm x^*)) \\\\ &\\leq A_{k+1}\\left(f(\\bm x_k)-f(\\bm x^*)-\\xi\\left(f(\\bm x_k)-f(\\bm x^*)\\right)+\\frac{\\xi^2C}2\\right) -A_k(f(\\bm x_k)-f(\\bm x^*)) \\\\ &=((1-\\xi)A_{k+1}-A_k)\\left(f(\\bm x_k)-f(\\bm x^*)\\right)+A_{k+1}\\frac{\\xi^2C}2\\qquad\\cdots\\cdots(**) \\end{align}

其中的第一个不等式来自定理 6,第二个不等式来自前文 (*) 式。接下来令 \\begin{align}A_k=\\frac{k(k+1)}{2}\\end{align}\\begin{align}\\xi=\\frac{2}{k+2}\\end{align} ,有

\\begin{align}(1-\\xi) A_{k+1}-A_k&=\\frac{k(k+1)}2-\\frac{k(k+1)}2=0\\\\A_{k+1}\\frac{\\xi^2}2&=\\frac{k+1}{k+2}\\leq1 \\end{align}

将这两个结果带入 (**) 式,有

\\begin{align}e_{k+1}-e_k\\le C \\end{align}

将该不等式从 0 累加到 k-1 ,结合 e_0=0 ,则对 \\forall\\,k>0

\\begin{align}e_k\\le kC\\enspace\\Rightarrow\\enspace f(\\bm x_k)-f(\\bm x^*)\\le \\frac{2C}{k+1}\\end{align}


最后,笔者这里给出两个来自他人的例子,这两个都是利用 FW 算法求解 lasso 问题,我们不再进行介绍:

笔者写这篇文章,也很大程度参考了上面两份材料。尤其 Fabian Pedregosa 的三篇博客对 FW 算法的讲解十分深入,这里一并贴出

关于 FW 更详尽的复杂度分析、以及改进的 FW 方法分析,读者可以参考 王奇超、文再文、蓝光辉、袁亚湘. (2020). 优化算法的复杂度分析. 中国科学:数学, 50(9), 66. faculty.bicmr.pku.edu.cn 的第四章。更多的相关文献亦可参考之。

| 首页 | 关于顺盈娱乐 | 顺盈新闻 | 顺盈注册 | 顺盈登录 | 顺盈平台 | 顺盈代理 | 顺盈APP下载 |

ICP备案:粤IP******** Copyright © 2002-2022 顺盈平台官方指定注册站 版权所有

平台注册入口