做机器学习绕不开的第一个算法几乎所有人入门时都会碰到线性回归。但很多人学完只会调sklearn.linear_model.LinearRegression真让写个公式推导或者用原生Python从零实现一遍就卡壳了。这篇文章我想把线性回归从数学原理到代码实现完整串一遍把“为什么这样做”讲清楚而不只是贴出一堆公式让读者死记硬背。如果你正准备学机器学习、或者学过但一直对数学部分似懂非懂这篇应该能帮你把最后那层窗户纸捅破。线性回归表面上是“画一条线拟合数据点”但它的内核是三个环环相扣的问题怎么定义“拟合得好”怎么算出让拟合效果最好的参数算不出来的时候怎么办这三个问题分别对应损失函数、正规方程和梯度下降。把这套逻辑吃透了后面学逻辑回归、岭回归、甚至神经网络都会顺很多因为它们骨子里是同一套思路——定义损失、优化参数。1. 从问题到数学表达线性回归到底在解什么1.1 先放下代码用一句话说清楚线性回归线性回归做的事情用最直白的话讲给定一堆数据点每个点有一些特征比如房子的面积、卧室数量对应一个要预测的值比如房价我们想找一条直线在特征多的时候是一个超平面让这条线尽量贴合所有数据点。贴合的意思是线上预测的值和真实值的差距整体最小。数学上假设有 ( m ) 个样本每个样本有 ( n ) 个特征我们把第 ( i ) 个样本写作 ( (x^{(i)}, y^{(i)}) )其中 ( x^{(i)} [x_1^{(i)}, x_2^{(i)}, ..., x_n^{(i)}] ) 是特征向量( y^{(i)} ) 是真实标签。线性回归的模型假设是[ h_\theta(x) \theta_0 \theta_1 x_1 \theta_2 x_2 ... \theta_n x_n ]这里的 ( \theta ) 就是我们要学习的参数( \theta_0 ) 是截距项bias其余 ( \theta_j ) 是每个特征对应的权重。为了方便矩阵运算通常会加一列全 1 的特征把 ( \theta_0 ) 也吸收进向量于是模型可以写成紧凑形式[ h_\theta(x) \theta^T x ]其中 ( \theta [\theta_0, \theta_1, ..., \theta_n]^T )( x [1, x_1, ..., x_n]^T )。这一步看起来简单但它是整个推导的基础——把求和写成矩阵乘法之后后面的求导和计算才能高效进行。注意很多人一开始搞不清“特征”和“参数”的区别。特征是数据里给你的、已知的参数是模型要学的、未知的。线性回归的任务就是找到一组参数让模型在所有样本上的预测误差最小。1.2 为什么线性模型值得先学简单但不简陋线性回归的假设是所有机器学习模型里最简单的特征和目标值之间是线性关系。但简单不等于没用。它有三个很实际的价值。第一它是理解复杂模型的脚手架。神经网络本质上就是多个线性变换叠加非线性激活函数逻辑回归则是在线性输出上套了一个 sigmoid。如果你连线性回归的损失函数和梯度都搞不明白后面这些东西只会越学越糊涂。第二线性模型有很强的可解释性。每个特征的权重 ( \theta_j ) 直接告诉我们这个特征每变化一个单位预测值变化多少。在很多需要向业务方解释模型的场景比如金融风控、医疗分析这种直接性非常宝贵黑盒模型反而没人敢用。第三线性回归的损失函数是凸函数意味着它只有一个全局最小值没有局部最优的困扰。这和我们后面会碰到的神经网络非常不一样。正因为凸优化性质好线性回归成了理解梯度下降的最佳试验田——算法跑崩了问题一定出在实现或调参上而不是数学上存在多个坑。2. 损失函数与最小二乘法如何量化“拟合得好不好”2.1 均方误差MSE是怎么来的以及为什么用它有了模型假设 ( h_\theta(x) )接下来要定义一个目标函数来衡量预测值与真实值的差距。常见的选择是均方误差Mean Squared Error, MSE也叫残差平方和RSS[ J(\theta) \frac{1}{2m} \sum_{i1}^{m} (h_\theta(x^{(i)}) - y^{(i)})^2 ]注意这里为什么有系数 ( \frac{1}{2m} )。( \frac{1}{m} ) 是求平均让损失不随样本量增长而膨胀方便不同数据集间比较。而 ( \frac{1}{2} ) 纯粹是为了后面求导时把平方项的 2 抵消掉属于“为了让数学更清爽”的约定不影响最优解。这是一个很典型的细节很多人背公式时不知道为什么有 1/2其实就是求导方便。为什么要用平方误差而不是绝对误差因为平方函数处处可导且对误差的惩罚是二次的——误差越大惩罚增长得越快。这带来的实际效果是模型会更努力地避开那些偏离很远的大误差样本。绝对误差在零点不可导求梯度时会有麻烦更关键的是一旦样本误差方向不同直接求和会把正负误差抵消平方之后则巧妙地让所有误差都变成对损失的“贡献”。一个更严谨的视角来自概率论。假设真实值与预测值之间的误差符合均值为零的高斯分布也就是噪声是随机的、无偏的那么用极大似然估计推导最后得到的优化目标就是最小化均方误差。这一点在关注“为什么选择这个损失函数”时很有分量不是我们拍脑袋选了 MSE而是在“噪声服从正态分布”这一合理假设下MSE 是最自然的产物。2.2 最小二乘法的核心思想让残差平方和最小“最小二乘”这个名字里的“二乘”指的就是平方。它的核心思想很朴素我们希望找一条线使所有样本点到这条线的竖直距离残差的平方和最小。这里有一个容易混淆的点是“竖直距离”而不是“垂直距离”。二维平面上点到直线的垂直距离是最短距离但线性回归用的是竖直距离预测误差方向。为什么因为我们的目标是要预测 ( y )误差只存在于 ( y ) 方向特征 ( x ) 被认为是精确观测到的。如果特征本身也有噪声那就变成另一个问题比如正交回归/总最小二乘了用到的数学也完全不同。这个区别在实际应用中会影响结果但在入门阶段先记住线性回归最小化的是竖直方向上的误差。有了损失函数 ( J(\theta) )问题转化为一个无约束优化问题[ \hat{\theta} \arg\min_{\theta} J(\theta) ]接下来的核心问题就是怎么求这个最优的 ( \theta )3. 求解参数的两条路正规方程与梯度下降3.1 正规方程一步到位的解析解既然 ( J(\theta) ) 是凸函数那么它的极小值点就是全局最小值点。求极小值的一个标准套路是求导并令导数等于零。线性回归的 MSE 是二次函数求导得到的是一组线性方程可以解析求解。先把损失函数写成矩阵形式。注意前面的 ( h_\theta(x) \theta^T x )对所有样本写成矩阵就是[ J(\theta) \frac{1}{2m} (X\theta - y)^T (X\theta - y) ]其中 ( X ) 是 ( m \times (n1) ) 的设计矩阵每一行是一个样本的特征向量含开头全为1的列( y ) 是 ( m \times 1 ) 的真实值向量。现在对 ( \theta ) 求梯度。这里会用到一个非常重要的矩阵求导结论[ \nabla_\theta (X\theta - y)^T (X\theta - y) 2 X^T (X\theta - y) ]展开来看由内积求导公式 ( \nabla_\theta (\theta^T X^T X \theta) 2 X^T X \theta ) 和 ( \nabla_\theta (-2 y^T X \theta) -2 X^T y ) 可推出。令梯度为零[ 2 X^T (X\theta - y) 0 ]即[ X^T X \theta X^T y ]如果 ( X^T X ) 可逆就得到正规方程[ \hat{\theta} (X^T X)^{-1} X^T y ]这就是最小二乘解的解析表达式。整个过程只有几步难点全在矩阵求导能不能写对。我自己学的时候会在草稿上把二维情况一个特征的导数先展开写一遍确认系数没问题再推广到矩阵形式这样不容易出错。3.2 正规方程的两个局限可逆性与计算复杂度正规方程很优雅但在实际工程中它有两个绕不开的问题。第一个问题是 ( X^T X ) 可能不可逆奇异。当特征之间存在完全共线性比如一个特征是另一个特征的倍数或者特征数量大于样本数量时( X^T X ) 就是奇异矩阵逆矩阵不存在公式用不了。解决办法包括删除冗余特征、使用主成分分析降维或者改用岭回归——岭回归在 ( X^T X ) 上加了一个 ( \lambda I ) 项既保证了可逆又对参数做了正则化约束。第二个问题是计算复杂度。( X^T X ) 是一个 ( (n1) \times (n1) ) 的矩阵求逆的计算复杂度大约是 ( O(n^3) )。当特征数量 n 是几百几千时计算量和内存占用都会变得很高。要人脸识别那种几万维的特征直接求逆基本不现实。更别提如果数据集特别大几百万样本构造 ( X^T X ) 本身就要消耗大量内存。所以正规方程适合特征维度低、样本量适中的场景一旦维度高或者数据量大我们通常会转向第二种方法——梯度下降。3.3 梯度下降一步步走到最低点梯度下降的思路可以这样理解你站在一座山坡上损失函数曲面眼睛看不见全局但能感觉到脚下哪个方向是下坡最陡的。每次沿着最陡的方向迈一小步走很多步之后就能到达山谷底部最小值点。数学上最陡下降方向就是梯度的反方向。梯度是函数上升最快的方向所以我们要沿着它的反方向更新参数。对于参数 ( \theta_j )更新规则是[ \theta_j : \theta_j - \alpha \frac{\partial J(\theta)}{\partial \theta_j} ]其中 ( \alpha ) 是学习率决定了每一步迈多大。把 MSE 的偏导数具体算出来。对单个样本 ( (x^{(i)}, y^{(i)}) ) 来说[ \frac{\partial}{\partial \theta_j} \frac{1}{2}(h_\theta(x^{(i)}) - y^{(i)})^2 (h_\theta(x^{(i)}) - y^{(i)}) x_j^{(i)} ]对所有样本求和再取平均得到批量梯度下降Batch Gradient Descent的参数更新公式[ \theta_j : \theta_j - \alpha \frac{1}{m} \sum_{i1}^{m} (h_\theta(x^{(i)}) - y^{(i)}) x_j^{(i)} ]这里有几个关键直觉值得说透每次更新的方向是误差 ( (h_\theta(x) - y) ) 与特征值 ( x_j ) 的乘积的均值。如果某个特征在误差大的样本上取值也大它对参数更新的“推动力”就大这符合直觉——大误差加上大特征值说明该参数对误差的“责任”更大。更新量的大小由学习率 ( \alpha ) 控制。学习率太大参数会在最小值附近来回震荡甚至发散学习率太小收敛慢得让人抓狂。每次更新都会用到全部样本因此批量梯度下降每次迭代的计算量是 ( O(mn) )。样本量大时每一轮迭代都很耗时这时候可以退而求其次用随机梯度下降每次只用一个样本更新或小批量梯度下降每次用一小批样本更新。提示梯度下降是一个通用优化框架它并不假设损失函数必须是凸函数。在线性回归里因为凸性梯度下降一定能收敛到全局最优但在神经网络里非凸的损失函数可能让梯度下降停在局部最优或鞍点这就是为什么深度学习里优化问题更复杂。3.4 学习率、迭代次数与特征缩放三个必须关心的调参项梯度下降有“三大命门”学习率、迭代次数、特征缩放。学习率 ( \alpha ) 的选择是最容易踩坑的地方。一个常用的办法是画损失函数随迭代次数变化的曲线如果损失在下降说明学习率合适如果损失震荡或上升说明学习率太大需要调小。实际中常见的取值是 0.01、0.03、0.1、0.3也可以先用非常小的值比如 0.001试跑观察曲线再逐步放大。没有一劳永逸的经验值只能多做实验。迭代次数决定算法跑多久。一种做法是设定一个固定次数比如 1000 次跑完后看损失是否收敛另一种是设定收敛条件——当相邻两次迭代的损失变化小于某个阈值比如 ( 10^{-5} )时就认为收敛并停止。我更推荐后者它可以避免白白浪费计算资源。特征缩放Feature Scaling是新手最容易忽略但影响极大的操作。如果某个特征的取值范围是 0~10000另一个是 0~1那么梯度下降在未缩放的损失曲面上会呈现出狭长的“碗形”导致参数更新时在一个方向上大幅度震荡、另一个方向上缓慢前进收敛非常慢。把每个特征缩放到相近的尺度比如均值归一化或标准化到 0 附近、方差为 1损失曲面会更接近正圆形收敛路径也直接得多。一句话总结正规方程适合特征少、数据量小的场景不用调参一步到位梯度下降适合特征多、数据量大的场景灵活可控但需要照顾好学习率、迭代次数和特征缩放这些细节。4. 从公式到代码用 NumPy 手写一个线性回归4.1 准备工作生成一份可复现的模拟数据为了让代码演示能直观地验证数学推导我先生成一份带噪声的线性数据。这里假设真实模型是 ( y 4 3x \epsilon )其中 ( \epsilon ) 是高斯噪声。使用numpy.random生成数据并设置随机种子保证每次运行结果一致方便调试复现。import numpy as np np.random.seed(42) # 生成 100 个样本特征 x 均匀分布在 [0, 2] X np.random.uniform(0, 2, size(100, 1)) y 4 3 * X np.random.randn(100, 1) * 0.5 # 可视化一下 import matplotlib.pyplot as plt plt.scatter(X, y, alpha0.7) plt.xlabel(x) plt.ylabel(y) plt.title(Simulated Linear Data) plt.show()这份数据只有 100 个样本、1 个特征非常适合用来同时跑通正规方程和梯度下降。加噪声是关键步骤——没有噪声的话数据完全落在一条直线上问题就退化成简单的解方程完全体现不出机器学习“从不确定中找规律”的乐趣。为了让矩阵运算方便需要给特征矩阵加一列全 1作为截距项对应的“哑特征”# 给 X 加一列 1作为截距项 X_b np.c_[np.ones((100, 1)), X]加了这一列之后模型就统一成了 ( h_\theta(x) \theta^T x )不需要再单独区分“截距”和“权重”这在代码实现里会省掉很多分支判断。4.2 用手写代码复现正规方程正规方程的代码简单到让人怀疑是不是真的有用# 正规方程θ (X^T X)^(-1) X^T y def normal_equation(X_b, y): theta np.linalg.inv(X_b.T X_b) X_b.T y return theta theta_ne normal_equation(X_b, y) print(正规方程结果: θ0 , theta_ne[0][0], , θ1 , theta_ne[1][0])跑完会发现结果非常接近真实的 ( \theta_04、\theta_13 )因为有噪声不会完全等于但已经很接近了。这个实现里最核心的一行就是np.linalg.inv(X_b.T X_b) X_b.T y和前面推导的正规方程公式一一对应。需要提醒的是np.linalg.inv在实际工程中并不推荐因为直接用inv求逆在数值上不如解线性方程组稳定。更专业的做法是用np.linalg.solve# 更稳定的做法直接解线性方程组 X^T X θ X^T y theta_ne np.linalg.solve(X_b.T X_b, X_b.T y)np.linalg.solve内部使用 LU 分解等方式求解线性方程组避免了直接求逆带来的数值误差和效率浪费。这也算是一个从“数学公式”到“工程代码”的典型差异公式里有逆矩阵但真正写代码时不一定要真的去算逆矩阵。4.3 用手写代码实现批量梯度下降接下来是梯度下降。需要做的事情有三件定义损失函数、计算梯度、迭代更新参数。写成如下代码def mse_loss(theta, X_b, y): m len(y) error X_b theta - y return float((1 / (2 * m)) * (error.T error)) def batch_gradient_descent(X_b, y, alpha0.1, n_iterations1000): m, n X_b.shape theta np.random.randn(n, 1) # 随机初始化 history [] for iteration in range(n_iterations): gradient (1 / m) * (X_b.T (X_b theta - y)) theta theta - alpha * gradient loss mse_loss(theta, X_b, y) history.append(loss) # 每 100 次迭代打印一次损失 if iteration % 100 0: print(fIteration {iteration}, Loss: {loss:.6f}) return theta, history # 标准化特征只标准化 x不标准化全为 1 的那列 def standardize(X): mean X.mean(axis0) std X.std(axis0) return (X - mean) / std, mean, std X_scaled, mean, std standardize(X) X_s np.c_[np.ones((100, 1)), X_scaled] theta_gd, history batch_gradient_descent(X_s, y, alpha0.1, n_iterations1000) print(梯度下降结果标准化后空间: θ0 , theta_gd[0][0], , θ1 , theta_gd[1][0])注意到这里我先把特征 ( x ) 做了标准化再在左边补一列 1。这一步不是可选项而是会让梯度下降收敛速度产生量级差距的关键操作。如果不做标准化学习率在 0.1 时可能直接导致损失变成 NaN或者收敛极慢标准化之后0.1 的学习率就能跑得很稳。训练完之后的损失曲线也值得看一眼它会呈现出典型的“快速下降、缓慢逼近”形态这是梯度下降收敛的正常模式。如果看到损失曲线震荡或上扬就得去检查学习率是不是太大。由于标准化的原因训练出的参数是在标准化特征空间里的。如果要把模型还原到原始特征空间需要做一个逆变换在标准化后的空间里模型是 ( \hat{y} \theta_0 \theta_1 x )其中 ( x (x - \mu) / \sigma )把它展开就能得到原始空间里的等效参数。# 把标准化空间中的参数还原到原始空间 theta1_orig theta_gd[1][0] / std[0] theta0_orig theta_gd[0][0] - theta1_orig * mean[0] print(f还原到原始空间: θ0 {theta0_orig:.4f}, θ1 {theta1_orig:.4f})跑完这段代码对比正规方程和还原后的梯度下降结果两个数字应该很接近。这其实是一个非常好的自检手段用正规方程的结果作为“标准答案”验证梯度下降实现是否正确。如果两者差距过大基本可以断定梯度下降的学习率、迭代次数或特征缩放处理出了问题。4.4 对比两种方法的效果与适用场景拿同一份数据分别跑正规方程和梯度下降数值结果几乎一样这验证了两种方法最终都在逼近同一个最优解。但在使用体验上两者区别明显对比维度正规方程批量梯度下降计算方式解析求解一步到位迭代逼近逐步优化是否需要调参不需要需要调学习率、迭代次数特征是否需要缩放不需要必须做特征缩放否则收敛慢特征数量很大时如 10 万维不可行矩阵求逆开销太大仍然可行每轮迭代 O(mn)样本数量很大时如 100 万样本构造 X^T X 可能耗尽内存可用随机/小批量梯度下降缓解数值稳定性直接求逆可能引入数值误差相对稳定取决于学习率选择表里的结论已经在实践中反复被验证小规模、特征少、一次性建模直接正规方程大规模、特征多、需要在线更新模型选择梯度下降或它的变体。当然现实工程中更多时候直接用sklearn的LinearRegression底层有 SVD 分解求解或SGDRegressor但理解背后的原理能帮你在模型跑出异常结果时快速定位问题。提示sklearn.linear_model.LinearRegression默认不是用正规方程求逆的它用的是最小二乘求解器底层通过 SVD 分解实现数值稳定性比直接求逆好很多。但它的数学本质和我们手写的正规方程是同一个解。5. 常见问题与排查技巧实录5.1 损失变成 NaN 或爆炸学习率太大还是特征没缩放这是梯度下降实现里最经典的问题。现象是训练过程中损失打印出来是nan或者第一次迭代的损失就巨大无比。先用排除法定位问题。第一步检查特征是否做了标准化——没标准化的特征量级可能从 0.001 到 100000梯度在某个方向会特别大更新一步就直接飞出去了。第二步检查学习率把学习率降到 0.001 或 0.0001 试跑如果问题消失说明是学习率太大的问题。第三步检查数据里有没有缺失值或无穷大这类脏数据会让矩阵运算直接产生 NaN。贴一个我自己的排查经验当时一个同事训练多项式回归损失前几轮在下降快到 100 轮时突然变成nan。查了半天发现不是学习率问题而是数据里有两行特征值特别大在梯度累积后数值溢出。最后做了标准化加裁剪clip问题立刻解决。梯度下降对数值范围非常敏感这是真实工程里最常见也最容易被忽略的坑。5.2 正规方程报错“Singular matrix”共线性或特征维度超样本量用np.linalg.inv或np.linalg.solve求解时如果数据存在完全共线性两个特征成比例或特征数量大于样本数量就会遇到奇异矩阵报错。解决思路有几个检查特征之间是否存在线性相关用np.corrcoef或画相关性矩阵热力图来排查删掉冗余的特征比如把“摄氏温度”和“华氏温度”同时放进模型这种低级错误改用岭回归它在 ( X^TX ) 上加了一个 ( \lambda I )保证矩阵可逆的同时还能约束参数大小是从数学上根治奇异性的方案。顺便提一句即使理论上矩阵可逆在数值计算中也可能因为同一个特征数量级差异过大导致“病态矩阵”计算出的参数极不稳定。这就得靠特征缩放和标准化来解决。5.3 多特征场景下参数含义发生变化先看尺度再看系数很多初学者在拿到多特征线性回归的参数后直接比较权重大小来判断特征重要性这其实很危险。当一个特征的单位是“平方米”数值几百上千另一个的单位是“房间数”数值单位个位数时前者的系数天然会很小后者的系数天然会很大但这不代表前者不重要。只有先把所有特征标准化到相同尺度再去比较系数的大小才有意义。所以在实际项目中如果要解释特征重要性一个标准流程是标准化→训练模型→比较标准化后系数的绝对值。要注意标准化的系数也要结合业务逻辑去验证不要盲目相信数字。5.4 拟合效果不好不全是模型的错先排查数据再升级模型如果你已经跑通了代码但发现模型在线性数据上的拟合效果很差先别急着换复杂模型。按这个顺序排查查看数据是否存在异常值一个离群点可能把回归线拉偏很多检查是否存在明显的非线性关系画个散点图看看确认特征是否真的和目标值相关数据里塞了完完全全的随机噪声特征模型学不到规律实属正常。线性回归拟合不好很多时候不是线性模型不行而是数据里本身就存在非线性的规律比如收入随年龄先升后降。这时候升级到多项式回归或引入交互特征才是正解。不过这是后话把线性回归的数学原理和代码实现吃透后面这些扩展都是顺手的事。写在后面的一点心得我自己刚学线性回归的时候也走过一段“背公式、调包、跑通就完事”的弯路。后来真正动手用 NumPy 从零实现一遍把每个公式展开、把每行代码和公式对应上才算是把最小二乘和梯度下降装进了脑子里。特别是矩阵求导那一步一开始总觉得“知道结论就行”但真正自己推一遍之后后面学到逻辑回归、多分类 Softmax 的梯度推导时明显顺畅了很多。再说一个建议不管你现在用的是 PyTorch、TensorFlow 还是 sklearn都值得在某个版本里用原生 NumPy 手写一遍线性回归。这个过程不需要很长时间但对你理解机器学习底层的优化机制、甚至后续排查模型不收敛问题帮助都是巨大的。数学公式看着枯燥但它们从代码里跑出结果的那一刻很多抽象的概念都会变得具体起来。这也是我把这篇文章的标题定为“从公式到代码实现”的原因——公式和代码本来就是一件事的两面。 线性回归从公式到代码：损失函数、正规方程与梯度下降全解析
做机器学习绕不开的第一个算法几乎所有人入门时都会碰到线性回归。但很多人学完只会调sklearn.linear_model.LinearRegression真让写个公式推导或者用原生Python从零实现一遍就卡壳了。这篇文章我想把线性回归从数学原理到代码实现完整串一遍把“为什么这样做”讲清楚而不只是贴出一堆公式让读者死记硬背。如果你正准备学机器学习、或者学过但一直对数学部分似懂非懂这篇应该能帮你把最后那层窗户纸捅破。线性回归表面上是“画一条线拟合数据点”但它的内核是三个环环相扣的问题怎么定义“拟合得好”怎么算出让拟合效果最好的参数算不出来的时候怎么办这三个问题分别对应损失函数、正规方程和梯度下降。把这套逻辑吃透了后面学逻辑回归、岭回归、甚至神经网络都会顺很多因为它们骨子里是同一套思路——定义损失、优化参数。1. 从问题到数学表达线性回归到底在解什么1.1 先放下代码用一句话说清楚线性回归线性回归做的事情用最直白的话讲给定一堆数据点每个点有一些特征比如房子的面积、卧室数量对应一个要预测的值比如房价我们想找一条直线在特征多的时候是一个超平面让这条线尽量贴合所有数据点。贴合的意思是线上预测的值和真实值的差距整体最小。数学上假设有 ( m ) 个样本每个样本有 ( n ) 个特征我们把第 ( i ) 个样本写作 ( (x^{(i)}, y^{(i)}) )其中 ( x^{(i)} [x_1^{(i)}, x_2^{(i)}, ..., x_n^{(i)}] ) 是特征向量( y^{(i)} ) 是真实标签。线性回归的模型假设是[ h_\theta(x) \theta_0 \theta_1 x_1 \theta_2 x_2 ... \theta_n x_n ]这里的 ( \theta ) 就是我们要学习的参数( \theta_0 ) 是截距项bias其余 ( \theta_j ) 是每个特征对应的权重。为了方便矩阵运算通常会加一列全 1 的特征把 ( \theta_0 ) 也吸收进向量于是模型可以写成紧凑形式[ h_\theta(x) \theta^T x ]其中 ( \theta [\theta_0, \theta_1, ..., \theta_n]^T )( x [1, x_1, ..., x_n]^T )。这一步看起来简单但它是整个推导的基础——把求和写成矩阵乘法之后后面的求导和计算才能高效进行。注意很多人一开始搞不清“特征”和“参数”的区别。特征是数据里给你的、已知的参数是模型要学的、未知的。线性回归的任务就是找到一组参数让模型在所有样本上的预测误差最小。1.2 为什么线性模型值得先学简单但不简陋线性回归的假设是所有机器学习模型里最简单的特征和目标值之间是线性关系。但简单不等于没用。它有三个很实际的价值。第一它是理解复杂模型的脚手架。神经网络本质上就是多个线性变换叠加非线性激活函数逻辑回归则是在线性输出上套了一个 sigmoid。如果你连线性回归的损失函数和梯度都搞不明白后面这些东西只会越学越糊涂。第二线性模型有很强的可解释性。每个特征的权重 ( \theta_j ) 直接告诉我们这个特征每变化一个单位预测值变化多少。在很多需要向业务方解释模型的场景比如金融风控、医疗分析这种直接性非常宝贵黑盒模型反而没人敢用。第三线性回归的损失函数是凸函数意味着它只有一个全局最小值没有局部最优的困扰。这和我们后面会碰到的神经网络非常不一样。正因为凸优化性质好线性回归成了理解梯度下降的最佳试验田——算法跑崩了问题一定出在实现或调参上而不是数学上存在多个坑。2. 损失函数与最小二乘法如何量化“拟合得好不好”2.1 均方误差MSE是怎么来的以及为什么用它有了模型假设 ( h_\theta(x) )接下来要定义一个目标函数来衡量预测值与真实值的差距。常见的选择是均方误差Mean Squared Error, MSE也叫残差平方和RSS[ J(\theta) \frac{1}{2m} \sum_{i1}^{m} (h_\theta(x^{(i)}) - y^{(i)})^2 ]注意这里为什么有系数 ( \frac{1}{2m} )。( \frac{1}{m} ) 是求平均让损失不随样本量增长而膨胀方便不同数据集间比较。而 ( \frac{1}{2} ) 纯粹是为了后面求导时把平方项的 2 抵消掉属于“为了让数学更清爽”的约定不影响最优解。这是一个很典型的细节很多人背公式时不知道为什么有 1/2其实就是求导方便。为什么要用平方误差而不是绝对误差因为平方函数处处可导且对误差的惩罚是二次的——误差越大惩罚增长得越快。这带来的实际效果是模型会更努力地避开那些偏离很远的大误差样本。绝对误差在零点不可导求梯度时会有麻烦更关键的是一旦样本误差方向不同直接求和会把正负误差抵消平方之后则巧妙地让所有误差都变成对损失的“贡献”。一个更严谨的视角来自概率论。假设真实值与预测值之间的误差符合均值为零的高斯分布也就是噪声是随机的、无偏的那么用极大似然估计推导最后得到的优化目标就是最小化均方误差。这一点在关注“为什么选择这个损失函数”时很有分量不是我们拍脑袋选了 MSE而是在“噪声服从正态分布”这一合理假设下MSE 是最自然的产物。2.2 最小二乘法的核心思想让残差平方和最小“最小二乘”这个名字里的“二乘”指的就是平方。它的核心思想很朴素我们希望找一条线使所有样本点到这条线的竖直距离残差的平方和最小。这里有一个容易混淆的点是“竖直距离”而不是“垂直距离”。二维平面上点到直线的垂直距离是最短距离但线性回归用的是竖直距离预测误差方向。为什么因为我们的目标是要预测 ( y )误差只存在于 ( y ) 方向特征 ( x ) 被认为是精确观测到的。如果特征本身也有噪声那就变成另一个问题比如正交回归/总最小二乘了用到的数学也完全不同。这个区别在实际应用中会影响结果但在入门阶段先记住线性回归最小化的是竖直方向上的误差。有了损失函数 ( J(\theta) )问题转化为一个无约束优化问题[ \hat{\theta} \arg\min_{\theta} J(\theta) ]接下来的核心问题就是怎么求这个最优的 ( \theta )3. 求解参数的两条路正规方程与梯度下降3.1 正规方程一步到位的解析解既然 ( J(\theta) ) 是凸函数那么它的极小值点就是全局最小值点。求极小值的一个标准套路是求导并令导数等于零。线性回归的 MSE 是二次函数求导得到的是一组线性方程可以解析求解。先把损失函数写成矩阵形式。注意前面的 ( h_\theta(x) \theta^T x )对所有样本写成矩阵就是[ J(\theta) \frac{1}{2m} (X\theta - y)^T (X\theta - y) ]其中 ( X ) 是 ( m \times (n1) ) 的设计矩阵每一行是一个样本的特征向量含开头全为1的列( y ) 是 ( m \times 1 ) 的真实值向量。现在对 ( \theta ) 求梯度。这里会用到一个非常重要的矩阵求导结论[ \nabla_\theta (X\theta - y)^T (X\theta - y) 2 X^T (X\theta - y) ]展开来看由内积求导公式 ( \nabla_\theta (\theta^T X^T X \theta) 2 X^T X \theta ) 和 ( \nabla_\theta (-2 y^T X \theta) -2 X^T y ) 可推出。令梯度为零[ 2 X^T (X\theta - y) 0 ]即[ X^T X \theta X^T y ]如果 ( X^T X ) 可逆就得到正规方程[ \hat{\theta} (X^T X)^{-1} X^T y ]这就是最小二乘解的解析表达式。整个过程只有几步难点全在矩阵求导能不能写对。我自己学的时候会在草稿上把二维情况一个特征的导数先展开写一遍确认系数没问题再推广到矩阵形式这样不容易出错。3.2 正规方程的两个局限可逆性与计算复杂度正规方程很优雅但在实际工程中它有两个绕不开的问题。第一个问题是 ( X^T X ) 可能不可逆奇异。当特征之间存在完全共线性比如一个特征是另一个特征的倍数或者特征数量大于样本数量时( X^T X ) 就是奇异矩阵逆矩阵不存在公式用不了。解决办法包括删除冗余特征、使用主成分分析降维或者改用岭回归——岭回归在 ( X^T X ) 上加了一个 ( \lambda I ) 项既保证了可逆又对参数做了正则化约束。第二个问题是计算复杂度。( X^T X ) 是一个 ( (n1) \times (n1) ) 的矩阵求逆的计算复杂度大约是 ( O(n^3) )。当特征数量 n 是几百几千时计算量和内存占用都会变得很高。要人脸识别那种几万维的特征直接求逆基本不现实。更别提如果数据集特别大几百万样本构造 ( X^T X ) 本身就要消耗大量内存。所以正规方程适合特征维度低、样本量适中的场景一旦维度高或者数据量大我们通常会转向第二种方法——梯度下降。3.3 梯度下降一步步走到最低点梯度下降的思路可以这样理解你站在一座山坡上损失函数曲面眼睛看不见全局但能感觉到脚下哪个方向是下坡最陡的。每次沿着最陡的方向迈一小步走很多步之后就能到达山谷底部最小值点。数学上最陡下降方向就是梯度的反方向。梯度是函数上升最快的方向所以我们要沿着它的反方向更新参数。对于参数 ( \theta_j )更新规则是[ \theta_j : \theta_j - \alpha \frac{\partial J(\theta)}{\partial \theta_j} ]其中 ( \alpha ) 是学习率决定了每一步迈多大。把 MSE 的偏导数具体算出来。对单个样本 ( (x^{(i)}, y^{(i)}) ) 来说[ \frac{\partial}{\partial \theta_j} \frac{1}{2}(h_\theta(x^{(i)}) - y^{(i)})^2 (h_\theta(x^{(i)}) - y^{(i)}) x_j^{(i)} ]对所有样本求和再取平均得到批量梯度下降Batch Gradient Descent的参数更新公式[ \theta_j : \theta_j - \alpha \frac{1}{m} \sum_{i1}^{m} (h_\theta(x^{(i)}) - y^{(i)}) x_j^{(i)} ]这里有几个关键直觉值得说透每次更新的方向是误差 ( (h_\theta(x) - y) ) 与特征值 ( x_j ) 的乘积的均值。如果某个特征在误差大的样本上取值也大它对参数更新的“推动力”就大这符合直觉——大误差加上大特征值说明该参数对误差的“责任”更大。更新量的大小由学习率 ( \alpha ) 控制。学习率太大参数会在最小值附近来回震荡甚至发散学习率太小收敛慢得让人抓狂。每次更新都会用到全部样本因此批量梯度下降每次迭代的计算量是 ( O(mn) )。样本量大时每一轮迭代都很耗时这时候可以退而求其次用随机梯度下降每次只用一个样本更新或小批量梯度下降每次用一小批样本更新。提示梯度下降是一个通用优化框架它并不假设损失函数必须是凸函数。在线性回归里因为凸性梯度下降一定能收敛到全局最优但在神经网络里非凸的损失函数可能让梯度下降停在局部最优或鞍点这就是为什么深度学习里优化问题更复杂。3.4 学习率、迭代次数与特征缩放三个必须关心的调参项梯度下降有“三大命门”学习率、迭代次数、特征缩放。学习率 ( \alpha ) 的选择是最容易踩坑的地方。一个常用的办法是画损失函数随迭代次数变化的曲线如果损失在下降说明学习率合适如果损失震荡或上升说明学习率太大需要调小。实际中常见的取值是 0.01、0.03、0.1、0.3也可以先用非常小的值比如 0.001试跑观察曲线再逐步放大。没有一劳永逸的经验值只能多做实验。迭代次数决定算法跑多久。一种做法是设定一个固定次数比如 1000 次跑完后看损失是否收敛另一种是设定收敛条件——当相邻两次迭代的损失变化小于某个阈值比如 ( 10^{-5} )时就认为收敛并停止。我更推荐后者它可以避免白白浪费计算资源。特征缩放Feature Scaling是新手最容易忽略但影响极大的操作。如果某个特征的取值范围是 0~10000另一个是 0~1那么梯度下降在未缩放的损失曲面上会呈现出狭长的“碗形”导致参数更新时在一个方向上大幅度震荡、另一个方向上缓慢前进收敛非常慢。把每个特征缩放到相近的尺度比如均值归一化或标准化到 0 附近、方差为 1损失曲面会更接近正圆形收敛路径也直接得多。一句话总结正规方程适合特征少、数据量小的场景不用调参一步到位梯度下降适合特征多、数据量大的场景灵活可控但需要照顾好学习率、迭代次数和特征缩放这些细节。4. 从公式到代码用 NumPy 手写一个线性回归4.1 准备工作生成一份可复现的模拟数据为了让代码演示能直观地验证数学推导我先生成一份带噪声的线性数据。这里假设真实模型是 ( y 4 3x \epsilon )其中 ( \epsilon ) 是高斯噪声。使用numpy.random生成数据并设置随机种子保证每次运行结果一致方便调试复现。import numpy as np np.random.seed(42) # 生成 100 个样本特征 x 均匀分布在 [0, 2] X np.random.uniform(0, 2, size(100, 1)) y 4 3 * X np.random.randn(100, 1) * 0.5 # 可视化一下 import matplotlib.pyplot as plt plt.scatter(X, y, alpha0.7) plt.xlabel(x) plt.ylabel(y) plt.title(Simulated Linear Data) plt.show()这份数据只有 100 个样本、1 个特征非常适合用来同时跑通正规方程和梯度下降。加噪声是关键步骤——没有噪声的话数据完全落在一条直线上问题就退化成简单的解方程完全体现不出机器学习“从不确定中找规律”的乐趣。为了让矩阵运算方便需要给特征矩阵加一列全 1作为截距项对应的“哑特征”# 给 X 加一列 1作为截距项 X_b np.c_[np.ones((100, 1)), X]加了这一列之后模型就统一成了 ( h_\theta(x) \theta^T x )不需要再单独区分“截距”和“权重”这在代码实现里会省掉很多分支判断。4.2 用手写代码复现正规方程正规方程的代码简单到让人怀疑是不是真的有用# 正规方程θ (X^T X)^(-1) X^T y def normal_equation(X_b, y): theta np.linalg.inv(X_b.T X_b) X_b.T y return theta theta_ne normal_equation(X_b, y) print(正规方程结果: θ0 , theta_ne[0][0], , θ1 , theta_ne[1][0])跑完会发现结果非常接近真实的 ( \theta_04、\theta_13 )因为有噪声不会完全等于但已经很接近了。这个实现里最核心的一行就是np.linalg.inv(X_b.T X_b) X_b.T y和前面推导的正规方程公式一一对应。需要提醒的是np.linalg.inv在实际工程中并不推荐因为直接用inv求逆在数值上不如解线性方程组稳定。更专业的做法是用np.linalg.solve# 更稳定的做法直接解线性方程组 X^T X θ X^T y theta_ne np.linalg.solve(X_b.T X_b, X_b.T y)np.linalg.solve内部使用 LU 分解等方式求解线性方程组避免了直接求逆带来的数值误差和效率浪费。这也算是一个从“数学公式”到“工程代码”的典型差异公式里有逆矩阵但真正写代码时不一定要真的去算逆矩阵。4.3 用手写代码实现批量梯度下降接下来是梯度下降。需要做的事情有三件定义损失函数、计算梯度、迭代更新参数。写成如下代码def mse_loss(theta, X_b, y): m len(y) error X_b theta - y return float((1 / (2 * m)) * (error.T error)) def batch_gradient_descent(X_b, y, alpha0.1, n_iterations1000): m, n X_b.shape theta np.random.randn(n, 1) # 随机初始化 history [] for iteration in range(n_iterations): gradient (1 / m) * (X_b.T (X_b theta - y)) theta theta - alpha * gradient loss mse_loss(theta, X_b, y) history.append(loss) # 每 100 次迭代打印一次损失 if iteration % 100 0: print(fIteration {iteration}, Loss: {loss:.6f}) return theta, history # 标准化特征只标准化 x不标准化全为 1 的那列 def standardize(X): mean X.mean(axis0) std X.std(axis0) return (X - mean) / std, mean, std X_scaled, mean, std standardize(X) X_s np.c_[np.ones((100, 1)), X_scaled] theta_gd, history batch_gradient_descent(X_s, y, alpha0.1, n_iterations1000) print(梯度下降结果标准化后空间: θ0 , theta_gd[0][0], , θ1 , theta_gd[1][0])注意到这里我先把特征 ( x ) 做了标准化再在左边补一列 1。这一步不是可选项而是会让梯度下降收敛速度产生量级差距的关键操作。如果不做标准化学习率在 0.1 时可能直接导致损失变成 NaN或者收敛极慢标准化之后0.1 的学习率就能跑得很稳。训练完之后的损失曲线也值得看一眼它会呈现出典型的“快速下降、缓慢逼近”形态这是梯度下降收敛的正常模式。如果看到损失曲线震荡或上扬就得去检查学习率是不是太大。由于标准化的原因训练出的参数是在标准化特征空间里的。如果要把模型还原到原始特征空间需要做一个逆变换在标准化后的空间里模型是 ( \hat{y} \theta_0 \theta_1 x )其中 ( x (x - \mu) / \sigma )把它展开就能得到原始空间里的等效参数。# 把标准化空间中的参数还原到原始空间 theta1_orig theta_gd[1][0] / std[0] theta0_orig theta_gd[0][0] - theta1_orig * mean[0] print(f还原到原始空间: θ0 {theta0_orig:.4f}, θ1 {theta1_orig:.4f})跑完这段代码对比正规方程和还原后的梯度下降结果两个数字应该很接近。这其实是一个非常好的自检手段用正规方程的结果作为“标准答案”验证梯度下降实现是否正确。如果两者差距过大基本可以断定梯度下降的学习率、迭代次数或特征缩放处理出了问题。4.4 对比两种方法的效果与适用场景拿同一份数据分别跑正规方程和梯度下降数值结果几乎一样这验证了两种方法最终都在逼近同一个最优解。但在使用体验上两者区别明显对比维度正规方程批量梯度下降计算方式解析求解一步到位迭代逼近逐步优化是否需要调参不需要需要调学习率、迭代次数特征是否需要缩放不需要必须做特征缩放否则收敛慢特征数量很大时如 10 万维不可行矩阵求逆开销太大仍然可行每轮迭代 O(mn)样本数量很大时如 100 万样本构造 X^T X 可能耗尽内存可用随机/小批量梯度下降缓解数值稳定性直接求逆可能引入数值误差相对稳定取决于学习率选择表里的结论已经在实践中反复被验证小规模、特征少、一次性建模直接正规方程大规模、特征多、需要在线更新模型选择梯度下降或它的变体。当然现实工程中更多时候直接用sklearn的LinearRegression底层有 SVD 分解求解或SGDRegressor但理解背后的原理能帮你在模型跑出异常结果时快速定位问题。提示sklearn.linear_model.LinearRegression默认不是用正规方程求逆的它用的是最小二乘求解器底层通过 SVD 分解实现数值稳定性比直接求逆好很多。但它的数学本质和我们手写的正规方程是同一个解。5. 常见问题与排查技巧实录5.1 损失变成 NaN 或爆炸学习率太大还是特征没缩放这是梯度下降实现里最经典的问题。现象是训练过程中损失打印出来是nan或者第一次迭代的损失就巨大无比。先用排除法定位问题。第一步检查特征是否做了标准化——没标准化的特征量级可能从 0.001 到 100000梯度在某个方向会特别大更新一步就直接飞出去了。第二步检查学习率把学习率降到 0.001 或 0.0001 试跑如果问题消失说明是学习率太大的问题。第三步检查数据里有没有缺失值或无穷大这类脏数据会让矩阵运算直接产生 NaN。贴一个我自己的排查经验当时一个同事训练多项式回归损失前几轮在下降快到 100 轮时突然变成nan。查了半天发现不是学习率问题而是数据里有两行特征值特别大在梯度累积后数值溢出。最后做了标准化加裁剪clip问题立刻解决。梯度下降对数值范围非常敏感这是真实工程里最常见也最容易被忽略的坑。5.2 正规方程报错“Singular matrix”共线性或特征维度超样本量用np.linalg.inv或np.linalg.solve求解时如果数据存在完全共线性两个特征成比例或特征数量大于样本数量就会遇到奇异矩阵报错。解决思路有几个检查特征之间是否存在线性相关用np.corrcoef或画相关性矩阵热力图来排查删掉冗余的特征比如把“摄氏温度”和“华氏温度”同时放进模型这种低级错误改用岭回归它在 ( X^TX ) 上加了一个 ( \lambda I )保证矩阵可逆的同时还能约束参数大小是从数学上根治奇异性的方案。顺便提一句即使理论上矩阵可逆在数值计算中也可能因为同一个特征数量级差异过大导致“病态矩阵”计算出的参数极不稳定。这就得靠特征缩放和标准化来解决。5.3 多特征场景下参数含义发生变化先看尺度再看系数很多初学者在拿到多特征线性回归的参数后直接比较权重大小来判断特征重要性这其实很危险。当一个特征的单位是“平方米”数值几百上千另一个的单位是“房间数”数值单位个位数时前者的系数天然会很小后者的系数天然会很大但这不代表前者不重要。只有先把所有特征标准化到相同尺度再去比较系数的大小才有意义。所以在实际项目中如果要解释特征重要性一个标准流程是标准化→训练模型→比较标准化后系数的绝对值。要注意标准化的系数也要结合业务逻辑去验证不要盲目相信数字。5.4 拟合效果不好不全是模型的错先排查数据再升级模型如果你已经跑通了代码但发现模型在线性数据上的拟合效果很差先别急着换复杂模型。按这个顺序排查查看数据是否存在异常值一个离群点可能把回归线拉偏很多检查是否存在明显的非线性关系画个散点图看看确认特征是否真的和目标值相关数据里塞了完完全全的随机噪声特征模型学不到规律实属正常。线性回归拟合不好很多时候不是线性模型不行而是数据里本身就存在非线性的规律比如收入随年龄先升后降。这时候升级到多项式回归或引入交互特征才是正解。不过这是后话把线性回归的数学原理和代码实现吃透后面这些扩展都是顺手的事。写在后面的一点心得我自己刚学线性回归的时候也走过一段“背公式、调包、跑通就完事”的弯路。后来真正动手用 NumPy 从零实现一遍把每个公式展开、把每行代码和公式对应上才算是把最小二乘和梯度下降装进了脑子里。特别是矩阵求导那一步一开始总觉得“知道结论就行”但真正自己推一遍之后后面学到逻辑回归、多分类 Softmax 的梯度推导时明显顺畅了很多。再说一个建议不管你现在用的是 PyTorch、TensorFlow 还是 sklearn都值得在某个版本里用原生 NumPy 手写一遍线性回归。这个过程不需要很长时间但对你理解机器学习底层的优化机制、甚至后续排查模型不收敛问题帮助都是巨大的。数学公式看着枯燥但它们从代码里跑出结果的那一刻很多抽象的概念都会变得具体起来。这也是我把这篇文章的标题定为“从公式到代码实现”的原因——公式和代码本来就是一件事的两面。