外梯度法在凸优化的最优二阶以及高阶加速
Published:
Paper reading: Optimal and Adaptive Monteiro-Svaiter Acceleration.
本文介绍最优的二阶乃至于高阶算法. 考虑凸优化问题$\min_{x \in R^d} f(x)$ ，该算法可以在 $\mathcal{O}(\epsilon^{-2/(3p+1)})$ 迭代内寻找到一个近似解 $\hat x$ 满足 $f(\hat x) - f^\ast \le \epsilon$, 其中每次迭代调用常数次 $p$ 阶Oracle $(\nabla f(x), \nabla^2 f(x), \cdots, \nabla^p f(x))$.
上述结果最早由 [1,2] 用不同的方法独立做出来，但原始的算法都稍显复杂。本文我们介绍 [1] 的算法，而且部分推导过程参照了后续工作 [3] 的推导。
Preliminaries
在本文中，我们假设目标函数满足
\[\begin{align*} \Vert \nabla^p f(x) -\nabla^p f(y) \Vert \le L \Vert x-y \Vert, \quad \forall x,y \in R^d. \end{align*}\]一个简单的推论是，定义 $p$ 阶Taylor展开
\[\begin{align*} T_{x,p}(z) = f(x) + \sum_{i=1}^p \frac{1}{i!} \nabla^i f(x) [z-x]^i, \quad \forall z \in R^d. \end{align*}\]那么我们有如下的逼近关系
\[\begin{align*} &\vert f(z) - T_{x,p}(z) \vert \le \frac{L}{(p+1)!} \Vert z - x \Vert^{p+1}; \\ &\Vert \nabla f(z) - \nabla T_{x,p}(z) \Vert \le \frac{L}{p!} \Vert z - x \Vert^{p}. \end{align*}\]Basic Tensor Step
本节，介绍基础的张量步，作为后续的oracle。
根据上一节介绍的Taylor展开的逼近关系，容易想到如下的操作，对于一个点 $x$, 利用Taylor展开定义一个局部的近似函数，然后求极小化，得到下一个点。用数学语言表达，我们有
\[\begin{align*} x^+ = \arg \min_{z \in R^d} T_{x,p}(z) + \frac{M}{(p+1)!} \Vert z - x \Vert^{p+1}, \end{align*}\]其中参数 $M \ge L$. 例如 $M = 2L$, 令 $\lambda = (M / p!) \Vert x^+ -x \Vert^{p-1}$, 可以验证上述迭代满足
\[\begin{align*} \Vert \nabla f(x^+) + \lambda (x^+ - x) \Vert=& \Vert \nabla f(x^+) - \nabla T_{x,p}(x^+) \Vert \le \frac{L}{p!} \Vert x^+ - x \Vert^p = \frac{\lambda}{2} \Vert x^+ - x \Vert. \end{align*}\]在后续算法的分析中，我们会用到这个条件。在这个章节，我们讨论上述张量步的具体实现。
- 对于 $p=1$ 的情况，可以简单发现这等价于步长为 $1/M$ 的梯度步。
- 对于 $p=2$ 的情况，可以发现这等价于之前介绍过的 三次正则牛顿步. 根据Nesterov 和 Polyak 著名的文章 [4], 我们知道这个二阶oracl可以在几乎和矩阵求逆等价的时间内实现。
- 对于 $p=3$ 的情况，这等价于实现一个三阶张量步。根据Nesterov的近期工作 [5], 根据我们在节末给出的矩阵不等式，事实上三阶方向导数 $\nabla^3 f(x)[h^3]$ 可以被二阶导数 $\nabla^2 f(x)[h]$ 以及四阶距离函数 $\Vert h \Vert^4$ 的和控制住。最后三阶导数对子问题的影响很小，仍然可以在几乎和矩阵求逆等价的时间内实现。这个矩阵不等式为，对于任意的 $\tau >0$, 成立
- 对于 $p=+\infty$ 的情况， 文章 [1] 同样也讨论了对于Hessian很稳定的函数，例如logistic function，也可以通过牛顿步进行实现。
Wamup: Basic Acceleration
本节，作为一个预热，我们给出一个相对简单的加速算法，可以取得 $\mathcal{O}(\epsilon^{-(p+1)})$ 的收敛率。这个方法完全可以由我们之前的内容得到。在外梯度法的最优一阶加速 中，我们介绍了使用anytime-online-to-batch conversion以及外梯度法可以得到一个最优的具有 $\mathcal{O}(\epsilon^{-1/2})$ 的一阶算法。下面我们将该算法推广到高阶。
给定算法访问的点 $w_1,\cdots, w_t$, 我们会定义一个加权平均 $x_t = \sum_{i=0}^{t} a_i w_i / A_t$, 其中 $A_t = \sum_{i=0}^{t} a_t$. 我们在之前文章中已经证明:
\[\begin{align*} f(x_T) - f(x^\ast) \le \frac{\langle a_t \nabla f(x_t), w_t - x^\ast \rangle}{A_T} := \frac{R_T(x^\ast) }{A_T}. \end{align*}\]我们下面设计一个算法保证遗憾界 $R_T(x^\ast) = \mathcal{O}(1)$。利用之前所介绍的外梯度法，为了达到该遗憾界，只需要构造出一个hint满足条件
\[\begin{align*} a_t \eta \Vert h_t - g_t \Vert \le \Vert w_t - v_t \Vert /2, \quad {\rm where} \quad g_t = \nabla f(x_t), \end{align*}\]其中 $v_t$ 为外梯度法中的对偶变量。注意到 $x_t = \sum_{i=0}^{t} a_i w_i / A_t$ 依赖于还没有访问的 $w_t$，我们考虑将其用 $v_t$ 来代替，得到如下的hint设计：
\[\begin{align*} z_t = \frac{\sum_{i=0}^{t-1} a_i w_i + a_t v_t}{A_t}, \quad h_t = \sum_{i=1}^{p-1} \frac{1}{i!} \nabla^i f(x_t)[z_t-x_t]^{i-1}. \end{align*}\]那么我们知道hint的近似误差满足
\[\begin{align*} \Vert h_t - g_t \Vert \le \frac{L}{p!} \Vert z_t - x_t \Vert^p = \frac{L}{p!} \left( \frac{a_t}{A_t} \right)^{p} \Vert w_t -v_t \Vert^p. \end{align*}\]通过将外梯度发的所有迭代都投影到一个固定大小的集合内，我们可以容易地保证 $\Vert w_t -v_t \Vert = \mathcal{O}(D)$, 那么就有
\[\begin{align*} \frac{a_t \eta \Vert h_t - g_t \Vert}{\Vert w_t - v_t \Vert} = \mathcal{O} \left( \frac{\eta L a_t^{p+1} D^{p-1}}{A_t^p} \right). \end{align*}\]选择步长 $\eta = \mathcal{O}(1/ (L D^{p-1}))$, 那么可以发现最大可选择的级数为 $a_t = \Omega(t^p)$。在这个设定下满足 $A_t = \Omega(t^{p+1})$. 然后，根据 外梯度法 的推导得到算法具有如下的收敛率:
\[\begin{align*} f(x_T) - f(x^\ast) \le \frac{L D^{p+1}}{T^{p+1}}. \end{align*}\]换句话说，找到一个 $\epsilon$-解的 $p$ 阶oracle复杂度为 $\mathcal{O}(\epsilon^{-(p+1)})$.
Optimal Acceleration
本节介绍 [1] 中所给出的算法。首先，初始化变量 $A_0= 0$, $v_0 = y_0 = x_0$. 然后进行如下迭代总共 $T$ 次。
\[\begin{align*} a'_{t+1} =& \frac{1 + \sqrt{1+4 \lambda'_{t+1} A_t}}{2 \lambda'_{t+1}}, \quad A'_{t+1} = A_t + a'_{t+1}, \quad y_t = \frac{A_t}{A'_{t+1}} x_t + \frac{a'_{t+1}}{A'_{t+1}} v_t \\ z_{t+1} =& \arg \min_{z \in R^d} T_{y_t,p}(z) + \frac{2L}{(p+1)!} \Vert z - y_t \Vert^{p+1}, \quad \lambda_{t+1} = 2L \Vert z_{t+1} - y_t \Vert^{p-1} \\ \gamma_{t+1} =& \min \left\{1, \frac{\lambda'_{t+1}}{\lambda_{t+1}} \right\}, \quad a_{t+1} = \gamma_{t+1} a'_{t+1}, \quad A_{t+1} = A_t + a_{t+1} \\ x_{t+1} =& \frac{(1-\gamma_{t+1}) A_t}{A_{t+1}} x_t + \frac{\gamma_{t+1} A'_{t+1}}{A_{t+1}} z_{t+1} \\ v_{t+1} = & v_t - a_{t+1} \nabla f(z_{t+1}). \end{align*}\]算法看起来有点复杂。我们进行一些简单的解释。框架上，算法大体上还是与 外梯度法 相似，用张量步更新变量原变量 $z$ 后，再紧接着更新对偶变量 $v$. 算法中的 $A_t = \sum_{i=0}^t a_t$ 也与之前加速算法中的量相同，我们希望证明算法的收敛率为 $\mathcal{O}(1/A_T)$ 并且希望 $A_t$ 是一个增长尽可能快的级数。$a’_{t+1}$ 的公式看着有点古怪，事实上它是一个二次函数的根，来源于 Nesterov加速梯度法 类似的步骤，我们很快会在后续的证明中看到这个公式的作用。算法中的 $\lambda’_{t+1}$ 会是一个动态选择的期望步长大小的倒数，而 $\lambda_{t+1}$ 则代表实际张量步得到的步长大小的倒数。当实际量 $\lambda_{t+1} \ge \lambda’_{t+1}$ 的时候，我们需要用系数 $\gamma_{t+1} = \lambda’_{t+1} / \lambda_{t+1}$ 减缓对偶变量的更新幅度。
Analysis
尽管算法稍显复杂，但结合证明可能更好理解。我们将证明拆解为以下几个步骤。
Potential Decrease
这一步中，我们定义势能函数并且证明其下降性。为了叙述方便，我们定义如下符号
\[\begin{align*} E_t: =f(x_t) - f(x^\ast), \quad D_t := \frac{1}{2} \Vert v_t - x^\ast \Vert^2, \quad N_{t+1}:= \frac{1}{2} \Vert z_{t+1} - y_t \Vert^2. \end{align*}\]与Nesterov加速梯度法类似，我们希望证明势能函数 $A_t E_t + D_t$ 是单调下降的。
根据 $A’_{t+1} = A_t + a’_{t+1}$. 我们有
\[\begin{align*} a'_{t+1} v_t = A'_{t+1} y_t - A_t x_t = a'_{t+1} z_{t+1} + A'_{t+1} (y_t - z_{t+1})- A_t (x_t - z_{t+1}). \end{align*}\]利用这个恒等式，我们将内积项写成如下的形式
\[\begin{align*} & a'_{t+1} \langle \nabla f(z_{t+1}), x^\ast - v_t \rangle \\ =& \langle \nabla f(z_{t+1}), A'_{t+1} y_t - A_t x_t = a'_{t+1} z_{t+1} + A'_{t+1} (y_t - z_{t+1})- A_t (x_t - z_{t+1}) \rangle \\ \le& a'_{t+1} (f(x^\ast) - f(z_{t+1})) + A'_{t+1} \langle \nabla f(z_{t+1}), z_{t+1} - y_t \rangle + A_t ( f(x_t) - f(z_{t+1}) ) \\ =& A_t E_t - A'_{t+1} ( f(z_{t+1}) - f(x^\ast) ) + A'_{t+1} \langle \nabla f(z_{t+1}), z_{t+1} - y_t \rangle \end{align*}\]其中不等式使用了目标函数 $f$ 的凸性。对于最后的内积项，使用恒等式 $\langle a,b\rangle=\frac{1}{2} \Vert a \Vert^2 - \frac{1}{2} \Vert a \Vert^2 - \frac{1}{2} \Vert b \Vert^2$ 可以进一步得到
\[\begin{align*} & A'_{t+1} \langle \nabla f(z_{t+1}), z_{t+1} - y_t \rangle = \frac{A'_{t+1}}{\lambda_{t+1}} \langle \lambda_{t+1} \nabla f(z_{t+1}), z_{t+1} - y_t \rangle \\ =& \frac{A'_{t+1}}{2 \lambda_{t+1}} \Vert \nabla f(z_{t+1}) + \lambda_{t+1} (z_{t+1} - y_t) \Vert^2 - \frac{A'_{t+1}}{2 \lambda_{t+1}} \Vert \nabla f(z_{t+1}) \Vert^2 - \frac{A'_{t+1}\lambda_{t+1}}{2} \Vert z_{t+1} - y_t \Vert^2 \\ \le&-\frac{3}{4} A'_{t+1}\lambda_{t+1} N_{t+1} - \frac{A'_{t+1}}{2 \lambda_{t+1}} \Vert \nabla f(z_{t+1}) \Vert^2, \end{align*}\]其中最后一步利用了之前所推到过的张量步的性质。因此，
\[\begin{align*} &A'_{t+1} (f(z_{t+1}) - f(x^\ast)) - A_t E_t \\ \le& a'_{t+1} \langle \nabla f(z_{t+1}), v_t - x^\ast \rangle -\frac{3}{4} A'_{t+1}\lambda_{t+1} N_{t+1} - \frac{A'_{t+1}}{2 \lambda_{t+1}} \Vert \nabla f(z_{t+1}) \Vert^2. \end{align*}\]利用恒等式 $A_{t+1} = A_t + a_{t+1}= A_t + \gamma_{t+1} a’_{t+1} = (1-\gamma_{t+1}) A_t +\gamma_{t+1} A’_{t+1}$ 以及函数的凸性，我们知道
\[\begin{align*} &A_{t+1} E_{t+1} \\ \le& (1-\gamma_{t+1}) A_t E_t + \gamma_{t+1} A'_{t+1} (f(z_{t+1} - f(x^\ast)) ) \\ \le& A_t E_t + a_{t+1} \langle \nabla f(z_{t+1}), v_t - x^\ast \rangle -\frac{3}{4} \gamma_{t+1} A'_{t+1} \lambda_{t+1} N_{t+1} - \frac{\gamma_{t+1} A'_{t+1}}{2 \lambda_{t+1}} \Vert \nabla f(z_{t+1}) \Vert^2. \end{align*}\]然后根据 $D_t$ 的定义以及对偶变量 $v_{t+1}$ 的外梯度更新公式，我们有
\[\begin{align*} D_{t+1}=& \frac{1}{2} \Vert v_{t+1} - x^\ast \Vert^2 = \frac{1}{2} \Vert v_t - a_{t+1} \nabla f(z_{t+1}) - x^\ast \Vert^2 \\ =& D_t + a_{t+1} \langle \nabla f(z_{t+1}), x^\ast - v_t \rangle + \frac{a_{t+1}^2}{2} \Vert \nabla f(z_{t+1}) \Vert^2. \end{align*}\]这里出现了相同（但是符号相反的）内积项，我们将两个式子相加后得到
\[\begin{align*} & (A_{t+1} E_{t+1} + D_{t+1}) - (A_t E_t + D_t) \\ \le& -\frac{3}{4} \gamma_{t+1} A'_{t+1} \lambda_{t+1} N_{t+1} - \left(\frac{\gamma_{t+1} A'_{t+1}}{2 \lambda_{t+1}} - \frac{a_{t+1}^2}{2} \right) \Vert \nabla f(z_{t+1}) \Vert^2. \end{align*}\]我们希望上式的最后一个一直是负数因此可以丢掉。这可以通过算法中定义 $a’_{t+1}$ 满足二次方程 $A’_{t+1} = A_t +a’_{t+1} = \lambda’_{t+1} (a’_{t+1})^2$ 来得到，这是因为在这个设定下，我们有
\[\begin{align*} \frac{\gamma_{t+1} A'_{t+1}}{2 \lambda_{t+1}} = \frac{\gamma_{t+1} \lambda'_{t+1} (a'_{t+1})^2}{2 \lambda_{t+1}} = \frac{\lambda'_{t+1} a_{t+1}^2}{2 \gamma_{t+1} \lambda_{t+1}} \ge \frac{a_{t+1}^2}{2}. \end{align*}\]因此，我们可以得到
\[\begin{align*} A_{t+1} E_{t+1} + D_{t+1}\le A_t E_t + D_t -\frac{3}{4} \gamma_{t+1} A'_{t+1} \lambda_{t+1} N_{t+1}. \end{align*}\]Estimation of the Growth Rate
下面，我们希望估计 $A_t$ 的增长速率。我们采用 [3] 中给出的一个简化版本的证明。
注意到算法中还有一个参数 $\lambda’_{t+1}$ 还没有选择，我们将根据上一小节得到的下降引理，选择最优的 $\lambda’_{t+1}$ 使得收敛最快。 在下面的证明中，我们分析会人为地将算法的运行过程分为多个epoch，其中每个epoch使得 $A_t$ 增长为2倍。我们分析每一个epoch应当有多长。我们记这样的epoch的开始时间为 $T_1$, 结束时间为 $T_2$, 那么 $T_2$ 也就是使得 $A_{t} \ge 2A_{T_1}$ 的最小时间戳 $t$。在每一个epoch内，我们会选择一个最优的，仅依赖于 $T_1$ 的参数 $\lambda’_{t+1}$.
在使用下降引理的时候，我们分别考虑 $\lambda_{t+1} \le \lambda’_{t+1}$ 以及 $\lambda_{t+1} > \lambda’_{t+1}$ 两种情况。它们分别对应着 $\gamma_{t+1} = 1$ 以及 $\gamma_{t+1}<1$ 的两种情况。我们将区间 $[T_1,T_2]$ 划分为这两种情况的不相交的并集，记作为 $S_{\rm down} \cup S_{\rm up}= [T_1,T_2]$.
在第一种情况下，我们有
\[\begin{align*} \sqrt{A_{t+1}} - \sqrt{A_t} = \frac{A_{t+1} - A_t}{\sqrt{A_{t+1}} + \sqrt{A_t}} = \frac{a_{t+1}}{\sqrt{A_{t+1}} + \sqrt{A_t}} = \frac{\sqrt{A_{t+1}/\lambda'_{t+1}}}{\sqrt{A_{t+1}} + \sqrt{A_t}} \ge \frac{1}{2 \sqrt{\lambda'_{t+1}}}. \end{align*}\]因此，我们可以推断出
\[\begin{align*} \vert S_{\rm down} \vert \le 2 \left( \sqrt{2} -1\right) \sqrt{A_{T_1} \lambda'_{t+1}}. \end{align*}\]在第二种情况下，根据关系式 $\lambda_{t+1} = 2 L \Vert z_{t+1} - y_t \Vert^{p-1} = 2L(2N_{t+1})^{(p-1)/2}$, 我们有
\[\begin{align*} D_0 \ge \frac{3}{4} \gamma_{t+1} A'_{t+1} \lambda_{t+1} N_{t+1} = \frac{3}{8} A'_{t+1} \lambda'_{t+1} \left( \frac{\lambda_{t+1}}{2L} \right)^{\frac{2}{p-1}} \ge \frac{3}{8} A'_{t+1} \lambda'_{t+1} \left( \frac{\lambda'_{t+1}}{2L} \right)^{\frac{2}{p-1}}. \end{align*}\]因此，这种情况下，我们可以推断出，
\[\begin{align*} \vert S_{\rm up} \vert \le \frac{8 D_0 (2L)^{2/(p-1)} }{3 A_{T_1} (\lambda'_{t+1})^{(p+1)/(p-1)}}. \end{align*}\]综合两种情况，我们在这个epoch内选择最优的 $\lambda’_{t+1}$ 之后，可以使得
\[\begin{align*} T_2 - T_1 = \vert S_{\rm down} \vert + \vert S_{\rm up} \vert = \mathcal{O} \left( \left( D_0^{\frac{p-1}{2}} L A_{T_1} \right)^{\frac{2}{3p+1}} \right). \end{align*}\]对所有epoch求和后，可以得到
\[\begin{align*} T = \mathcal{O} \left( \left( D_0^{\frac{p-1}{2}} L A_{T} \right)^{\frac{2}{3p+1}} \right). \end{align*}\]移项后就得到了
\[\begin{align*} A_T = \Omega \left( \frac{T^{(3p+1)/2}}{D_0^{(p-1)/2} L } \right ). \end{align*}\]根据 $E_T \le D_0 / A_T$, 我们知道最终的收敛率满足
\[\begin{align*} E_T = \mathcal{O} \left( \frac{D_0^{(p+1)/2} L}{T^{(3p+1)/2}} \right). \end{align*}\]换句话说，想要找到一个 $\epsilon$-解的复杂度为 $\mathcal{O}(\epsilon^{-2/(3p+1)})$.
Reference
[1] Carmon Y, Hausler D, Jambulapati A, Jin Y, Sidford A. Optimal and adaptive Monteiro-Svaiter acceleration. In NeurIPS, 2022.
[2] Kovalev D, Gasnikov A. The first optimal acceleration of high-order methods in smooth convex optimization. In NeurIPS, 2022.
[3] Adil D, Bullins B, Jambulapati A, Sidford A. Convex optimization with $p$-norm oracles. arXiv preprint arXiv:2410.24158. 2024.
[4] Nesterov Y, Polyak BT. Cubic regularization of Newton method and its global performance. Mathematical programming. 2006.
[5] Nesterov Y. Implementable tensor methods in unconstrained convex optimization. Mathematical Programming. 2021.