# 写在前面

介绍L1/L2L_1/L_2最小化的加速算法 [1],我在本科的时候就开始跟着我的导师看娄易非老师的文章(ℓ1−ℓ2\ell_1-\ell_2),现在虽然方向不一致,但是看文献总能学到新的技巧,我想以后会一直留意娄老师的最新文献的。附上娄老师主页[2]。

# 优化模型

# L0L_0最小化

min⁡x∈Rn∥x∥0s.t.Ax=b\min_{\mathbf x \in \mathbb R^n} \|\mathbf x\|_0 \quad \text{s.t.} \quad A \mathbf x = \mathbf b

NP 难问题

  • 贪婪算法(当nn很大时,缺乏求解精度)
    • 正交匹配追踪(OMP)
    • 正交最小二乘(OLS)
    • 采样匹配追踪(CoSaMp)
  • 近似算法()
    • 凸松弛:L1L_1替换L0L_0的基追踪(BP)
    • 非凸松弛:多种改进

# 比率优化问题L1L2\frac{L_1}{L_2}

min⁡x∈Rn∥x∥1∥x∥2s.t.Ax=b\min_{\mathbf x \in \mathbb R^n} \frac{\|\mathbf x\|_1}{\|\mathbf x\|_2} \quad \text{s.t.} \quad A \mathbf x = \mathbf b

  • 尺度不变性
  • 不含参数
  • 非凸的稀疏度度量
  • 使用ADMM算法求解

引入两个辅助变量,构造增广的拉格朗日函数

L(x,y,z,v,w)=∥z∥1∥y∥2+I(Ax−b)+ρ12∥x−y+1ρ1v∥22+ρ22∥x−z+1ρ2w∥22\begin{aligned} L(x,y,z,v,w) =& \frac{\|z\|_1}{\|y\|_2} + I(Ax-b) \\ &+ \frac{\rho_1}{2}\|x-y +\frac{1}{\rho_1}v\|_2^2 \\ &+ \frac{\rho_2}{2}\|x-z+\frac{1}{\rho_2}w\|_2^2 \end{aligned}

其中I(⋅)I(\cdot)为示性函数

