弹性网络

弹性网络（elastic net，常称弹性网回归）是一种同时使用 $L^1$ 和平方 $L^2$ 惩罚的回归方法。前者可把部分系数压到零，后者使系数收缩，并在正惩罚下保证解的唯一性。这里的网络指惩罚结构，不是神经网络。

目标函数与参数约定

设 $X\in {ℝ}^{n\times p}$、$y\in {ℝ}^n$。截距不受惩罚；将各列及响应中心化后，本文使用 $$\underset{\beta }{min}\frac{1}{2n}\Vert y-X\beta {\Vert }_2^2+{\lambda }_1\Vert \beta {\Vert }_1+\frac{{\lambda }_2}{2}\Vert \beta {\Vert }_2^2,{\lambda }_1\ge 0,{\lambda }_2\ge 0.$$ 常见软件用 ${\lambda }_1=\alpha r$、${\lambda }_2=\alpha (1-r)$，$r\in [0,1]$。原始论文还区分直接惩罚估计与另行缩放的版本；本文讨论上面明确写出的、现代实现常用的惩罚目标，不再额外缩放系数。

${\lambda }_2>0$ 时，二次惩罚使目标严格凸且强制随 $\Vert \beta \Vert \to \mathrm{\infty }$ 增大，所以最优解存在且唯一，即使预测变量多于观测。若 ${\lambda }_2=0$，退化为Lasso，系数则不总唯一；若 ${\lambda }_1=0$，就是此尺度约定下的岭回归。

一个坐标的软阈值推导

固定其余系数，将去掉第 $j$ 项后的残差记为 $r^{(-j)}$。令 $c_j=X_j^{𝖳}{r}^{(-j)}/n$、$z_j=X_j^{𝖳}{X}_j/n$，则与 ${\beta }_j=b$ 有关的目标为 $$\frac{1}{2}(z_j+{\lambda }_2)b^2-c_{j}b+{\lambda }_1|b|+\text{常数}.$$ 在 $b>0$ 区间令导数为零，得 $b=(c_j-{\lambda }_1)/(z_j+{\lambda }_2)$，只有 $c_j>{\lambda }_1$ 才落在该区间；负区间同理。若 $|c_j|\le {\lambda }_1$，零点的左右导数夹住零，因此最优值为零。合并成 $${\beta }_j\leftarrow \frac{\mathrm{sign}(c_j)(|c_j|-{\lambda }_1{)}_{+}}{z_j+{\lambda }_2},$$ 其中 $a_{+}=max(a,0)$，分母要求为正。依次更新各坐标便得到一种坐标下降计算。

正交设计的完整算例

取四个教学观测，两个预测列为 $X_1=(-1,-1,1,1{)}^{𝖳}$、$X_2=(-1,1,-1,1{)}^{𝖳}$，响应为 $y=3X_1-0.4X_2$。各列均值零，且 $X^{𝖳}X/4=I$，两个坐标可以分别求解。取 ${\lambda }_1=0.5,{\lambda }_2=0.25$，则 $${\hat{\beta }}_1=\frac{3-0.5}{1.25}=2,{\hat{\beta }}_2=0.$$ 预测为 $(-2,-2,2,2)$。第一项先去掉阈值的零点五，再除以一点二五；只做一次普通软阈值而不处理平方惩罚，会得到不同答案。

阈值区间由L1项决定，区间外斜率由L2项调整。标点为主例的c一等于三、估计系数二。

变式。 保持 ${\lambda }_1=0.5$，把 ${\lambda }_2$ 改为一，系数变为 $(1.25,0)$。若两列完全相同且 ${\lambda }_2>0$，交换两列系数不改变目标，而唯一性迫使两个系数相等；这给出完全相关特征的分组效应特例。

带最优性残差的计算

下面实现中心化后的坐标下降，返回原尺度截距和系数。代码不自动标准化，输入各列的单位应由建模者事先处理。要求 ${\lambda }_2>0$，通过次梯度最优性条件检查停止，而不是把迭代次数用完当成成功；输入不被改写，依赖NumPy。

Python 3

import math
import numpy as np

def elastic_net_coordinate(x, y, l1, l2, tolerance=1e-9, max_sweeps=10000):
    x, y = np.asarray(x, dtype=float), np.asarray(y, dtype=float)
    if not (x.ndim == 2 and y.ndim == 1 and x.shape[0] == len(y)
            and min(x.shape) > 0 and np.all(np.isfinite(x)) and np.all(np.isfinite(y))):
        raise ValueError("nonempty finite paired arrays required")
    if not (all(math.isfinite(v) for v in (l1, l2, tolerance))
            and l1 >= 0 and l2 > 0 and tolerance > 0
            and type(max_sweeps) is int and max_sweeps > 0):
        raise ValueError("invalid penalties or stopping rule")
    xm, ym = x.mean(axis=0), y.mean()
    xc, yc = x - xm, y - ym
    n, p = xc.shape
    z = np.sum(xc * xc, axis=0) / n
    beta, residual = np.zeros(p), yc.copy()
    for _ in range(max_sweeps):
        for j in range(p):
            residual += xc[:, j] * beta[j]
            c = xc[:, j] @ residual / n
            beta[j] = np.sign(c) * max(abs(c) - l1, 0) / (z[j] + l2)
            residual -= xc[:, j] * beta[j]
        gradient = -xc.T @ residual / n + l2 * beta
        violation = np.where(beta == 0, np.maximum(abs(gradient) - l1, 0),
                             abs(gradient + l1 * np.sign(beta)))
        if np.max(violation) <= tolerance:
            return float(ym - xm @ beta), beta
    raise RuntimeError("optimality residual did not reach tolerance")

x = np.array([[-1, -1], [-1, 1], [1, -1], [1, 1]], dtype=float)
y = 3 * x[:, 0] - 0.4 * x[:, 1]
print(elastic_net_coordinate(x, y, 0.5, 0.25))

结果为截距零、系数 [2. 0.]。每轮时间 $O(np)$，含中心化副本的存储为 $O(np+n+p)$；总时间还取决于达到容差所需轮数。交叉验证中的中心化、标准化和惩罚选择须在各训练折内完成。某系数被压到零不等于已证明该变量没有因果作用。