线性回归及其代码实现

前言

目前所流行的机器学习与深度学习算法，可以被总结为是一类利用数据驱动的算法。而回归算法可以说是最早的数据驱动的算法，其简洁有效，且最基础的线性回归算法具有及其稀有的显示解。虽然其模型复杂度不高，且对于数据有一定要求，但是对于满足要求的数据，线性回归具有如今一众大模型都达不到的效果（毕竟线性回归模型小，还有显示解），因此至今在工业界也被广泛使用，有一种大巧不工的美感。而其后来者，如在机器学习时代大名鼎鼎的SVM，其也是基于线性回归的，而现在所流行的深度学习，其中最重要的算法神经网络，其本质也是多个线性模型，而反过来，线性回归其实就是一个单层的神经网络。

综上，线性回归在机器学习领域具有举足轻重的地位。但是对于回归相关的内容，由于其始于统计，而现于机器学习学习中大放异彩，因此侧重统计与侧重机器学习的讲法会完全不同。权衡再三，本文将侧重机器学习讲解回归相关的内容，而侧重统计的内容待我另起炉灶。

注：本文会涉及到梯度下降算法，对此不了解的读者可以移步本站的梯度下降博客

多元线性回归

在高中数学课本中，我们就接触过最简单的一元线性回归，也就是$y = ax + b$，其常常与如下可视化图一同出现

而在实际问题中，往往并不止一个自变量$x$，这也就引出了多元线性回归，其数学方程为

$$y = w_0 + w_1 x_1 + \cdots + w_p x_p$$

其中$p$为自变量的数量。

损失函数

确定了模型，接下来就是如何计算出模型的参数，对于数据驱动的回归模型，那自然是借助于样本数据。这里我们记第$i$个样本为$(x_{i1}, x_{i2}, \cdots, x_{ip}, y_i)$，样本数量为$n$。

确定好符号后我们继续，我们的目的肯定是希望最终的模型越准确越好，那至少，这个模型在我们现有的样本数据上的表现要好吧。那现在，我们随便取一个数据样本$(x_{i1}, x_{i2}, \cdots, x_{ip}, y_i)$，再随便取一组模型参数，也就是$(w_0, w_1, w_2, \cdots, w_p)$，这个样本的真实因变量为$y_i$，模型的预测值为$\hat{y} = w_0 + w_1 x_{i1} + \cdots + w_p x_{ip}$，我们希望真实值与预测值之间的差距尽可能小，不过，差距有正负之分，而我们这里并不关注正负，只关注差距大小，因此再数学上我们可以使用绝对值或者平方来忽略掉正负，这里我们就选用平方了，具体原因后续再说。那也就是说，我们希望$(y_i - \hat{y})^2$要尽可能的小。

好，现在停下来看看我们本来在做什么，我们想要计算模型的参数，那么什么样的参数是好的呢，使$(y_i - \hat{y})^2$越小的参数越好。当然，$(y_i - \hat{y})^2$只是参考了一个数据样本，我们当然是希望每一个样本数据的$(y_i - \hat{y})^2$都要小，那我们把所有样本的$(y_i - \hat{y})^2$加起来不就好了，再将这个和记为一个函数，就是

$$L(\boldsymbol{w}) = \sum_{i=1}^{n} (y_i - \hat{y})^2 = \sum_{i=1}^{n} (y_i - w_0 - w_1 x_{i1} - \cdots - w_p x_{ip})^2 = ||\boldsymbol{y} - X\boldsymbol{w}||^2$$

其中$\boldsymbol{w} = (w_0, w_1, w_2, \cdots, w_p)^T$。不难看出这个函数里，样本数据都是给定的，那其自变量就只有$w_0, w_1, w_2, \cdots, w_p$了，为了简洁，这里我们将其写为矩阵形式。而$\boldsymbol{y}$与$X$分别为

$$\boldsymbol{y} = (y_1, y_2, \cdots, y_n)^T \\ X = \begin{bmatrix} 1&x_{11}&x_{12}&\ldots&x_{1p}\\ 1&x_{21}&x_{22}&\ldots&x_{2p}\\ 1&x_{31}&x_{32}&\ldots&x_{3p}\\ &&\cdots\\ 1&x_{n1}&x_{n2}&\ldots&x_{np} \end{bmatrix}$$

其实就是样本数据的自变量与因变量，而矩阵$X$的第一列$1$则是为了与$w_0$相乘。

注1：这里的矩阵写法只是为了简洁，本质上与前面的式子是一样的，对矩阵不熟悉的同学看不懂也没有关系，看懂前面的式子也是一样的。

注2：这里的$||\ \ ||$符号表示$L_2$范数，在高中数学中的名字是向量的模。这里有同学可能就会问了，向量的模的符号不是$|\ \ |$吗。这是因为在线性代数中，行列式的符号也是$|\ \ |$，为了避免歧义，范数就改用$||\ \ ||$符号了。

还有一个小点，在有些地方，$L_2$范数也会写成$||\ \ ||_2$，这里的下标$2$其实就是$L_2$中的$2$，同理$||\ \ ||_1$表示$L_1$范数，而没有下标一般默认为$L_2$范数，感兴趣的同学可以查$L_p$范数进一步了解。

最后还有范数的计算，对于向量$\boldsymbol{a} = (a_1, a_2, \cdots, a_p)$， $|| \boldsymbol{a} ||_2= \sqrt{a_1^2 + a_2^2 + \cdots + a_p^2} $， $|| \boldsymbol{a} ||_1= |a_1| + |a_2| + \cdots + |a_p| $

到这里，计算模型的参数的问题已经被我们转化为了求解函数$L(\boldsymbol{w})$的最小值问题。也就是说，多元线性回归模型的计算参数计算公式为

$$\hat{\boldsymbol{w}} = \mathop{min} \limits_{\boldsymbol{w}} L(\boldsymbol{w}) = \mathop{min} \limits_{\boldsymbol{w}} ||\boldsymbol{y} - X\boldsymbol{w}||^2$$

如$L(\boldsymbol{w})$这类函数，在机器学习中被称为损失函数，这里我们采用的该损失函数被称为平方误差损失函数(MSE)。

解析解与数值解

解析解

最后，就剩下一个问题，如何求解损失函数。这里我们需要先介绍一下解析解与数值解。解析解是指用解析表达式来表达的解，例如我们初中就接触过的二次函数$y = ax^2 + bx +c$，对其求导并令导数为$0$，即可得到其最值一定在$x = - \frac{b}{2a}$处取得，这就是解析解。而对于多元线性回归来说，其损失函数的导数为

$$\frac{ \partial L }{ \partial \boldsymbol{w} } = 2 X^{T} (X \boldsymbol{w} - y)$$

令损失函数的导数等于$0$，可以解得多元线性回归的解析解

$$2 X^{T} (X \boldsymbol{w} - y) = 0 \quad \Rightarrow \quad \hat{\boldsymbol{w}}= (X^T X)^{-1} X^T \boldsymbol{y}$$

数值解

而什么是数值解呢，就是只关注$x$的具体取值，而并不关注其表达式，通常通过一些优化算法搜索求得，比如大名鼎鼎的梯度下降，这是机器学习与深度学习中最常用的优化算法，例如函数$y = x^2$，这个函数的最小值在$x=0$处，而使用优化算法求解数值解可能会得到的结果可能为$x = 0.0001$或者$x = -0.0002$等等，你应该发现了这并不是最优解，没错，数值解通常会接近最优解，但是几乎不会达到最优解。

对比与区别

好，现在让我们来对比一下解析解与数值解。毫无疑问解析解是优秀的，因为解析解相当于一个通解，只要损失函数是确定的，最优解也就出来了，而使用计算机进行计算时这会节省大量的计算时间，在以大量计算著称的机器学习领域这无疑是巨大的优势。而相比之下，数值解相当于是一个特解，具体问题具体分析，每拿到一个新的数据都要算出具体的损失函数，然后用优化算法搜索最小值，这非常的浪费算力。

现在让我们来思考一下，既然解析解全都是优点，那为什么还要有数值解呢。首先，在纷繁复杂的损失函数中，解析解是非常稀有的，我们上面刚刚讲过用$(y_i - \hat{y})^2$表示真实值与预测值之间的差距，只要把其替换为$|y_i - \hat{y}|$，那么损失函数就会变为

$$L(\boldsymbol{w}) = \sum_{i=1}^{n} |y_i - \hat{y}| = \sum_{i=1}^{n} |y_i - w_0 - w_1 x_{i1} - \cdots - w_p x_{ip}| = ||\boldsymbol{y} - X\boldsymbol{w}||_1$$

这个损失函数马上就没有解析解了，这也就是为什么上面我们要用平方而不是绝对值。至少以人类现在的数学发展水平来说，很多函数都是没有解析解的。

其次，让我们观察一下上面给出的多元线性回归的解析解，不难看到其需要计算逆矩阵，相信学习过线性代数的读者大多跟笔者一样觉得计算逆矩阵是一件非常麻烦的事情，这对于计算机也是一样的。而在深度学习领域，由于模型的层层堆叠，即便存在解析解，多个庞大矩阵的逆矩阵计算也会将计算量拉到一个不可接受的程度。

现在我们来整理一下，解析解是最优雅的解决方案，但是其可遇不可求，并且在深度学习领域，随着模型的层层堆叠，其计算量可能会达到一个不可接受的程度。而数值解虽然有很多缺点，但是其不挑损失函数，且计算量相对可控，因此，在多个领域中，数值解都是非常泛用的存在。最后，之所以要在这里提到解析解与数值解，是因为我们上面所使用的损失函数是少数拥有解析解的函数，多数教材也是在线性回归处提到解析解，毕竟后续算法基本都是数值解。

岭回归与Lasso回归

损失函数

岭回归与Lasso回归是对多元线性回归的一种改进，其在损失函数中加入了正则化项，以此解决了多元线性回归可能存在的问题。所谓的正则化项其实就是参数$\boldsymbol{w}$的范数，$L_1$正则化项就是参数$\boldsymbol{w}$的$L_1$范数$||\boldsymbol{w}||_1$，而$L_2$正则化项是参数$\boldsymbol{w}$的$L_2$范数的平方$||\boldsymbol{w}||_2^2$。到这里，我们可以给出岭回归与Lasso回归的损失函数，岭回归的损失函数为

$$L(\boldsymbol{w}) = ||\boldsymbol{y} - X\boldsymbol{w}||_2^2 + \alpha ||\boldsymbol{w}||_2^2$$

