原文地址: https://www.csie.ntu.edu.tw/~cjlin/papers/l1.pdf.
前篇: https://welts.xyz/2022/01/27/l1/.
我们继续对《A Comparison of Optimization Methods and Software for Large-scale L1-regularized Linear Classification》这篇论文进行解读，主要是作者列举的以及他们所提出的带L1正则化优化问题的求解算法。这些算法还是可以分成三组:
- 分解方法；
- 转换成带约束问题的求解方法；
- 其他方法.
因为篇幅过大，笔者只会选择自认为重要的一部分讲解。
分解方法
先这里介绍几种序列选择和随机选择变量的方法。
循环坐标下降法
当算法处于第$k$轮大迭代时，设此时的解向量为$\pmb w$。此时算法依次更新各坐标分量，生成$\pmb w^{k,j}$，$k=1,\cdots,n+1$:
\[\pmb w^{k,j}=[w_1^{k+1},\cdots,w_{j-1}^{k+1},w_j^k,\cdots,w_n^k]^T,j=2,\cdots,n\]要将$\pmb w^{k,j}$更新到$\pmb w^{k,j+1}$，需要求解下面的优化问题:
\[\min_z\quad g_j(z)=\vert w_j^{k,j}+z\vert-\vert w_j^{k,j}\vert+L_j(z;\pmb w^{k,j})-L_j(0;\pmb w^{k,j})\tag1\]其中
\[L_j(z;\pmb w^{k,j})=L(\pmb w^{k,j}+z\pmb e_j)\]即在第$j$个分量更新后的损失函数项。循环坐标下降的流程如下图所示:
伪代码中的(17)对应本文的(1).
BBR算法
BBR算法用于求解带L1正则的对率回归问题。在每轮外循环中，BBR使用带置信域的牛顿法求解子问题(1)，带置信域的牛顿法可参考https://welts.xyz/2021/12/19/tron/，该方法需要求解一个模型子问题来调整置信域。尽管$g_j(z)$不可微，但还是能得到下面的式子:
\[g_j(z)=g_j(0)+g'_j(0)z+\frac12g_j''(\eta z)z^2\]其中$0<\eta<1$，同时定义一阶和二阶导数
\[g_j'(0)=\begin{cases} L_j'(0)+1&\text{if }w_j^{k,j}>0\\ L_j'(0)-1&\text{if }w_j^{k,j}<0\\ \end{cases}\\ g_j''(z)=L_j''(\eta z)\]当$w_j^{k,j}=0$时，$g_j(z)$不可导，我们后面讨论这一特殊情况。BBR算法会去寻找二阶导数在置信域内的上界:
\[U_j\geq g_j''(z),\forall\vert z\vert\leq\Delta_j\]同理也就得到了子问题函数的上界:
\[\hat{g}_j(z)=g_j(0)+g_j'(0)z+\frac12U_jz^2\]然后就可以有：让$\hat{g}_j(z)$下降的$z$必然也可以让$g_j(z)$下降:
\[g_j(z)-g_j(0)=g_j(z)-\hat{g}_j(0)\leq\hat{g}_j(z)-\hat{g}_j(0)<0\]如果损失函数是Logistics loss，那么BBR算法建议的$U_j$值设置为
\[U_j=C\sum_{i=1}^lx_{ij}^2F(y_i(\pmb w^{k,j})^T\pmb x_j,\Delta_j\vert x_{ij}\vert)\]其中$\Delta_j$为置信域大小，
\[F(r,\delta)=\begin{cases} 0.25, &\text{if }\vert r\vert\leq\delta\\ \dfrac1{2+e^{|r|-\delta}+e^{\delta-|r|}}&\text{otherwise} \end{cases}\]当$w_j^{k,j}\neq0$，BBR在置信域限制下求解子问题，得到闭式解
\[d=\min(\max(P(-\frac{g'_j(0)}{U_j},w_j^{k,j}),-\Delta_j), \Delta_j)\tag{2}\]其中
\[P(z,w)=\begin{cases} z&\text{if }\text{sgn}(w+z)=\text{sgn}(w)\\ -w&\text{otherwise} \end{cases}\]现在考虑$w_{j}^{k,j}=0$，也就是不可微的情景，当$L_j’(0)+1<0$，我们定义$g_j’(0)=L_{j}’(0)+1$，此时任意的$0<z\leq-g’_j(0)/U_j$都能让$\hat g_j(z)<\hat g_j(0)$，此时用(2)式将该点投影到$[-\Delta_j,\Delta_j]$中。$L’_j(0)-1>0$的情况与上述类似。如果$-1\leq L_j’(0)\leq1$，显然$z=0$已经是最优解了，因此不需要考虑。
BBR算法求解子问题(1)的伪代码如下:
BBR算法没有理论上的收敛证明，笔者也觉得该方法启发式意味太重了。
CDN算法
BBR利用上界对二阶导数进行替换。如果我们坚持使用二阶导数并在$z=0$处获取牛顿方向，则收敛速度会更快。这种方法在求解L2正则的L2 loss-SVM中已经被提及(参考https://welts.xyz/2022/01/18/cdm/)，此处作者打算将该方法拓展到L1正则化的对率回归问题上。
CDN方法不考虑L1正则项的不可导问题，只考虑损失项，也就是
\[\min_z\quad\vert w_{j}^{k,j}-z\vert-\vert w_j^{k,j}\vert+L'_j(0)z+\frac12L_j''(0)z^2\tag{3}\]显然这是一个软阈值函数(参考https://welts.xyz/2022/01/22/soft_thresholding/)。可以直接给出闭式解
\[d=\begin{cases} -\frac{L_j'(0)+1}{L_j''(0)}&\text{if }L_j'(0)+1\leq L_j''(0)w_j^{k,j}\\ -\frac{L_j'(0)-1}{L_j''(0)}&\text{if }L_j'(0)-1\geq L_j''(0)w_j^{k,j}\\ -w_j^{k,j}&\text{otherwise} \end{cases}\tag4\]因为(3)只是对原子问题的二次近似，所以$d$无法确保目标函数下降。为了确保收敛性，作者加上线搜索步骤来寻找保证下降的牛顿步径$\lambda d,\lambda\in(0,1)$，也就是
\[f(\pmb w^{k,j}+\lambda d\pmb e_j)-f(\pmb w^{k,j})=g_j(\lambda d)-g_j(0)\leq \sigma(\lambda d)^2\]其中$\sigma\in(0,1)$是常数。但Tseng和Yun在2007年指出上述条件在不平滑正则项$\Vert w\Vert_1$下效果不好。作者听取了他们的意见，改用下述条件
\[g_j(\lambda d)-g_j(0)\leq\sigma\lambda(L'_j(0)d+|w_j^{k,j}+d|-|w_{j}^{k,j}|)\tag5\]为了找到合适的$\lambda$，CDN算法使用回溯直线搜索，依次查询$\lambda=1,\beta,\beta^2,\cdots$，直到$\lambda$满足上述的下降条件。CDN算法流程如下:
伪代码中的(27)对应本文的(4)，(28)对应(5)
作者证明了上述线搜索步骤会在有限步终止，同时，在对率回归问题下，$\{\pmb w^k\}$会在有限步内到达最优点。
约束问题求解
置信域牛顿法(TRON)
考虑L1正则化问题的等价问题
\[\begin{aligned} \min_{\pmb w^+,\pmb w^-}\quad&\sum_{j=1}^nw_j^++\sum_{j=1}^nw_j^-+C\sum_{i=1}^n\xi(\pmb w^+-\pmb w^-;\pmb x_i,y_i)\\ \text{s.t.}\quad& w_j^+\geq0,w_j^-\geq0,\forall j. \end{aligned}\tag{6}\]然后我们令
\[\pmb w=\begin{bmatrix} \pmb w^+\\ \pmb w^-\\ \end{bmatrix}\in\mathbb{R}^{2n}\]将问题(6)转换成
\[\begin{aligned} \min_{\pmb w}\quad&\bar{f}(\pmb w)\\ \text{s.t.}\quad&\pmb w\in\Omega=\{l_j\leq w_j\leq u_j,\forall j\} \end{aligned}\]其中$\bar{f}$是(6)中的目标函数，和原问题目标函数作区分。$\pmb l$和$\pmb u$是下界与上界。可以用带置信域的牛顿法(TRON)解决该问题，仍然可以参考https://welts.xyz/2021/12/19/tron/。
内点法(IPM)
原理和方法可参考https://welts.xyz/2022/01/30/ipm/。内点法将下面的不等式约束问题
\[\begin{aligned} \min_{\pmb w,\pmb u}\quad&\sum_{j=1}^nu_j+C\sum_{i=1}^l\xi(\pmb w;\pmb x_i,y_i)\\ \text{s.t.}\quad&-u_j\leq w_j\leq u_j,j=1,\cdots,n \end{aligned}\tag{7}\]转换成下面的无约束优化问题
\[\phi(b,\pmb w,\pmb u)=t(\sum_{j=1}^nu_j+C\sum_{i=1}^j\xi(\pmb w,b;\pmb x_i,y_i))-\sum_{i=1}^n\log(u_j^2-w_j^2)\tag{8}\]$t$是参数，用来衡量问题(7)和(8)的最优解间的差距，当$t\to\infty$，(8)的解将收敛到$(7)$的解。
L1正则化的L2-loss SVM
前面讨论的都是L1正则化的对率回归问题，这里讨论L2-loss SVM。优化问题可写作
\[\min_{\pmb w}\quad f(\pmb w)=\Vert w\Vert_1+C\sum_{i\in I(\pmb w)}b_i(\pmb w)^2\]其中
\[b_i(\pmb w)=1-y_i\pmb w^T\pmb x_i,\quad I(\pmb w)=\{i\vert b_i(\pmb w)>0\}\]因此，损失函数项显然就是
\[L(\pmb w)=C\sum_{i\in I(\pmb w)}b_i(\pmb w)^2\]损失函数项一阶可导，但二阶不可导，因此前面很多方法不一定适用。但是损失函数在大部分位置都是二阶可导的，此外，$\nabla L(\pmb w)$是Lipschitz连续的，因而可以广义黑塞矩阵处处存在。在下面的部分中，作者对CDN和TRON算法进行拓展，应用到L1正则的L2-loss SVM中。
CDN求解
在问题(1)中，我们有
\[L_j(z)=C\sum_{i\in I(\pmb w^{k,j}+z\pmb e_j)}b_i(\pmb w^{k,j}+z\pmb e_j)^2\]其梯度
\[L_j'(0)=-2C\sum_{i\in I(\pmb w^{k,j})}y_ix_{ij}b_i(\pmb w^{k,j})\]不幸的是如果出现$b_i(\pmb w^{k,j})=0$的样本，$L’‘_j(0)$不是良定义的。因而定义广义的二阶导数
\[2C\sum_{i\in I(\pmb w^{k,j})}x_{ij}^2\]通过这样的替换，CDN算法就可以用于求解L2-loss SVM了。有一个问题是，当$x_{ij}=0$，上述广义二阶导数为0，由于在CDN中常常将其作为分母，所以此时需要保证二阶导数恒正:
\[\max(2C\sum_{i\in I(\pmb w^{k,j})}x_{ij}^2,\varepsilon)\]其中$\varepsilon$是小正数。
TRON求解
和前面(6)一样，带置信域的牛顿法是用来求解原问题的等价带约束优化问题。我们要求$f(\bar{\pmb w})$的梯度和Hessian矩阵。梯度不难求:
\[\nabla\bar{f}(\pmb w)=\pmb e+2C(\bar{X}_{I,:}^T\bar{X}_{I,:}\pmb w-\bar{X}_{I,:}^T\pmb y_{I})\]其中$e\in\mathbb{R}^{2n}$是全1矩阵。$\bar{X}=[X\quad-X]$。至于Hessian矩阵，我们使用广义Hessian
\[2C\bar{X}^TD\bar{X}\]其中$D\in\mathbb{R}^{l\times l}$是对角矩阵，对角元素满足
\[D_{ii}=\begin{cases} 1&\text{if }b_i(\pmb w)>0\\ 0&\text{if }b_i(\pmb w)\leq0 \end{cases}\]然后就可以像前面那样执行TRON算法。
