4.1 什么是多元线性回归

从一个特征到多个特征

多元线性回归就是把一元线性回归推广到多个特征的情形，代价函数和优化逻辑不变。

模型从$\hat{y}=wx+b$ 变成：

$\hat{y}=w_1{x}_1+w_2{x}_2+\cdots +w_n{x}_n+b$

其中 n 是特征数量，
$x1,x2,\dots ,xn$ 是每个特征的值，$w_1,w_2,\dots ,w_n$ 是对应的权重，b 依然是偏置。每个权重 $w_j$
表示"在其他特征不变的情况下，第 j 个特征增加一个单位，预测值会增加多少"。

几何含义

一元线性回归在二维平面上是一条直线。多元线性回归在高维空间里是一个超平面。两个特征时是三维空间里的一张平面，三个以上特征时已经无法可视化，但数学结构完全一样。这也是为什么第4章要引入矩阵——当特征数量很多时，逐个写出每个 $w_j$ 非常繁琐，矩阵提供了一种紧凑优雅的表示方式。

参数数量的变化

一元线性回归有2个参数（w 和 b）。n 个特征的多元线性回归有n+1 个参数（n 个权重加一个偏置）。参数多了，梯度下降就需要同时对每个参数计算偏导并更新，手写会很麻烦。矩阵运算可以把这些操作打包成简洁的公式，一次性处理所有参数，这是下一节之后要做的事情。

多元线性回归的代价函数

代价函数的形式和一元时完全一样，依然是MSE：

$J(w_1,\dots ,w_n,b)=\frac{1}{m}{\sum }_{i=1}^m({\hat{y}}^{(i)}-y^{(i)}{)}^2$

只是现在 ${\hat{y}}^{(i)}$ 是由多个特征共同计算出来的。损失函数的形式没变，优化的逻辑也没变，变的只是参数的数量和计算的复杂度。

多元回归带来的新挑战

特征多了之后会出现一个一元回归没有的问题：不同特征的量纲和数值范围差异很大。比如房屋面积可能在50到300之间，而房龄可能在1到50之间，楼层可能在1到30之间。这种差异会让损失曲面变成一个极度扁平的椭圆形，梯度下降收敛极慢，就像第10节里看到的那种震荡情形。解决方案是特征缩放，后续会涉及到。

举例如下：

数据集预览（前5条样本）：
      面积       楼层       房龄         价格
   106.2      1.9     32.5      324.6
   192.6     19.5      5.1      646.0
   159.8     10.1      8.9      579.1
   139.8     15.7     45.0      446.0
    73.4     27.3     30.7      238.3

面积范围：50.8 ~ 198.0
楼层范围：1.2 ~ 29.6
房龄范围：1.2 ~ 49.5

注意三个特征的数值范围差异很大，这正是需要特征缩放的原因

4.2 特征缩放和归一化

特征缩放

参数更新公式是

$w_j:=w_j-\alpha \cdot \frac{\partial J}{\partial w_j}$

由于梯度大小直接受特征数值大小影响，面积特征数值大，对应的梯度会很大，则参数$w_1$每一步的更新幅度很大，而房龄更新幅度小。结果损失曲面在不同方向上坡度不均衡，变成极度拉伸的椭圆，梯度下降会剧烈震荡，收敛慢。

特征缩放是在训练前把所有特征的数值范围统一到同一个尺度，让损失曲面变得更“圆”，梯度下降可以平稳快速的收敛。

两种常用的缩放方法

方法一：Min-Max归一化，把每个特征的值线性压缩到[0, 1] 区间：

$x^{\prime }=\frac{x-x_{min}}{x_{max}-x_{min}}$

直觉上就是：原来最小值变成0，最大值变成1，其他值按比例落在中间。优点是结果范围固定，缺点是对异常值（极端大或极端小的值）很敏感，一个异常值会把所有数据都压缩到一个很窄的范围里。

方法二：标准化（Z-score Normalization），把每个特征变换成均值为0、标准差为1的分布：