Lasso回归的损失函数为

$$L(\boldsymbol{w}) = ||\boldsymbol{y} - X\boldsymbol{w}||_2^2 + \alpha ||\boldsymbol{w}||_1$$

其中的$\alpha$被称为正则化系数，我们后面再讲。

正则化项

我们先来解决一个问题，就是为什么要加入正则化项。如果是在统计学教材中，则一定会讲到解决多重共线性问题，这也是该算法的本意，且有严格的数学推导。不过既然是侧重与机器学习的文章，这里我就用机器学习的角度来解释一下正则项。

在机器学习中，正则项常常被用于控制模型复杂度，防止模型过拟合。简单来说，在机器学习中，模型的复杂度过大或者过小都不是我们所期望的，复杂度过大会导致过拟合，复杂度过小会导致欠拟合，而我们需要在这之间找到一个平衡点。而解决欠拟合的手段非常的单一，就是用更复杂的模型，上更多的参数。所以，如今一般的策略是先做到过拟合，再通过限制模型的手段将其复杂度往下降，例如决策树的剪枝就是一种防止过拟合的手段。

现在让我们来看看回归模型，我们会发现一个问题，回归模型的参数数量是固定的，那我们该如何降低复杂度呢。还有一种手段，就是限制参数的取值范围，这也可以达到降低复杂度的作用。而此时我们会发现，参数$\boldsymbol{w}$的范数正好就是一个很好的度量。以$L_2$范数为例，其计算公式为

$$|| \boldsymbol{w} ||_2= \sqrt{w_1^2 + w_2^2 + \cdots + w_p^2}$$

可以直观的看出，各参数的数值部分（也就是不管正负）越大，范数越大。而现在损失函数中加入了范数，我们希望损失函数尽可能的小，这就要求在兼顾模型的拟合效果（也就是将真实值与预测值之间的差距缩小）的同时，还需要控制模型的复杂度（就是范数不能太大）。这就需要在两者之间找到一个平衡点。而上述损失函数中的正则化系数$\alpha$就是用来控制正则化的强度的（也就是给我们调的），$\alpha$越大，模型复杂度的影响就越大，训练出的模型复杂度会更低，当然其精度也可能出现下降。

注：在统计视角的回归分析教材中，对于参数$\alpha$常常是使用如岭迹图等方法进行选择，而在机器学习中，选择参数$\alpha$的方法与其他机器学习算法并无二致，即对参数$\alpha$绘制参数曲线，选择分数最高的参数。

解析解

最后，上述两个损失函数其实也是存在解析解的，这里我们直接给出，岭回归的解析解为

$$\hat{\boldsymbol{w}} = (X^T X + \alpha I)^{-1} X^T \boldsymbol{y}$$

Lasso回归的解析解为

$$\hat{\boldsymbol{w}} = (X^T X)^{-1} (X^T \boldsymbol{y} - \frac{\alpha}{2} I)$$

logistic回归与softmax回归

logistic回归与softmax回归是使用多元线性回归模型进行分类算法，logistic回归用于二分类问题，而softmax回归用于多分类问题。

logistic回归

首先我们来看二分类问题，也就是logistic回归。在多元线性回归模型的基础上，要将其用于分类，我们需要解决两个问题。第一，回归模型的输出是一个连续值，而分类问题的输出是一个离散值，如何让回归模型输出离散值。第二，既然要解决的问题变成了分类问题，那原来的损失函数也就用不了了，应该使用何种损失函数。

sigmoid函数

那让我们先来解决第一个问题，如何让回归模型输出离散值。那我们先来看看，我们希望的输出值是什么样的，这里我们在处理的是二分类问题，用$0,1$分别代表两个类别，也就是我希望我的模型最后输出为$0$或$1$，这想法很自然。不过还能更进一步，那就是让模型的输出为一个范围在$(0,1)$之间的实数，而这个数可以被解释为样本属于类别$1$的概率，这会有很多好处，我们后面再讲。

回到让回归模型输出离散值这个问题，我们给出的答案是套娃，在原线性模型上再套一层函数。我们看，多元线性回归的模型是

$$f(\boldsymbol{x}) = w_0 + w_1 x_1 + \cdots + w_p x_p = \boldsymbol{w}^{T} \widetilde{\boldsymbol{x}}$$

其中$\boldsymbol{w} = (w_0, w_1, w_2, \cdots, w_p)^T$，$\boldsymbol{x} = (x_1, x_2, \cdots, x_p)^T$，$\widetilde{\boldsymbol{x}} = (1, x_1, x_2, \cdots, x_p)^T$。

我们需要认识到，数学上的函数与代码中的函数一样，都是有固定形状的输入与输出的。例如这个函数，其输入为$p$维向量，输出为一个实数，当然范围是$(-\infty,+\infty)$。而现在，我们希望输出为一个范围在$(0,1)$之间的实数，那我们只需要找到一个函数，其输入为$(-\infty,+\infty)$上的实数，而输出为$(0,1)$上的实数即可。满足这个要求的函数有很多，而logistic回归中使用的为sigmoid函数，函数表达式为

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

函数图像如下

可以看到，这个函数完美的符合我们的要求。那么，我们的logistic回归模型就可以写为

$$\hat{y} = F(\boldsymbol{x}) = \sigma(f(\boldsymbol{x})) = \frac{1}{1 + e^{-f(\boldsymbol{x})}} = \frac{1}{1 + e^{-\boldsymbol{w}^{T} \widetilde{\boldsymbol{x}}}}$$

那现在，我们来看看输出概率而非直接输出类别有何好处。其主要的好处就是使得我们的模型更精细，就比如说，在模型的实际使用中，有两个样本的logistic回归模型输出分别为0.7和0.8，那我们可以判断出输出为0.8的样本被归于类别1的概率是更大的，而如果直接输出类别，那这两个样本在我们看来就是一样的。并且，直接输出概率还可以自由的设定分类阈值，这也能更加满足实际应用中的需求，例如像银行贷款这种需要高置信度的场景，我们可能需要设定其概率达到0.7，0.8甚至更高，我们才会允许其进行贷款，或者说这个周期内概率最高的100人可以获得贷款等。在这种场景下，输出概率值会更加灵活。

注：虽然我们前面一直说概率值，但是严格来说logistic回归的输出并不是概率值，只是其数学性质与概率值基本相同，并且称其为概率值会更好理解，因此一般教材中也会直接称其为概率值，对此不必太钻牛角尖。

对数似然损失函数

然后我们来解决第二个问题，用什么损失函数。那我们会想，用原来的损失函数不行吗，仔细想想也不是说用不了啊。没错，确实是可以用，但是前人的研究告诉我们效果并不好，不过好用的损失函数前人也帮我们研究好了。

首先，我们前面说过logistic回归的输出$F(\boldsymbol{x})$代表样本属于类别1的概率，反之，$1 - F(\boldsymbol{x})$就代表样本属于类别0的概率，那么可以构造函数

$$p(\boldsymbol{x}) = [F(\boldsymbol{x})]^{y} [1 - F(\boldsymbol{x})]^{1 - y}$$

我们看，当样本的真实值为类别1，也就是$y = 1$时，$p(\boldsymbol{x}) = F(\boldsymbol{x})$，也就是样本被分别类别1的概率。反之，当样本的真实值为类别0，也就是$y = 0$时，$p(\boldsymbol{x}) = 1 - F(\boldsymbol{x})$，也就是样本被分别类别0的概率。综上，函数$p(\boldsymbol{x})$代表样本被分类到正确类别的概率。

将所有样本的$p(\boldsymbol{x})$相乘，也就是

$$\prod_{i=1}^n p(\boldsymbol{x_i}) = \prod_{i=1}^n [F(\boldsymbol{x_i})]^{y_i} [1 - F(\boldsymbol{x_i})]^{1 - y_i}$$

其中$\boldsymbol{x_i} = (x_{i1}, x_{i2}, \cdots, x_{ip})^T$。这就代表了将所有样本预测正确的概率，我们希望将其最大化，也就是说，这就是我们的损失函数了。为了方便运算，一般对其取对数再取负，那么我们最终的损失函数为

$$L(\boldsymbol{w}) = - \mathrm{ln} \left( \prod_{i=1}^n p(\boldsymbol{x_i}) \right) = - \mathrm{ln} \left( \prod_{i=1}^n [F(\boldsymbol{x_i})]^{y_i} [1 - F(\boldsymbol{x_i})]^{1 - y_i} \right) = \\ - \sum_{i=1}^{n} y_i \mathrm{ln} \left[ F(\boldsymbol{x_i}) \right] + (1 - y_i)\mathrm{ln} \left[ 1 - F(\boldsymbol{x_i}) \right] = - \sum_{i=1}^{n} y_i \mathrm{ln} \left[ \frac{1}{1 + e^{-\boldsymbol{w}^{T} \widetilde{\boldsymbol{x}}_i }} \right] + (1 - y_i)\mathrm{ln} \left[ \frac{e^{-\boldsymbol{w}^{T} \widetilde{\boldsymbol{x}}_i}}{1 + e^{-\boldsymbol{w}^{T} \widetilde{\boldsymbol{x}}_i}} \right]$$

也可以写作矩阵形式

$$L(\boldsymbol{w}) = - \sum_{i=1}^{n} y_i \mathrm{ln} \hat{y}_i + (1 - y_i)\mathrm{ln} (1 - \hat{y}_i) = - \boldsymbol{y}^T \ln \hat{\boldsymbol{y}} - (\boldsymbol{1} - \boldsymbol{y})^T \ln (\boldsymbol{1} - \hat{\boldsymbol{y}})$$

该损失函数被称为对数似然损失函数。

注：这里取对数是为了方便计算，而取负则是为了将最大化转换为最小化，因为在机器学习中损失函数默认都是最小化，如果你的损失函数要求最大化，那就加一个负号就可以了。

数值解与正则化

最后我们来做一些收尾工作，首先是损失函数的求解，在logistic回归中是没有解析解的，因此损失函数的求解依靠优化算法来实现，例如使用鼎鼎大名的梯度下降，这里我们不对优化算法做过多的介绍。

最后是正则化，logistic回归也可以通过正则化项来防止过拟合，这与上文的岭回归与Lasso回归是一样的，即

