手撕统计机器学习 - 岭回归、感知机

本章所需知识点：

• 线性回归 (Linear Regression)
• 梯度下降 (Gradient Descent)

---

岭回归 (Ridge Regression)

线性回归有什么问题？

线性回归简单易用，但是除了不够高级，它还有什么缺点？为了回答这个问题我们需要回到如下正则方程的解：
$$\mathbf{\beta }=({\mathbf{X}}^{\mathrm{T}}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}$$
有个关键前提就是：${\mathbf{X}}^{\mathrm{T}}\mathbf{X}$ 一定可逆吗？这个条件虽然不苛刻，但是也可能出现较为病态的情况（行列式过于接近 $0$ 导致难以求逆），例如多重共线性的情况.

---

多重共线性

多重共线性表示特征之间具有较高的相关性关系，一个极端的例子就是：
$$\mathbf{X}=\left(\begin{array}{cc}1& 2\\ 2.4& 4.8\\ 0.3& 0.6\end{array}\right),\mathbf{y}=\left(\begin{array}{c}3\\ 7.2\\ 0.9\end{array}\right)$$
我们可以发现，两个特征 $p_1,p_2$ 之间的取值满足 $2p_1=p_2$ ，这就说明二者具有相关关系，计算 ${\mathbf{X}}^{\mathrm{T}}\mathbf{X}$ 可以得到
$${\mathbf{X}}^{\mathrm{T}}\mathbf{X}=\left(\begin{array}{cc}6.85& 13.7\\ 13.7& 27.4\end{array}\right)$$
此时行列式为 $0$ ，这就不可逆了，多重共线性的这个特性导致有很多种线性模型都可以适配，例如
$$y=x_1+x_2,y=3x_1,\cdots$$
等等，哪个线性模型有比较好的预测性能？这点无从得知. 越是接近多重共线性，就越有着出上述问题的风险.

---

岭回归

为了解决上述问题，岭回归 (Ridge Regression) 应运而生，其更改非常简洁，就是在损失函数上加一个 ${\ell }_2$ 范数惩罚项：
$$L(\mathbf{\beta })=\Vert \mathbf{X\beta }-\mathbf{y}{\Vert }_2^2+\lambda \Vert \mathbf{\beta }{\Vert }_2^2$$
这里 $\lambda >0$ 为惩罚参数，多重共线性越严重，就越需要将其调大，正是这个惩罚项，它就减轻了多重共线性带来的影响.

---

为什么它能避免多重共线性？我们下面来直接求解这个损失函数，首先先展开：
$$\begin{array}{rl}L(\mathbf{\beta })& =(\mathbf{X}\mathbf{\beta }-\mathbf{y}{)}^{\mathrm{T}}(\mathbf{X\beta }-\mathbf{y})+\lambda {\mathbf{\beta }}^{\mathrm{T}}\mathbf{\beta }\\ & ={\mathbf{\beta }}^{\mathrm{T}}({\mathbf{X}}^{\mathrm{T}}\mathbf{X})\mathbf{\beta }-{\mathbf{y}}^{\mathrm{T}}\mathbf{X\beta }-{\mathbf{\beta }}^{\mathrm{T}}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}+{\mathbf{y}}^{\mathrm{T}}\mathbf{y}+\lambda {\mathbf{\beta }}^{\mathrm{T}}\mathbf{\beta }\end{array}$$
对 $\mathbf{\beta }$ 求导有
$$\frac{\partial L}{\partial \mathbf{\beta }}=2{\mathbf{X}}^{\mathrm{T}}\mathbf{X}\mathbf{\beta }-2{\mathbf{X}}^{\mathrm{T}}\mathbf{y}+2\lambda \mathbf{\beta }$$
此时的正规方程为
$$({\mathbf{X}}^{\mathrm{T}}\mathbf{X}+\lambda \mathbf{I})\mathbf{\beta }={\mathbf{X}}^{\mathrm{T}}\mathbf{y}\Rightarrow \mathbf{\beta }=({\mathbf{X}}^{\mathrm{T}}\mathbf{X}+\lambda \mathbf{I}{)}^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}$$
只要 $\lambda$ 足够大，就可以使矩阵达到对角占优，此时必定是可逆的. 当然，考虑到计算成本，选择参数时采用适当大的参数就可以了 （利用交叉验证进行选择）.

