回归与迭代优化

这是一篇关于回归与迭代优化的文章，主要内容包括线性回归、逻辑回归、梯度下降、牛顿法等。

Linear Regression

线性回归 （Linear Regression）是一种用于建立输入变量与输出变量之间关系的线性模型。线性回归的目标是找到一条直线，使得所有数据点到直线的距离之和最小。这个距离可以是点到直线的垂直距离的平方和，也可以是点到直线的水平距离的平方和。

线性回归的目标函数 / 估计函数 （Hypothesis）可以表示为：

$$h_{\theta}(x) = \theta_0 + \theta_1 x_1 + \theta_2 x_2 + \cdots + \theta_n x_n$$

简单起见，我们把截距项 $\theta_0$ 也加入到 $\theta$ 中，这样可以统一表示为：

$$h_{\theta}(x) = \theta^T x$$

其中

$$x = \begin{bmatrix} 1 \\ x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}, \quad \theta = \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \\ \vdots \\ \theta_n \end{bmatrix}$$

Idea ：最小化损失

我们常用均方误差（Mean Squared Error）平方误差作为损失函数（Loss Function）：

$$L(h_{\theta}(x), y) = (h_{\theta}(x) - y)^2$$

我们的目标是找到一组参数 $\theta$，使得所有样本的损失函数之和最小：

$$E(\theta) = \sum_{i=1}^{m} (h_{\theta}(x^{(i)}) - y^{(i)})^2$$

改为矩阵表示：

$$E(\theta) = ||X\theta - y||^2 = (X\theta - y)^T(X\theta - y) = \theta^T X^T X \theta - 2 \theta^T X^T y + y^T y$$

让我们把 $X^T X$ 记为 $A$，$X^T y$ 记为 $b$，$y^T y$ 记为 $c$，则：

$$E(\theta) = \theta^T A \theta - 2 \theta^T b + c$$

我们对目标函数 $E(\theta)$ 求导，得到：

$$\frac{\partial E(\theta)}{\partial \theta} = 2 A \theta - 2 b$$

即

$$\frac{\partial E(\theta)}{\partial \theta} = 2 X^T X \theta - 2 X^T y$$

所以，现在我们的目标是找到 $\theta$，使得 $\frac{\partial E(\theta)}{\partial \theta} = 0$

正规方程法（Normal Equation）

这可以分为两种情况：

1. 当 $X^T X$ 可逆时，即 $X^T X$ 是满秩矩阵，这时 $\theta = (X^T X)^{-1} X^T y$，这里的 $(X^T X)^{-1} X^T$ 称为伪逆矩阵（Pseudo Inverse Matrix）。
2. 当 $X^T X$ 不可逆时，可以使用正则化（Regularization）来使得 $X^T X$ 可逆，或者使用梯度下降（Gradient Descent）等迭代优化方法。

总之，线性回归可以分为三步：

1. 构造特征矩阵 $X$ 和标签 $y$
2. 求解伪逆矩阵 $(X^T X)^{-1} X^T$
3. 计算 $\theta$ 并进行预测即， $\theta = (X^T X)^{-1} X^T y$

这种方法有时候也称为正规方程法。在简单情况下，即损失函数是一个凸函数时，我们可以直接解出导数为 0 的点，因此可以使用这种方法。

当然，由于要求逆矩阵，这种方法的时间复杂度会比较高

在实现中，主要的复杂度在于求逆矩阵，如何求逆矩阵是一个线性代数课程中的重要内容，这里不再赘述，因为，如果我是一个大一新生，或许我还会，但现在我是一个大二老东西，我只会调用 numpy.linalg.inv 函数。

代码实现

import numpy as np
from sklearn.metrics import r2_score
from typing import Optional
from sklearn.model_selection import train_test_split
from sklearn import datasets

class LinearRegression_normal:
    def __init__(self):
        self.coef_: Optional[np.ndarray] = None
        self.intercept_: Optional[float] = None
        self._theta: Optional[np.ndarray] = None

    def fit(self, X_train: np.ndarray, y_train: np.ndarray) -> None:
        assert X_train.shape[0] == y_train.shape[0], "the size of X_train must be equal to the size of y_train"
        X_b = np.hstack([np.ones((X_train.shape[0], 1)), X_train])

        # 如果矩阵是满秩的，那么伪逆矩阵等于：X_b.T.dot(X_b).I.dot(X_b.T)
        if np.linalg.matrix_rank(X_b) == X_b.shape[1]:
            pseudo_inverse = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T)
        else:
            raise ValueError("the matrix is not full rank")
        
        self._theta = pseudo_inverse.dot(y_train)

        self.intercept_ = self._theta[0] # 截距
        self.coef_ = self._theta[1:] # 系数

    def predict(self, X_predict: np.ndarray) -> np.ndarray:
        def _hypothese(theta: np.ndarray, X: np.ndarray) -> np.ndarray:
            return X.dot(theta)
        
        assert self._theta is not None, "must fit before predict"
        assert X_predict.shape[1] == len(self.coef_), "the feature number of X_predict must be equal to X_train"
        # 添加一列全为1的列
        X_b = np.hstack([np.ones((X_predict.shape[0], 1)), X_predict])
        return _hypothese(self._theta, X_b)

    def score(self, X_test: np.ndarray, y_test: np.ndarray) -> float:
        y_predict = self.predict(X_test)
        return r2_score(y_test, y_predict)

    def __repr__(self) -> str:
        return "LinearRegression_normal()"