$$L(\boldsymbol{w}) = - \boldsymbol{y}^T \ln \hat{\boldsymbol{y}} - (\boldsymbol{1} - \boldsymbol{y})^T \ln (\boldsymbol{1} - \hat{\boldsymbol{y}}) + \alpha ||\boldsymbol{w}||_2^2 \\ L(\boldsymbol{w}) = - \boldsymbol{y}^T \ln \hat{\boldsymbol{y}} - (\boldsymbol{1} - \boldsymbol{y})^T \ln (\boldsymbol{1} - \hat{\boldsymbol{y}}) + \alpha ||\boldsymbol{w}||_1$$

softmax回归

然后我们来看多分类问题，处理多分类问题的方法有很多，这里我们会简单介绍一下常用的方法，但是主要介绍softmax回归。这是因为，多元线性可以看作一个处理回归问题的单层神经网络，而softmax回归则可以看作一个处理分类问题的单层神经网络，可以说softmax回归是深度学习的重要基础。

多分类学习策略

经过logistic回归的学习，我们可以发现logistic回归只能处理二分类问题，那么遇到多分类问题应该如何处理呢。答案就是将其处理为多个二分类问题，而其处理方式就被称为多分类学习策略。

注1：有同学可能会想，为什么不能把多个分类进行编码，比如把三个类编码为1，2，3，然后训练一个输出1，2，3其中之一的模型呢。这是因为，将多个分类进行编码一定会出现类别之间关系的紊乱。就比如说三个类编码为1，2，3，本来这三个类的地位应该是平起平坐的，但是这么一编码，就凭空给他们创造了一个数学上的关系，比如第二个类比第一个类大1，第三个类比第一个类大2。我们当然可以看得出来这个关系是错误的，但是模型没有那么聪明，模型是分辨不出来的，这必然会导致最终模型的效果不佳。

注2：对比决策树算法我们就会发现，决策树就没有这样的问题，这是因为决策树不需要把类别的编码用于计算，而线性回归是需要的（损失函数）。

OvO策略

OvO(One vs One)策略是指，将所有类别（共$n$类）两两组合，训练出$n(n-1)/2$个二分类器，最终结果由所有分类器投票决定。光说可能太抽象了，我们举个例子，比如说三个类别1，2，3，两两组合可以得到$(1,2),(1,3),(2,3)$，我们分别训练出三个分类器。比如说$(1,2)$，我们先把所有属于第三类的数据都筛掉，剩下的数据不就只有两个类了，我们就拿剩下的数据做logistic回归。就这样得到三个二分类器，现在拿到一个新样本，第一个分类器说他属于第一类，第二个分类器说他属于第一类，第三个分类器说他属于第三类，投票得到这个样本属于第一类。OvO策略就是这样一个流程。

OvR策略

OvR(One vs Rest)策略是指，将所有类别（共$n$类）各取一次作为正类，其余类别作为反类，训练出$n$个二分类器，最终取最高置信度的类别作为结果。听着是不是更抽象了，我们跟OvO对比着看，OvO训练出来的二分类器，比如说$(1,2)$，告诉我们的结果是这个样本属于第一类还是第二类，而OvR训练出来的分类器告诉我们的结果是这个样本属于第一类还是不属于第一类，或者说这个样本属于第一类还是属于其他类。还是拿三个类别1，2，3举例，三个类别分别做一次老大，比如说现在第一类做老大，跟OvO不一样，这里我们不用筛数据，而是说，现在的老大是第一类，我们给第一类编号为1（正类），剩下的第二类第三类统统编号为0（反类），然后做logistic回归。以此类推得到三个二分类器，经过上文的学习我们知道logistic回归会输出样本属于正类的概率，现在拿到一个新样本，第一个分类器输出0.9，第二个分类器输出0.3，第三个分类器输出0.4，第一个分类器的输出最高，也就是置信度最高，故这个样本属于第一类。OvR策略就是这样一个流程。下面我们用一个图再对比一下两种策略

注：从上述讲解中，我们可以发现，OvO策略的复杂度是比OvR策略要高的，但是在模型表现上，大多数情况下两者是差不多的，因此目前OvR策略可以说是主流。

MvM策略

当然在上述两种策略之上还有更为高级的MvM(Many vs Many)策略，就是将部分类别作为正类，剩余类别作为反类。可以看出OvO与OvR都是MvM的一种特殊情况。由于MvM策略较为复杂，因为如何进行划分是个大问题，目前也有多种MvM技术，因此这里只做简介，让大家知道有MvM策略这么个东西，而不去深究。

One-Hot编码

了解完多分类学习策略之后，我们正式开始讲softmax回归，这里可以先告诉大家，softmax回归可以被归类为使用OvR策略。首先，上文我们提到过对于多分类问题编码是一个大问题，不能用1，2，3这种编码。那既然一个数解决不了，那我们就用多个数不就好了，这也就是本节的主角，One-Hot编码(独热编码)。One-Hot编码使用一个$m \times 1$向量对类别进行编码，其中$m$为类别数。其编码规则为，对属于第$i$类的样本，则向量的第$i$行为$1$，其余行均为$0$。比如说三个类别，属于第一类的被编码为$(1, 0, 0)^T$，属于第二类的被编码为$(0, 1, 0)^T$，属于第三类的被编码为$(0, 0, 1)^T$。因为这种编码方式一个向量中只会有一个数为$1$，因此被称为One-Hot编码(独热编码)。

在改变编码方式之后，则会引出下一个问题，之前的模型输出是一个数，现在输出要改为一个向量，那模型结构连带着也会发生改变。我们先不管sigmoid函数，就看内层的模型

$$y = f(\boldsymbol{x}) = w_0 + w_1 x_1 + \cdots + w_p x_p = \boldsymbol{w}^{T} \widetilde{\boldsymbol{x}}$$

其中$\boldsymbol{w} = (w_0, w_1, w_2, \cdots, w_p)^T$，$\boldsymbol{x} = (x_1, x_2, \cdots, x_p)^T$，$\widetilde{\boldsymbol{x}} = (1, x_1, x_2, \cdots, x_p)^T$。

现在将改为

$$\boldsymbol{y} = f(\boldsymbol{x}) = W^{T} \widetilde{\boldsymbol{x}}$$

其中$\boldsymbol{y} = (y_1, y_2, \cdots, y_m)^T$，$\boldsymbol{x} = (x_1, x_2, \cdots, x_p)^T$，$\widetilde{\boldsymbol{x}} = (1, x_1, x_2, \cdots, x_p)^T$，而

$$W = \begin{bmatrix} w_{01}&w_{02}&w_{03}&\ldots&w_{0m}\\ w_{11}&w_{12}&w_{13}&\ldots&w_{1m}\\ w_{21}&w_{22}&w_{23}&\ldots&w_{2m}\\ &&\cdots\\ w_{p1}&w_{p2}&w_{p3}&\ldots&w_{pm}\\ \end{bmatrix}$$

softmax函数

内层的模型发生了变化，那logistic回归中使用的sigmoid函数也得有其替代品，这就是softmax函数。我们来回顾一下，原来的编码方式，输出被解释为属于正类的概率，因此我们用sigmoid函数把输出压到了$(0,1)$。那么One-Hot编码的解释方式是什么呢。答案是，输出向量的第$i$行代表样本属于第$i$类的概率，比如说$(0.8,0.05,0.15)$，代表样本属于第一类的概率是0.8，样本属于第二类的概率是0.05，样本属于第三类的概率是0.15。而由于是概率，这个向量需要满足各行之和为1，因此我们需要一个函数来完成这一需求，即softmax函数，其函数表达式为

$$\overline{\boldsymbol{y}} = \mathrm{softmax}(\boldsymbol{y}) , \quad \overline{y}_i = \frac{e^{y_i}}{ \displaystyle \sum_{j=1}^{m} e^{y_j} }$$

其中$\boldsymbol{y} = (y_1, y_2, \cdots, y_m)^T$，$\overline{\boldsymbol{y}} = (\overline{y}_1, \overline{y}_2, \cdots, \overline{y}_m)^T$。

可能会有同学看着有点懵，这是因为softmax函数是一个向量值函数，大家可能对此接触不多。其实很简单，就是说softmax函数输入是一个向量，输出也是一个向量。输入为向量$\boldsymbol{y}$，而输出为向量$\overline{\boldsymbol{y}}$，而$\overline{\boldsymbol{y}}$的每一个分量为

$$\overline{y}_i = \frac{e^{y_i}}{ \displaystyle \sum_{j=1}^{m} e^{y_j} }$$

我们也可以把$\overline{\boldsymbol{y}}$完全写开，就是

$$\overline{\boldsymbol{y}} = (\overline{y}_1, \overline{y}_2, \cdots, \overline{y}_m)^T = ( \frac{e^{y_1}}{ \displaystyle \sum_{j=1}^{m} e^{y_j} }, \frac{e^{y_2}}{ \displaystyle \sum_{j=1}^{m} e^{y_j} }, \cdots, \frac{e^{y_m}}{ \displaystyle \sum_{j=1}^{m} e^{y_j} })^T$$

这样应该比较清晰了。不难看出，softmax函数的想法是，先用指数函数将所有值均映射为正，然后再对数据进行概率归一化，以此达到目的。当然，softmax函数也可以写作矩阵形式

$$\mathrm{softmax}(\boldsymbol{y}) = \frac{ e^{\boldsymbol{y}} }{\boldsymbol{1}^T e^{\boldsymbol{y}}}$$

注：只有当输入为向量时，上述矩阵形式才是正确的

至此，我们就可以得到softmax回归的模型为

$$F(\boldsymbol{x}) = \mathrm{softmax}(f(\boldsymbol{x})) = \mathrm{softmax}(W^{T} \widetilde{\boldsymbol{x}})$$

交叉熵损失函数

最后就是解决损失函数了，其推导过程与logistic回归的对数似然损失函数是类似的，我们可以对比着来看。在上述推导中，我们构造了一个函数$p(\boldsymbol{x})$用于计算样本被分类到正确类别的概率，这里也是一样的。首先，我们来约定一下符号，模型的输出$F(\boldsymbol{x})$与真实值$\boldsymbol{y}$均为$m \times 1$向量，故这里我们约定用$F^{(j)}(\boldsymbol{x})$，$\boldsymbol{y}^{(j)}$表示向量的第$j$个分量，也就是第$j$行。那么，我们可以构造函数

$$p(\boldsymbol{x}) = \prod_{j=1}^m [F^{(j)}(\boldsymbol{x})]^{\boldsymbol{y}^{(j)}}$$