$x^{\prime }=\frac{x-\mu }{\sigma }$

其中 $\mu$ 是该特征的均值，$\sigma$ 是标准差。直觉上就是：把特征的中心移到0，然后用标准差作为单位来衡量每个值偏离中心多远。它对异常值的鲁棒性更好，在深度学习里使用更广泛。

需要注意的是：计算$\mu$ 和$\sigma$（或者$x_{min}$ 和$x_{max}$）时，只能用训练集的数据来计算，然后把同样的数值应用到验证集和测试集上，不能用测试集自己的统计量来缩放自己。
原因是：测试集模拟的是真实世界中模型从未见过的新数据。如果用测试集自己的统计量，相当于提前"偷看"了测试数据的信息，评估结果就不再真实可靠了。这个原则在机器学习里叫做避免数据泄露（Data Leakage），违反它会导致你以为模型表现很好，但部署到真实场景后效果很差。

import numpy as np

np.random.seed(42)
m = 100

# 三个特征：面积、楼层、房龄
x1 = np.random.uniform(50, 200, m)
x2 = np.random.uniform(1, 30, m)
x3 = np.random.uniform(1, 50, m)
X = np.column_stack([x1, x2, x3])  # 把三列合并成一个矩阵，shape: (100, 3)
y = 3*x1 + 2*x2 - 1.5*x3 + 50 + np.random.randn(m) * 20

# 划分训练集和测试集
split = int(0.8 * m)
X_train, X_test = X[:split], X[split:]
y_train, y_test = y[:split], y[split:]

# ===== 标准化：只用训练集的均值和标准差 =====
mu    = X_train.mean(axis=0)   # 每个特征的均值，shape: (3,)
sigma = X_train.std(axis=0)    # 每个特征的标准差，shape: (3,)

# 用训练集统计量缩放训练集和测试集
X_train_scaled = (X_train - mu) / sigma
X_test_scaled  = (X_test  - mu) / sigma   # 注意：用的是训练集的mu和sigma

print("缩放前，训练集各特征的范围：")
for j, name in enumerate(['面积', '楼层', '房龄']):
    print(f"  {name}：{X_train[:,j].min():.1f} ~ {X_train[:,j].max():.1f}")

print("\n缩放后，训练集各特征的范围：")
for j, name in enumerate(['面积', '楼层', '房龄']):
    print(f"  {name}：{X_train_scaled[:,j].min():.2f} ~ {X_train_scaled[:,j].max():.2f}")

print(f"\n缩放后各特征均值（应接近0）：{X_train_scaled.mean(axis=0).round(4)}")
print(f"缩放后各特征标准差（应接近1）：{X_train_scaled.std(axis=0).round(4)}")

训练结果为

缩放前，训练集各特征的范围：
  面积：50.8 ~ 198.0
  楼层：1.2 ~ 29.6
  房龄：1.2 ~ 49.5

缩放后，训练集各特征的范围：
  面积：-1.52 ~ 1.72
  楼层：-1.61 ~ 1.76
  房龄：-1.88 ~ 1.52

缩放后各特征均值（应接近0）：[-0. -0.  0.]
缩放后各特征标准差（应接近1）：[1. 1. 1.]

从参差不齐的状态，统一变成了大约-2 到+2 的相近范围，均值接近0，标准差接近1。这就是标准化的效果。且处理后模型的预测能力完全不受影响，但$w_j$的数值含义会发生变化，不再是原始特征的直接参数，如果需要解读现实焊锡，可以把缩放操作，反推还原参数。一般情况下只关心预测结果是否准确，不解读参数具体含义。
总体来说，特征缩放，让损失曲面更圆，梯度下降更稳；计算缩放统计量只能用训练集，测试集用同样的统计量来变换，避免数据泄露。

4.3 Pytorch实现线性回归

