原文链接: https://www.csie.ntu.edu.tw/~cjlin/papers/cdl2.pdf.
对于分类数据集$(\pmb x_j,y_j),j=1,\cdots,l,x_j\in\mathbb{R}^n,y_i\in\{-1,+1\}$，SVM求解的是下面的无约束优化问题:
\[\min_{\pmb w}\quad\frac12\pmb w^T\pmb w+C\sum_{j=1}^l\xi(\pmb w;\pmb x_j,y_j)\]其中$\xi$是损失函数，$C$是惩罚参数。不同的损失函数对应不同的模型，比如Hinge loss对应的就是L1-SVM模型:
\[f(\pmb w)=\frac12\pmb w^T\pmb w+C\sum_{j=1}^l\max(1-y_j\pmb w^T\pmb x_j,0)\]Hinge loss损失的平方作为损失时则是L2-SVM模型:
\[f(\pmb w)=\frac12\pmb w^T\pmb w+C\sum_{j=1}^l\max(1-y_j\pmb w^T\pmb x_j,0)^2\]如果采用对数损失，模型变成了对率回归模型:
\[f(\pmb w)=\frac12\pmb w^T\pmb w+C\sum_{j=1}^l\log(1+e^{-y_j\pmb w^T\pmb x_j})\]我们可以通过对数据进行增广，已实现简化问题的目的:
\[\pmb x_j^T\gets[\pmb x_j^T\quad1],\pmb w^T\gets[\pmb w^T\quad b]\]这样，决策函数的偏置项$b$便会被消去:
\[\pmb w^T\pmb x\gets\pmb w^T\pmb x+b\]这篇文章讨论的是在大规模数据下用坐标下降法求解线性L2-SVM原问题。对于特征多，样本少的场景，求解对偶问题十分有效，我们在https://welts.xyz/2021/12/02/dcdm/中讨论了用坐标下降法求解SVM对偶问题的算法；而在数据多特征少的场景下，求解原问题的效率是更高的。
我们设$x_{ji}$是数据集的第$j$个数据的第$i$个特征，而数据集的样本数为$l$，特征数为$n$。坐标下降法必然是从初始点$\pmb w^0$开始，形成一个序列$\{\pmb w^k\}_{k=0}^\infty$。从$\pmb w^k$更新到$\pmb w^{k+1}$，需要经历$n$次单坐标迭代。因此我们设$\pmb w^{k,i}\in\mathbb{R}^n,i=1,\cdots,n$:
\[\pmb w^{k,i}=[w_1^{k+1}, \cdots,w_{i-1}^{k+1},w_i^k,\cdots,w_n^k]^T\text{ for }i=1,\cdots,n\]也就是第$k$轮迭代中，对前$i$个分量进行更新后的结果。将$\pmb w^{k,i}$更新到$\pmb w^{k,i+1}$的过程中，我们需要解决单变量子问题:
\[\min_{z}f(w_1^{k+1},\cdots,w_{i-1}^{k+1},w_{i}^k+z,w_{i+1}^{k},\cdots.w_n^k)\equiv\min_{z}f(\pmb w^{k,i}+z\pmb e_i)\]其中$\pmb e_i$是长度为$n$，第$i$个元素为1，其他元素都为0的向量。对子问题进行改写:
\[\begin{aligned} D_i(z)&=f(\pmb w^{k,i}+z\pmb e_i)\\ &=\frac12(\pmb w^{k,i}+z\pmb e_i)^T(\pmb w^{k,i}+z\pmb e_i)+C\sum_{j\in I(\pmb w^{k,i}+z\pmb e_i)}(b_j(\pmb w^{k,i}+z\pmb e_i))^2 \end{aligned}\]其中$b_j(\pmb w)=1-y_j\pmb w^T\pmb x_j$，$I(\pmb w)=\{j\vert b_j(\pmb w)>0\}$。考虑损失函数的性质，只有$b_j>0$对应的数据$\pmb x_j$才会对目标函数值产生影响。对于$I(\pmb w^{k,i}+z\pmb e_i)$不变的区间中，$D_i(z)$是一个二次函数。因此$D_i(z),z\in\mathbb R$是一个分段(piecewise)二次函数。由此文章选择使用牛顿法来优化该问题。如果$D_i(z)$是二次可微的，在给定$\bar{z}$处的牛顿步径应该是
\[-\frac{D_i'(\bar z)}{D_i''(\bar z)}\]函数的一阶微分:
\[D'_i(z)=w_i^{k,i}+z-2C\sum_{j\in I(\pmb w^{k,i}+z\pmb e_i)}y_jx_{ji}(b_j(\pmb w^{k,i}+z\pmb e_i))\]但函数无法进行二阶微分，这是因为$\max$函数不可微。面对这种情况，Mangasarian在2002年定义了泛化的二阶微分:
\[\begin{aligned} D_i''(z)&=1+2C\sum_{j\in I(\pmb w^{k,i}+z\pmb e_i)}y_{j}^2x^2_{ji}\\ &=1+2C\sum_{j\in I(\pmb w^{k,i}+z\pmb e_i)}x_{ji}^2 \end{aligned}\]所以设$z^0=0$，使用朴素的牛顿法，也就是按照下式迭代$z$直到$D_i’(z)=0$:
\[z^{t+1}=z^{t}-\frac{D_i'(z^t)}{D_i''(z^t)}\text{ for }t=0,1,\cdots\]所以L2-SVM的坐标下降算法:
Algorithm 1 L2-SVM的坐标下降算法
for $k=0,1,\cdots$
$\quad$for $i=1,2,\cdots$
$\qquad$固定$\pmb w^{k,i}$并求解上面所描述的子问题，对$w_{i}^{k+1}$进行更新.
Mangasarian证明了在某种假设下，上面的牛顿法可以在有限步结束并且正确求解子问题。但这种假设在真实情况下不会恒成立，这就导致取完整的牛顿步径进行下降不一定会让目标函数下降；此外，这种求解子问题的方法的求解代价很昂贵。
早期研究坐标下降求解L2-SVM的文章不打算精确求解子问题。在2001年的文章中，Zhang和Oles提出了CMLS算法，它求出来的近似解被限制在一个区域中，通过评估函数在这个区域内的二阶导数上界，然后用这个上界去替代二阶导数，也就是牛顿步径的分母。这种设置保证了函数下降。但仍存在两个问题:
而在这篇文章中，作者选取了保证坐标下降法收敛的充分条件:
\[D_i(z)-D_i(0)\leq-\sigma z^2\]其中$z$是当前所在的位置，$\sigma\in(0,\frac12)$是任意常数。考虑第一次迭代
\[d=-\frac{D_i'(0)}{D_i''(0)}\]第一次迭代是否恒满足收敛条件是很重要的，也就是判断
\[D_i(d)-D_i(0)\leq -\sigma d^2\]因为$d$很靠近0，$D_i(z)$也是一个二次函数，因此有
\[D_i(d)-D_i(0)=D_i'(0)d+\frac12D_{i}''(0)d^2\]因为$D_i’‘(z)>1$，所以
\[\begin{aligned} D_i'(0)d+\frac12D''_i(0)d^2+\sigma d^2 &=-\frac{D_i'(0)^2}{D_i''(0)}+(\frac12D''_i(0)+\sigma)\frac{D'_i(0)^2}{D''_i(0)^2}\\ &=-\frac12\frac{D_i'(0)^2}{D_i''(0)}+\sigma\frac{D_i'(0)^2}{D_i''(0)}\\ &=(\sigma-\frac12)\frac{D_i'(0)^2}{D_i''(0)}\leq 0\\ D_i(d)-D_i(0)&\leq-\sigma d^2 \end{aligned}\]也就是满足收敛条件。考虑到$D_i(z)$是分段二次函数，所以对于$z=d$，收敛条件不一定成立，但作者使用了线搜索。作者也提出了下面的定理:
Theorem 1 给定上述的牛顿步径$d$，$\forall\lambda\in[0,\bar\lambda],z=\lambda d$始终满足上面的收敛条件。其中
\[\bar{\lambda}=\frac{D_i''(0)}{H_i/2+\sigma},H_i=1+2C\sum_{j=1}^l x_{ji}^2\]
这样，作者希望通过线搜索找到满足上述条件的$\lambda$中的最大值，将$\lambda d$作为牛顿步径，也就是下面的线搜索算法:
Algorithm 2 用带线搜索的牛顿步径求解子问题
在算法2中，计算最耗时的一步是计算$D_i(\lambda d)$，现在考虑如何减少计算$D_i(\lambda d)$的次数。理论上，我们需要依次计算$D_i(d),D_i(\beta d),D_i(\beta^2d)\cdots$，逐一判断是否满足收敛条件。实际上，考虑到我们有定理一，所以在计算更小的$D_i(\lambda d)$之前，先判断$\lambda$是否满足定理一中的收敛条件，如果满足则直接返回，不需要进行$D_i(\lambda d)$的计算。由于$D_i’‘(0)$和$H_i$都只需要计算一次，因此相比于计算$D_i(\lambda d)$，与$\bar\lambda$进行比较更节省时间。此外，如果$\lambda=1$满足定理1，我们也没必要进行后面的计算了，可以直接退出。由此，算法2的计算时间可以大大缩短。
算法中有两个超参数：$\beta$和$\sigma$。因为$\lambda=1$时常常收敛，因此算法对$\beta$是不敏感的，设它为0.5即可；$\sigma$通常选择0.01。
算法1所生成的$\{\pmb w^k\}$序列是线性收敛的，也就是说，存在常数$\mu\in(0,1)$使得
\[f(\pmb w^{k+1})-f(\pmb w^*)\leq(1-\mu)(f(\pmb w^k)-f(\pmb w^*)),\forall k\]此外，该解序列可以收敛到全局最优$\pmb w^*$，它在$O(nC^3P^6(\#nz)^3\log(1/\epsilon))$次循环内能够到达$\epsilon$精确解，也就是
\[f(\hat{\pmb w})\leq f(\pmb w^*)+\epsilon\]其中$\#nz$是训练数据中非零值的数目，显然在稀疏数据中，这个值会很小，可以看出该算法适合稀疏数据集分类；$P$是数据集元素的绝对值上界:
\[P=\max_{ji}\vert x_{ji}\vert\]由于分类数据预处理中会进行归一化，所以$P\leq1$常常成立。具体的证明过程参考原论文附录。
现在考察算法1的一个外层循环，也就是从$\pmb w^k$到$\pmb w^{k+1}$的更新的计算复杂度。算法2的第二步主要是计算$D_i’(0)$和$D’‘_i(0)$:
\[\begin{aligned} D'_i(0)&=w_i^{k,i}-2C\sum_{j\in I(\pmb w^{k,i})}y_jx_{ji}(b_j(\pmb w^{k,i}))\\ D_i''(0)&=1+2C\sum_{j\in I(\pmb w^{k,i})}x_{ji}^2 \end{aligned}\]考虑数据集是稀疏的，计算$b_j(\pmb w),j=1,\cdots,l$的计算复杂度是$O(\#nz)$的，这样计算量是很大的。这里作者用到了一个trick:
\[b_j(\pmb w+z\pmb e_i)=b_j(\pmb w)-zy_jx_{ji}\]也就是说，我们需要保存并维护$b_j(\pmb w)$：求解上一个子问题得到$z$后，按照上面的公式对$b_j(\pmb w)$进行更新，注意到一次更新只需要考虑数据集的第$i$个特征，所以执行一次算法2的第二步的计算复杂度为$O(\frac{\#nz}{n})$，也就是每个特征中非零数据数目的平均值。类似的，在算法2的第三步，我们需要指定多次线搜索，每次线搜索需要计算$D_i(\lambda d)$，计算复杂度也是$O(\frac{\#nz}{n})$。实验表明算法2中线搜索的次数均不多，有时甚至不需要，因此算法2的复杂度，也就是算法1每轮迭代的复杂度约为$nO(\frac{\#nz}{n})=O(\#nz)$。
论文在这一节列举了3种实现上述算法的技巧，可以提高程序速度。第一种方法是从数据表示入手，由于该算法适合稀疏数据，因此数据存储也应该是稀疏矩阵。常用的存储稀疏矩阵的方法有“行存储”和“列存储“，这样存储能够快速读取稀疏数据集元素。由于该算法是以特征为单位求解子问题，因此按列存储是好主意。
算法1中求解子问题是按照特征顺序进行，特征的顺序有时会影响收敛速度，因此作者提出每一轮可以按照不同的特征顺序求解子问题，也就是对原数据集特征顺序进行一个排列。作者也证明了这样打乱顺序后，收敛效果和原算法相同。
如果数据的特征数极大，我们常常不需要在每个循环中遍历全部特征。相反，我们可以随机选择一个特征进行下降，也就是从$\pmb w^k$到$\pmb w^{k+1}$的过程中，只有一个分量被更新:
Algorithm 3 算法1的Online版本
初始化解向量$\pmb w^0$.
for $k=0,1,\cdots$
(a) 随机选取$i_k\in\{1,2,\cdots,n\}$.
(b) 固定解向量除$w_{i_k}$以外的分量，求解子问题。