可以看出，当样本属于第$j$类时，$\boldsymbol{y}^{(j)} = 1$，而$\boldsymbol{y}$的其余分量均为$0$，因此$p(\boldsymbol{x}) = F^{(j)}(\boldsymbol{x})$，即函数$p(\boldsymbol{x})$代表样本被分类到正确类别的概率，完美完成任务。

接下来的流程与对数似然损失函数的推导是相同的，将所有样本的$p(\boldsymbol{x})$相乘，即

$$\prod_{i=1}^n p(\boldsymbol{x_i}) = \prod_{i=1}^n \prod_{j=1}^m [F^{(j)}(\boldsymbol{x}_i)]^{\boldsymbol{y}_{i}^{(j)}}$$

最后取对数再取负，得到损失函数为

$$L(W) = - \mathrm{ln} \left( \prod_{i=1}^n p(\boldsymbol{x_i}) \right) = - \mathrm{ln} \left( \prod_{i=1}^n \prod_{j=1}^m [F^{(j)}(\boldsymbol{x}_i)]^{\boldsymbol{y}_{i}^{(j)}} \right) = - \sum_{i=1}^{n} \sum_{j=1}^{m} {y}_{i}^{(j)} \mathrm{ln} [F^{(j)}(\boldsymbol{x}_i)]$$

也可以写作矩阵形式

$$L(W) = - \sum_{i=1}^{n} \sum_{j=1}^{m} {y}_{i}^{(j)} \ln \hat{y}_{i}^{(j)} = - \sum_{i=1}^{n} \boldsymbol{y}_{i} \ln \hat{\boldsymbol{y}}_{i} = - \mathrm{tr} (Y^T \ln \hat{Y})$$

这个损失函数被称为交叉熵损失函数，与对数似然函数相同，交叉熵损失函数也是使用优化算法进行求解，并且也可以加入正则化项。

注：交叉熵是一个信息论中的概念，该损失函数的推导也可从信息论的角度出发。

softmax回归与OvR策略

刚才我们有提到，softmax回归可以被归类为使用OvR策略，这里我们展开讲讲。首先，我们来思考一个问题，One-Hot编码的每一个分量，以及softmax回归的输出的每一个分量代表什么，比如说样本真实值$(1, 0, 0)^T$，softmax回归输出的预测值为$(0.9, 0.07, 0.03)^T$，放在一起看我们会说样本属于第一类的概率是0.9，样本属于第二类的概率是0.07，样本属于第三类的概率是0.03，那单独看0.9这一个分量呢。这个0.9代表样本属于第一类的概率是0.9，属于其他类的概率为0.1。其他类，这与OvR策略的思想不谋而合，我们再来想想，如果一个三分类问题用OvR策略做，那对于一个属于第一类的样本，第一个模型的真实值是不是也应该是1，后面两个模型的真实值是不是也应该是0。到这里我们就可以发现，softmax回归本质上就是把OvR策略构建的多个二分类模型合并成了一个用矩阵表示的模型，我们这里反过来，把softmax回归的模型写开，上文说过softmax回归模型为

$$\boldsymbol{y} = f(\boldsymbol{x}) = W^{T} \widetilde{\boldsymbol{x}}$$

把矩阵全部展开就是

$$\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ \cdots\\ y_{m} \end{bmatrix} = \begin{bmatrix} w_{01}&w_{11}&w_{21}&\ldots&w_{p1}\\ w_{02}&w_{12}&w_{22}&\ldots&w_{p2}\\ w_{03}&w_{13}&w_{23}&\ldots&w_{p3}\\ &&\cdots\\ w_{0m}&w_{1m}&w_{2m}&\ldots&w_{pm}\\ \end{bmatrix} \begin{bmatrix} 1\\ x_{1}\\ x_{2}\\ \cdots\\ x_{p} \end{bmatrix}$$

把$\boldsymbol{y}$的每个分量全部算出来，写成方程组为

$$\begin{cases} y_1 = w_{01} + w_{11} x_1 + w_{21} x_2 + \cdots + w_{p1} x_p \\ y_2 = w_{02} + w_{12} x_1 + w_{22} x_2 + \cdots + w_{p2} x_p \\ y_3 = w_{03} + w_{13} x_1 + w_{23} x_2 + \cdots + w_{p3} x_p \\ \qquad \qquad \qquad \qquad \cdots\\ y_m = w_{0m} + w_{1m} x_1 + w_{2m} x_2 + \cdots + w_{pm} x_p \\ \end{cases}$$

我们可以发现整个softmax回归的内层模型被拆分成了$m$个模型，这与OvR策略构建的$m$个模型是相同的。至此，softmax回归与OvR策略的关系就清楚了，两者本质是相同的，其唯一区别就是最后一步压缩成概率的时候，softmax回归使用softmax函数，而OvR则是用sigmoid函数对单独处理，下面我们用两张图对比一下两者的区别

sklearn代码实现

最后我们来介绍一下线性回归的sklearn代码实现，运行环境为jupyter notebook

注1：sklearn是机器学习中常用的python库，其封装了大多数常见的机器学习算法，使得复杂的机器学习算法可以在短短几行代码中被实现，俗称“调包”（bushi。下面的内容会建立在对sklearn有简单了解的基础上，没有接触过sklearn的读者可移步sklearn快速入门作简单了解。

注2：以下内容会介绍实现算法的函数及其重要参数，并且用sklearn中自带的小型数据集跑一个简单的demo。主要目的是将sklearn中的参数跟理论部分做一个对照，知道理论部分中讲到的参数在sklearn中是哪个，怎么调。还有就是简单把这个算法跑一遍，不过由于机器学习是一门以实践为主的学科，我们用sklearn中自带的小型数据集做的demo的作用也就仅此而已了，其对实际工程中的调参没有任何参考意义，仅仅是告诉你这个算法怎么跑起来。

注3：本文并不会把函数的细节讲得面面俱到，甚至参数都只是调重要或常用的讲，想要详细了解的可以参考sklearn官方文档或者sklearn中文文档。

我们先做一些前期工作，导入依赖库

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from sklearn.linear_model import LinearRegression  # 线性回归
from sklearn.linear_model import Ridge  # 岭回归
from sklearn.linear_model import Lasso  # Lasso回归
from sklearn.linear_model import LogisticRegression  # Logistic回归 

from sklearn.datasets import load_wine  # 红酒数据集
from sklearn.datasets import load_diabetes  # 糖尿病数据集

from sklearn.model_selection import train_test_split  # 切分数据
from sklearn.model_selection import cross_val_score  # 交叉验证

import warnings
warnings.filterwarnings("ignore")

我们使用红酒数据集与糖尿病数据集作为示例数据，这里我们使用sklearn自带的数据模块，读取数据代码为

# 读取数据
wine = load_wine()  # 红酒数据集，用于分类
diabetes = load_diabetes()  # 糖尿病数据集，用于回归

# 提取特征与标签（红酒数据集）
data_c = pd.DataFrame(wine.data)
data_c.columns = wine.feature_names
target_c = pd.DataFrame(wine.target)
target_c.columns = ['class']

# 提取特征与标签（糖尿病数据集）
data_r = pd.DataFrame(diabetes.data)
data_r.columns = diabetes.feature_names
target_r = pd.DataFrame(diabetes.target)
target_r.columns = ['Y']

好，做完准备工作我们正式开始讲算法实现

多元线性回归

首先是多元线性回归，该算法用于解决回归问题，因此我们使用糖尿病数据集作为示例数据集，对糖尿病数据集进行训练集与测试集的划分

# 划分训练集与测试集
Xtrain, Xtest, Ytrain, Ytest = train_test_split(data_r, target_r, test_size=0.2)

在sklearn中，实现多元线性回归的函数为LinearRegression()，由于模型简单且sklearn中使用解析解，所以其参数非常少，只有四个参数，并且大多数时候全默认就可以，参数列表如下

参数	说明
fit_intercept	是否计算此模型的截距。默认为True，如果设置为False，则在计算中将不使用截距（即，数据应中心化）。
normalize	默认为False，如果为True，则在回归之前通过减去均值并除以L2-范数来对回归变量X进行归一化。
copy_X	默认为True，将复制X，否则X可能会被覆盖。
n_jobs	用于计算的核心数。

可以看到，线性回归的参数简单，而且是真的没什么可以调的，下面是sklearn经典建模步骤

# 线性回归
clf = LinearRegression()
clf.fit(Xtrain, Ytrain)
clf.score(Xtest, Ytest)

如果想要获知确切的参数，线性回归算法还有两个属性可以用

clf.coef_  # 特征系数
clf.intercept_  # 截距

岭回归与Lasso回归

在sklearn中，实现岭回归与Lasso回归的函数为Ridge()与Lasso()，其重要的参数只有一个，就是正则化系数alpha

# 岭回归
clf = Ridge(alpha=0.02)
clf.fit(Xtrain, Ytrain)
clf.score(Xtest, Ytest)

# Lasso回归
clf = Lasso(alpha=0.02)
clf.fit(Xtrain, Ytrain)
clf.score(Xtest, Ytest)

而该参数的调整方法与众多机器学习算法的调整方法是一样的，就是暴力搜索，找到效果最好的，一般会画个超参数图像看看模型效果在不同参数下是怎么样的，下面是绘制超参数图像的兼用卡代码

# 调参
test = []
for i in range(100):
    clf = Lasso(alpha=i*0.01)
    score = cross_val_score(clf, Xtrain, Ytrain, cv=10).mean()  # K折交叉验证
    test.append(score)

plt.plot(np.arange(0.01, 1.01, 0.01),test,color="red",label="alpha")
plt.legend()

test_max = max(test)
print("最高准确率：", test_max)
print("最佳参数", test.index(test_max)*0.02)

注：如果是一些调整范围比较大的参数，例如需要从1到200进行暴力搜索，为了节省算力与时间一般会使用变步长搜索。很简单，就是先设置一个比较大的步长，比如5，然后假如搜索出来最佳参数是45，那我们就缩小范围，例如30到60，再缩小步长，例如1，做更高精度的搜索，这就是所谓“变步长”，一般来个一次也就差不多了。

logistic回归与softmax回归

最后是logistic回归与softmax回归，这两个算法均用于解决分类问题，因此我们使用红酒数据集作为示例数据集，对红酒数据集划分训练集与测试集

# 划分训练集与测试集
Xtrain, Xtest, Ytrain, Ytest = train_test_split(data_c, target_c, test_size=0.2)

