线性模型（线性回归，logistic回归）公式推导+numpy实现

用 NumPy 手写所有主流 ML 模型，普林斯顿博士后 David Bourgin 开源了一个非常剽悍的项目
矩阵求导公式

线性模型

线性回归

Code

class LinearRegression:
    def __init__(self, fit_intercept=True):
        '''
        通过法线方程拟合的普通最小二乘回归模型
        '''
        self.beta = None
        self.fit_intercept = fit_intercept

    def fit(self, X: np.ndarray, y):
        '''
        X np.ndarray <N, M>
        y np.ndarray <N, k>
        '''
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]  # 拓展第一列，用来和截距做内积

        pseudo_inverse = np.dot(np.linalg.inv(np.dot(X.T, X)), X.T)  # 
        self.beta = np.dot(pseudo_inverse, y)

    def predict(self, X:np.ndarray):
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]  # 拓展第一列，用来和截距做内积

        return np.dot(X, self.beta)

分析

假设我们有数据集$X^{n\times m}$，$n$是样本数量，$m$是特征数，$x_i^{1\times m}$。
线性回归的损失函数就是$E_{(w,b)}=\frac{1}{2}{\sum }_{i=1}^n(x_i\cdot w+b-y_i{)}^2$，目标是求出$E_{(w,b)}$最小时的$(w,b)=\hat{w}$，这个目标就最小二乘法，不需要梯度下降，梯度下降的结果也一样。
代码

if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]  # 拓展第一列，用来和截距做内积

就是为了满足$(w,b)=\hat{w}$这一步，将$X$在列方向拓展一列，变成${\hat{X}}^{n\times (m+1)}$，结果${\hat{w}}^{(m+1)\times 1}$。
接着将损失函数$E_{(w,b)}=\frac{1}{2}{\sum }_{i=1}^n(x_i\cdot w+b-y_i{)}^2$写成$$E_{\hat{w}}=\frac{1}{2}{(\hat{X}\cdot \hat{w}-Y)}^T\cdot (\hat{X}\cdot \hat{w}-Y)$$
然后就是矩阵求导，这里使用分母布局矩阵求导公式
$$\frac{\partial E_{\hat{w}}}{\partial \hat{w}}=\frac{1}{2}\times 2\times (\hat{X}\cdot \hat{w}-Y)\frac{\partial \hat{X}\cdot \hat{w}}{\partial \hat{w}}$$
$$\frac{\partial E_{\hat{w}}}{\partial \hat{w}}={\hat{X}}^T\cdot (\hat{X}\cdot \hat{w}-Y)$$
让一阶导数为0
$$\frac{\partial E_{\hat{w}}}{\partial \hat{w}}={\hat{X}}^T\cdot \hat{X}\cdot \hat{w}-\hat{X}\cdot Y=0$$
$${\hat{X}}^T\cdot \hat{X}\cdot \hat{w}={\hat{X}}^T\cdot Y$$
求得
$$\hat{w}={({\hat{X}}^T\cdot \hat{X})}^{-1}\cdot {\hat{X}}^T\cdot Y$$

几何解释

有$N$个样本点，每个样本点有$m$维，这样就有一个$m$维空间。我们的目标$Y$一般不在这个$m$维空间中（噪声，随机性等），如果$Y$在$m$维空间中，那么就是求方程组
$$XW=Y$$
$$W=X^{-1}\cdot Y$$
不过现实数据中$X$一般都不可逆，上述公式无解
最小二乘法是$Y$投影到$m$维空间上，在$m$维空间表示为$\acute{Y}=XW$，值得注意的是$Y-\acute{Y}$和$m$维空间正交，也就是
$$X^T\cdot (Y-XW)=0$$
MIT线性代数矩阵投影

概率视角

概率视角来看，最小二乘隐含假设假设噪声服从高斯分布
对应的代码
求${({\hat{X}}^T\cdot \hat{X})}^{-1}\cdot {\hat{X}}^T$

pseudo_inverse = np.dot(np.linalg.inv(np.dot(X.T, X)), X.T)

算$\hat{w}={({\hat{X}}^T\cdot \hat{X})}^{-1}\cdot {\hat{X}}^T\cdot Y$

self.beta = np.dot(pseudo_inverse, y)

Ridge 回归

加了$L2$正则化
目标函数
$E_{(w,b)}=\frac{1}{2}{\sum }_{i=1}^n(x_i\cdot w+b-y_i{)}^2+\frac{\lambda }{2}||\hat{w}|{|}^2$写成$$E_{\hat{w}}=\frac{1}{2}{(\hat{X}\cdot \hat{w}-Y)}^T\cdot (\hat{X}\cdot \hat{w}-Y)+\frac{\lambda }{2}{\hat{w}}^T\cdot \hat{w}$$

矩阵求导环节
$$\frac{\partial E_{\hat{w}}}{\partial \hat{w}}={\hat{X}}^T\cdot \hat{X}\cdot \hat{w}-\hat{X}\cdot Y+\lambda \hat{w}=0$$
$I$是单位矩阵
$$\frac{\partial E_{\hat{w}}}{\partial \hat{w}}=({\hat{X}}^T\cdot \hat{X}+\lambda I)\cdot \hat{w}-\hat{X}\cdot Y=0$$
求得
$$\hat{w}={({\hat{X}}^T\cdot \hat{X}+\lambda I)}^{-1}\cdot {\hat{X}}^T\cdot Y$$
这就是岭回归的解析解，抑制过拟合，提一下，截距部分$b$不参与正则化，只抑制$w_i$

code

def fit(self, X, y):
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]
		# 在训练部分，加上单位矩阵就行了
        A = self.alpha * np.eye(X.shape[1])

        pseudp_inverse = np.dot(np.linalg.inv(X.T @ X + A), X.T)
        self.beta = np.dot(pseudp_inverse, y)