I(t)={0,t=0+∞,otherwiseI(t) = \begin{cases} 0, & t = 0 \\ +\infty, & \text{otherwise} \end{cases}

# 混合优化问题L1−αL2L_1 - \alpha L_2

min⁡x∈Rn∥x∥1−α∥x∥2s.t.Ax=b\min_{\mathbf x \in \mathbb R^n} \|\mathbf x\|_1 - \alpha \|\mathbf x\|_2 \quad \text{s.t.} \quad A \mathbf x = \mathbf b

娄老师在这一问题上有很多成果,可以参考主页[2:1]。

# 联系

下面的性质给出比率模型L1L2\frac{L_1}{L_2}与混合模型L1−αL2L_1 - \alpha L_2之间的关系:

α∗=inf⁡x{∥x∥1∥x∥2s.t.Ax=b}T(α)=inf⁡x{∥x∥1−α∥x∥2s.t.Ax=b}\begin{aligned} \alpha^* &= \inf_{\mathbf x} \left\{\frac{\|\mathbf x\|_1}{\|\mathbf x\|_2} \quad \text{s.t.} \quad A \mathbf x = \mathbf b \right\}\\ T(\alpha) &= \inf_{\mathbf x} \left\{\|\mathbf x\|_1 - \alpha \|\mathbf x\|_2 \quad \text{s.t.} \quad A \mathbf x = \mathbf b \right\} \end{aligned}

则

  • 若T(α)<0T(\alpha) < 0,则α>α∗\alpha > \alpha^*。
  • 若T(α)≥0T(\alpha) \geq 0,则α≤α∗\alpha \leq \alpha^*。
  • 若T(α)=0T(\alpha) = 0,则α=α∗\alpha = \alpha^*。

这表明比率模型的最优值是T(α)T(\alpha)的根,因此可通过分半搜索(bisection search)来求得。

由向量范数的等价关系

∥x∥2≤∥x∥1≤n∥x∥2,∀x∈Rn\|\mathbf x\|_2 \leq \|\mathbf x\|_1 \leq \sqrt{n}\|\mathbf x\|_2,\quad\forall \mathbf x \in \mathbb R^n

可以设置初始参数α(0)∈[1,n]\alpha^{(0)} \in [1,\sqrt n],并根据T(α(0))T(\alpha^{(0)})进行分半搜索

  • 若T(α)=0T(\alpha) = 0,则混合模型的最优值是比率模型的最优值。

  • 若T(α)<0T(\alpha) < 0,则更新α\alpha的范围[α(0),n][\alpha^{(0)},\sqrt n]。

  • 若T(α)≥0T(\alpha) \geq 0,则更新α\alpha的范围[1,α(0)][1,\alpha^{(0)}]。因为在下一次迭代时,混合模型L1−∥x(k+1)∥1∥x(k+1)∥2L2L_1 - \frac{\|x^{(k+1)}\|_1}{\|x^{(k+1)}\|_2}L_2的目标函数不超过00,可进一步考虑区间[1,∥x(k+1)∥1∥x(k+1)∥2][1,\frac{\|x^{(k+1)}\|_1}{\|x^{(k+1)}\|_2}]。

  • 更新范围后,选择区间断点的均值作为参数α(1)\alpha^{(1)},并迭代下去。

# 凸函数之差算法

对于非凸函数f(x)f(x),若存在两个凸函数g(x),h(x)g(x),h(x)满足f(x)=g(x)−h(x)f(x) = g(x) - h(x),那么如下优化问题

min⁡xf(x)\min_x f(x)

可通过对h(x)h(x)线性化,得到迭代DCA序列

x(k+1)=arg⁡min⁡xg(x)−⟨x,∇h(x(k)⟩x^{(k+1)} = \arg\min_x g(x) - \left \langle x, \nabla h(x^{(k)} \right \rangle

# 求解算法

由于L1−αL2L_1 - \alpha L_2模型的目标函数可化成两个凸函数之差

g(x)=∥x∥1+I(Ax−b),h(x)=α∥x∥2g(x) = \|x\|_1 + I(Ax-b), h(x) = \alpha\|x\|_2

因此可通过如下DCA迭代序列求解L1−αL2L_1 - \alpha L_2模型

x(k+1)=arg⁡min⁡xg(x)−⟨x,αx(k)∥x(k)∥2⟩x^{(k+1)} = \arg\min_x g(x) - \left \langle x, \frac{\alpha x^{(k)}}{\|x^{(k)}\|_2} \right \rangle

# 分半搜索算法流程

解释:

  • line 2:lb,ublb,ub为参数α\alpha的区间范围
  • line 4:通过DCA序列求解混合模型
  • line 11:通过区间的均值作为参数α\alpha的更新
  • 整个算法是基于比率模型和混合模型的解关系来设计的分半搜索

# 参数自适应迭代

用当前解的比率作为参数α\alpha的更新,得到如下模型(A1)(A1)

{x(k+1)=arg⁡min⁡x{g(x)−⟨x,α(k)x(k)∥x(k)∥2⟩}α(k+1)=∥x(k+1)∥1/∥x(k+1)∥2\left\{\begin{array}{l} \mathbf{x}^{(k+1)}=\arg \min _{\mathbf{x}}\left\{g(\mathbf{x})-\left\langle\mathbf{x}, \frac{\alpha^{(k)} \mathbf{x}^{(k)}}{\left\|\mathbf{x}^{(k)}\right\|_{2}}\right\rangle\right\} \\ \alpha^{(k+1)}=\left\|\mathbf{x}^{(k+1)}\right\|_{1} /\left\|\mathbf{x}^{(k+1)}\right\|_{2} \end{array}\right.

求解x\mathbf{x}的子问题是一个线性规划问题,但由于问题的无界性,不能保证最优解存在。下面整合一个二次项来增加模型的鲁棒性(A2)(A2)

{x(k+1)=arg⁡min⁡x{g(x)−⟨x,α(k)x(k)∥x(k)∥2⟩+β2∥x−x(k)∥22}α(k+1)=∥x(k+1)∥1/∥x(k+1)∥2\left\{\begin{array}{l} \mathbf{x}^{(k+1)}=\arg \min _{\mathbf{x}}\left\{g(\mathbf{x})-\left\langle\mathbf{x}, \frac{\alpha^{(k)} \mathbf{x}^{(k)}}{\left\|\mathbf{x}^{(k)}\right\|_{2}}\right\rangle+\frac{\beta}{2}\left\|\mathbf{x}-\mathbf{x}^{(k)}\right\|_{2}^{2}\right\} \\ \alpha^{(k+1)}=\left\|\mathbf{x}^{(k+1)}\right\|_{1} /\left\|\mathbf{x}^{(k+1)}\right\|_{2} \end{array}\right.

  • 模型(A1)(A1)

arg⁡min⁡x∥x∥1+I(Ax−b)−⟨x,α(k)x(k)∥x(k)∥2⟩\arg \min _{\mathbf{x}} \|\mathbf{x}\|_1 + I(A\mathbf{x}-b)-\left\langle\mathbf{x}, \frac{\alpha^{(k)} \mathbf{x}^{(k)}}{\left\|\mathbf{x}^{(k)}\right\|_{2}}\right\rangle

任意实数(向量)可表示为两正数(向量)之差,可假设x=x+−x−\mathbf{x} = \mathbf{x}^+ - \mathbf{x}^-,其中x+,x−>0\mathbf{x}^+,\mathbf{x}^- >0,记xˉ=[x+;x−]T,Aˉ=[A−A]\bar{\mathbf{x}} = [\mathbf{x}^+ ;\mathbf{x}^-]^T,\bar A = [A \quad -A],则模型(A1)(A1)转化为

min⁡xˉ≥0[1+α(k)x(k)∥x(k)∥2;1−α(k)x(k)∥x(k)∥2]Txˉs.t.Aˉxˉ=b\min_{\bar{\mathbf{x}} \geq 0} \left[1+\frac{\alpha^{(k)} \mathbf{x}^{(k)}}{\left\|\mathbf{x}^{(k)}\right\|_{2}};1-\frac{\alpha^{(k)} \mathbf{x}^{(k)}}{\left\|\mathbf{x}^{(k)}\right\|_{2}}\right]^T \bar{\mathbf{x}} \quad \text{s.t.} \quad \bar A \bar{\mathbf{x}} = \mathbf{b}

该线性规划问题可通过Gurobi求解器解决。

  • 模型(A2)(A2)

arg⁡min⁡x∥x∥1+I(Ax−b)−⟨x,α(k)x(k)∥x(k)∥2⟩+β2∥x−x(k)∥22\arg \min _{\mathbf{x}}\|\mathbf{x}\|_1 + I(A\mathbf{x}-b)-\left\langle\mathbf{x}, \frac{\alpha^{(k)} \mathbf{x}^{(k)}}{\left\|\mathbf{x}^{(k)}\right\|_{2}}\right\rangle+\frac{\beta}{2}\left\|\mathbf{x}-\mathbf{x}^{(k)}\right\|_{2}^{2}

引入二次项后,模型是一个二次规划问题,可通过ADMM算法求解。引入辅助变量并构造拉格朗日函数

Lρ(x,y,u)=∥y∥1+I(Ax−b)−⟨x,α(k)x(k)∥x(k)∥2⟩+β2∥x−x(k)∥22+uT(x−y)+ρ2∥x−y∥22\begin{aligned} L_\rho (\mathbf x, \mathbf y, \mathbf u) =& \|\mathbf{y}\|_1 + I(A\mathbf{x}-b)-\left\langle\mathbf{x}, \frac{\alpha^{(k)} \mathbf{x}^{(k)}}{\left\|\mathbf{x}^{(k)}\right\|_{2}}\right\rangle\\ &+\frac{\beta}{2}\left\|\mathbf{x}-\mathbf{x}^{(k)}\right\|_{2}^{2} + \mathbf u^T(\mathbf x-\mathbf y) + \frac{\rho}{2}\left\|\mathbf{x}-\mathbf{y}\right\|_{2}^{2} \end{aligned}

对x\mathbf x的子问题可视为Ax=bAx=b约束下的最小二乘问题,通过投影算子写成显式解;对y\mathbf y的子问题是一类经典的ℓ1\ell_1范数的优化问题,可通过软阈值算子写成其显式解。这里就不一一细说了。

# 参数自适应迭代的算法流程

# 与其他工作的关联

文章算法与参数选择、广义逆幂和基于梯度的方法均产生关联,下面介绍后两种的联系

# 广义逆幂

逆幂方法是一种求解半正定对称矩阵的最小特征值方法,其通过迭代求解下面线性方程

Bx(k+1)=x(k)B\mathbf x^{(k+1)} = \mathbf x^{(k)}

记迭代收敛结果(即最小特征值对应的特征向量)为x∗\mathbf x^*,则最小特征值估计为λ=q(x∗)\lambda = q(\mathbf x^*),其中q(⋅)q(\cdot)为瑞利商

q(x)=⟨x,Bx⟩∥x∥22q(\mathbf x) = \frac{\left\langle\mathbf x,B\mathbf x\right\rangle}{\|\mathbf x\|_2^2}

迭代求解线性方程等价为如下优化问题

x(k+1)=arg⁡min⁡x12⟨x,Bx⟩−⟨x(k),x⟩\mathbf x^{(k+1)} = \arg\min_{\mathbf x}\frac{1}{2}\left\langle\mathbf x,B\mathbf x\right\rangle-\left\langle\mathbf x^{(k)},\mathbf x\right\rangle

以上思路可以推广至非线性情形,得到广义逆幂方法。对任意函数r(⋅),s(⋅)r(\cdot), s(\cdot),考虑一个广义的商函数

q(x)=r(x)s(x)q(\mathbf x) = \frac{r(\mathbf x)}{s(\mathbf x)}

以及对应优化问题和特征值的更新

{x(k+1)=arg⁡min⁡xr(x)−λ(k)⟨∇s(x(k)),x⟩λ(k+1)=r(x(k+1))/s(x(k+1))\left\{\begin{array}{l} \mathbf{x}^{(k+1)}=\arg\min_{\mathbf x}r(\mathbf x)- \lambda^{(k)}\left\langle \nabla s(\mathbf x^{(k)}),\mathbf x\right\rangle \\ \lambda^{(k+1)}=r(\mathbf{x}^{(k+1)})/s(\mathbf{x}^{(k+1)}) \end{array}\right.

显然,当r(x)=g(x),s(x)=∥x∥2,λ=αr(\mathbf{x}) = g(\mathbf{x}),s(\mathbf{x}) = \|\mathbf{x}\|_2,\lambda = \alpha时,以上广义逆幂与模型(A1)(A1)一致。通过最速下降流(steepest descent flow)整合二次项则得到模型(A2)(A2)。

# 梯度方法

回到比率优化问题L1L2\frac{L_1}{L_2}

min⁡x∈Rn∥x∥1∥x∥2s.t.Ax=b\min_{\mathbf x \in \mathbb R^n} \frac{\|\mathbf x\|_1}{\|\mathbf x\|_2} \quad \text{s.t.} \quad A \mathbf x = \mathbf b

根据KKT条件,x∗≠0\mathbf x^*\neq \mathbf 0是优化问题的临界点当且仅当存在向量s\mathbf s满足

{0∈∂∥x∗∥1∥x∗∥2−∥x∗∥1∥x∗∥22x∗∥x∗∥2+ATs0=Ax∗−b\left\{\begin{array}{l} 0 \in \frac{\partial\left\|\mathbf{x}^{*}\right\|_{1}}{\left\|\mathbf{x}^{*}\right\|_{2}}-\frac{\left\|\mathbf{x}^{*}\right\|_{1}}{\left\|\mathbf{x}^{*}\right\|_{2}^{2}} \frac{\mathbf{x}^{*}}{\left\|\mathbf{x}^{*}\right\|_{2}}+A^{T} \mathbf{s} \\ 0=A \mathbf{x}^{*}-\mathbf{b} \end{array}\right.

对向量缩放s^=∥x∗∥2s\hat{\mathbf s} = \|\mathbf{x}^{*}\|_2 \mathbf s,第一个次梯度条件为

0∈∂∥x∗∥1−∥x∗∥1∥x∗∥2x∗∥x∗∥2+ATs^0 \in \partial\left\|\mathbf{x}^{*}\right\|_{1}-\frac{\left\|\mathbf{x}^{*}\right\|_{1}}{\left\|\mathbf{x}^{*}\right\|_{2}} \frac{\mathbf{x}^{*}}{\left\|\mathbf{x}^{*}\right\|_{2}}+A^{T} \hat{\mathbf{s}}

该最优性条件往往也对应符合此梯度相同的其他优化问题

min⁡xg(x)+w(x)\min_{\mathbf x} g(\mathbf x) + w(\mathbf x)

其中w(⋅)w(\cdot)的形式无法明确确定,但仅需要其梯度信息

g(x)=∥x∥1+I(Ax−b)∇w(x)=−∥x∥1∥x∥2x∥x∥2\begin{aligned} g(\mathbf{x}) &= \|\mathbf{x}\|_1 + I(A\mathbf{x}-\mathbf{b})\\ \nabla w(\mathbf{x}) &= -\frac{\left\|\mathbf{x}\right\|_{1}}{\left\|\mathbf{x}\right\|_{2}} \frac{\mathbf{x}}{\left\|\mathbf{x}\right\|_{2}} \end{aligned}

对于该优化问题,采用近邻梯度法(proximal gradient)得到迭代

x(k+1)=proxgβ(x(k)−1β∇w(x(k)))\mathbf x^{(k+1)} = \textbf{prox}_{\frac{g}{\beta}} (\mathbf x^{(k)} - \frac{1}{\beta}\nabla w(\mathbf x^{(k)}))

其中近邻算子为

proxg(y)=arg⁡min⁡zg(z)+12∥z−y∥22\textbf{prox}_g (\mathbf y) = \arg\min_{\mathbf z} g(\mathbf z) + \frac{1}{2}\|\mathbf z-\mathbf y\|_2^2

这与模型(A2)(A2)相同,因此模型(A2)(A2)可以被理解为一种梯度下降方法。

这些一致性为算法的收敛性提供了良好的性质,至于收敛性就不细说了。

# 参考文献


  1. C. Wang, M. Yan, Y. Rahimi and Y. Lou. Accelerated Schemes for the L1/L2 Minimization, IEEE Transactions on Signal Processing, vol. 68, pp. 2660-2669, 2020. ↩︎

  2. 娄易非老师主页: https://sites.google.com/site/louyifei (opens new window) ↩︎ ↩︎