在sklearn中，实现logistic回归与softmax回归的函数为LogisticRegression()，这个函数参数就比较多了，我们分类来看

正则化

参数	说明
penalty	正则化方法，可填{‘L1’, ‘L2’, ‘elasticnet’, ‘none’}，分别代表L1正则化，L2正则化，同时使用L1正则化与L2正则化，以及不使用正则化。默认为使用L2正则化。
C	正则化系数的倒数，注意是倒数，也就是数值越小，正则化强度越高。

多分类

参数	说明
multi_class	多分类策略，可填{‘auto’, ‘ovr’, ‘multinomial’}，分别表示自动选择，OvR策略与MvM策略。

优化算法

参数	说明
solver	优化算法，可填{‘newton-cg’, ‘lbfgs’, ‘liblinear’, ‘sag’, ‘saga’}，默认为‘lbfgs’。 - 对于小型数据集，“ liblinear”是一个不错的选择，而对于大型数据集，“ sag”和“ saga”更快。 - 对于多类分类问题，只有’newton-cg’ ，‘sag’，‘saga’ 和 ‘lbfgs’ 可以使用 ‘multinomial’。“ liblinear”仅能用‘ovr’。 - ‘newton-cg’，‘lbfgs’，'sag’和’saga’只能处理L2正则化 - 'liblinear’和’saga’也可以处理L1正则化 - ‘saga’还支持“ elasticnet”正则化
max_iter	优化算法的最大迭代轮次，可以用参数曲线调出效果最好的迭代轮次。
random_state	随机数种子，用于确保优化算法每次运行结果一致。

虽然参数挺多的，但是需要调的并不是很多，大多数参数还是根据数据集的情况做选择就好了，真正需要调的也就是正则化系数，调整方法跟上面的线性回归一样，照例放一组兼用卡代码

# Logistic回归
clf = LogisticRegression(multi_class='ovr', penalty='l1', solver='saga', C=1.0, random_state=1008, max_iter=10000)
clf.fit(Xtrain, Ytrain)
clf.score(Xtest, Ytest)

补充：从零实现线性回归

随着机器学习与深度学习的发展，许多常用的机器学习算法都有非常好用的依赖库可以进行调用，这对于我们运用机器学习算法解决实际问题自然是十分方便的，但是这也导致长时间以来，对于在学习机器学习算法时是否需要从零实现算法这个问题，一直有着两派观点。一派认为学习机器学习算法只需要理解原理即可，而另一派则认为从零实现算法有助于更好的理解算法本身，因此，笔者选择将这部分内容放入补充，以上。

准备工作

虽然说我们要从零实现算法，但这并不代表我们不会使用任何的第三方库，我认为实现算法就是一个将数学语言翻译为计算机编程语言的过程，在上文中，我们应该可以感受到线性回归中存在大量的矩阵运算，而python对矩阵运算的原生支持是很差的，因此，我们会大量依赖numpy这个第三方库。

注：如R语言，MATLAB等对矩阵运算的支持就非常好而且自然，写起来非常舒服，而python毕竟还是计算机语言，即使用上了numpy也有一些在学数学的人看来很奇怪的地方，在后续实现中我会提到，不过这都不是什么大问题。

numpy是一个科学计算库，其对矩阵运算有着很好的支持，对于从零实现算法是合适的选择，同时我们还会使用sklearn库中的数据生成器，生成数据用于验证算法是否正确实现。下面是导入依赖库的代码，运行环境为jupyter notebook

import numpy as np

# 生成数据集
from sklearn.datasets import make_regression  # 回归数据
from sklearn.datasets import make_classification  # 分类数据

from sklearn.model_selection import train_test_split  # 划分数据集

注：对numpy不熟悉的读者可以移步NumPy官方快速入门，也有对应的NumPy官方快速入门中文版。在后文中，所有实现过程都将使用numpy中的array数据类型。

梯度下降

然后是梯度下降的实现，在本站的梯度下降博客中已经给出了小批量梯度下降的numpy实现，这里就直接挪用，代码如下

# 线性回归训练
def LR_train(features, labels, criterion, epochs_num, batch_size=10, alpha=0.01):
    # features：特征矩阵，labels：标签向量，criterion：损失函数，epochs_num：数据集遍历次数，batch_size：批量大小，alpha：学习率
    
    # 小批量数据读取
    def data_iter(batch_size, features, labels):
        # batch_size：批量大小，features：特征，labels：标签
        
        data_size = np.size(features, 0)  # 数据集数量
        index = np.random.permutation(data_size)  # 生成随机索引
        
        for i in range(0, data_size, batch_size):
            batch_index = index[i: min(i + batch_size, data_size)]  # batch的索引列表
            X, y = features[batch_index], labels[batch_index]
            X, y = X.reshape(np.size(batch_index), -1), y.reshape(np.size(batch_index), -1)  # 固定形状
            yield X, y
    
    # 梯度下降
    def gradient_descent(w, gradient, alpha):
        # w：更新的参数，gradient：梯度
        w = w - alpha * gradient  # 更新参数
        return w
    
    w = np.random.normal(loc=0, scale=0.1, size=(np.size(features, 1)+1, 1))  # 生成初始化参数
    
    for epoch in range(epochs_num):
        for X, y in data_iter(batch_size=batch_size, features=features, labels=labels):
            gradient = criterion.gradient(w, X, y)  # 计算梯度
            w = gradient_descent(w, gradient, alpha)  # 梯度下降
            
        train_l = criterion(w, features, labels)  # 本轮epoch的训练误差
        print(f'epoch {epoch + 1}, loss {float(train_l):f}')
    
    return w

当然也进行了一些小改动，由于本章节后续需要实现多个算法，每一个算法都是一个不同的损失函数，因此将损失函数改为了参数。

多元线性回归numpy实现

数据集生成

使用sklearn.datasets中的make_regression函数生成用于回归问题的样本数据，生成数据的代码如下

bias = 8
X, y, coef = make_regression(
    n_samples=1008,  # 样本数量
    n_features=3,  # 特征维度
    n_informative=3,  # 有效特征
    n_targets=1,  # 输出维度
    bias=bias,  # 偏差
    coef=True,  # 返回权重
    noise=1,  # 噪声标准差
    random_state=1008  # 随机数种子
)

具体各参数的含义已经写在注释中了，这里重点说几个参数，首先是bias与coef=True，这两个吃饭时代表了该数据集的“正确答案”，后续实现的算法计算得到的权重越接近返回的coef，偏置越接近bias，则说明其效果越好。然后是n_informative有效特征数，其代表参与计算y的特征数量，如果其小于特征数n_features，则代表有特征与y无关，是无效特征，其权重的“正确答案”应该为0，这里为简单起见，设定所有特征均有效。

解析解实现

现在，我们正式进入多元线性回归的实现，我们先从解析解开始。解析解的实现及其的简单，简单到其核心只有一行代码，就是

w = np.linalg.inv(X.T @ X) @ X.T @ y

回顾上述所讲，虽然我们说求解线性回归等价于求解损失函数的最小值，但是解析解已经把这个最小值表示出来了，这代表着我们根本就不需要管什么损失函数，只需要把解析解这个式子算出来就行了，就这么简单。线性回归的解析解数学上写为

$$\hat{\boldsymbol{w}} = (X^T X)^{-1} X^T \boldsymbol{y}$$

numpy中可以用运算符@实现矩阵乘法，用.T代表转置，逆矩阵可以使用np.linalg.inv()函数实现，现在我们就可以把上面的式子翻译成代码了，也就是刚才那一行核心代码，将其封装成函数。需要注意，上述公式中的$X$矩阵是做了增广的，而输入的样本矩阵是没有经过增广的，因此，这里需要对$X$矩阵做增广

# 线性回归解析解计算
def LR_compute(X, y):
    X = np.hstack((np.ones((np.size(X, 0), 1)), X))  # 增广X矩阵
    
    w = np.linalg.inv(X.T @ X) @ X.T @ y
    w = w.reshape(-1, 1)
    return w

再写一个函数用于预测，核心就是计算$\boldsymbol{y} = X\boldsymbol{w}$，完整代码如下

def LR_predict(X, w):
    X = X.reshape(-1, len(w)-1)  # 保证形状
    X = np.hstack((np.ones((np.size(X, 0), 1)), X))  # 增广X矩阵
    
    y_hat = X @ w
    return y_hat

注：在线性代数中，向量与矩阵都是二维的，但是在numpy中常常会使用一维的array来储存向量，例如上述生成数据集中生成的y就是一个一维array。在使用numpy实现算法的过程中，一维的array就像一个幽灵，容易在不经意间让你的代码报错。

因此，我们会常常使用.reshape方法来保证array有正确的形状。例如上述代码中的X = X.reshape(-1, len(w)-1)，其中-1代表由程序指定，例如X本身的形状是(504, 2)，那就会被转化为(1, 1008)，如果X的形状是(1008,)，那就会被转化为(1, 1008)，这也是reshape方法常用的写法。

数值解与正则化实现

损失函数及其梯度

然后来到了多元线性回归的数值解，梯度下降上面已经实现好了，现在需要做的就是两件事，第一，实现损失函数的计算，第二，实现损失函数梯度的计算。对于多元线性回归来说，损失函数为

$$L(\boldsymbol{w}) = ||\boldsymbol{y} - X\boldsymbol{w}||^2 = (\boldsymbol{y} - X\boldsymbol{w})^{T} (\boldsymbol{y} - X\boldsymbol{w})$$

损失函数的梯度为

$$\frac{ \partial L }{ \partial \boldsymbol{w} } = 2 X^{T} (X \boldsymbol{w} - y)$$

将上述两个计算公式翻译为代码，我们就得到了两行核心代码

# 计算损失函数
l = (y - X @ w).T @ (y - X @ w)
# 计算损失函数梯度
lg = 2 * (X.T @ (X @ w - y))

正则化

由于岭回归与Lasso回归，即正则化的实现，就是给损失函数添加正则项，故在此一并做实现。正则化的实现其实也非常简单，由于正则项是加在损失函数后面的，因此，其梯度也是直接加在损失函数的梯度上，故我们只需要实现两种正则项以及其梯度的计算，再根据需要加在损失值与梯度值上即可。$L_1$正则项与$L_2$正则项的计算公式分别为

$$|| \boldsymbol{w} ||_1 = |w_1| + |w_2| + \cdots + |w_p|$$

$$|| \boldsymbol{w} ||_2^2 = w_1^2 + w_2^2 + \cdots + w_p^2 = \boldsymbol{w}^T \boldsymbol{w}$$