logistics回归

在二分逻辑回归中，标签是 { 0 , 1 } \{0,1\} {0,1}，而且逻辑回归不仅可以预测类别，还可以得到概率。后面直接用$w^T\cdot x_i$代替$w^T\cdot x_i+b$，为了书写方便。
逻辑回归中使用激活函数 $sigmoid$ 处理线性回归的计算结果。
$$p(y_i=1|x_i)=sigmoid(w^T\cdot x_i)=\frac{1}{1+e^{-w^T\cdot x_i}}$$ $$p(y_i=0|x_i)=1-p(y_i=1|x_i)=\frac{e^{-w^T\cdot x_i}}{1+e^{-w^T\cdot x_i}}$$
那么将上述两个结合在一起就可以得到$$p(y_i|x_i)=p(y_i=1|x_i{)}^{y_i}\times p(y_i=0|x_i{)}^{1-y_i}$$，然后做极大似然估计$MLE$，乘法符号又不好处理，使用$log$把目标变成加法形式。我们对$$log(p(y_i|x))$$做极大似然估计，得到$$y_{i}ln(p(y_i=1|x_i))+(1-y_i)ln(p(y_i=0|x_i))$$，会发现这就是交叉熵，然后再拓展到全部的数据并且加上正则化得到如下的公式：
$$\text{NLL}=-\frac{1}{N}\left[\left(\sum _{i=0}^N{y}_i\log (\widehat{y_i})+(1-y_i)\log (1-\widehat{y_i})\right)-\frac{\gamma }{2}||\mathbf{w}|{|}_2\right]$$

梯度下降

回到原来的公式$$y_{i}ln(p(y_i=1|x_i))+(1-y_i)ln(p(y_i=0|x_i))$$等于$$y_{i}ln(p(y_i=1|x_i))+(1-y_i)ln(1-p(y_i=1|x_i))$$再化简 $$ln(1-p(y_i=1|x_i))+y_{i}ln(\frac{p(y_i=1|x_i)}{1-p(y_i=1|x_i)})$$ 把$p(y_i=1|x_i)=\frac{1}{1+e^{-w^T\cdot x_i}}$代入得到$$-ln(1+e^{w^T\cdot x_i})+y_i(w^T\cdot x_i)$$这就是化简的最终结果，这是负的交叉熵，前面还要加一个负号最后得到$$L(w)=ln(1+e^{w^T\cdot x_i})-y_i(w^T\cdot x_i)$$

下面开始求导

$$\frac{\partial L(w)}{\partial w}=\frac{\partial ln(1+e^{w^T\cdot x_i})}{\partial w}-\frac{\partial y_i(w^T\cdot x_i)}{\partial w}=\frac{e^{w^T\cdot x_i}}{1+e^{w^T\cdot x_i}}{x}_i-y_i{x}_i$$ 注意$\frac{e^{w^T\cdot x_i}}{1+e^{w^T\cdot x_i}}$就是$p(y_i=1|x_i)$所以最后求得导数形式就是$$\frac{\partial L(w)}{\partial w}=(p(y_i=1|x_i)-y_i)x_i$$

Code

numpy写主流 ML 模型，普林斯顿博士后 David Bourgin的项目
我自己对照实现了一遍，简单解读一下，关键的三个部分就是损失函数NLL、损失函数求导NLL_grad、梯度下降fit。这里使用的是$l2$正则化。

class LogisticRegression:
    def __init__(self, penalty="l2", gamma=0, fit_intercept=True):
        '''
        penalty 正则化系数
        gamma   正则化的权重
        '''
        err_msg = "penalty must be 'l1' or 'l2', but got: {}".format(penalty) 
        assert penalty in ["l2", "l1"], err_msg
        self.beta = None
        self.gamma = gamma
        self.penalty = penalty
        self.fit_intercept = fit_intercept

    def fit(self, X: np.ndarray, y, lr=0.01, tol=1e-7, max_iter=1e7):
        '''
        max_iter 最大的迭代轮数
        '''
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]

        l_prev = np.inf
        self.beta = np.random.rand(X.shape[1])  
        for _ in range(int(max_iter)):
            y_pred = sigmoid(np.dot(X, self.beta))  # 预测值，就是y为1的概率
            loss = self._NLL(X, y, y_pred)
            if l_prev - loss < tol:
                return
            l_prev = loss  # 更新损失，当损失值基本不再变动，就返回结果
            self.beta -= lr * self._NLL_grad(X, y, y_pred)

    def _NLL(self, X, y, y_pred):
        '''
        计算损失并添加惩罚
        '''
        N,M = X.shape
        order = 2 if self.penalty == "l2" else 1
        # 这里y的真实值就是 1，0，交叉损失函数结构如下
        nll = -np.log(y_pred[y==1]).sum() - np.log(1 - y_pred[y==0]).sum()
        penalty = 0.5 * self.gamma * np.linalg.norm(self.beta, ord=order) ** 2  # 求l2范数

        return (nll + penalty) / N

    def _NLL_grad(self, X: np.ndarray, y, y_pred):
        '''
        计算梯度
        '''
        N, M = X.shape
        p = self.penalty
        beta = self.beta
        gamma = self.gamma
        l1norm = lambda X: np.linalg.norm(X, 1)
        d_penalty = gamma * beta  if p == "l2" else gamma * l1norm(beta) * np.sign(beta)
        return -(np.dot(y-y_pred, X) + d_penalty) / N  # 这一步涉及矩阵求导

    def predict(self, X):
        if self.fit_intercept:
            X = np.c_[np.ones(X.shape[0]), X]
        return sigmoid(np.dot(X, self.beta))

欢迎一起来参与leetcode刷题项目

刷题的GitHub: github链接.