breast = datasets.load_breast_cancer()
X = breast.data
y = breast.target

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=666)
lr = LinearRegression_normal()
lr.fit(X_train, y_train)
print(lr.score(X_test, y_test))

Logistic Regression

逻辑回归（Logistic Regression）是一种用于解决分类问题的线性模型。逻辑回归的目标是找到一个函数，使得输入变量的线性组合映射到一个概率值，显然概率的值域是 0 到 1 之间。

首先引入 Sigmoid 函数（Logistic Function）：

$$\sigma(z) = \frac{1}{1 + e^{-z}}$$

这个函数的图像是一条 S 型曲线，值域在 0 到 1 之间。

我们知道，线性回归可以将输入变量映射到输出变量，通常输出的值域是实数域。然而，在分类问题中，我们需要的是一个离散值，即 0 或 1。

逻辑回归的思想是将线性回归的输出通过 Sigmoid 函数映射到 0 到 1 之间，然后根据阈值来判断类别。

回顾线性回归的估计函数：

$$h_{\theta}(x) = \theta^T x$$

我们将其通过 Sigmoid 函数映射到 0 到 1 之间，即为逻辑回归的估计函数 ：

$$y = \frac{1}{1 + e^{-\theta^T x}}$$

可以变化为：

$$ln \frac{y}{1-y} = \theta^T x$$

若将 y 视为正类的概率，1-y 视为负类的概率，那么两者的比值 $y / (1-y)$ 称为几率（Odds）。几率的对数称为对数几率（Log Odds），也称为对数几率比（Log Odds Ratio）。对数几率比可以看作是输入变量的线性组合。

则，我们可以继续重写为

$$ln \frac{p(y = 1|x)}{1 - p(y = 1|x)} = \theta^T x$$

显然有

$$p(y = 1|x) = \frac{e^{\theta^T x}}{1 + e^{\theta^T x}}$$

$$p(y = 0|x) = \frac{1}{1 + e^{\theta^T x}}$$

于是，我们可以得到逻辑回归的似然函数（Likelihood Function）：

$$L(\theta) = \prod_{i=1}^{m} p(y^{(i)}|x^{(i)};\theta)$$

我们可以使用最大似然估计（Maximum Likelihood Estimation）来求解 $\theta$，即：

$$\theta = \arg \max_{\theta} L(\theta)$$

通过 Cross-entropy loss 的定义，我们可以得到逻辑回归的损失函数：

$$E(\theta) = - \sum_{i=1}^{m} [ -y^{(i)} \theta^T x^{(i)} + ln(1 + e^{\theta^T x^{(i)}})]$$

为了使得损失函数最小，我们对损失函数求导：

$$\frac{\partial E(\theta)}{\partial \theta} = - \sum_{i=1}^{m} [ -y^{(i)} x^{(i)} + \frac{e^{\theta^T x^{(i)}}}{1 + e^{\theta^T x^{(i)}}} x^{(i)}]$$

观察到我们可以利用 Sigmoid 函数的性质来简化上式：

设 Sigmoid 函数为 $z = \frac{1}{1 + e^{-\theta^T x}}$，则有：

$$\frac{\partial E(\theta)}{\partial \theta} = - \sum_{i=1}^{m} [ z({\theta^T x^{(i)}}) - y^{(i)}] x^{(i)} = 0$$

注意到导数为 0 的条件是 $z({\theta^T x^{(i)}}) = y^{(i)}$，这正是我们的目标。

迭代优化（Iterative Optimization）

然而，逻辑回归的损失函数是非凸函数，无法直接求解最优解。我们可以使用迭代优化 （Iterative Optimization）的方法

迭代优化分为两种

1. 确定性方法 (Deterministic Methods)
2. 随机性方法 (Stochastic Methods)

对于 Deterministic Methods，分为：

1. First Order Methods：methods that use only the gradient.
2. Second Order Methods：methods that use the gradient and the Hessian matrix.
   其中 Hessian 矩阵是损失函数的二阶导数矩阵，即：

   $$H (f)_i,_j = \frac{\partial^2}{\partial x_i \partial x_j} f(x)$$