---

具体实现时要注意的点

• 截距项将会是我们的实现与 sklearn 是否能保持一致的重要因素.
• 考虑截距项的方法是考虑中心化数据，也就是说每列数据都减去其平均值.
• 如果不需要截距项，那么就不需要进行中心化.

---

感知机

线性回归做分类可行吗？

线性回归用来做分类可以吗？这个问题对二分类而言就非常简单，我给出一条直线或者超平面 (对 $n$ 维空间，$n-1$ 维的子空间就是其超平面)，对于 ${\mathbb{R}}^p$ 空间我相当于将空间划为两个部分，这两个部分作为我的分类依据就很简单了.

直觉上这么做如此，但是线性回归输出的值必然是连续的，为了使其变为分类器，不妨直接使用符号函数，对于单条样本数据 $\mathbf{x}$：
$$f(\mathbf{x})=\mathrm{sign}({\mathbf{\beta }}^{\mathrm{T}}\mathbf{x}+{\beta }_0)=\mathrm{sign}\left(\sum _{i=1}^p{\beta }_i{x}_i+{\beta }_0\right)$$
此时 f:Rp→{−1,+1}f: \mathbb{R}^{p}\to \left\lbrace -1,+1 \right\rbracef:Rp→{−1,+1} 就成为了一个分类器.

---

损失的设置

我们之前就说过，机器学习通常就需要在已经给定的数据集上用某种函数或者模型进行拟合，优化训练的损失（误差），意图达到良好的泛化性能.

现在的问题就是，如何设计损失函数 (loss)？一个很自然的想法就是，误分类点的个数：
Loss1=♯{yi∣yi≠yi∗} \mathrm{Loss}_{1}= \sharp \left\lbrace y_{i}\mid y_{i} \neq y_{i}^{*} \right\rbrace Loss1​=♯{yi​∣yi​=yi∗​}
其中 $y_i$ 是预测标签，$y_i^{\ast }$ 是正确标签，这种损失符合直觉但是不适合计算，并且会具有不唯一性，如下图所示，虚线和实现似乎都能准确的划分，但是哪个泛化性能更好呢？这就是刚刚这个损失函数无法确定的地方.

---

损失函数

误分类点的个数不行，那么就考虑另一个连续值：误分类点与超平面之间的距离之和. 首先引入点到超平面距离公式，设 ${\mathbf{x}}_0$ 为超平面外一点，那么：

$$\frac{1}{\Vert \mathbf{\beta }{\Vert }_2}|{\mathbf{\beta }}^{\mathrm{T}}{\mathbf{x}}_0+{\beta }_0|$$

设误分类点全体为 $M$ ，则

$$\mathrm{Loss}=\frac{1}{\Vert \mathbf{\beta }{\Vert }_2}\sum _{{\mathbf{x}}_i\in M}|{\mathbf{\beta }}^{\mathrm{T}}{\mathbf{x}}_i+{\beta }_0|$$

接下来就是考虑绝对值如何去掉（优化需要），注意到的是对于误分类的点而言， $-y_i\cdot ({\mathbf{\beta }}^{\mathrm{T}}{\mathbf{x}}_0+{\beta }_0)⩾0$ ，因此可以将损失改为

$$\mathrm{Loss}=-\frac{1}{\Vert \mathbf{\beta }{\Vert }_2}\sum _{{\mathbf{x}}_i\in M}({\mathbf{\beta }}^{\mathrm{T}}{\mathbf{x}}_i+{\beta }_0)y_i$$

最终，令 $\Vert \mathbf{\beta }{\Vert }_2=1$ ，就得到了最终的损失函数：

$$L(\mathbf{\beta },{\beta }_0)=-\sum _{{\mathbf{x}}_i\in M}({\mathbf{\beta }}^{\mathrm{T}}{\mathbf{x}}_i+{\beta }_0)y_i$$