可以看到$L_1$正则项是将参数向量全部取绝对值相加，而$L_2$正则项是全部取平方再相加，当然也可以使用向量形式进行计算，不仅如此numpy中也提供了linalg.norm()函数用于计算范数，因此，实现$L_1$正则项与$L_2$正则项的方法非常多，下面给出几种

# L1正则项
L1 = np.sum(np.abs(w))  # 取绝对值再相加
L1 = np.linalg.norm(w, ord=1)  # 计算L1范数
# L2正则项
L2 = np.sum(w ** 2)  # 取平方再相加
L2 = w.T @ w  # 向量形式计算
L2 = np.linalg.norm(w, ord=2) ** 2  # 计算L2范数的平方

而$L_1$正则项与$L_2$正则项的梯度分别为

$$\frac{ \partial || \boldsymbol{w} ||_1 }{ \partial \boldsymbol{w} } = \mathrm{sign}(\boldsymbol{w})$$

$$\frac{ \partial || \boldsymbol{w} ||_2^2 }{ \partial \boldsymbol{w} } = 2\boldsymbol{w}$$

翻译为代码则为

# L1正则项梯度
L1_g = np.sign(w)
# L2正则项梯度
L2_g = 2 * w

封装

最后，我们选择将上面的成果封装为类，代码有点长，我们先看完整代码再来讲，完整代码如下

class MSELoss():
    def __init__(self, regularization='none', weight_decay=1):
        # regularization：正则化类型, weight_decay：正则化系数
        
        # 判断regularization值的合法性
        valid_regularizations = ['none', 'L1', 'L2']
        if regularization not in valid_regularizations:
            raise ValueError(f"Invalid regularization type. Expected one of {valid_regularizations}, but got '{regularization}'")
        
        self.regularization = regularization
        self.weight_decay = weight_decay
    
    # 计算损失值
    def __call__(self, w, X, y):
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        
        # 确保w, y为array类型，并控制其形状
        w = np.array(w)
        w = w.reshape(-1, 1)  # w:(num_features+1, 1)
        y = np.array(y)
        y = y.reshape(-1, 1)
        
        l = (y - X @ w).T @ (y - X @ w)
        l = l / np.size(X, 0) # 损失值除以样本数
        
        # 正则化
        if self.regularization == 'L1':
            coef = np.delete(w, 0, axis=0)  # 偏置不参与正则化
            l = l + self.weight_decay * np.linalg.norm(coef, ord=1)
        elif self.regularization == 'L2':
            coef = np.delete(w, 0, axis=0)
            l = l + self.weight_decay * (np.linalg.norm(coef, ord=2) ** 2)
        
        return l
    
    # 计算梯度值
    def gradient(self, w, X, y):
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        
        # 确保w, y为array类型，并控制其形状
        w = np.array(w)
        w = w.reshape(-1, 1)  # w:(num_features+1, 1)
        y = np.array(y)
        y = y.reshape(-1, 1)
        
        lg = 2 * (X.T @ (X @ w - y))
        lg = lg / np.size(X, 0)  # 梯度值除以样本数
        
        # 正则化
        if self.regularization == 'L1':
            coef = w
            coef[0] = 0
            lg = lg + self.weight_decay * np.sign(coef)
        if self.regularization == 'L2':
            coef = w
            coef[0] = 0
            lg = lg + self.weight_decay * 2 * coef
        
        return lg

首先，由于考虑了正则化，我们需要定义变量来控制是否使用正则化，以及正则化系数是多少，这部分在类的初始化中实现。然后，计算损失值与计算梯度被写成了类的两个方法，计算损失值的部分使用了__call__方法，这是python类中的一个魔法方法，其作用是让类实例可以像调用函数一样被调用，计算梯度的部分则被写入了gradient方法。接着，在计算损失值与梯度前的代码，实现了对$X$矩阵的增广，以及保证$\boldsymbol{w}$与$y$的形状正确。

最后，是正则化的实现部分，这里有两个需要注意的点，第一，损失值与梯度不包含正则项的部分，需要先除以样本数，再与正则项相加，以损失值的计算为例，使用$L_2$正则化，其损失函数为

$$L(\boldsymbol{w}) = ||\boldsymbol{y} - X\boldsymbol{w}||^2 + \alpha ||\boldsymbol{w}||^2$$

不难发现，前面一项是与样本量有关的，因为$\boldsymbol{y}$跟$X$都是跟样本有关的嘛，样本量越大其数值越大，而正则项部分是与样本量无关的，因此，将前面一项除以样本量可以保证训练时调整batch_size参数对正则化强度无影响，也就是

$$L(\boldsymbol{w}) = \frac{1}{n} ||\boldsymbol{y} - X\boldsymbol{w}||^2 + \alpha ||\boldsymbol{w}||^2$$

如果不这样处理，那么在训练调参时，将batch_size调大，则前一项会变大，这其实相当于将weight_decay调小（在损失函数中，数值越大就越有话语权，想想为什么正则化系数越大，正则化强度越高，应该不难理解吧）。

第二，偏置项$w_0$不参与正则化，因此需要对$\boldsymbol{w}$做一些处理，在计算损失值部分，由于正则项是一个标量，直接删除$w_0$再做正则项的计算即可，而在计算梯度部分，由于梯度是一个向量，故将$w_0$改为$0$再计算正则项梯度，这样可以保证其计算得到的梯度也为$0$，相加后不会对原梯度造成影响。

解析解与数值解测试

最后，我们需要测试一下上述实现代码，由于回归问题是有所谓的”正确答案“的，对于两种求解方式，我们的测试方法就是看求解出来的参数结果是否接近”正确答案“，这里我们先把”正确答案“定义出来，以便后续使用

w_true = np.insert(coef, 0, bias).reshape(-1, 1)  # ”正确答案“

先从解析解开始测试，使用上述实现的函数进行求解出w，并与w_true相减

w = LR_compute(X, y)
w - w_true

输出结果为

array([[ 0.00714635],
       [-0.04096568],
       [ 0.03152433],
       [ 0.01598949]])

然后是数值解，代码与输出结果如下

criterion = MSELoss()
w = LR_train(X, y, criterion, epochs_num=10, batch_size=64, alpha=0.05)

epoch 1, loss 474.756124
epoch 2, loss 16.673691
epoch 3, loss 1.520759
epoch 4, loss 0.959252
epoch 5, loss 0.949348
epoch 6, loss 0.946817
epoch 7, loss 0.946365
epoch 8, loss 0.947011
epoch 9, loss 0.946535
epoch 10, loss 0.946679

w - w_true

array([[-0.00812155],
       [-0.03824893],
       [ 0.02039796],
       [ 0.02474028]])

可以看出两种求解方式的效果都非常理想。

正则化测试

然后是正则化的测试，上文就提到，正则化是用于防止过拟合的一项技术，因此，我们先将数据集变得复杂一些，主要是加入一些无效特征，使得数据集更容易过拟合，然后将数据集划分为训练集与测试集

X, y, coef = make_regression(
    n_samples=1008,  # 样本数量
    n_features=30,  # 特征维度
    n_informative=3,  # 有效特征
    n_targets=1,  # 输出维度
    bias=8,  # 偏差
    coef=True,  # 返回权重
    noise=10,  # 噪声标准差
    random_state=1008  # 随机数种子
)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=1008)

使用多元线性回归与岭回归，在训练集上分别计算出结果

criterion = MSELoss()
w = LR_train(X_train, y_train, criterion, epochs_num=10, batch_size=64, alpha=0.05)

epoch 1, loss 366.210005
epoch 2, loss 121.428784
epoch 3, loss 99.172641
epoch 4, loss 97.034523
epoch 5, loss 96.523046
epoch 6, loss 96.505926
epoch 7, loss 96.980217
epoch 8, loss 96.564764
epoch 9, loss 96.713952
epoch 10, loss 96.512866

criterion_r = MSELoss(regularization='L2', weight_decay=0.01)
w_r = LR_train(X_train, y_train, criterion_r, epochs_num=10, batch_size=64, alpha=0.05)

epoch 1, loss 401.374109
epoch 2, loss 155.507344
epoch 3, loss 135.353584
epoch 4, loss 133.308615
epoch 5, loss 132.809783
epoch 6, loss 132.812304
epoch 7, loss 132.908994
epoch 8, loss 132.649182
epoch 9, loss 132.956535
epoch 10, loss 132.915263

再分别计算两者在测试集上的损失值

criterion(w, X_test, y_test)  # 输出为array([[145.26178254]])

criterion(w_r, X_test, y_test)  # 输出为array([[139.72114463]])

可以看到，在复杂的数据集上，使用了正则化的模型，在测试集上的表现要优于不使用正则化的模型。

logistic回归numpy实现

数据集生成

然后，我们来到logistic回归的实现，作为处理分类问题的模型，我们自然不能继续使用上面的数据集了。使用sklearn.datasets中的make_classification函数生成用于分类问题的样本数据，生成数据的代码如下

X, y = make_classification(
    n_samples=1008,
    n_features=3,  # 特征维度
    n_informative=3,  # 所有特征均与标签相关
    n_redundant=0,  # 无冗余特征
    n_classes=2,  # 类别
    class_sep=1,  # 样本分离度
    flip_y=0.1,  # 标签噪声
    random_state=1008  # 随机数种子
)

这里重点说两个参数，首先是class_sep样本分离度，由于生成分类数据的原理是，先生成几个类别中心，再在类别中心周围生成数据，而class_sep代表了类别中心的分离度，class_sep越大不同的类别中心距离越远，则数据的划分会变得更容易。然后是flip_y标签噪声，其代表有一定比例的样本会被随机赋予标签，同理，噪声越大，数据集越难。

sigmoid函数实现

搞定了数据集，按理说可以进入模型的实现了，但是在此之前，我们需要先实现sigmoid函数，后面的实现过程会多次用上，sigmoid函数的表达式为

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

翻译为代码就是

1 / (1 + np.exp(-x))

但是事情并没有那么简单，这涉及到计算机中的算术溢出问题，这里我们不详细解释，仅作简单说明。

算术溢出问题

简单来说，计算机中的所有变量都是有表示范围的，如果学习过如C语言这类比较传统的计算机语言的读者，对此应该并不陌生，例如32位有符号整数的表示范围为-2147483648至2147483647，当变量的数值超过其表示范围时，就会出现各种各样的bug，例如著名的2147483647 + 1 = -2147483648。