3.7节使用纯Numpy实现线性回归，这一节用PyTorch再实现一遍。
Pytorch的实现方式是神经网络训练的标准写法，Numpy版本可以帮助理解原理，Pytorch把许多细节封装好了，代码更加简洁，也更容易拓展。

Pytorch训练的标准结构

PyTorch把训练流程拆成了几个固定的组件，每个组件都有对应的模块：用 nn.Linear 定义模型，用 nn.MSELoss 定义损失函数，用 torch.optim 里的优化器来更新参数。这三个组件组合在一起，构成了PyTorch训练的标准骨架，后面每一章都是在这个骨架上填入不同的内容。
nn.Linear

nn.Linear(in_features, out_features) 是PyTorch内置的线性层，它封装了$\hat{y}=xw+b$ 的计算，自动创建并管理w 和b 两个参数，也自动把它们标记为 requires_grad=True，不需要你手动处理。括号里的两个数字分别是输入特征数和输出数量，对于一元线性回归就是nn.Linear(1, 1)，意思是1个输入、1个输出。

import torch
import torch.nn as nn
import numpy as np

# ===== 模块一：准备数据 =====
np.random.seed(42)
m = 50
x_np = np.linspace(0, 10, m)
y_np = 2.5 * x_np + 5 + np.random.randn(m) * 3

# 转换成PyTorch张量，并调整形状为(m,1),这里-1是自动计算，前面讲过
# reshape(-1,1)把一维数组变成列向量，因为nn.Linear要求输入是二维的
x = torch.tensor(x_np, dtype=torch.float32).reshape(-1, 1)
y = torch.tensor(y_np, dtype=torch.float32).reshape(-1, 1)

# ===== 模块二：定义模型 =====
# nn.Linear(1,1)：1个输入特征，1个输出
# 它自动创建w和b，并追踪梯度
model = nn.Linear(1, 1)

# ===== 模块三：定义损失函数和优化器 =====
criterion = nn.MSELoss()                           # MSE损失函数
optimizer = torch.optim.Adam(model.parameters(),   # 用Adam优化器
                             lr=0.1)               # 学习率

# ===== 模块四：训练循环 =====
loss_history = []

for epoch in range(500):
    # 第一步：前向传播，计算预测值
    y_pred = model(x)

    # 第二步：计算损失
    loss = criterion(y_pred, y)
    loss_history.append(loss.item())   # .item()把张量转成普通Python数字

    # 第三步：清零梯度（每次反向传播前必须清零，否则梯度会累加）
    optimizer.zero_grad()

    # 第四步：反向传播，自动计算所有梯度
    loss.backward()

    # 第五步：优化器更新参数
    optimizer.step()

    if epoch % 100 == 0:
        print(f"Epoch {epoch:4d} | Loss={loss.item():.4f}")

# ===== 模块五：查看结果 =====
# model.weight和model.bias就是学到的w和b
w_learned = model.weight.item()
b_learned = model.bias.item()
print(f"\n训练结果：w={w_learned:.4f}，b={b_learned:.4f}")
print(f"真实参数：w=2.5000，b=5.0000")

运行结果如下：

Epoch    0 | Loss=231.7257
Epoch  100 | Loss=8.5531
Epoch  200 | Loss=7.6832
Epoch  300 | Loss=7.4608
Epoch  400 | Loss=7.4291
Epoch  500 | Loss=7.4264

训练结果：w=2.3296，b=5.1703
真实参数：w=2.5000，b=5.0000

训练循环的五个步骤：前向传播—计算损失—清零梯度—反向传播—更新参数。
NumPy版本需要手动写梯度公式 grad_w = np.mean(error * x_train)，PyTorch版本只需要 loss.backward() 一行，梯度自动算好。NumPy版本需要手动写参数更新 w = w - lr * grad_w，PyTorch版本只需要 optimizer.step() 一行。模型越复杂，这种自动化带来的优势就越大

4.4 矩阵和矩阵的运算

省略矩阵和矩阵乘法、矩阵转置等基本含义