---

为什么可以令 $\Vert \mathbf{\beta }{\Vert }_2=1$ ？
因为这样我们就确保了唯一性，常数倍的 $\mathbf{\beta }$ 在去掉 $\Vert \mathbf{\beta }\Vert$ 项之前的损失函数里是等价的，令 $\Vert \mathbf{\beta }{\Vert }_2=1$ 可以简化损失函数并确定最终的 $\mathbf{\beta }$ .

到此，我们最终就是想求得
$$\hat{\mathbf{\beta }}=\arg \underset{\mathbf{\beta }}{\min }-\sum _{{\mathbf{x}}_i\in M}({\mathbf{\beta }}^{\mathrm{T}}{\mathbf{x}}_i+{\beta }_0)y_i$$

---

感知机的优化算法

损失函数确定后，关键就是如何优化，如果直接进行梯度下降，误分类点随时可能会发生变化，那么这个时候的收敛性不好保证. 那么一个想法就是每次都只选取一个误分类点进行梯度下降并反复进行.

首先求导有
$$\frac{\partial L}{\partial \mathbf{\beta }}=\sum _{{\mathbf{x}}_{\mathbf{i}}\in M}{y}_i\cdot {\mathbf{x}}_i,\frac{\partial L}{\partial {\beta }_0}=\sum _{{\mathbf{x}}_i\in M}{y}_i$$
每次随机抽取一个误分类点 ${\mathbf{x}}_i$ ，此时的参数更新方法为
$$\mathbf{\beta }\leftarrow \mathbf{\beta }+\eta y_i{\mathbf{x}}_i,{\beta }_0\leftarrow {\beta }_0+\eta y_i$$
其中 $\eta \in (0,1]$ 称为学习率 (learning rate)，用于控制更新的步长. 反复进行优化之后，最终将会收敛到没有误分类点（此时损失为 $0$）.

---

最终形式

参数更新的过程中，我们可以发现每次实际上都加的是一个恒定值，只是迭代更新的次数不同而已，设对样本点 $({\mathbf{x}}_i,y_i)$ ，学习过程当中更新了 $n_i$ 次，那么最终的参数为
$$\mathbf{\beta }=\sum _{i=1}^n\eta \cdot n_i\cdot y_i{\mathbf{x}}_i,{\beta }_0=\sum _{i=1}^n\eta n_i{y}_i$$
这个最终形式不免引起我们的遐想，$n_i$ 越大说明什么？说明对 ${\mathbf{x}}_i$ 这个样本点，它距离超平面过近，所以很容易误分类.

---

感知机的对偶形式

现在考虑一个新的场景，维度 $p$ 很高，但是样本数 $n$ 却没有那么多，这个时候原始算法的计算成本在哪里体现？事实上将会体现在 $({\mathbf{\beta }}^{\mathrm{T}}\mathbf{x}+{\beta }_0)y$ 的正负性判断上，每一次判断时，都需要进行 $O(p)$ 的计算，这在维度较高的情况下会使计算变得困难.

为了适应这个场景，不妨认为 $\mathbf{\beta }$ 已经是刚刚讨论的最终形式，只是每个 $n_i=0$ ，然后迭代更新时只考虑更新 $n_i$ . 对样本点 ${\mathbf{x}}_j$ 判断时判断的是
$$\left(\sum _{i=1}^n\eta n_i{y}_i{\mathbf{x}}_i^{\mathrm{T}}\cdot {\mathbf{x}}_j+\sum _{i=1}^n\eta n_i{y}_i\right)y_j$$
的正负性.

它如何降低计算量的？我们是把 ${\mathbf{\beta }}^{\mathrm{T}}\mathbf{x}$ 转换为了 ${\mathbf{x}}_i^{\mathrm{T}}{\mathbf{x}}_j$ ，似乎还是要做内积运算，但是这个时候，内积的运算就可以提前进行了，只需要提前计算下面的 Gram 矩阵
$$G={\left(\begin{array}{c}{\mathbf{x}}_i^{\mathrm{T}}{\mathbf{x}}_j\end{array}\right)}_{n\times n}$$
就能通过查表来简化计算.

