线性回归及其代码实现

前言

目前所流行的机器学习与深度学习算法，可以被总结为是一类利用数据驱动的算法。而回归算法可以说是最早的数据驱动的算法，其简洁有效，且最基础的线性回归算法具有及其稀有的显示解。虽然其模型复杂度不高，且对于数据有一定要求，但是对于满足要求的数据，线性回归具有如今一众大模型都达不到的效果（毕竟线性回归模型小，还有显示解），因此至今在工业界也被广泛使用，有一种大巧不工的美感。而其后来者，如在机器学习时代大名鼎鼎的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})$这类函数，在机器学习中被称为损失函数，这里我们采用的该损失函数被称为平方误差损失函数。

解析解与数值解

最后，就剩下一个问题，如何求解损失函数。这里我们需要先介绍一下解析解与数值解。解析解是指用解析表达式来表达的解，例如我们初中就接触过的二次函数，其最大或最小值一定在$x = - \frac{b}{2a}$处取得，这就是解析解。而什么是数值解呢，就是只关注$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$$

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

现在我们来整理一下，解析解是最优雅的解决方案，但是其可遇不可求，而数值解虽然有很多缺点，但是其不挑函数，因此数值解非常泛用。最后，之所以要在这里提到解析解与数值解，是因为我们上面所使用的损失函数是少数拥有解析解的函数，多数教材也是在线性回归处提到解析解，毕竟后续算法基本都是数值解。最后，我们给出多元线性回归的解析解

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

岭回归与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回归模型就可以写为

$$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]$$

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

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

数值解与正则化

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

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

$$L(\boldsymbol{w}) = - \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] + \alpha ||\boldsymbol{w}||_2^2 \\ L(\boldsymbol{w}) = - \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] + \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回归的模型为

$$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)]$$

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

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

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是一个科学计算库，其对矩阵运算有着很好的支持，对于从零实现算法是合适的选择，最后我们还会使用matplotlib库用于实现一些可视化的部分，下面是导入依赖库的代码，运行环境为jupyter notebook

import matplotlib.pyplot as plt
import numpy as np

import warnings
warnings.filterwarnings("ignore")

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

数据集生成

然后是数据集的选择，由于是自己从零实现，我们自然会希望可以验证算法的效果，所以这里我们选择自己生成数据集，这样我们就可以提前知道“正确答案”。生成数据集的思路就是，首先给定“正确答案”，即参数向量$\boldsymbol{w} = (w_0, w_1, w_2, \cdots, w_p)^T$，再随机生成$n$个自变量$X$，最后根据线性回归模型

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

计算出对应的$n$个因变量$y$，记为向量$\boldsymbol{y} = (y_1, y_2, \cdots, y_n)^T$，为方便计算，上述计算公式可以写为矩阵形式

$$\boldsymbol{y} = X\boldsymbol{w}$$

最后，我们还需要给生成的因变量$\boldsymbol{y}$添加随机噪声，当然这个噪声需要服从均值为0的正态分布（在统计中这是线性回归的假设之一，在模型拟合完成后需要对残差进行检验来验证数据是否满足该假设，由于本文是机器学习视角的，故没有提及，但是这里为了生成数据集的严谨，还是得用上）。记噪声向量$\boldsymbol{\varepsilon} = (\varepsilon_1, \varepsilon_2, \cdots, \varepsilon_n)^T$，其中$\varepsilon_i \stackrel{\mathrm{i.i.d}}{\sim} N(0, 1), \ i = 1,2,\cdots,n$，$\sim$上面的$\mathrm{i.i.d}$代表独立同分布（Independent and identically distributed），意思是$\varepsilon_i$相互独立且均服从标准正态分布，那现在我们得到了因变量的完整计算公式

$$\boldsymbol{y} = X\boldsymbol{w} + \boldsymbol{\varepsilon}$$

理清了生成数据集的思路，接下来就是写代码了。首先，我们要给定数据集的大小，还有参数向量$\boldsymbol{w}$，这里我们先随便填，注意参数向量$\boldsymbol{w}$要使用numpy中的array数据类型。

n = 1008  # 数据集数量
w = np.array([8, 3, -4])  # 权重向量

然后就是随机生成自变量$X$了，一个简单的思路是，先确定定义域$[a,b]$，再生成服从均匀分布$U(a,b)$的随机变量即可。我们的定义域是$R$，那就随便取一个比较大的范围就行，这里我们先取$[-1000,1000]$。至于随机生成，可以使用np.random模块中的函数实现，代码如下

X = np.random.uniform(-1000, 1000, [n, np.size(w)-1])  # 随机生成自变量

需要注意的是，我们生成的自变量矩阵$X$大小为$n \times p$，而非$n \times (p+1)$，因此代码中写为[n, np.size(w)-1]。接下来就是代进公式里计算了