用Pytorch来做矩阵运算

import torch

A = torch.tensor([[1., 2., 3.],
                  [4., 5., 6.]])          # shape: (2, 3)

B = torch.tensor([[1., 2.],
                  [3., 4.],
                  [5., 6.]])              # shape: (3, 2)

# 矩阵乘法：用@符号，A是2×3，B是3×2，结果是2×2
C = A @ B
print(f"A的形状：{A.shape}")             # torch.Size([2, 3])
print(f"B的形状：{B.shape}")             # torch.Size([3, 2])
print(f"A@B的形状：{C.shape}")           # torch.Size([2, 2])
print(f"A@B =\n{C}")

# 转置
print(f"\nA的转置形状：{A.T.shape}")     # torch.Size([3, 2])
print(f"A^T =\n{A.T}")

# 逐元素乘法（不是矩阵乘法，要求形状相同）
A2 = torch.tensor([[1., 2.], [3., 4.]])
B2 = torch.tensor([[5., 6.], [7., 8.]])
print(f"\n逐元素乘法：\n{A2 * B2}")     # 不是矩阵乘法，对应位置相乘

# 验证矩阵乘法不满足交换律
print(f"\nA@B =\n{A @ B}")
# B@A的形状是3×3，和A@B的2×2完全不同，说明顺序不能交换
print(f"B@A =\n{B @ A}")

运行结果如下：

A的形状：torch.Size([2, 3])
B的形状：torch.Size([3, 2])
A@B的形状：torch.Size([2, 2])
A@B =
tensor([[22., 28.],
        [49., 64.]])

A的转置形状：torch.Size([3, 2])
A^T =
tensor([[1., 4.],
        [2., 5.],
        [3., 6.]])

逐元素乘法：
tensor([[ 5., 12.],
        [21., 32.]])

A@B =
tensor([[22., 28.],
        [49., 64.]])
B@A =
tensor([[ 9., 12., 15.],
        [19., 26., 33.],
        [29., 40., 51.]])

4.5 矩阵表示线性回归

设计矩阵X和参数向量w

把所有样本的特征组织成一个矩阵 X，每一行是一个样本，每一列是一个特征：

$X=\left(\begin{array}{cccc}{x}_1^{(1)}& x_2^{(1)}& \cdots & x_n^{(1)}\\ x_1^{(2)}& x_2^{(2)}& \cdots & x_n^{(2)}\\ ⋮& ⋮& \ddots & ⋮\\ x_1^{(m)}& x_2^{(m)}& \cdots & x_n^{(m)}\end{array}\right)\text{形状：}m\times n$

把所有权重组织成一个列向量 $\mathbf{w}$
$\mathbf{w}=\left(\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_n\end{array}\right)\text{形状：}n\times 1$

做一次矩阵乘法 $X\mathbf{w}$，形状是$(m\times n)\times (n\times 1)=m\times 1$，正好得到所有 m 个样本的预测值组成的列向量。这一步矩阵乘法同时完成了所有样本的计算，比逐个循环高效得多。

把偏置b融入矩阵

现在还有一个偏置b 需要处理。常见的做法是在 X 的最左边加一列全为1的列，同时把 b 也并入参数向量 $\mathbf{w}$，这样偏置就被"吸收"进了矩阵乘法：

$\tilde{X}=\left(\begin{array}{cccc}1& x_1^{(1)}& \cdots & x_n^{(1)}\\ 1& x_1^{(2)}& \cdots & x_n^{(2)}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_1^{(m)}& \cdots & x_n^{(m)}\end{array}\right),\tilde{\mathbf{w}}=\left(\begin{array}{c}b\\ w_1\\ w_2\\ ⋮\\ w_n\end{array}\right)$

于是整个预测公式变成一行：

$\hat{\mathbf{y}}=\tilde{X}\tilde{\mathbf{w}}$