---

---

总结

本次手撕的两个新模型为岭回归和感知机，都是基于线性回归的改进模型：

• 岭回归解决了线性回归容易受到多重共线性影响的问题；
• 感知机则使用线性回归解决了二分类的问题.

使用 numpy 手撕的过程中，上述的算法都需要注意数值稳定性.

---

实现

岭回归

class RidgeRegression:
    def __init__(self, penalty=1.0, fit_intercept=True):
        self.beta: np.array | None = None
        self._is_fit = False
        self.penalty = penalty
        self.fit_intercept = fit_intercept
        self.intercept_ = None
    
    def fit(self, X: np.array, y: np.array) -> None:
        if self.fit_intercept:
            X_mean = X.mean(axis=0)
            y_mean = y.mean()
            X_centered = X - X_mean
            y_centered = y - y_mean
        else:
            X_centered = X
            y_centered = y

        # 惩罚矩阵：对所有的特征系数进行惩罚，因为没有截距项了
        penalty_matrix = self.penalty * np.eye(X_centered.shape[1])

        A = X_centered.T @ X_centered + penalty_matrix
        b = X_centered.T @ y_centered
        
        try:
            self.coef_ = np.linalg.solve(A, b)
        except np.linalg.LinAlgError:
            self.coef_ = np.linalg.pinv(A) @ b

        if self.fit_intercept:
            self.intercept_ = y_mean - X_mean @ self.coef_
        else:
            self.intercept_ = 0.0

        self._is_fit = True

    
    def predict(self, X: np.array) -> np.array:
        if not self._is_fit:
            raise ValueError('请先进行模型拟合！')

        return X @ self.coef_ + self.intercept_

感知机

class PerceptronClassifier:
    def __init__(self, learning_rate=0.1, max_iter=1000) -> None:
        self.beta = None
        self.beta0 = None
        self._is_fit = False
        self.learning_rate = learning_rate
        self.mapping = None    # 编码映射
        self.max_iter = max_iter  # 防止数据集线性不可分，无限循环
    

    def fit(self, X: np.array, y: np.array) -> None:
        """感知机训练算法"""
        n, p = X.shape

        if np.unique(y).shape[0] != 2:
            raise Exception("[ERROR] 仅支持二分类，请检查标签数组！")
        
        # 编码
        self.mapping = {
            np.unique(y)[0]: -1,
            np.unique(y)[1]: 1,
        }
        lfunc = np.vectorize(lambda e: self.mapping[e])
        y = lfunc(y)


        # 接下来考虑感知机的训练，先初始化
        self.beta = np.random.randn(p)
        self.beta0 = np.random.randn()

        # 迭代更新，梯度下降
        for epoch in range(self.max_iter):
            misclassified = False
            for i in range(n):
                if y[i] * (self.beta @ X[i] + self.beta0) <= 0:
                    self.beta += self.learning_rate * y[i] * X[i]
                    self.beta0 += self.learning_rate * y[i]
                    misclassified = True
            if not misclassified:
                break
        

    def predict(self, X: np.array) -> np.array:
        """预测函数"""
        if not self._is_fit:
            raise Exception("[ERROR] 模型未训练，请先调用 fit 方法！")
        
        # 计算预测值并映射回原始标签
        scores = X @ self.beta + self.beta0
        y_pred = np.sign(scores)
        
        # 创建反向映射
        reverse_mapping = {v: k for k, v in self.mapping.items()}
        lfunc = np.vectorize(lambda e: reverse_mapping[e])
        return lfunc(y_pred)

---

感知机可以通过绘图来查看效果，生成线性可分数据集可以用：

from sklearn.datasets import make_blobs

X, y = make_blobs(n_samples=100,centers=2, n_features=2, random_state=42)

生成出来应该是

根据随机数初始化，效果可能有略微差别，最终训练结果为：