# 对矩阵X做增广
I = np.ones((n, 1))  # 全1列
X_hat = np.hstack((I, X))  # 拼接得到增广X矩阵
# 计算因变量
y = X_hat @ w.T
y = y.reshape(n, 1)  # 调整形状
y = y + np.random.randn(n, 1)  # 添加噪声

上述代码可以被拆分为两个部分，首先是对矩阵$X$的增广，回顾上述讲解原理的部分，为了方便计算，我们使用的矩阵$X$都是第一列为全$1$列的增广$X$矩阵，但是实际情况中，原始的矩阵$X$与我们上面随机生成的一样，是没有全$1$列的，因此在计算因变量之前我们需要对矩阵$X$进行增广，并且不能覆盖原来的矩阵$X$，毕竟后面还要用的。

全$1$列可以使用np.one()函数生成，大小为$n \times 1$，即(n, 1)，再将全$1$列与矩阵$X$拼接即可，注意全$1$列在左，矩阵$X$在右，矩阵拼接可以使用np.hstack()函数实现，该函数用于水平拼接矩阵，对应的还有一个np.vstack()函数用于垂直拼接矩阵。

最后，就是计算因变量了，numpy中可以用运算符@实现矩阵乘法，所以数学上的$X\boldsymbol{w}$直接写成X_hat @ w.T就可以了。这里的.T代表转置，因为我们前面定义的w是一行的，我们要将其转置变成我们数学上定义的一列。

本来事情到这里就差不多了，只需要再加一个噪声向量就行，但是事情并没有这么简单，我上面就提到过numpy有一些在学数学的人看来很奇怪的地方，这里就出现了，我们运行一下下面的代码

y = X_hat @ w.T
y.shape

这段代码会输出X_hat @ w.T的形状，数学上的来说应该是$n \times 1$，也就是这段代码应该输出(1008, 1)，但实际输出的是(1008,)，这表示y只有一个维度，大小为$n$。在数学上只有单个数字我们才会认为是只有一个维度，而一行或一列数字我们会认为其有两个维度，只是其中一个维度为$1$，但是numpy里面是存在一行或一列数字只有一个维度的，使用numpy实现算法可以说一定会遇到这个问题。

解决这个问题的方法也很简单，用reshape方法强制改变形状，把形状从$n$改写为$n \times 1$，例如，上述代码y.reshape(n, 1)的作用就是将y的形状从(1008,)改为(1008, 1)。

当然，也不止这一种解决办法，y之所以维度为$1$是追根溯源因为w的维度为$1$，我们也可以先改变w的维度，再计算y，则y的维度也会是二维

# 写法一
w = np.array([[8, 3, -4]])  # 权重向量

# 写法二
# w = np.array([8, 3, -4])  # 权重向量
# w = w.reshape(1, 3)  # 调整形状

# 计算y
y = X_hat @ w.T
y.shape

注：这也提醒我们，在使用numpy进行矩阵运算时，要注意各计算值的维度，只要维度正确，基本就相当于把数学公式往代码里写了，还是非常方便的。

最后就剩下添加噪声了，依旧是使用np.random模块中的函数，生成一个由服从标准正态分布的随机变量组成的随机向量，形状为$n \times 1$，再进行矩阵加法即可，也就是上述代码中的

y = y + np.random.randn(n, 1)  # 添加噪声

现在，将我们上面的成果封装为函数

# 数据集生成函数
def synthetic_data(n, w):
    # n：数据集数量，w：权重向量
    w = np.array(w)  # 确保w为array类型
    
    X = np.random.uniform(-1000, 1000, [n, np.size(w)-1])  # 随机生成自变量
    
    # 对矩阵X做增广
    I = np.ones((n, 1))  # 全1列
    X_hat = np.hstack((I, X))  # 拼接得到增广X矩阵
    
    # 计算因变量
    y = X_hat @ w.T
    y = y.reshape(n, 1)  # 调整形状
    y = y + np.random.randn(n, 1)  # 添加噪声
    return X, y

并且填入参数，生成数据集

# 数据集生成
n = 1008  # 数据集数量
w = [8, 3, -4]  # 权重

X, y = synthetic_data(n, w)

多元线性回归：解析解

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

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

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

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

矩阵乘法跟转置在刚才生成数据集的时候我们都见过了，求逆矩阵可以使用np.linalg.inv()函数实现，现在我们就可以把上面的式子翻译成代码了，也就是刚才那一行核心代码，将其封装成函数，别忘了增广$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
    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

其中的X = X.reshape(-1, len(w)-1)是用来防止报错的，还是上面提到的老问题，如果我们输入的X是一维的，那么就会报错，我们需要将其转化为二维，而-1代表不指定，例如X.reshape(1, -1)表示将X转换为一行，至于列数交由程序判断，例如X本身的形状是(504, 2)，那就会被转化为(1, 1008)，这也是reshape方法常用的写法。

总结

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