Gradient Descent

Motivation: to minimize the local first-order Taylor approximation of 𝑓

梯度下降 （Gradient Descent，Cauchy 1847） 是一种 First Order Methods，运用了 Taylor 展开式，其想法就是最小化损失函数的局部一阶泰勒展开式，即：

$$min_{x} f(x) \approx min_{x} f(x_0) + \nabla f(x_0)^T (x - x_0)$$

那么更新规则即为：

$$x_{t+1} = x_t - \eta _t \nabla f(x_t)$$

其中 $\eta _t > 0$ 为步长 (Step Size)，也称为学习率 (Learning Rate)。

通过不断迭代，最终我们可以找到损失函数的临界点（Critical Points），即梯度为 0 的点，又称为稳定点（Stationary Points），即驻点。

当然，如果我们可以直接解出导数为 0 的点，比如在线性回归 （矩阵满秩）中，我们可以直接求解出最优解为 $\theta = (X^T X)^{-1} X^T y$，这时候我们可以直接求解出最优解。

Newton’s Methods

牛顿法 （Newton’s Methods ，Newton 1685）是一种 Second Order Methods，其想法是最小化损失函数的局部二阶泰勒展开式，即：

$$min_{x} f(x) \approx min_{x} f(x_0) + \nabla f(x_t)^T (x - x_t) + \frac{1}{2} (x - x_t)^T \nabla^2 f(x_t) (x - x_t)$$

对两边求导得到

$$\frac{\partial f(x)}{\partial x} = \nabla f(x_t) + \nabla^2 f(x_t) (x - x_t) = 0$$

更新规则为

$$x = x_t - [\nabla^2 f(x_t)]^{-1} \nabla f(x_t)$$

Go back to Logistic Regression

现在我们学习了 迭代优化的方法，让我们回到难住我们的那个问题：

我这里将 \theta 换成 w，和 ppt 保持一致

$$\nabla E(w) = - \sum_{i=1}^{m} [ -y^{(i)} x^{(i)} + \frac{e^{w^T x^{(i)}}}{1 + e^{w^T x^{(i)}}} x^{(i)}]$$

$$\nabla E(w) = - \sum_{i=1}^{m} [ z({w^T x^{(i)}}) - y^{(i)}] x^{(i)} = X^T (z(Xw) - y) = 0$$

让我们用 牛顿法 来解决这个问题

还记得我们的更新规则吗？

$$x = x_t - [\nabla^2 f(x_t)]^{-1} \nabla f(x_t)$$

（这里引入 Hessian 矩阵以使得式子更简洁）

回顾梯度 （Gradient）的定义： 梯度是一个向量，其每个元素是函数在对应坐标轴上的偏导数。

而海森矩阵 （Hessian Matrix）是一个矩阵，其每个元素是函数在对应坐标轴上的二阶偏导数。

我们接下来来求这几个东西

数学演算

损失函数

$$E(w) = - \sum_{i=1}^{m} [ -y^{(i)} w^T x^{(i)} + ln(1 + e^{w^T x^{(i)}})]$$

梯度（损失函数的一阶导数）

$$\nabla E(w) = X^T (z(Xw) - y)$$

Hessian 矩阵（损失函数的二阶导数）

这里不会写就不解了，直接给出结果：

$$H(w) = \sum_{i=1}^{m} x_i z(w^T x_i) z(-w^T x_i) x_i^T = X^T R X$$

其中 R 是一个对角矩阵，对角线上的元素为 $R_{ii} = z(w^T x_i) z(-w^T x_i)$

更新规则

综上，我们就可以得到逻辑回归的更新规则：

$$w = w_t - H(w_t)^{-1} \nabla E(w_t)$$

其中 H(w) 是损失函数的 Hessian 矩阵，$\nabla E(w)$ 是损失函数的梯度。

与 Linear Regression 的比较

对于线性回归

损失函数

$$E(w) = ||Xw - y||^2 = (Xw - y)^T(Xw - y)$$

Gradient

$$\nabla E(w) = X^T X w - X^T y$$

Hessian

$$H = \nabla^2 E(w) = X^T X$$

可以看到，线性回归的 Hessian 矩阵是一个常数矩阵，如果我们使用牛顿法

$$w = w_t - (X^T X)^{-1} (X^T X w_t - X^T y) = (X^T X)^{-1} X^T y$$

牛顿法直接求解出了最优解，这正好与我们之前提到的正规方程法一致。

而逻辑回归的 Hessian 矩阵是一个变化的矩阵，因此我们需要不断迭代求解。

参考资料

• 辛鑫老师的 PPT (非常羞愧，上课睡大觉，考试周破防才发觉老师做的 PPT 很好很用心)
• 《统计学习方法》李航
• 《机器学习》周志华