形状是$(m\times (n+1))\times ((n+1)\times 1)=m\times 1$，干净利落。这种把偏置融入矩阵的技巧在很多论文和教材里都会用到，看到"在特征矩阵加一列1"你就知道它在做什么了。

代价函数的矩阵形式

代价函数也可以用矩阵写成一行。残差向量是 $\hat{\mathbf{y}}-\mathbf{y}$，MSE就是残差向量的每个元素平方后求平均：

$J(\mathbf{w})=\frac{1}{m}(\hat{\mathbf{y}}-\mathbf{y}{)}^T(\hat{\mathbf{y}}-\mathbf{y})$

其中 $(\hat{\mathbf{y}}-\mathbf{y}{)}^T(\hat{\mathbf{y}}-\mathbf{y})$是向量和自身的点积，结果正好等于每个元素的平方和。这个写法在论文里非常常见，看到它你就知道是在算MSE。

import numpy as np

np.random.seed(42)
m = 100   # 样本数
n = 3     # 特征数

# 生成数据：三个特征
x1 = np.random.uniform(50, 200, m)
x2 = np.random.uniform(1, 30, m)
x3 = np.random.uniform(1, 50, m)
y  = 3*x1 + 2*x2 - 1.5*x3 + 50 + np.random.randn(m) * 20

# 构建特征矩阵X，形状(m, n) = (100, 3)
X = np.column_stack([x1, x2, x3])

# 在X左边加一列全1，把偏置融入矩阵，形状变成(100, 4)
ones = np.ones((m, 1))
X_tilde = np.hstack([ones, X])           # shape: (100, 4)

# 参数向量：[b, w1, w2, w3]，随机初始化
w_tilde = np.array([0.0, 0.0, 0.0, 0.0])  # shape: (4,)

# 矩阵形式的前向传播：一行代替m个循环
y_pred = X_tilde @ w_tilde               # shape: (100,)

# 矩阵形式的MSE
residual = y_pred - y                    # 残差向量，shape: (100,)
loss = residual @ residual / m           # 等价于np.mean(residual**2)

print(f"X_tilde的形状：{X_tilde.shape}")   # (100, 4)
print(f"w_tilde的形状：{w_tilde.shape}")   # (4,)
print(f"y_pred的形状： {y_pred.shape}")    # (100,)
print(f"初始MSE损失：{loss:.2f}")

运行结果为：

X_tilde的形状：(100, 4)
w_tilde的形状：(4,)
y_pred的形状： (100,)
初始MSE损失：183086.90

这一节的矩阵表示方法是理解神经网络的直接铺垫。神经网络的每一层本质上就是一次$X\mathbf{w}+b$的矩阵运算，只是后面还跟着一个激活函数。你在论文里看到的模型公式几乎全是矩阵形式，能熟练读懂$\hat{\mathbf{y}}=X\mathbf{w}+b$ 这种写法，后面的推导就会顺畅很多。

4.6 最小二乘法及其推导

最小二乘法

最小二乘法（Least Squares）就是第3.6节提到的"数学方法"在多元线性回归下的完整版本。它的目标和梯度下降完全一样——最小化MSE损失——但它不是迭代逼近，而是直接用矩阵运算一步算出最优参数。"最小二乘"这个名字来自它的核心思想：最小化残差的平方和（二乘就是平方的意思）。

推导过程

有了上一节的矩阵表示，推导就非常自然。代价函数的矩阵形式是：

$J(\mathbf{w})=\frac{1}{m}(\tilde{X}\mathbf{w}-\mathbf{y}{)}^T(\tilde{X}\mathbf{w}-\mathbf{y})$

要找到让 J 最小的 $\mathbf{w}$，令梯度等于零：

$\frac{\partial J}{\partial \mathbf{w}}=\frac{2}{m}{\tilde{X}}^T(\tilde{X}\mathbf{w}-\mathbf{y})=0$

化简这个方程：

${\tilde{X}}^T\tilde{X}\mathbf{w}={\tilde{X}}^T\mathbf{y}$