而浮点数的情况会更复杂一点，例如，深度学习常用32位有符号浮点数，其表示范围约为$1.17 \times 10^{-38}$至$3.4 \times 10^{38}$，以及$-3.4 \times 10^{38}$至$-1.17 \times 10^{-38}$，不难发现，浮点数的表示范围在$0$附近有一片空白，无论正负，如果浮点数的数值过于靠近$0$，也是会超出表示范围的。

最后，在python中，情况又会有一些不同，毕竟python在设计之初就是面向不熟悉计算机底层的使用者，其设计时也考虑了算术溢出问题。python的处理是对发生了算术溢出的变量赋予零或无穷，python中使用符号inf代表数学中的无穷符号，当变量数值超过表示范围的最大值时，变量会被赋值inf，同理，变量数值小于表示范围时则赋值-inf，而当变量数值落入$0$附近无法表示的范围时，则赋值0。

现在，让我们回到sigmoid函数，不难发现，由于指数函数的存在，不需要很大的x就会引发算术溢出的问题，例如，np.exp(1000) = inf，np.exp(-1000) = 0。而解决办法也很简单，就是截断，例如设定输入范围为$[-15, 15]$，则超过这个范围的输入均截断为$15$或$-15$，并且截断之后的函数值与真实函数值不会有过大差距，误差是可以接受的，具体代码实现为

def fun_sigmoid(x, clip_threshold=15):
    # clip_threshold: 输入上下限
    x = np.clip(x, -clip_threshold, clip_threshold)  # 对输入进行截断，防止输出为0或1
    return 1 / (1 + np.exp(-x))

注：有读者可能会疑惑，即便产生了数值问题，sigmoid函数的函数值也就是为0或1，虽然对sigmoid函数来说这是不可能的输出，但是也不会造成什么问题吧。非也，问题出在后面的对数似然损失函数，如果sigmoid函数的函数值为0或1，则会出现$\ln0$这种非法输入，这可是实打实会出bug的。

损失函数与梯度计算

现在正式进入logistic回归的实现，经过前面多元线性回归的铺垫，我们只需要实现logistic回归损失函数及其梯度的计算即可，剩下的东西都是一样的，损失函数及其梯度的表达式为

$$L(\boldsymbol{w}) = - \sum_{i=1}^{n} y_i \mathrm{ln} \hat{y}_i + (1 - y_i)\mathrm{ln} (1 - \hat{y}_i) = - \boldsymbol{y}^T \ln \hat{\boldsymbol{y}} - (\boldsymbol{1} - \boldsymbol{y})^T \ln (\boldsymbol{1} - \hat{\boldsymbol{y}})$$

$$\frac{ \partial L }{ \partial \boldsymbol{w} } = - X^{T} [ (\boldsymbol{y} \odot (\boldsymbol{1} - \hat{\boldsymbol{y}})) + ((\boldsymbol{1} - \boldsymbol{y}) \odot \hat{\boldsymbol{y}}) ]$$

其中

$$\hat{\boldsymbol{y}} = \sigma(X \boldsymbol{w})$$

注：$\odot$代表Hadamard积，其表示两个形状相同的向量或矩阵逐元素相乘，在numpy中直接使用*运算符即可实现

将上述表达式翻译为代码

y_hat = fun_sigmoid(X @ w).reshape(-1, 1)
# 计算损失值
l = - (y.T @ np.log(y_hat)) - ((1 - y).T @ np.log(1 - y_hat))
# 计算损失值梯度
lg = X.T @ ((1 - y) * y_hat - y * (1 - y_hat))

最后封装为类，完整代码如下

class NLLLoss():
    def __init__(self, regularization='none', weight_decay=1, clip_threshold=15):
        # regularization：正则化类型, weight_decay：正则化系数, clip_threshold：sigmoid函数截断值
        
        # 判断regularization值的合法性
        valid_regularizations = ['none', 'L1', 'L2']
        if regularization not in valid_regularizations:
            raise ValueError(f"Invalid regularization type. Expected one of {valid_regularizations}, but got '{regularization}'")
        
        self.regularization = regularization
        self.weight_decay = weight_decay
        self.clip_threshold = clip_threshold
    
    def __call__(self, w, X, y):
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        
        # 确保w, y为array类型，并控制其形状
        w = np.array(w)
        w = w.reshape(-1, 1)  # w:(num_features+1, 1)
        y = np.array(y)
        y = y.reshape(-1, 1)
        
        y_hat = self.__fun_sigmoid(X @ w).reshape(-1, 1)
        l = - (y.T @ np.log(y_hat)) - ((1 - y).T @ np.log(1 - y_hat))
        l = l / np.size(X, 0)
        
        # 正则化
        if self.regularization == 'L1':
            coef = np.delete(w, 0, axis=0)  # 偏置不参与正则化
            l = l + self.weight_decay * np.linalg.norm(coef, ord=1)
        elif self.regularization == 'L2':
            coef = np.delete(w, 0, axis=0)
            l = l + self.weight_decay * (np.linalg.norm(coef, ord=2) ** 2)
        
        return l
    
    def gradient(self, w, X, y):
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        
        # 确保w, y为array类型，并控制其形状
        w = np.array(w)
        w = w.reshape(-1, 1)  # w:(num_features+1, 1)
        y = np.array(y)
        y = y.reshape(-1, 1)
        
        y_hat = self.__fun_sigmoid(X @ w).reshape(-1, 1)
        lg = X.T @ ((1 - y) * y_hat - y * (1 - y_hat))
        lg = lg / np.size(X, 0)
        
        # 正则化
        if self.regularization == 'L1':
            coef = w
            coef[0] = 0
            lg = lg + self.weight_decay * np.sign(coef)
        if self.regularization == 'L2':
            coef = w
            coef[0] = 0
            lg = lg + self.weight_decay * 2 * coef
        
        return lg
    
    def __fun_sigmoid(self, x):
        # clip_threshold: 输入上下限
        x = np.clip(x, -self.clip_threshold, self.clip_threshold)  # 对输入进行截断，防止输出为0或1
        return 1 / (1 + np.exp(-x))

代码就是上面多元线性回归的代码，改了损失函数跟梯度的计算公式，再加了个sigmoid函数，这里就不多赘述了。

测试

最后，我们对上面实现的代码简单做个测试，由于分类问题的数据集没有所谓的“正确答案”，因此我们直接看准确率，写一个简单的函数用于测试

def logistic_regression_test(w, X, y):
    # sigmoid函数
    def fun_sigmoid(x, clip_threshold=15):
        # clip_threshold: 输入上下限
        x = np.clip(x, -clip_threshold, clip_threshold)  # 对输入进行截断，防止输出为0或1
        return 1 / (1 + np.exp(-x))

    # 对矩阵X做增广
    I = np.ones((np.size(X, 0), 1))  # 全1列
    X = np.hstack((I, X))  # 拼接得到增广X矩阵

    y_hat = fun_sigmoid(X @ w).reshape(-1)
    y_hat = (y_hat >= 0.5).astype(int)  # 根据阈值划分类别，需要将bool转化为int
    return np.sum(y == y_hat) / np.size(X, 0)

criterion = NLLLoss()
w = LR_train(X, y, criterion, epochs_num=10, batch_size=64, alpha=0.01)

epoch 1, loss 0.735997
epoch 2, loss 0.693674
epoch 3, loss 0.657632
epoch 4, loss 0.626926
epoch 5, loss 0.600610
epoch 6, loss 0.578014
epoch 7, loss 0.558276
epoch 8, loss 0.541183
epoch 9, loss 0.526265
epoch 10, loss 0.513075

logistic_regression_test(w, X, y)  # 输出为0.873015873015873

softmax回归实现

数据集生成

最后，我们来到了softmax回归的实现，经过了前面的铺垫，到这里基本上没有什么特别要讲的了，分类数据集的生成同logistic回归，使用make_classification函数，只需将类别数改为多类即可

X, y = make_classification(
    n_samples=1008,
    n_features=3,  # 特征维度
    n_informative=3,  # 所有特征均与标签相关
    n_redundant=0,  # 无冗余特征
    n_classes=3,  # 类别
    class_sep=1.1,  # 样本分离度
    flip_y=0.15,  # 标签噪声
    random_state=1008  # 随机数种子
)

softmax函数实现

然后是softmax函数的实现，其相比于sigmoid函数还是要复杂一点，在此之前我们还是先来看看表达式

$$\mathrm{softmax}(\boldsymbol{x}) = \frac{ e^{\boldsymbol{x}} }{\boldsymbol{1}^T e^{\boldsymbol{x}}}$$

在numpy中，由于广播机制的存在，我们并不需要严格按照数学表达式来，代码实现如下

np.exp(x) / np.sum(np.exp(x), 0)

好，现在我们来解决softmax函数的问题，首先，由于softmax函数是一个向量值函数，而我们需要考虑小批量的情况，也就是说，输入可能为一个向量，也可能是多个向量组成的矩阵，无论何种情况都需要保证正确的输出，对于向量值函数，这两种输入是有区别的，而对于sigmoid函数来说就没有这个问题。

其次，sigmoid函数中提到的算术溢出问题，在softmax函数中依然存在，不过解决方法稍有不同，softmax函数中的解决方法是，将输入向量均减去向量最大值，也就是对向量数值作一个平移，这并不会改变softmax函数的函数值，非常优雅，读者可以自行验证。综上，完整实现代码如下

def fun_softmax(x):
    # x:(num_classes)或(batch_size, num_classes)
    if np.ndim(x) == 1:
        axis = 0
    else:
        axis = 1
    
    x = x - np.max(x, axis=axis, keepdims=True)  # 平移解决数值问题    
    
    return np.exp(x) / np.sum(np.exp(x), axis).reshape(-1, 1)  # 调整形状，使得广播机制可以正常运作

损失函数与梯度计算

最后是损失函数及其梯度计算，交叉熵损失函数的表达式及其梯度表达式为

$$L(W) = - \mathrm{tr} (Y^T \ln \hat{Y})$$

$$\frac{ \partial L }{ \partial W } = - X^{T} (\hat{Y} - Y)$$

其中

$$\hat{Y} = \mathrm{softmax}(XW)$$

将上述公式翻译为代码

Y_hat = fun_softmax(X @ W)  # Y_hat:(n, num_classes)
# 计算损失函数
l = - np.trace(Y @ np.log(Y_hat.T))
# 计算梯度
lg = X.T @ (Y_hat - Y)

