先是一个简单的一元线性回归，理解它在求解什么？
\[ \begin{aligned} x_i &= i,\quad i = 1,2,\ldots,100, \\ y_i &= x_i + \varepsilon_i,\quad \varepsilon_i \sim \mathcal{N}(0,4). \end{aligned} \]
对应的代码和可视化代码：
set.seed(123)
x <- 1:100
y <- x + rnorm(100, 0, 4)
library(ggplot2)
ggplot(data.frame(x = x, y = y), aes(x, y)) +
geom_point(color = "steelblue", size = 2) +
geom_smooth(method = "lm", formula = y ~ x, color = "red", se = FALSE) +
labs( title = "Scatter Plot of x vs y")
对于熟悉 R 语言的数据科学工作者来说，下一步的操作几乎是肌肉记忆——直接调用 lm(y ~ x)，瞬间就能得到最优的 \(\beta_0\) 和 \(\beta_1\)。
fit <- lm(y ~ x)
library(stargazer)
stargazer(fit, type = 'text')
===============================================
Dependent variable:
---------------------------
y
-----------------------------------------------
x 1.010***
(0.013)
Constant -0.146
(0.737)
-----------------------------------------------
Observations 100
R2 0.985
Adjusted R2 0.985
Residual Std. Error 3.658 (df = 98)
F Statistic 6,352.370*** (df = 1; 98)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01然而，当我们跨入可微编程（Differentiable Programming）和深度学习的大门时，这个熟悉且优雅的“一行代码”突然消失了。取而代之的，是繁琐的循环、学习率设定以及梯度更新。
这种视角的转换，往往是许多传统统计建模者在初学深度学习时遇到的第一个“认知门槛”。为了弄清楚我们为什么要“舍近求远”，我们需要深入模型底层的求解逻辑。
在深入学习 R torch 的自动微分机制之前，我们需要先解决一个 R 语言老手最常见的困惑：我们为什么要通过编写循环、设置学习率来“训练”模型，而不是像调用 lm() 那样直接算出一个结果？
要回答这个问题，我们需要重新审视统计建模与深度学习在求解方法上的本质差异。
无论是经典的线性回归，还是深度神经网络，在回归任务中，我们的核心目标往往是一致的：最小化残差平方和（Residual Sum of Squares, RSS）。
也就是找到一组参数 \(w\)，使得预测值与真实值之间的误差最小 (Hastie 等 2009)：
\[ \text{Minimize } L(w) = \sum_{i=1}^{N} (y_i - \hat{y}_i)^2 = \sum_{i=1}^{N} (y_i - X_i w)^2 \]
既然目标一致，为什么解决问题的方法会发生如此剧烈的变化？这取决于我们的目标函数有多复杂。
在 lm() 的世界里，我们可以直接通过数学推导找到最优解。
回忆一下微积分：函数的极值点通常出现在导数为 0 的位置。对于线性回归 \(L(w)\) 这样一个完美的凸函数，我们可以直接对 \(w\) 求导，并令其等于 0：
\[ \frac{\partial L}{\partial w} = -2X^T(y - Xw) = 0 \]
解这个方程，我们不仅不需要迭代，反而能直接得出一个完美的公式，也就是正规方程 (Normal Equations)：
\[ w = (X^T X)^{-1} X^T y \]
这就是解析解 (analytical solution)，也被称为闭式解 (Closed-form solution) 。它就像是拥有了“瞬间移动”的能力——只要数据 \(X\) 和 \(y\) 给定，我们不需要试错，直接把数字代入公式，一步就能算出那个让 Loss 最小的完美的 \(w\)。这正是 R 语言中 lm()、solve() 背后的逻辑1：精确、直接、一步到位。
我们以 R 内置的 mtcars 数据集为例，通过汽车的车重 (wt) 和 马力 (hp) 来预测 油耗 (mpg)。
\[ \text{mpg} = w_1 \cdot \text{wt} + w_2 \cdot \text{hp} + b \]
为了确保后续手动实现的梯度下降算法能够快速收敛，数据标准化（Scaling） 是至关重要的一步。如果特征之间的尺度差异过大（例如 wt 约为 3，而 hp 约为 100），梯度下降会在寻找最优解时走弯路（Zigzag），导致训练极其困难。
# 选取特征与目标
x_raw <- mtcars[, c("wt", "hp")]
y_raw <- mtcars$mpg
# 标准化特征 (Z-score normalization)
x_scaled <- scale(x_raw)
# 设计矩阵：加上截距列
X <- cbind(1, x_scaled) # 1 是截距项
y <- y_raw
# OLS 解析解：beta = (X^T X)^(-1) X^T y
beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y
print(beta_hat)
# [,1]
# 20.090625
# wt -3.794292
# hp -2.178444
# 等价于实现了
# fit_lm <- lm(y_raw ~ x_scaled)
# print(coef(fit_lm))运行上述代码，我们可以得到截距项约为 20.09，wt 系数约为 -3.79，hp 系数约为 -2.18）。
既然解析解这么好，为什么深度学习不用？主要有两个原因阻断了这条路：
因此，深度学习被迫放弃了“瞬间移动”，选择了数值优化（Numerical Optimization）。
我们不再试图一步算出 \(w\)，而是采用迭代更新的策略。这就是著名的梯度下降 (Gradient Descent) 算法：
\[ w_{t+1} = w_t - \eta \cdot \nabla_w L(w_t) \]
这个公式描述了一个动态的修正过程：
整个训练过程变成了一个循环：
这种方法，需要我们手动推导均方误差（MSE）关于参数的导数公式：
\[ \frac{\partial L}{\partial w} = \frac{2}{N} X^T (Xw + b - y) \]
\[ \frac{\partial L}{\partial b} = \frac{2}{N} \sum (Xw + b - y) \]
写成 R 代码：
set.seed(123)
# w: 对应 wt 和 hp 的权重 (2行1列)
w <- matrix(rnorm(2), nrow = 2)
b <- 0
# 超参数
learning_rate <- 0.1
epochs <- 100
n_samples <- nrow(x_scaled) # 样本数量 N
# 开始训练循环
for (i in 1:epochs) {
# 1. 前向传播 (Forward Pass)
y_pred <- x_scaled %*% w + b # 使用 %*% 进行矩阵乘法
# 2. 计算误差 (Residuals)
error <- y_pred - y_raw
# 计算 Loss (仅用于观察，不参与运算)
loss <- mean(error^2)
# 3. 手动计算梯度 (The Hard Part!)
# 对应公式：(2/N) * X^T * error
w_grad <- (2 / n_samples) * (t(x_scaled) %*% error)
b_grad <- (2 / n_samples) * sum(error)
# 4. 更新参数 (Gradient Descent)
# 向梯度的反方向移动
w <- w - learning_rate * w_grad
b <- b - learning_rate * b_grad
if (i %% 10 == 0) {
cat(sprintf("Epoch %d: Loss = %.4f\n", i, loss))
}
}
cat(b, t(w))
# 20.09062 -3.793605 -2.179132运行这段代码，你会惊讶地发现，经过 100 轮简单的加减乘除，我们得到的系数与 lm() 的结果惊人地相似（w 约为 -3.7934 和 -2.1793）。
既然 Base R 也能做，为什么我们需要 torch？
请注意上面代码中的第 3 步——手动计算梯度。对于简单的线性回归，导数公式还是高中数学水平。但想象一下，如果我们要构建一个拥有 100 层、包含各种非线性激活函数（Sigmoid, ReLU, Tanh）的复杂网络，手动推导并编写每一个参数的导数公式将是一场噩梦。只要有一个公式推导错误，整个模型就会崩溃。
这就是深度学习框架存在的意义。接下来，我们将看到 Torch 最核心的魔法——自动微分（Automatic Differentiation）。它将把我们从繁琐的导数推导中彻底解放出来。
很多初学者最大的困惑在于：Torch 到底是怎么在“不知道”函数全貌的情况下算出导数的？
其实，自动微分（Autograd）并不是对着一个长长的公式 \(y = f(g(h(x)))\) 去推导全局导数。它采用的是“盲人摸象”策略：Torch 的 C++ 底层预先写好了所有基本算子（加法、乘法、指数、对数等）的“局部导数规则”。
Torch 不需要知道公式，它只需要把这一路上的“局部导数”连乘起来，最终结果就是我们要的全局梯度。回忆一下高中数学的定理——链式法则 (Chain Rule)。
对于复合函数 \(y = f(g(x))\)，如果我们令 \(u = g(x)\)，那么：
\[ \frac{dy}{dx} = \frac{dy}{du} \cdot \frac{du}{dx} \]
Torch 的 Autograd 系统本质上就是在这个链路中自动执行乘法运算。在深度学习框架中，每一次数学运算不仅是数值的计算，更是图结构的构建。R torch 采用的是动态计算图 (Dynamic Computational Graph)，即“运行时定义 (Define-by-Run)”模式。
从图论的角度看，神经网络是一个有向无环图 \(G = (V, E)\)：
当我们执行一段 R 代码时，torch 在后台同时维护着两套逻辑：
为了实现这一机制，torch_tensor 对象在 C++ 底层维护了一个关键指针：grad_fn。
grad_fn 为 NULL。让我们通过代码，像调试器一样逐层检查这个结构。
library(torch)
# 1. 定义叶子节点
# requires_grad = TRUE 是构建图的开关
# 只有开启它，torch 才会分配内存去记录后续的操作
w <- torch_tensor(2.0, requires_grad = TRUE)
b <- torch_tensor(5.0, requires_grad = TRUE)
x <- torch_tensor(3.0) # 输入数据通常不需要梯度
# 2. 执行运算（构建图）
# 步骤 A: 乘法
u <- w * x
# 步骤 B: 加法
v <- u + b
# 步骤 C: 平方
loss <- v ^ 2
# 3. 计算图结构分析：追踪 grad_fn
# Level 0: 输出层
print(loss$grad_fn)
# 输出: <PowBackward0>
# 解释: 这是一个“幂运算”的反向函数对象
# Level 1: 中间层
# 我们可以通过 next_functions 访问图的上一层
# 注意：这是 R torch 暴露出来的底层接口，用于调试图结构
# loss 是由 v 算出来的，所以它的上一级应该是 AddBackward (v 的生成逻辑)
print(loss$grad_fn$next_functions)
# 输出 list 中包含: <AddBackward0>
# Level 2: 更深层
# AddBackward0 (v) 是由 u 和 b 算出来的
# 所以它的 next_functions 应该包含 MulBackward0 (u) 和 AccumulateGrad (b)
v_grad_fn <- loss$grad_fn$next_functions[[1]]$fn
print(v_grad_fn$next_functions)在执行前向计算 loss <- v^2 的瞬间，torch 并不是只算出了 121 这个数，而是在内存中生成了一个 PowBackward0 对象，并让这个对象持有了一个指向 AddBackward0 的指针。
这一系列对象构成了一个单向链表。当我们调用 loss$backward() 时，torch 的引擎只是简单地遍历这个链表，依次调用每个节点的 Apply 方法计算偏导数。
因为是动态图，所以 R 语言的 if 或 for 可以用来改变网络结构的拓扑结构：
dynamic_op <- function(x, w) {
if (as.numeric(x$mean()) > 0) {
# 路径 A: 包含 Log 和 Mul 运算
return(torch_log(x) * w)
} else {
# 路径 B: 包含 Abs 和 Add 运算
return(torch_abs(x) + w)
}
}
w <- torch_tensor(1.0, requires_grad = TRUE)
y1 <- dynamic_op(torch_tensor(2.0), w) # 输入为正
y1$grad_fn$print()
# MulBackward0
y2 <- dynamic_op(torch_tensor(-2.0), w) # 输入为负
y2$grad_fn$print()
# AddBackward0在 R torch 中，不存在一个“固定的”全局图。每次前向传播，都是在创建一张全新的图。这虽然带来了一定的 R 运行时开销，但换来了极致的灵活性，使得处理变长序列（RNN）和图神经网络（GNN）变得符合直觉。
在深度学习的训练循环中，梯度的生命周期遵循严格的“三部曲”：开启追踪 -> 反向传播 -> 梯度清零。
一、开启追踪 (requires_grad)
在 R torch 中，张量被创建时默认是不分配梯度内存的。这是为了最大化工程效率——输入数据（Features）和标签（Labels）通常不需要更新，如果对它们也分配梯度内存，显存占用将翻倍。
我们必须显式地声明哪些张量是“参数（Parameters）”。
# 普通数据张量 (Default)
x <- torch_tensor(3.0)
# 模型参数张量
w <- torch_tensor(2.0, requires_grad = TRUE)
# 或者后期原地修改 (In-place)
b <- torch_tensor(5.0)
b$requires_grad_(TRUE) 当你使用 nn_linear 或 nn_conv2d 等标准层时，torch 已经自动将其内部的 weight 和 bias 设置为了 requires_grad = TRUE。
二、反向传播(backward)
loss$backward() 是 Autograd 引擎的启动键。需要特别注意的是，这个函数没有返回值。
在 R 的常规函数中，我们习惯 grad <- calculate_grad(loss)。但在 torch 中，梯度是被“写入”到计算图叶子节点的 $grad 属性中的。我们可以看一个例子,设函数 \(y = 3x^2 + 2\)，当 \(x=2\) 时：
x <- torch_tensor(2.0, requires_grad = TRUE)
(y <- 3 * x^2 + 2) # 前向传播，打印 y
# torch_tensor
# 14
# [ CPUFloatType{1} ][ grad_fn = <AddBackward1> ]
y$backward() # 遍历计算图，计算梯度，并将结果写入 x$grad
print(x$grad) # 检查反向传播后的结果
# torch_tensor
# 12
# [ CPUFloatType{1} ]为了让你彻底放心 torch 的数学能力，我们来看一个多元微积分例子。假设函数 \(z = x^2 + y^3\)。我们要计算在点 \((x=2, y=3)\) 处的梯度。
x <- torch_tensor(2.0, requires_grad = TRUE)
y <- torch_tensor(3.0, requires_grad = TRUE)
z <- x^2 + y^3
z$backward()
cat("x 的梯度 (2x -> 4): ", x$grad$item(), "\n")
cat("y 的梯度 (3y^2 -> 27): ", y$grad$item(), "\n")看到代码输出精确地吻合数学推导，是不是有一种掌控感？我们完全不需要像第一节一样，将每个数学推导都来一遍，这正是“可微编程（Differentiable Programming）”的魅力所在。
三、梯度的累加机制(Accumulation)
这是 R torch 设计中最反直觉的地方。当你再次调用 backward() 时，新计算出的梯度不会覆盖旧的梯度，而是会“加”在旧梯度上。即：
\[ \text{grad}_{\text{new}} = \text{grad}_{\text{old}} + \frac{\partial L}{\partial w} \]
可以做个小实验：
library(torch)
# 定义一个变量 w = 2
w <- torch_tensor(2.0, requires_grad = TRUE)
# 假设函数 y = 3 * w
# 导数 dy/dw = 3
# 理论上，无论我们算多少次，导数都应该是 3
y1 <- 3 * w # 第 1 次反向传播
y1$backward()
cat("第 1 次 backward 后，w$grad =", w$grad$item(), "\n")
y2 <- 3 * w # 第 2 次反向传播，我们没有清空梯度，直接再算一次
y2$backward()
cat("第 2 次 backward 后，w$grad =", w$grad$item(), "\n")
y3 <- 3 * w # 第 3 次反向传播
y3$backward()
cat("第 3 次 backward 后，w$grad =", w$grad$item(), "\n")返回的结果：
# 第 1 次 backward 后，w$grad = 3
# 第 2 次 backward 后，w$grad = 6
# 第 3 次 backward 后，w$grad = 9 但这并非设计缺陷，而是为了工程上的灵活性。假如 GPU 显存塞不下 Batch Size = 64 的数据，我们可以把数据切成 4 个 Batch Size = 16 的小块，分别 Backward，梯度会自动累加，最后 Update 一次。这等价于训练了 Batch Size = 64。
这种机制意味着，在常规的训练循环中，我们必须手动干预，显式地把梯度归零。否则，第 100 轮训练时的梯度，将包含前 99 轮所有梯度的总和，这会导致更新步长过大，模型瞬间发散。
在底层 API 中，我们需要调用张量的原地操作 zero_()：
# 手动清空梯度
w$grad$zero_()
cat("清空后，w$grad =", w$grad$item(), "\n")
# 输出: 0在 R torch 中，只要输入张量的 requires_grad = TRUE，系统就会默认你正在进行模型训练，从而不惜一切代价记录操作历史。但在验证（Validation）或预测（Inference）阶段，我们需要暂停这个机制，这就是 with_no_grad 上下文管理器的作用。 在 with_no_grad 代码块内，即便声明了张量的 requires_grad = TRUE，torch 也不会为它们构建计算图。
w <- torch_randn(100, 100, requires_grad = TRUE)
x <- torch_randn(100, 100)
# 训练模式
start_time <- Sys.time()
for (i in 1:5000) {
y <- torch_matmul(w, x)
}
Sys.time() - start_time
# Time difference of 0.2006299 secs
# 推断模式
start_time <- Sys.time()
with_no_grad({
for (i in 1:5000) {
y <- torch_matmul(w, x)
}
})
Sys.time() - start_time
# Time difference of 0.06900597 secs推断模式下不构建计算图，所以速度会更快。同 with_no_grad() 经常一起出现的还有 model$eval()，它的主要作用是改变层的数学行为：比如训练的时候有 dropout，但 eval 时就不要做 dropout 动作了。
另外，R 用户非常习惯使用 list 或者 data.frame 来存储循环中的结果，但在 torch 中会造成内存泄露。
history <- list()
# 假设我们在训练循环中
for (i in 1:1000) {
# Forward ...
loss <- (w * x)^2 # 假设这是一个复杂的计算图终点
# [错误做法]
# history[[i]] <- loss
# loss 对象身上挂着整个计算图（grad_fn 链表）。
# 只要 list 引用了 loss，整个图及其所有中间变量就无法被垃圾回收（GC）。
# 跑几千轮后，内存将被撑爆。
# [正确做法 A]：只存数值
# item() 将 0-dim 张量转换为 R 的原生数值 (numeric)
history[[i]] <- loss$item()
# [正确做法 B]：剥离张量
# detach() 创建一个新的张量，数据共享，但切断了与图的联系
detached_loss <- loss$detach()
# history[[i]] <- detached_loss
}如果你需要将 Tensor 转化为 R 的原生数据（如用于 ggplot2 绘图），或者需要将 Tensor 存起来供下一轮使用（如 RNN 的隐藏状态），务必使用 detach() 或 item()。
回到前面的案例，我们重新用张量实现一遍。
首先我们需要定义模型参数 \(W\) 和 \(b\)，并告诉 Torch：“请帮我盯着这几个变量，我要对它们求导。”这就是 requires_grad = TRUE 的用武之地。
torch_manual_seed(123) # torch 的随机数种子
# 数据转换
x_tensor <- torch_tensor(x_scaled, dtype = torch_float())
y_tensor <- torch_tensor(y_raw, dtype = torch_float())
y_tensor <- y_tensor$view(c(-1, 1)) # 转为 (N, 1) 形状
# 权重 W (形状: 2行1列)
# 开启 requires_grad = TRUE，让 Autograd 记录对 W 的所有操作
w <- torch_randn(c(2, 1), dtype = torch_float(), requires_grad = TRUE)
# 偏置 b (形状: 1)
b <- torch_zeros(c(1, 1), dtype = torch_float(), requires_grad = TRUE)
print(w) # 看下 w 长什么样
# torch_tensor
# -0.1115
# 0.1204
# [ CPUFloatType{2,1} ][ requires_grad = TRUE ]设置学习率和迭代次数
learning_rate <- 0.1
epochs <- 100前面开启了梯度追踪，紧跟着是反向传播 loss$backward() 和梯度清零：
for (epoch in 1:epochs) {
# Forward，建立计算图：y = XW + b
y_pred <- torch_matmul(x_tensor, w) + b
# 计算均方误差 (MSE)
loss <- (y_pred - y_tensor)$pow(2)$mean()
# 反向传播 (Backward Pass) - Autograd
# 自动计算 loss 关于 w 和 b 的梯度，并存储在 w$grad 和 b$grad 中。
loss$backward()
# 参数更新 (Update Weights)
# 必须在 with_no_grad() 中进行，防止更新操作被记录到计算图中
with_no_grad({
# 使用原地减法 sub_ (substract in-place)
w$sub_(learning_rate * w$grad)
b$sub_(learning_rate * b$grad)
# 手动清空梯度！否则下一轮的梯度会叠加在这一轮的梯度上
w$grad$zero_()
b$grad$zero_()
})
if (epoch %% 10 == 0) {
cat(sprintf("Epoch: %d, Loss: %.4f\n", epoch, loss$item()))
}
}对比解析解、R 原生梯度下降、Torch_Autograd 三个版本的结果。
| Term | 解析解 | R 原生梯度下降 | Torch_Autograd |
|---|---|---|---|
| Intercept | 20.090625 | 20.09062 | 20.0906 |
| wt | -3.794292 | -3.793605 | -3.7936 |
| hp | -2.178444 | -2.179132 | -2.1792 |
三种方法的计算结果基本没什么差异。
仔细观察上一节的代码，你会发现参数更新的步骤略显繁琐：我们需要手动通过 with_no_grad 暂停追踪，手动做减法更新，还要手动清空梯度。
# 手动更新的繁琐写法
with_no_grad({
w$sub_(learning_rate * w$grad)
b$sub_(learning_rate * b$grad)
w$grad$zero_()
b$grad$zero_()
})如果模型有 100 层，每层都有权重和偏置，手动写这几行代码将是灾难。
R torch 提供了一个更高层的抽象——优化器 (Optimizer)。优化器本质上是一个通过“托管”参数来自动化更新流程的对象。它将上述的样板代码封装成了两个标准方法：
optimizer$step(): 替代了手动的减法更新。optimizer$zero_grad(): 替代了手动的 grad$zero_()。让我们用 optim_sgd（随机梯度下降优化器）来重构之前的线性回归代码。你会发现，无论模型多么复杂，训练循环的代码结构永远保持不变。
# 定义优化器
optimizer <- optim_sgd(params = list(w, b), lr = learning_rate)
# 训练循环
for (epoch in 1:epochs) {
y_pred <- torch_matmul(x_tensor, w) + b
loss <- (y_pred - y_tensor)$pow(2)$mean()
# 1. 清空过往梯度
# 必须在 backward 之前调用，否则梯度会累加
optimizer$zero_grad()
# 2. 反向传播
# 计算当前梯度，填入 w$grad 和 b$grad
loss$backward()
# 3. 执行更新
# 根据 w$grad 和学习率，自动修正 w 的值
optimizer$step()
if (epoch %% 10 == 0) {
cat(sprintf("Epoch: %d, Loss: %.4f\n", epoch, loss$item()))
}
}这段代码展示了深度学习工程中最通用的范式。此时，optimizer 就像是一个仅仅执行命令的工兵，我们告诉它用最简单的 SGD 算法（仅仅是减去梯度）。
在传统的统计建模中，最基础的假设是可加性（Additivity）。但现实数据往往充满复杂的交互效应（Interaction Effects）。
考虑一个经典的异或（XOR）风格的数据分布：
这种关系在数学上可以近似为乘法关系 \(y \propto x_1 \cdot x_2\)。此时，\(x_1\) 对 \(y\) 的边际效应取决于 \(x_2\) 的正负——这就是典型的交互效应。
看一下 R 代码的实现：
library(tidyverse)
set.seed(42)
# 1. 生成模拟数据 (N = 1000)
n_samples <- 1000
x1 <- runif(n_samples, min = -1, max = 1)
x2 <- runif(n_samples, min = -1, max = 1)
# 定义真实的数据生成机制：y = x1 * x2 + 噪声
# 注意：这是一个纯粹的交互关系，没有截距，也没有单独的线性项
y_true <- x1 * x2
noise <- rnorm(n_samples, mean = 0, sd = 0.05)
y <- y_true + noise
# 整理为 Tibble 以便分析
df <- tibble(
x1 = x1,
x2 = x2,
y = y
)可视化：颜色代表 y 的值
ggplot(df, aes(x = x1, y = x2, color = y)) +
geom_point(alpha = 0.6) +
scale_color_gradient2(low = "blue", mid = "white", high = "red") +
theme_minimal()
如果我们尝试使用标准的线性回归模型的加法模型和交互模型来拟合这份数据，结果会如何？
# 线性回归 (加法模型)
model_linear <- lm(y ~ x1 + x2, data = df)
summary(model_linear)$r.squared
# [1] 0.00160... (完全失败，因为平面无法拟合马鞍面)
# 交互回归 (手动特征工程)
model_interaction <- lm(y ~ x1 + x2 + x1:x2, data = df)
# 也可写作 x1*x2, 等价于 x1:x2
summary(model_interaction)$r.squared
# [1] 0.9234... (完美拟合)虽然 x1:x2 或 x1*x2 能解决问题，但在深度学习面临的高维数据（如图像像素）中，比如输入是一张 \(28 \times 28\) 的图像（784 个像素），我们要构建交互项，理论上有 \(C_{784}^2\) 种两两组合，甚至更多的高阶组合。手动构造特征在维数灾难面前是不可行的。
我们希望构建一个神经网络，只输入原始的 \(x_1, x_2, ...\)，让网络自动学会“相乘”这个非线性关系。
初学者常有一个误区：既然单层网络能力有限，那我多堆叠几层（Deep Neural Network）不就变强了吗？
比如这个两层网络(没有激活函数)：
# 这是一个没有激活函数的两层网络
net <- nn_sequential(
nn_linear(2, 10), # 第一层：把 2 维变成 10 维
nn_linear(10, 1) # 第二层：把 10 维变成 1 维
)如果你训练这个网络，你会发现 Loss 永远降不下去，预测结果和 lm(y ~ x1 + x2) 一样糟糕。为什么？
让我们做一个简单的数学推导：
我们将第一层的公式代入第二层： \[ y = W_2 (W_1 x + b_1) + b_2 = (W_2 W_1) x + (W_2 b_1 + b_2) \]
如果我们令新的权重矩阵 \(W' = W_2 W_1\)，新的偏置 \(b' = W_2 b_1 + b_2\)，那么整个网络就变成了：
\[ y = W' x + b' \]
结论令人震惊： 两个线性层的叠加，在数学上严格等价于一个线性层。
这意味着，无论你堆叠了 100 层还是 1000 层线性网络，它的表达能力（Capacity）并不比第 1 节里的那个单层线性回归强。这就是所谓的“线性塌缩”。为了打破这种诅咒，我们需要在每一层的输出后面插入一个非线性变换 \(\sigma\)：
\[ y = W_2 \cdot \sigma(W_1 x + b_1) + b_2 \]
正是这个非线性的 \(\sigma\)，让神经网络拥有了“弯曲空间”的能力，从而能够逼近交互项。
我们在代码中加入 nn_relu() 再试一次：
net <- nn_sequential(
nn_linear(2, 10),
nn_relu(), # <--- 关键的非线性层
nn_linear(10, 1)
)训练后，这个网络能完美拟合交互数据（400 轮迭代后 R-squared 为 0.9232）。网络在没有被显式告知“相乘”规则的情况下，自动学会了模拟 \(x_1 \times x_2\) 的关系。
虽然它本质上是把 \(x_1 \times x_2\) 的曲面看作是由无数个小平面（ReLU 的折线）拼凑出来的，但在效果上，它自动学会了交互。
一、Sigmoid，非常漂亮的 S 型曲线
\[ \sigma(x) = \frac{1}{1 + e^{-x}} \]
能将任意实数压缩到 \((0, 1)\) 区间。它模拟了生物神经元的“饱和”特性——输入极小是 0（抑制），极大是 1（激活），中间平滑过渡。
但它最大的问题是，当输入 \(x\) 很大或很小时，Sigmoid 曲线非常平坦，导数趋近于 0。这会导致反向传播时，梯度信号在传递几层后就消失了，深层网络无法更新。而且输出恒为正，非零中心化。所以它几乎不在隐藏层使用。主要用于二分类任务的输出层（将 Logits 转化为概率）。
二、Tanh，零中心化的改进版
\[ \tanh(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}} \]
能将输入压缩到 \((-1, 1)\) 区间。它是零中心化的（Zero-centered），收敛速度通常更快。但依然存在梯度消失问题。在循环神经网络（RNN/LSTM）中依然是标配。
三、ReLU (Rectified Linear Unit)，简单粗暴的王者 \[ \text{ReLU}(x) = \max(0, x) \]
正数保持不变，负数直接置零。优点是计算极快，只有比较和赋值，没有昂贵的指数运算。
在 \(x > 0\) 的区域，导数恒为 1。这意味着无论网络多深，只要神经元处于激活状态，梯度信号就可以无损地传回前端。负输入输出为 0（神经元“死掉”），这在某种程度上引入了稀疏表达，模拟了生物大脑的节能特性。是目前深度学习中隐藏层默认的首选激活函数。
library(torch)
library(tidyverse)
# 定义辅助函数
get_data <- function(name, act_fn) {
x <- torch_linspace(-5, 5, 200, requires_grad = TRUE)
y <- act_fn(x) # 向前传播
y$sum()$backward() # 反向传播
# 提取数据
data.frame(
x = as_array(x$detach()),
y = as_array(y$detach()),
grad = as_array(x$grad),
func = name
)
}
df_sigmoid <- get_data("Sigmoid", torch_sigmoid)
df_tanh <- get_data("Tanh", torch_tanh)
df_relu <- get_data("ReLU", torch_relu)
# 合并数据
df_all <- bind_rows(df_sigmoid, df_tanh, df_relu) %>%
pivot_longer(cols = c("y", "grad"), names_to = "type", values_to = "value")
ggplot(df_all, aes(x = x, y = value, color = func)) +
geom_line(linewidth = 0.9) +
facet_grid(type ~ func, scales = "free_y") +
scale_color_manual(values = c(
"Sigmoid" = "#1f77b4",
"Tanh" = "#2ca02c",
"ReLU" = "#d62728"
)) +
theme_bw() +
labs( y = "Value / Gradient")
运行上述代码，重点观察第二行（grad）的图像：
在编写 R torch 代码时，你会发现激活函数有两种调用方式。初学者容易混淆，但它们各有用途。
方式一：函数式 (Functional API)
位于 torch_relu(), torch_sigmoid() 等。 它们是纯粹的函数，输入一个张量，输出一个张量，不保存任何状态。这个场景用的最多。
library(torch)
x <- torch_tensor(c(-2, -1, 0, 1, 2))
out <- torch_relu(x) # 直接调用函数
# 输出: [0, 0, 0, 1, 2]适用场景：当你自己在 forward 方法中手写复杂的运算逻辑时。
方式二：模块式 (Modular API)
位于 nn_relu(), nn_sigmoid() 等。 它们是 R6 类（对象）。你需要先实例化它，把它变成模型层的一部分2。
relu_layer <- nn_relu() # 实例化一个层对象
out <- relu_layer(x) # 调用它 (内部会自动调用 forward)当我们把线性层和非线性激活函数交替堆叠，并且至少引入一个中间层（隐藏层）时，我们就构建了深度学习中最基础、最经典的模型——多层感知机（Multilayer Perceptron, MLP）。
从数学结构上看，MLP 的本质就是复合函数。如果说单层线性网络只能拟合直线（或超平面），而激活函数提供了非线性的弯曲能力，那么 MLP 就是通过层层嵌套，获得逼近任意复杂函数的能力。假设有一个 13 个输入特征的数据集，一个“输入层 -> 隐藏层 -> 输出层”的 MLP 数学表达如下：
\[ \text{Price} = W_2 \cdot \underbrace{f(\underbrace{W_1 \cdot x + b_1}_{\text{线性变换}})}_{\text{隐藏特征 } h} + b_2 \]
其中：
构建 MLP 最直观的方法，就是把这些数学运算像搭积木一样按顺序拼起来。R torch 提供了 nn_sequential 容器，它完全对应了数学上的复合函数逻辑。我们需要：
代码如下3：
model <- nn_sequential(
nn_linear(13, 32), # 输入13个特征 -> 32个神经元
nn_relu(), # 非线性关系的层
nn_linear(32, 16), # 32 -> 16
nn_relu(), # 激活函数
nn_linear(16, 1) # 16 -> 输出1个预测值
)
print(model) # 打印模型结构在回归任务中，目标变量是一个连续的实数，我们希望输出范围不受限制，因此最后一层通常保持线性输出。
nn_module在后续章节的实战中，我们将面对更复杂的训练需求。为了更好地管理代码和参数，我们需要掌握 torch 的核心范式：nn_module。
在 R torch 中，这通过 R6 类系统实现。定义一个模型类需要遵循两个核心步骤：
让我们把上面的房价预测模型重构为一个标准的 R6 类：
net <- nn_module(
classname = "net",
# 1. 构造函数：定义模型组件
initialize = function(in_dim = 13, h1 = 32, h2 = 16) {
self$fc1 <- nn_linear(in_dim, h1) # 输入层
self$fc2 <- nn_linear(h1, h2) # 隐藏层
self$fc3 <- nn_linear(h2, 1) # 输出层
},
# 2. 前向传播：定义数据流向
forward = function(x) {
x <- self$fc1(x) # 第一个线性层
x <- torch_relu(x) # relu 层
x <- self$fc2(x) # 第二个线性层
x <- torch_relu(x) # relu 层
x <- self$fc3(x) # 回归任务保持线性输出
x
}
)
# 实例化模型
model <- net(in_dim = 13, h1 = 32, h2 = 16)最后一步实例化结束后，就可以在 update 模块里调用 net 这个模型了。
虽然 nn_sequential 很好用，但作为工程师，我们需要更灵活的控制权，类的写法为未来留出了扩展空间，比如我们可能需要：
在前面的章节中，我们反复提到“损失函数”的概念——它是衡量模型预测值与真实值之间差异的函数。如果说模型结构定义了如何计算预测，那么损失函数则定义了什么是好的预测。
对于回归问题（预测连续值），最常用的损失函数是均方误差（Mean Squared Error, MSE）：
\[ L_{\text{MSE}} = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2 \]
其中 \(y_i\) 是真实值，\(\hat{y}_i\) 是预测值。MSE 的数学特性非常友好：
很多情况下会把 MSE 开平方，得到 RMSE（Root Mean Squared Error，均方根误差）作为评价指标。它与原始目标变量具有相同的量纲，更便于解释模型的平均预测偏差。
但 MSE 对异常值（Outliers）比较敏感，因为平方会放大异常值的影响。作为替代，平均绝对误差（MAE）对异常值更鲁棒：
\[ L_{\text{MAE}} = \frac{1}{n} \sum_{i=1}^n |y_i - \hat{y}_i| \]
对于分类问题（预测离散类别），最常用的是交叉熵损失（Cross-Entropy Loss）。它衡量的是预测概率分布与真实分布之间的差异。
二分类交叉熵损失（Binary Cross-Entropy）： \[ L_{\text{BCE}} = -\frac{1}{n} \sum_{i=1}^n \left[ y_i \log(\hat{y}_i) + (1 - y_i) \log(1 - \hat{y}_i) \right] \] 其中 \(y_i \in \{0, 1\}\) 是真实标签，\(\hat{y}_i \in (0, 1)\) 是预测为正类的概率（通常经过 Sigmoid 激活）。
多分类交叉熵损失（Categorical Cross-Entropy）： \[ L_{\text{CCE}} = -\frac{1}{n} \sum_{i=1}^n \sum_{j=1}^C y_{ij} \log(\hat{y}_{ij}) \] 其中 \(C\) 是类别数，\(y_{ij}\) 是样本 \(i\) 属于类别 \(j\) 的 one-hot 编码，\(\hat{y}_{ij}\) 是预测概率（通常经过 Softmax 激活）。
交叉熵损失有几个重要特性：
R torch 实现：
| 任务类型 | 输出层激活 | 损失函数 | 函数 |
|---|---|---|---|
| 回归 | 无（或线性） | MSE | nn_mse_loss() |
| 回归 | 无 | MAE | nn_l1_loss() |
| 二分类 | Sigmoid | 二元交叉熵 | nn_bce_with_logits_loss() |
| 多分类 | Softmax | 分类交叉熵 | nn_cross_entropy_loss() |
不同的损失函数会产生不同的梯度特性，直接影响优化过程：
MSE 的梯度： \[ \frac{\partial L_{\text{MSE}}}{\partial \hat{y}_i} = \frac{2}{n} (\hat{y}_i - y_i) \] 梯度与误差成正比，线性关系。
交叉熵的梯度（对于 Sigmoid 输出）： \[
\frac{\partial L_{\text{BCE}}}{\partial \hat{y}_i} = \frac{\hat{y}_i - y_i}{\hat{y}_i(1 - \hat{y}_i)}
\] 当预测接近 0 或 1 时，分母很小，梯度会变得非常大。这可能导致训练不稳定，这也是为什么我们常用 nn_bce_with_logits_loss 的原因——它在数值上更稳定。
常见业务环境下“排序有关的损失函数”也很重要，例如 nn_margin_ranking_loss 和 Bayesian Personalized Ranking 损失等，篇幅所限，不做详细讨论。
本节我们将 modeldata(Kuhn 2025) 中的 Ames Housing 数据集为基础，探讨如何利用 torch 的模块化设计构建一个回归模型。
改数据集有：
modeldata::ames 中）Sale_Price（房屋最终成交价）变量（不完整）信息：
| 类别 | 示例变量 (英文) | 主要说明 |
|---|---|---|
| 标识与分区 | PID, MS_SubClass, MS_Zoning |
地块标识、建筑类型、区域规划 |
| 地块与位置 | Lot_Area, Neighborhood, Condition1 |
占地面积、所属社区、邻近条件 |
| 建筑与质量 | Overall_Qual, Year_Built, Year_Remod |
整体质量评分、建成年份、改建年份 |
| 外观与结构 | Exter_Qual, Roof_Style, Foundation |
外墙质量、屋顶样式、地基类型 |
| 生活空间 | Gr_Liv_Area, Full_Bath, TotRmsAbvGrd |
地上居住面积、卫生间数量、总房间数 |
| 地下室 | Total_Bsmt_SF, Bsmt_Qual |
地下室总面积、地下室质量 |
| 车库 | Garage_Type, Garage_Cars, Garage_Area |
车库类型、车库容量、车库面积 |
| 设施与功能 | Fireplaces, Heating_QC, Central_Air |
壁炉数量与质量、供暖质量、中央空调 |
| 户外设施 | Wood_Deck_SF, Open_Porch_SF, Pool_Area |
木平台、开放式门廊、泳池面积 |
| 销售信息 | Sale_Type, Sale_Condition, Yr_Sold |
销售类型、销售条件、售出年份/月份 |
| 目标变量 | SalePrice |
房屋最终售价 |
首先是数据预处理，利用 tidymodels 包中的函数将数据分为训练集和测试集
library(tidymodels)
library(torch)
library(tidyverse)
data(ames)
set.seed(123)
torch_manual_seed(123)
ames_split <- initial_split(ames, prop = 0.80)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)定义数据处理的 pipeline，依次是
Sale_Price 以外的全部，数据源为 ames_train。以上 10 个操作步骤对应的 R 代码：
rec <- recipe(Sale_Price ~ ., data = ames_train) %>%
step_normalize(all_outcomes()) %>%
step_rm(matches("Id|Longitude|Latitude")) %>%
step_log(Gr_Liv_Area, Lot_Area, base = 10) %>%
step_other(all_nominal_predictors(), threshold = 0.01) %>%
step_dummy(all_nominal_predictors()) %>%
step_nzv(all_predictors()) %>%
step_normalize(all_predictors()) %>%
step_impute_median(all_numeric_predictors()) %>%
step_impute_mode(all_nominal_predictors())这里的目标变量执行了 Normalization 操作，在预测结束之后还需要做逆向操作。即预测值加上 ames_train 数据集的目标变量均值，再乘以 ames_train 数据集的目标变量的方差。
上面的函数定义了数据处理流程，但并不执行。接下来开始执行：
prep_rec <- prep(rec)
train_df <- bake(prep_rec, new_data = NULL)
test_df <- bake(prep_rec, new_data = ames_test)train_df 就是自己，但 test_df 的均值方差、中位数、众数要从 ames_train 继承，这样才能保证信息没有泄露给 test_df 数据集。
接着将训练集和测试集数据转为 tonsor：
x_train <- torch_tensor(
as.matrix(train_df %>% select(-Sale_Price)), dtype = torch_float())
y_train <- torch_tensor(
matrix(train_df$Sale_Price, ncol=1), dtype = torch_float())
x_test <- torch_tensor(
as.matrix(test_df %>% select(-Sale_Price)), dtype = torch_float())
y_test <- torch_tensor(
matrix(test_df$Sale_Price, ncol=1), dtype = torch_float())构建一个多层感知机，它的结构如下：
输入层 (ncol(x_train)) -> 隐藏层 (128, ReLU) -> 隐藏层 (64, ReLU) -> 输出层 (1)
结构并不复杂，我们使用 nn_sequential 来实现该模型（读者也可以尝试改成 nn_module 来实现）：
model <- nn_sequential(
nn_linear(ncol(x_train), 128), # 输入特征 -> 128个神经元
nn_relu(), # 激活函数
nn_linear(128, 64), # 128 -> 64
nn_relu(), # 激活函数
nn_linear(64, 1) # 64 -> 输出1个预测值
)定义优化器和损失函数：
learning_rate <- 0.0025
optimizer <- optim_adam(model$parameters, lr = learning_rate)
loss_fn <- nn_mse_loss()接着开始循环训练：
epochs <- 1000
# 用于记录 loss 变化以便画图
train_losses <- numeric(epochs)
for (t in 1:epochs) {
# 前向传播
y_pred <- model(x_train)
loss <- loss_fn(y_pred, y_train)
# 梯度求解三部曲
optimizer$zero_grad();loss$backward();optimizer$step()
train_losses[t] <- loss$item() # 用于记录与打印
}可以看到模型的 loss 下降很快：
plot(train_losses, type = "l", col = "blue", lwd = 2,
xlab = "Epoch", ylab = "MSE Loss")
手工计算 RMSE
# 切换到评估模式 (虽然本例没用 Dropout，但养成好习惯)
model$eval()
with_no_grad({
preds_nn <- as.numeric(model(x_test))
})
# 计算 RMSE
truth <- as.numeric(y_test) # 真实标签
rmse_nn <- sqrt(mean((truth - preds_nn)^2))对比原生 lm 方法
truth <- as.numeric(y_test) # 真实标签
rmse_nn <- sqrt(mean((truth - preds_nn)^2))
cat(sprintf("Neural Network RMSE: %.4f\n", rmse_nn))Neural Network RMSE: 0.3635# 测试 lm 使用与神经网络完全相同的数据集
fit_lm <- lm(Sale_Price ~ ., data = train_df)
preds_lm <- predict(fit_lm, newdata = test_df)
# 计算 RMSE (均方根误差)
rmse_lm <- sqrt(mean((test_df$Sale_Price - preds_lm)^2))
cat(sprintf("Linear Regression RMSE: %.4f\n", rmse_lm))Linear Regression RMSE: 0.3914即便是在我们盲猜参数没有任何优化的情况下，MLP 因为引入了非线性变化，对于 RMSE 这个指标而言，还是略胜一筹。
打印真实值和预测值的散点图做两个方法的对比：
results <- tibble(
Truth = truth,
Linear_Model = preds_lm,
Neural_Network = preds_nn
) %>%
pivot_longer(cols = c(Linear_Model, Neural_Network),
names_to = "Model", values_to = "Prediction")
library(ggplot2)
ggplot(results, aes(x = Truth, y = Prediction, color = Model)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(
values = c("Linear_Model" = "black", "Neural_Network" = "red")) +
labs(
x = "True price", y = "Pred price"
) +
theme_minimal() +
coord_obs_pred() # 保证 x 和 y 轴比例一致