两边左乘 $({\tilde{X}}^T\tilde{X}{)}^{-1}$，得到正规方程（Normal Equation）：

${\mathbf{w}}^{\ast }=({\tilde{X}}^T\tilde{X}{)}^{-1}{\tilde{X}}^T\mathbf{y}$

这就是目录里标注的那个公式 $(X^{T}X{)}^{-1}{X}^{T}y$。代入数据矩阵，直接算出最优参数，不需要任何迭代。

逐步理解这个公式的结构

这个公式乍看很复杂，但拆开来每一部分都有清晰的含义。${\tilde{X}}^T\tilde{X}$的形状是$(n+1)\times (n+1)$
，它捕捉了特征之间的相关性结构。$({\tilde{X}}^T\tilde{X}{)}^{-1}$ 是对这个结构求逆，相当于"消除特征相关性的影响"。${\tilde{X}}^T\mathbf{y}$的形状是 $(n+1)\times 1$，它捕捉了每个特征和标签之间的关联程度。把两者相乘，就得到了在当前数据结构下，让预测误差最小的参数向量。

什么时候正规方程会失效

正规方程并不是万能的，有两种情况会出问题。

第一种是 ${\tilde{X}}^T\tilde{X}$不可逆的情况。这发生在特征之间存在完全线性相关时，比如你有"面积（平方米）"和"面积（平方英尺）"两个特征，它们完全成比例，这时${\tilde{X}}^T\tilde{X}$ 是奇异矩阵，无法求逆。解决方法是删掉冗余特征，或者加入正则化项（后面章节会涉及）。

第二种是特征数量太多。求矩阵逆的计算复杂度是 $O(n^3)$，当特征数 n 很大时（比如图像识别里一张图片有几万个像素作为特征），计算量大得无法接受。这时只能用梯度下降。

import numpy as np

np.random.seed(42)
m = 100   # 样本数

# 生成三个特征的数据
x1 = np.random.uniform(50, 200, m)
x2 = np.random.uniform(1,  30,  m)
x3 = np.random.uniform(1,  50,  m)
y  = 3*x1 + 2*x2 - 1.5*x3 + 50 + np.random.randn(m) * 20

# 构建加了一列1的特征矩阵，shape: (100, 4)
X_tilde = np.hstack([np.ones((m, 1)),
                     np.column_stack([x1, x2, x3])])

# 正规方程：w* = (X^T X)^{-1} X^T y
# 逐步计算，方便理解每一步的形状
XtX     = X_tilde.T @ X_tilde          # shape: (4, 4)
XtX_inv = np.linalg.inv(XtX)           # shape: (4, 4)，求逆
Xty     = X_tilde.T @ y                # shape: (4,)
w_star  = XtX_inv @ Xty                # shape: (4,)，最优参数

print("正规方程求解结果：")
print(f"  b  = {w_star[0]:.4f}，真实值 = 50.0000")
print(f"  w1 = {w_star[1]:.4f}，真实值 =  3.0000")
print(f"  w2 = {w_star[2]:.4f}，真实值 =  2.0000")
print(f"  w3 = {w_star[3]:.4f}，真实值 = -1.5000")

# 验证：计算训练集上的MSE
y_pred = X_tilde @ w_star
mse    = np.mean((y_pred - y)**2)
print(f"\n训练集MSE：{mse:.4f}")

运行结果为：

正规方程求解结果：
  b  = 38.9483，真实值 = 50.0000
  w1 = 3.0555，真实值 =  3.0000
  w2 = 2.1438，真实值 =  2.0000
  w3 = -1.3656，真实值 = -1.5000

训练集MSE：371.3162

误差来自数据噪声

最小二乘法是线性回归的"专属捷径"，在特征数不多、特征不冗余的情况下，一步给出精确解，不需要调任何超参数。梯度下降是通用优化工具，适用于任何模型，是深度学习的基石。理解最小二乘法的推导，最大的价值不是"以后都用它"，而是让你深刻理解"令梯度为零求解"这个思路，以及为什么在神经网络里这条路走不通。