再将其封装为类

class CrossEntropyLoss():
    def __init__(self, regularization='none', weight_decay=1):
        # regularization：正则化类型, weight_decay：正则化系数
        
        # 判断regularization值的合法性
        valid_regularizations = ['none', 'L1', 'L2']
        if regularization not in valid_regularizations:
            raise ValueError(f"Invalid regularization type. Expected one of {valid_regularizations}, but got '{regularization}'")
        
        self.regularization = regularization
        self.weight_decay = weight_decay
    
    def __call__(self, W, X, Y):
        # W:(num_features+1, num_classes), X:(n, num_features), Y:(n, 1)
    
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        # X:(n, num_features+1)
        
        # 对Y做one-hot编码
        num_classes = np.max(Y) + 1  # 获取类别数量
        Y = np.eye(num_classes)[Y].reshape(-1, num_classes)  # one-hot编码
        # Y:(n, num_classes)
        
        Y_hat = self.__fun_softmax(X @ W)  # Y_hat:(n, num_classes)
        l = - np.trace(Y @ np.log(Y_hat.T))
        l = l / np.size(X, 0)
        
        # 正则化
        if self.regularization == 'L1':
            coef = np.delete(W, 0, axis=0)  # 偏置不参与正则化
            l = l + self.weight_decay * np.sum(np.abs(coef))  # 矩阵范数的定于与向量范数不一致
        elif self.regularization == 'L2':
            coef = np.delete(W, 0, axis=0)
            l = l + self.weight_decay * np.sum(coef ** 2)  # 矩阵范数的定于与向量范数不一致
        
        return l
    
    def gradient(self, W, X, Y):
        # W:(num_features+1, num_classes), X:(n, num_features), Y:(n, 1)
    
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        # X:(n, num_features+1)
        
        # 对Y做one-hot编码
        num_classes = np.max(Y) + 1  # 获取类别数量
        Y = np.eye(num_classes)[Y].reshape(-1, num_classes)  # one-hot编码
        # Y:(n, num_classes)
        
        Y_hat = self.__fun_softmax(X @ W)  # Y_hat:(n, num_classes)
        
        lg = X.T @ (Y_hat - Y)
        lg = lg / np.size(X, 0)
        
        # 正则化
        if self.regularization == 'L1':
            coef = np.delete(W, 0, axis=0)
            J = np.zeros((1, np.size(W, 1)))
            coef = np.vstack((J, coef))
            lg = lg + self.weight_decay * np.sign(coef)
        if self.regularization == 'L2':
            coef = np.delete(W, 0, axis=0)
            J = np.zeros((1, np.size(W, 1)))
            coef = np.vstack((J, coef))
            lg = lg + self.weight_decay * 2 * coef
        
        return lg
    
    def __fun_softmax(self, x):
        # x:(num_classes)或(batch_size, num_classes)
        if np.ndim(x) == 1:
            axis = 0
        else:
            axis = 1
            
        x = x - np.max(x, axis=axis, keepdims=True)  # 平移解决数值问题
        
        return np.exp(x) / np.sum(np.exp(x), axis).reshape(-1, 1)  # 调整形状，使得广播机制可以正常运作

补充

关于上述softmax回归的实现，这里再补充两点，首先是关于正则化，在上文中我们讲到，$L_1$正则项就是参数的$L_1$范数，$L_2$正则项就是参数的$L_2$范数平方，但是在softmax回归中，这种说法是错的。这是因为softmax回归中的使用的是参数矩阵，矩阵的范数计算跟向量是不一样的，而参数矩阵计算正则项与向量矩阵一样，$L_1$正则项就是所有参数取绝对值相加，$L_2$正则项就是所有参数取平方相加，实现的时候就按照这两种计算方式去实现，千万不能用numpy提供的linalg.norm()函数，这样计算出来的是矩阵范数，是错误的。

然后是关于计算性能，虽然将损失函数与梯度写作矩阵形式从数学上看非常的优雅，并且用代码写出来也非常简单，但是对于softmax回归来说，这里面其实有很多的性能优化空间。首先就是计算损失函数时，矩阵形式的表达式中使用了迹（trace）运算，这相当于$Y^T \ln \hat{Y}$这个矩阵除了对角线元素外，其余元素的计算都是多余的。再者，按照数学公式进行计算，意味着要先对标签做One-Hot编码，由于One-Hot编码产生的矩阵是一个稀疏矩阵，其本身就会浪费一定的内存空间，再者，由于矩阵中多为$0$，后续产生的许多运算也是无效运算，例如损失函数的计算，假设一个样本的正确类别是第$r$类，那这个样本的损失值其实等价于算$\ln \hat{y}^{(r)}$，因为错误类别的$y^{(j)}$都是$0$，算出了$\ln \hat{y}^{(j)}$最后也是乘以$0$，等于没算，而正确类别的$y^{(r)}$是$1$，也不用乘，总结来说就是算$\ln \hat{y}^{(r)}$就完事了。下面给出经过优化的损失函数与梯度的实现

n = X.shape[0]  # 样本数量
Y_hat = fun_softmax(X @ W)  # Y_hat:(n, num_classes)
# 计算损失函数
correct_probs = Y_hat[np.arange(n), Y.ravel()]  # 直接索引正确类别概率
l = -np.sum(np.log(correct_probs))
# 计算梯度
error = Y_hat.copy()
error[np.arange(n), Y.ravel()] -= 1  # 直接修改正确类别的误差
lg = X.T @ error

用上述实现替代原实现，完整代码为

class CrossEntropyLossOptimization():
    def __init__(self, regularization='none', weight_decay=1):
        # regularization：正则化类型, weight_decay：正则化系数
        
        # 判断regularization值的合法性
        valid_regularizations = ['none', 'L1', 'L2']
        if regularization not in valid_regularizations:
            raise ValueError(f"Invalid regularization type. Expected one of {valid_regularizations}, but got '{regularization}'")
        
        self.regularization = regularization
        self.weight_decay = weight_decay
    
    def __call__(self, W, X, Y):
        # W:(num_features+1, num_classes), X:(n, num_features), Y:(n, 1)
    
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        # X:(n, num_features+1)
        
        n = X.shape[0]  # 样本数量
        
        Y_hat = self.__fun_softmax(X @ W)  # Y_hat:(n, num_classes)
        correct_probs = Y_hat[np.arange(n), Y.ravel()]  # 直接索引正确类别概率
        l = -np.sum(np.log(correct_probs))
        l = l / np.size(X, 0)
        
        # 正则化
        if self.regularization == 'L1':
            coef = np.delete(W, 0, axis=0)  # 偏置不参与正则化
            l = l + self.weight_decay * np.sum(np.abs(coef))  # 矩阵范数的定于与向量范数不一致
        elif self.regularization == 'L2':
            coef = np.delete(W, 0, axis=0)
            l = l + self.weight_decay * np.sum(coef ** 2)  # 矩阵范数的定于与向量范数不一致
        
        return l
    
    def gradient(self, W, X, Y):
        # W:(num_features+1, num_classes), X:(n, num_features), Y:(n, 1)
    
        # 对矩阵X做增广
        I = np.ones((np.size(X, 0), 1))  # 全1列
        X = np.hstack((I, X))  # 拼接得到增广X矩阵
        # X:(n, num_features+1)
        
        n = X.shape[0]  # 样本数量

        Y_hat = self.__fun_softmax(X @ W)  # Y_hat:(n, num_classes)
        
        error = Y_hat.copy()
        error[np.arange(n), Y.ravel()] -= 1  # 直接修改正确类别的误差
        lg = X.T @ error
        lg = lg / np.size(X, 0)
        
        # 正则化
        if self.regularization == 'L1':
            coef = np.delete(W, 0, axis=0)
            J = np.zeros((1, np.size(W, 1)))
            coef = np.vstack((J, coef))
            lg = lg + self.weight_decay * np.sign(coef)
        if self.regularization == 'L2':
            coef = np.delete(W, 0, axis=0)
            J = np.zeros((1, np.size(W, 1)))
            coef = np.vstack((J, coef))
            lg = lg + self.weight_decay * 2 * coef
        
        return lg
    
    def __fun_softmax(self, x):
        # x:(num_classes)或(batch_size, num_classes)
        if np.ndim(x) == 1:
            axis = 0
        else:
            axis = 1
            
        x = x - np.max(x, axis=axis, keepdims=True)  # 平移解决数值问题
        
        return np.exp(x) / np.sum(np.exp(x), axis).reshape(-1, 1)  # 调整形状，使得广播机制可以正常运作

测试

最后我们来测试一下上述代码

def softmax_regressor_test(W, X, y):
    # softmax函数
    def fun_softmax(x):
        # x:(num_classes)或(batch_size, num_classes)
        if np.ndim(x) == 1:
            axis = 0
        else:
            axis = 1
            
        x = x - np.max(x, axis=axis, keepdims=True)  # 平移解决数值问题    
        
        return np.exp(x) / np.sum(np.exp(x), axis).reshape(-1, 1)  # 调整形状，使得广播机制可以正常运作
    
    # 对矩阵X做增广
    I = np.ones((np.size(X, 0), 1))  # 全1列
    X = np.hstack((I, X))  # 拼接得到增广X矩阵
    
    Y_hat = fun_softmax(X @ W)  # Y_hat:(n, num_classes)
    Y_hat = np.argmax(Y_hat, axis=1).reshape(-1)
    
    return np.sum(y == Y_hat) / np.size(X, 0)  # 正确率

criterion = CrossEntropyLoss()
W = LR_train(X, y, criterion, epochs_num=10, batch_size=64, alpha=0.1)

epoch 1, loss 0.773139
epoch 2, loss 0.714879
epoch 3, loss 0.694765
epoch 4, loss 0.686554
epoch 5, loss 0.682150
epoch 6, loss 0.679736
epoch 7, loss 0.678368
epoch 8, loss 0.677144
epoch 9, loss 0.676270
epoch 10, loss 0.675807

softmax_regressor_test(W, X, y)

总结

本文从机器学习与深度学习的角度介绍了常用的线性回归算法，从最简单的多元线性回归出发，到引入了正则化的岭回归与Lasso回归，再到结合sigmoid函数与softmax算子实现使用线性回归模型解决分类问题。虽然，在现如今，线性回归模型只能在少数及其适合他的环境中发光发热，但是其深刻的算法思想，以及其作为神经网络基础的重要地位，使得其仍旧是广大学习者无法绕过的一环。