案例 1：从 NumPy 广播到能源负荷线性回归#
学习目标与记号#
本案例根据建筑的相对紧凑度、表面积、墙体面积、屋顶面积、高度、朝向、玻璃面积和玻璃分布预测供暖负荷，并把“读取数据—核对数组—定义损失—估计参数—分析误差”组织成一个完整实验。这样既能复习 NumPy 的基本用法，也能理解这些命令在线性回归中的实际作用。
本案例用普通小写字母表示标量，用粗体小写字母表示向量，用粗体大写字母表示矩阵；\(m\) 表示样本量，\(n\) 表示特征数量，批量样本按行存放。主要学习目标如下：
区分样本轴和特征轴，主动检查设计矩阵 \(\boldsymbol{X}\in\mathbb{R}^{m\times n}\) 与观测标签向量 \(\boldsymbol{y}\in\mathbb{R}^{m\times1}\) 的维度；
解释广播为何能把偏置 \(b\) 加到整批预测上，并识别维度不合适却仍能运行的错误广播；
在同一个训练集上比较最小二乘解析解与批量梯度下降的结果；
从均方误差推导向量化梯度，并用中心差分做数值检查；
结合测试集上的均方根误差（RMSE）、决定系数 \(R^2\) 和残差图判断模型的预测表现与局限。
对应正文： “Python 基本命令”“NumPy 数组的广播”“一元变量的线性回归”以及“线性回归与逻辑回归”。
数据来源、许可与任务边界#
数据说明页：https://archive.ics.uci.edu/dataset/242/energy%2Befficiency
匿名下载直链：https://archive.ics.uci.edu/static/public/242/energy+efficiency.zip
许可：CC BY 4.0
数据规模：768 个建筑设计、8 个输入特征和 2 个连续观测标签。本案例只预测供暖负荷。
下载日期基准：2026-08-12。下载函数会显示文件的 SHA-256 摘要；正式发布时可以把该摘要记录在
catalog.json中，以便核对文件是否发生变化。
这是教学用回归案例。它说明变量之间的统计关系，不应被直接当成实际建筑设计规范或能耗认证工具。
# 统一环境、随机种子与下载缓存
from pathlib import Path
import hashlib
import os
import shutil
import tempfile
import urllib.request
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
SEED = 42
rng = np.random.default_rng(SEED)
CACHE_ROOT = Path(os.environ.get(
"AI_COURSE_DATA_DIR",
Path.home() / ".cache" / "ai-course-cases",
)).expanduser()
CACHE_ROOT.mkdir(parents=True, exist_ok=True)
def download(url, filename, sha256=None):
"""匿名下载到独立缓存；已有文件仍会计算摘要，避免在没有提示的情况下使用错误数据。"""
destination = CACHE_ROOT / filename
if not destination.exists():
request = urllib.request.Request(url, headers={"User-Agent": "ai-course-case/1.0"})
with urllib.request.urlopen(request, timeout=120) as response:
with tempfile.NamedTemporaryFile(dir=CACHE_ROOT, delete=False) as tmp:
shutil.copyfileobj(response, tmp)
temporary = Path(tmp.name)
temporary.replace(destination)
digest = hashlib.sha256(destination.read_bytes()).hexdigest()
if sha256 is not None and digest != sha256:
destination.unlink(missing_ok=True)
raise ValueError(f"SHA-256 不匹配：{digest}")
print(f"{destination.name}：{destination.stat().st_size / 1024**2:.2f} 兆字节，SHA-256={digest}")
return destination
plt.rcParams["figure.dpi"] = 120
print("缓存目录：", CACHE_ROOT)
缓存目录： /private/tmp/ai-course-case-data
第 1 步：下载并读取数据#
压缩包中的 Excel 文件是数据的原始发布格式。程序先查看压缩包中的文件名，再读取其中唯一的 .xlsx 文件，因此不依赖某台计算机上的绝对路径。重命名各列后，立即检查数据类型、缺失值和重复行。是否删除重复设计应由研究问题决定：本案例保留原始记录，因为不同朝向和玻璃分布本来就是实验设计的一部分。
import io
import zipfile
URL = "https://archive.ics.uci.edu/static/public/242/energy+efficiency.zip"
archive = download(URL, "uci_energy_efficiency.zip")
with zipfile.ZipFile(archive) as zf:
members = [name for name in zf.namelist() if name.lower().endswith(".xlsx")]
if len(members) != 1:
raise RuntimeError(f"预期 1 个 xlsx，实际为 {members}")
raw = pd.read_excel(io.BytesIO(zf.read(members[0])))
columns = [
"relative_compactness", "surface_area", "wall_area", "roof_area",
"overall_height", "orientation", "glazing_area",
"glazing_area_distribution", "heating_load", "cooling_load",
]
column_labels = {
"relative_compactness": "相对紧凑度", "surface_area": "表面积",
"wall_area": "墙体面积", "roof_area": "屋顶面积",
"overall_height": "建筑高度", "orientation": "朝向",
"glazing_area": "玻璃面积", "glazing_area_distribution": "玻璃分布",
"heating_load": "供暖负荷", "cooling_load": "制冷负荷",
}
raw = raw.iloc[:, :10].copy()
raw.columns = columns
raw = raw.apply(pd.to_numeric, errors="coerce").dropna()
print("数据维度：", raw.shape)
display(raw.head().rename(columns=column_labels))
display(raw.isna().sum().rename(index=column_labels).to_frame("缺失数量").T)
assert raw.shape[1] == 10 and raw["heating_load"].notna().all()
uci_energy_efficiency.zip：0.07 兆字节，SHA-256=499441eee27929a4b00417f58fd8c63c9cc14b8a71520cd0dd27fcb626738351
数据维度： (768, 10)
| 相对紧凑度 | 表面积 | 墙体面积 | 屋顶面积 | 建筑高度 | 朝向 | 玻璃面积 | 玻璃分布 | 供暖负荷 | 制冷负荷 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.98 | 514.5 | 294.0 | 110.25 | 7.0 | 2 | 0.0 | 0 | 15.55 | 21.33 |
| 1 | 0.98 | 514.5 | 294.0 | 110.25 | 7.0 | 3 | 0.0 | 0 | 15.55 | 21.33 |
| 2 | 0.98 | 514.5 | 294.0 | 110.25 | 7.0 | 4 | 0.0 | 0 | 15.55 | 21.33 |
| 3 | 0.98 | 514.5 | 294.0 | 110.25 | 7.0 | 5 | 0.0 | 0 | 15.55 | 21.33 |
| 4 | 0.90 | 563.5 | 318.5 | 122.50 | 7.0 | 2 | 0.0 | 0 | 20.84 | 28.28 |
| 相对紧凑度 | 表面积 | 墙体面积 | 屋顶面积 | 建筑高度 | 朝向 | 玻璃面积 | 玻璃分布 | 供暖负荷 | 制冷负荷 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 缺失数量 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
第 2 步：先理解变量，再训练模型#
描述统计可以帮助比较各个特征的数值尺度。例如，面积可能达到数百，而相对紧凑度小于 1；若直接进行梯度下降，数值较大的特征可能对参数更新产生过强影响。相关矩阵只反映两个变量的线性相关程度，不能据此判断因果关系。下面同时绘制供暖负荷的分布，检查是否存在明显的异常值、取值截断或多个峰。
feature_names = columns[:8]
description = raw[feature_names + ["heating_load"]].describe().T
description = description.rename(index=column_labels, columns={
"count": "样本数", "mean": "均值", "std": "标准差",
"min": "最小值", "25%": "下四分位数", "50%": "中位数",
"75%": "上四分位数", "max": "最大值",
})
display(description)
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
raw["heating_load"].hist(bins=25, ax=axes[0])
axes[0].set(title="供暖负荷分布", xlabel="供暖负荷", ylabel="频数")
corr = raw[feature_names + ["heating_load"]].corr(numeric_only=True)
im = axes[1].imshow(corr, vmin=-1, vmax=1, cmap="coolwarm")
corr_labels = [column_labels[name] for name in corr.columns]
axes[1].set_xticks(range(len(corr)), corr_labels, rotation=90, fontsize=7)
axes[1].set_yticks(range(len(corr)), corr_labels, fontsize=7)
axes[1].set_title("皮尔逊相关矩阵（非因果）")
fig.colorbar(im, ax=axes[1], shrink=.75)
plt.tight_layout()
| 样本数 | 均值 | 标准差 | 最小值 | 下四分位数 | 中位数 | 上四分位数 | 最大值 | |
|---|---|---|---|---|---|---|---|---|
| 相对紧凑度 | 768.0 | 0.764167 | 0.105777 | 0.62 | 0.6825 | 0.75 | 0.8300 | 0.98 |
| 表面积 | 768.0 | 671.708333 | 88.086116 | 514.50 | 606.3750 | 673.75 | 741.1250 | 808.50 |
| 墙体面积 | 768.0 | 318.500000 | 43.626481 | 245.00 | 294.0000 | 318.50 | 343.0000 | 416.50 |
| 屋顶面积 | 768.0 | 176.604167 | 45.165950 | 110.25 | 140.8750 | 183.75 | 220.5000 | 220.50 |
| 建筑高度 | 768.0 | 5.250000 | 1.751140 | 3.50 | 3.5000 | 5.25 | 7.0000 | 7.00 |
| 朝向 | 768.0 | 3.500000 | 1.118763 | 2.00 | 2.7500 | 3.50 | 4.2500 | 5.00 |
| 玻璃面积 | 768.0 | 0.234375 | 0.133221 | 0.00 | 0.1000 | 0.25 | 0.4000 | 0.40 |
| 玻璃分布 | 768.0 | 2.812500 | 1.550960 | 0.00 | 1.7500 | 3.00 | 4.0000 | 5.00 |
| 供暖负荷 | 768.0 | 22.307195 | 10.090204 | 6.01 | 12.9925 | 18.95 | 31.6675 | 43.10 |
第 3 步：先划分，再用训练集统计量标准化#
随机划分必须早于标准化。若先用全体数据计算均值和标准差，测试集的信息会进入训练过程，形成轻微泄漏。这里固定随机种子，并保留 20% 作为最终测试集。
设训练集设计矩阵为 \(\boldsymbol{X}\in\mathbb{R}^{m\times n}\)，各列均值组成的向量为 \(\boldsymbol{\mu}\in\mathbb{R}^{n}\)，各列标准差组成的向量为 \(\boldsymbol{\sigma}\in\mathbb{R}^{n}\)。表达式 \((\boldsymbol{X}-\boldsymbol{\mu})/\boldsymbol{\sigma}\) 利用广播，把两个长度为 \(n\) 的向量逐行用于 \(\boldsymbol{X}\)。观测标签 \(\boldsymbol{y}\) 明确保留为 \((m,1)\) 的二维列向量，避免 \((m,)\) 数组与 \((m,1)\) 数组相减后被广播成 \((m,m)\) 数组。
from sklearn.model_selection import train_test_split
X = raw[feature_names].to_numpy(dtype=np.float64)
y = raw[["heating_load"]].to_numpy(dtype=np.float64) # 保留二维列向量
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.20, random_state=SEED,
)
mu = X_train.mean(axis=0, keepdims=True)
sigma = X_train.std(axis=0, keepdims=True)
sigma[sigma == 0] = 1.0
X_train_z = (X_train - mu) / sigma
X_test_z = (X_test - mu) / sigma
assert X_train_z.shape == (len(X_train), len(feature_names))
assert y_train.shape == (len(X_train), 1)
print("训练特征、训练标签、测试特征、测试标签的维度：",
X_train_z.shape, y_train.shape, X_test_z.shape, y_test.shape)
print("训练集标准化后列均值最大绝对值：", np.abs(X_train_z.mean(0)).max())
训练特征、训练标签、测试特征、测试标签的维度： (614, 8) (614, 1) (154, 8) (154, 1)
训练集标准化后列均值最大绝对值： 1.9217345774782027e-15
第 4 步：广播、向量化与维度核对#
线性模型为 \(\widehat{\boldsymbol{y}}=\boldsymbol{X}\cdot\boldsymbol{w}+b\)，其中 \(\boldsymbol{w}\in\mathbb{R}^{n\times1}\)，\(b\) 是标量。矩阵乘法一次得到 \(m\) 个样本的预测值，\(b\) 通过广播加到每一行。均方误差为
下面先用一个小批量核对循环计算与向量化计算是否一致，再展示一个常见的维度错误。维度断言不是装饰：它可以尽早发现“程序能够运行，但计算含义不正确”的广播。
w_demo = rng.normal(size=(X_train_z.shape[1], 1))
b_demo = 0.3
small_X = X_train_z[:5]
pred_loop = np.array([
sum(row[j] * w_demo[j, 0] for j in range(row.size)) + b_demo
for row in small_X
]).reshape(-1, 1)
pred_vectorized = small_X @ w_demo + b_demo
np.testing.assert_allclose(pred_loop, pred_vectorized)
print("循环与向量化最大差：", np.max(np.abs(pred_loop - pred_vectorized)))
dangerous_y = y_train.ravel()
print("(m,1) - (m,) 的错误广播结果：", (y_train - dangerous_y).shape)
assert y_train.ndim == 2 and y_train.shape[1] == 1
循环与向量化最大差： 4.440892098500626e-16
(m,1) - (m,) 的错误广播结果： (614, 614)
第 5 步：最小二乘解析解#
给 \(\boldsymbol{X}\) 增加一列全 1 后，可以把偏置也写入参数向量。最小二乘估计量满足正规方程 \((\boldsymbol{X}^{\mathsf T}\cdot\boldsymbol{X})\cdot\boldsymbol{\theta}=\boldsymbol{X}^{\mathsf T}\cdot\boldsymbol{y}\)。数值计算时不显式计算矩阵的逆，而使用 np.linalg.lstsq 求解最小二乘问题；当特征之间高度相关或矩阵接近奇异时，这种方法通常更加稳定。该解析解可以作为基准，用于检查后面的梯度下降程序。
Xb_train = np.c_[X_train_z, np.ones((len(X_train_z), 1))]
Xb_test = np.c_[X_test_z, np.ones((len(X_test_z), 1))]
theta_closed, *_ = np.linalg.lstsq(Xb_train, y_train, rcond=None)
pred_closed = Xb_test @ theta_closed
print("参数维度（8 个权重 + 1 个偏置）：", theta_closed.shape)
参数维度（8 个权重 + 1 个偏置）： (9, 1)
第 6 步：推导并检查梯度#
令 \(\boldsymbol{e}=\boldsymbol{X}\cdot\boldsymbol{w}+b-\boldsymbol{y}\)，则
中心差分分别把某个参数增加和减少一个很小的正数 \(\epsilon\)，再用 \([\mathcal{L}(\theta+\epsilon)-\mathcal{L}(\theta-\epsilon)]/(2\epsilon)\) 近似该参数的导数。由于每检查一个参数都要额外计算两次损失，它只适合在小型例子中抽查梯度。解析梯度与数值梯度足够接近后，再继续训练模型。
def mse_and_grad(X_batch, y_batch, w, b):
residual = X_batch @ w + b - y_batch
loss = float(np.mean(residual ** 2))
grad_w = (2.0 / len(X_batch)) * X_batch.T @ residual
grad_b = float(2.0 * residual.mean())
return loss, grad_w, grad_b
w0 = rng.normal(scale=.1, size=(X_train_z.shape[1], 1))
b0 = 0.0
loss0, grad_w0, grad_b0 = mse_and_grad(X_train_z[:32], y_train[:32], w0, b0)
epsilon = 1e-6
j = 2
w_plus, w_minus = w0.copy(), w0.copy()
w_plus[j, 0] += epsilon
w_minus[j, 0] -= epsilon
loss_plus = mse_and_grad(X_train_z[:32], y_train[:32], w_plus, b0)[0]
loss_minus = mse_and_grad(X_train_z[:32], y_train[:32], w_minus, b0)[0]
numeric = (loss_plus - loss_minus) / (2 * epsilon)
relative_error = abs(numeric - grad_w0[j, 0]) / max(1.0, abs(numeric), abs(grad_w0[j, 0]))
print({"解析梯度": grad_w0[j, 0], "数值梯度": numeric, "相对误差": relative_error})
assert relative_error < 1e-6
{'解析梯度': np.float64(-23.41832222459604), '数值梯度': -23.418322257384716, '相对误差': np.float64(1.400129193812618e-09)}
第 7 步：批量梯度下降#
每轮都使用整个训练集，因此曲线应相对平滑。学习率过小会收敛慢，过大则损失振荡甚至溢出。这里保存损失但不偷看测试集；测试集只能在模型确定后使用。
w = np.zeros((X_train_z.shape[1], 1))
b = 0.0
history = []
learning_rate = 0.03
for step in range(3000):
loss, grad_w, grad_b = mse_and_grad(X_train_z, y_train, w, b)
w -= learning_rate * grad_w
b -= learning_rate * grad_b
if step % 25 == 0:
history.append((step, loss))
history = np.asarray(history)
plt.semilogy(history[:, 0], history[:, 1])
plt.xlabel("更新步")
plt.ylabel("训练 MSE（对数轴）")
plt.title("批量梯度下降收敛曲线")
plt.show()
print("末次记录的训练 MSE：", history[-1, 1])
末次记录的训练 MSE： 8.40006291044383
第 8 步：在模型确定后评价测试集#
RMSE 与供暖负荷使用相同的单位，便于理解预测误差的大致大小；\(R^2\) 衡量当前模型相对于“把所有测试样本都预测为训练集观测标签均值”这一简单方法的改进程度。下面还绘制“真实值—预测值”图和“预测值—残差”图。若残差呈现明显弧线或随预测值增大而扩散，线性关系或等方差假设可能不合适。梯度下降与最小二乘解析解的结果接近，也能为程序实现正确提供证据。
from sklearn.metrics import mean_squared_error, r2_score
pred_gd = X_test_z @ w + b
metrics = pd.DataFrame({
"模型": ["最小二乘解析解", "梯度下降"],
"均方根误差（RMSE）": [
mean_squared_error(y_test, pred_closed) ** .5,
mean_squared_error(y_test, pred_gd) ** .5,
],
"决定系数（R²）": [r2_score(y_test, pred_closed), r2_score(y_test, pred_gd)],
})
display(metrics)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].scatter(y_test, pred_gd, s=15, alpha=.7)
low, high = float(y_test.min()), float(y_test.max())
axes[0].plot([low, high], [low, high], "k--")
axes[0].set(xlabel="真实值", ylabel="预测值", title="预测—真实")
residual = y_test - pred_gd
axes[1].scatter(pred_gd, residual, s=15, alpha=.7)
axes[1].axhline(0, color="black", linestyle="--")
axes[1].set(xlabel="预测值", ylabel="残差", title="残差诊断")
plt.tight_layout()
| 模型 | 均方根误差（RMSE） | 决定系数（R²） | |
|---|---|---|---|
| 0 | 最小二乘解析解 | 3.025427 | 0.912185 |
| 1 | 梯度下降 | 3.036574 | 0.911536 |
第 9 步：比较标准化前后的训练结果#
若把同一个学习率直接用于未经标准化的面积等特征，不同参数的更新大小会很不均衡。下面只运行少量更新，并在损失出现无穷大或无效数值时停止。这个比较并不是要说明“所有模型都必须标准化”，而是帮助理解为什么采用梯度下降的线性模型容易受到特征数值尺度的影响；其他类型的模型可能具有不同表现。
w_raw = np.zeros((X_train.shape[1], 1))
b_raw = 0.0
raw_losses = []
for step in range(100):
loss, gw, gb = mse_and_grad(X_train, y_train, w_raw, b_raw)
raw_losses.append(loss)
if not np.isfinite(loss):
break
w_raw -= learning_rate * gw
b_raw -= learning_rate * gb
print("原始量纲前 5 个损失：", raw_losses[:5])
print("是否保持有限：", bool(np.isfinite(raw_losses[-1])))
原始量纲前 5 个损失： [591.7693806188926, 579997563904.2013, 7.405944828480892e+20, 9.456659314631366e+29, 1.2075218960306289e+39]
是否保持有限： False
阅读结果时的检查顺序#
运行程序后，不要只记录最后一行指标。建议依次回答：第一，训练损失是否稳定下降，还是因学习率不合适而明显波动；第二，最小二乘解析解与梯度下降在同一测试集上的差异是否足够小；第三，残差是否围绕 0 随机分布，是否随着预测值增大而扩散；第四，最大残差对应的建筑特征是否超出训练集中的常见范围；第五，改变随机划分后结论是否稳定。
这个顺序分别检查“优化是否成功”“程序是否正确”“模型假设是否合理”和“模型能否用于不同样本”。这些问题不能由一个 \(R^2\) 完全回答。若优化失败，应先检查程序或学习率；若优化正常但残差仍有明显规律，应考虑调整特征或模型；若只有少量边缘设计预测较差，应如实说明模型的适用范围，而不应随意删除这些样本。
可复现记录#
固定随机种子只控制本案例明确设置的数据划分与参数初始化，并不能保证不同计算机上的每一位数都完全相同。正式记录还应包括 Python、NumPy、pandas、scikit-learn 和 openpyxl 的版本，原始 ZIP 文件的 SHA-256 摘要，训练集与测试集的样本编号，快速模式设置以及运行时间。数据文件摘要和软件版本有助于查找重复实验结果不一致的原因。
结论、局限与常见错误#
广播按照数组轴和维度工作，并不会自动理解编程者的意图。预测值与观测标签都应保持为 \((m,1)\) 的二维列向量，并用断言检查维度。
最小二乘解析解和梯度下降求解的是同一个线性回归目标；若二者差异很大，应先检查标准化步骤、梯度中的 \(2/m\) 因子和学习率。
高 R² 不等于因果关系。样本是模拟建筑组合，真实地区、材料价格、气候和施工约束均未进入模型。
随机划分适用于本数据的独立设计组合；带时间或空间结构的数据必须改用相应划分。
Excel 解析依赖 openpyxl。服务器应把数据缓存挂载到源码树之外，不要把解压文件打包进 HTML。
综合练习#
同时预测供暖负荷与制冷负荷，使 \(\boldsymbol{y}\) 变为 \((m,2)\) 的矩阵。写出 \(\boldsymbol{w}\) 与 \(\boldsymbol{b}\) 的新维度，并验证广播仍然正确。
用
timeit比较 10,000 次循环预测与矩阵乘法。重复计时多次，并解释为什么小数组的单次计时结果可能波动。加入二次项或岭回归，仍只用训练集选择正则化强度。比较残差图，而不只比较 R²。
把学习率改成 0.003、0.03、0.3，画在一张图上；解释“收敛快”和“泛化好”为什么不是同一概念。
将本次下载显示的 SHA-256 摘要写入
catalog.json，再改变下载文件中的一个字节，验证下载函数能否识别文件已经发生变化。