4.7 Python实现最小二乘法

三种实现方式

求解线性回归最优参数有三种代码层面的路径：手动实现正规方程、用NumPy内置的最小二乘求解器、用PyTorch实现。三者结果应该完全一致，但抽象层次不同。

import numpy as np
import torch
import torch.nn as nn

np.random.seed(42)
m = 100

# 生成数据
x1 = np.random.uniform(50, 200, m)
x2 = np.random.uniform(1,  30,  m)
x3 = np.random.uniform(1,  50,  m)
y  = 3*x1 + 2*x2 - 1.5*x3 + 50 + np.random.randn(m) * 20

X_tilde = np.hstack([np.ones((m, 1)),
                     np.column_stack([x1, x2, x3])])  # shape: (100, 4)

# ===== 方式一：手动实现正规方程 =====
# 直接翻译公式 w* = (X^T X)^{-1} X^T y
w_manual = np.linalg.inv(X_tilde.T @ X_tilde) @ X_tilde.T @ y

print("===== 方式一：手动正规方程 =====")
print(f"  b={w_manual[0]:.4f}  w1={w_manual[1]:.4f}  w2={w_manual[2]:.4f}  w3={w_manual[3]:.4f}")

# ===== 方式二：NumPy内置求解器 =====
# np.linalg.lstsq是"least squares"的缩写，专门求解最小二乘问题
# 比手动求逆更数值稳定（即使X^TX接近奇异也能处理）
w_numpy, _, _, _ = np.linalg.lstsq(X_tilde, y, rcond=None)

print("\n===== 方式二：NumPy lstsq =====")
print(f"  b={w_numpy[0]:.4f}  w1={w_numpy[1]:.4f}  w2={w_numpy[2]:.4f}  w3={w_numpy[3]:.4f}")

# ===== 方式三：PyTorch训练=====
import torch
import torch.nn as nn

X_torch = torch.tensor(np.column_stack([x1, x2, x3]), dtype=torch.float32)
y_torch = torch.tensor(y, dtype=torch.float32).reshape(-1, 1)

model     = nn.Linear(3, 1)
criterion = nn.MSELoss()
optimizer = torch.optim.Adam(model.parameters(), lr=1.0)  # 原始数据量纲大，学习率调大

for epoch in range(5000):  # 多跑几步让它充分收敛
    optimizer.zero_grad()
    loss = criterion(model(X_torch), y_torch)
    loss.backward()
    optimizer.step()

w_result = model.weight.detach().numpy().flatten()
b_result = model.bias.item()

print("===== 方式三：PyTorch训练=====")
print(f"  b={b_result:.4f}  w1={w_result[0]:.4f}  w2={w_result[1]:.4f}  w3={w_result[2]:.4f}")

print("\n===== 真实参数 =====")
print(f"  b=50.0000  w1=3.0000  w2=2.0000  w3=-1.5000")

运行结果

===== 方式一：手动正规方程 =====
  b=38.9483  w1=3.0555  w2=2.1438  w3=-1.3656

===== 方式二：NumPy lstsq =====
  b=38.9483  w1=3.0555  w2=2.1438  w3=-1.3656
===== 方式三：PyTorch训练=====
  b=38.9484  w1=3.0555  w2=2.1439  w3=-1.3655

===== 真实参数 =====
  b=50.0000  w1=3.0000  w2=2.0000  w3=-1.5000

方式一手动求逆np.linalg.inv(X_tilde.T @ X_tilde)在理论上完全正确，但在工程上有一个隐患：当

${\tilde{X}}^T\tilde{X}$ 接近奇异（行列式接近零）时，直接求逆会产生很大的数值误差，结果不稳定。np.linalg.lstsq内部用的是SVD分解（奇异值分解）来求解，即使矩阵接近奇异也能得到稳定的结果。在实际工程中应该优先用方式二，方式一主要用于理解原理。