关于为什么要写这个课程 lab
其实暑假的时候一直在推的时 Prof.Li Mu的 d2l ，推着推着发现了一些问题。d2l本身的讲解其实还算细致，但实际上更偏向于从宏观角度进行把握。也就是说，跟着 d2l 能够搭建起对基本模型和概念的认知，但由于 pytorch(用的是这个) 库的“过于”成熟，搭建起的模型在底层计算上常常是通过 torch 暴露的接口进行直接调用。这就导致不读源码就会对底层实现是一知半解甚至浑然不知，进而在尝试进行 debug或者优化时无从下手。然鹅经过多年发展，torch 库的架构现在变得及其庞大，即使有 Document 辅助查询也耗时繁多(更不用提 torch 文档的加载问题)，所以考虑跟进 CMU 的这门，尝试从0搭建出类 numpy 库和 torch 的底层架构
ok 那么进入相关项目，个人的项目仓库放在 这了 ，我争取对Jupter Notebook 完成全中翻(当然是机翻(❍ᴥ❍ʋ))，个人配置是跑在 x86 下的，所以有些实现会出现不同。我尽量在仓库的 .ipynb 和当前的笔记中追加说明。
HW0
hw0 的目的只是确认当前你的知识储备是否足够应对这个项目，从 hw0 整体来看，它至少要求你要掌握：
- 基础的 python 语法， numpy 库的使用
- 基本的 C++ 语法
- Bash 的基本使用
- 由于参与课程的人员提交的项目可以在 Colab 环境下，个人完成项目时建议设备有基础的深度学习配置并且建议有 CUDA 支持
问题一：基本的 add 函数和测试/自动提交
这一部分没啥好说的，simple_ml.add()就是一个很基本的加法函数， return a + b 即可
注意一下 x86 架构下调用 python 用的是 python 而不是 python3
自动提交略(后续皆同)，因为用不到
问题二：加载 MNIST 数据
这一部分要求我们完成 src/simple_ml.py 下的 parse_mnist_data() 函数，在 pytorch 中相关工作是由 Dataloader 完成的，我们只需要实例化一个 Dataloader 并指定它包含的对应数据集即可，但此处需要我们自己读取数据。
由 YANN 关于MNIST 格式的介绍可知，对于 MNIST 数据集，其通常为四个特殊格式化的二进制文件，具体来说
train-images-idx3-ubyte为训练集图像train-label-idx1-ubyte为训练集标签t10k-images-idx3-utype为测试集图像t10k-labels-idx1-ubyte测试集标签
所有文件都为一种简单的 IDX 文件格式，用于储存向量和多维矩阵
对于不同文件，我们通过前 8 个字节的值来进行类型分辨
比如对 train-images-idx3-ubyte :
| 偏移量 | 数据类型/值 | 描述 |
| 0 | 32-Integer/2051 | Magic Number |
| 4 | 32-Integer/* | 图像数量 |
以及 train-labels-idx1-ubyte:
| 偏移量 | 数据类型/值 | 描述 |
| 0 | 32-Integer/2049 | Magic Number |
| 4 | 32-Integer/* | 标签数量 |
通过 Magic Number 的校正，我们能够辨别当前读入文件的类型并进行相应操作。
所以思路很明确了，对当前文件名进行二进制读入，并对前 8 个字节先后进行魔数校验和总数提取
其中二进制读入我们依赖 gzip 包，字节校验依照 struct 包，处理结果返回 numpy array. 给出 pass 的参考代码如下：
import struct
import numpy as np
import gzip
# 加载图像和标签
with gzip.open(label_filename, 'rb') as f:
magic, num_images = struct.unpack('>II', f.read(8)) # 当前架构为 x86 64 大端序 且前8个字节为魔数和图像数量
if magic != 2049:
raise ValueError('非法的标签魔数: {}'.format(magic))
labels = np.frombuffer(f.read(), dtype=np.uint8)
with gzip.open(image_filename, 'rb') as f:
magic, num_images, num_rows, num_cols = struct.unpack('>IIII', f.read(16))
if magic != 2051:
raise ValueError('非法的图像魔数: {}'.format(magic))
images = np.frombuffer(f.read(), dtype=np.uint8).reshape(num_images, num_rows*num_cols)
# 归一化图像至 [0,1]
images = images.astype(np.float32) / 255.0
return images, labels
注意我们当前的架构为 x86，读入顺序为大端序， Linux/MacOS 下为小端序读入， struct 的 unpack 需要注意参数。
问题三：Softmax loss
这一部分要求实现实现 src/simple_ml.py 中的 softmax_loss() 函数中定义的 softmax(又称交叉熵)损失,，交叉熵损失定义为：
对于可以取值 $y\in\{1,\dots, k\}$ 的多类输出，$\text{Softmax loss}$ 将 $\text{logits}$ 向量 $z\in\mathbb{R}^k$ 作为输入，真类 $y\in\{1,\dots, k\}$ 的损失表述为：
$$\ell _{Softmax}(z,y)=\log\sum\limits_{i=1}^k\exp(z_i-z_y)$$
所以实现过程相当简单，照着公式拼就行：
def softmax_loss(Z, y):
# 数值稳定性：减去每行的最大值
max_Z = np.max(Z, axis=1, keepdims=True)
exp_Z = np.exp(Z - max_Z) # 形状: (batch_size, num_classes)
# 计算softmax分母
sum_exp_Z = np.sum(exp_Z, axis=1, keepdims=True) # 形状: (batch_size, 1)
# 获取每个样本真实类别的概率
n = Z.shape[0]
true_class_probs = exp_Z[np.arange(n), y] # 形状: (batch_size,)
# 计算log softmax概率 for true class
log_probs = np.log(true_class_probs / sum_exp_Z.flatten()) # 除法后形状匹配
# 返回平均损失
return -np.mean(log_probs)
注意一下为了防止上溢，代码中我们对输入进行了最大值归一化，但没有对分母进行 $\epsilon$ 加和以避免下溢。
讲义中还提示一行代码，在此一并给出：
return -np.mean( # 平均值
np.log( # 对数
np.exp( # 指数
(Z-np.max(Z,axis=1,keepdims=True))[np.arange(Z.shape[0]),y] # 取最大值后再取对应标签的logit
)/np.sum( # 求和
np.exp(Z-np.max(Z,axis=1,keepdims=True)),axis=1,keepdims=True
) # 求softmax
)
)
问题四：Softmax 回归的随机梯度下降法
这一问要求我们实现 Softmax 的单步梯度下降(Stochastic Gradient Descent, SGD)，严格来说，应该是基于Softmax loss 的单步梯度下降，所以流程其实相对还是比较清楚的：
$$\text{Forward}\rightarrow\text{Calculate Loss}\rightarrow \text{Calculate Gradient}\rightarrow\text{Update Parameters} $$
所以我们可以写出相应的伪代码：
def sgd():
probs = softmax(logits)
loss = cross_entropy(logits, y_one_hot)
gradient = compute_gradient(loss)
但讲义同时提示了，Softmax 的梯度函数实际上我们是能精确求得的：
$$L = -\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{C}y_{ij}\log(p_{ij})$$
其中 $p_{ij} = \frac{e^{z_{ij}}}{\sum\limits_{k = 1}^{C}}e^{z_{ik}$，为 Softmax 概率
那么由链式法则，我们可以得到：
$$\frac{\partial L}{\partial z_{ij}} = p_{ij} – y_{ij}$$
所以我们不需要显式计算 loss，直接计算偏导即可
def softmax_regression_epoch(X, y, theta, lr = 0.1, batch=100):
"""
对数据执行单轮SGD训练实现softmax回归，使用步长lr和指定的批处理大小。该函数应原地修改theta矩阵，并按顺序遍历X中的批次（不进行随机打乱）。
参数说明：
X (np.ndarray[np.float32]): 二维输入数组，尺寸为
(样本数 x 输入维度)
y (np.ndarray[np.uint8]): 一维类别标签数组，尺寸为 (样本数,)
theta (np.ndarray[np.float32]): softmax回归参数的二维数组，形状为
(输入维度, 类别数)
lr (float): SGD的步长（学习率）
batch (int): SGD小批量样本数量
返回值：
无
"""
n_samples = X.shape[0]
n_classes = theta.shape[1]
for i in range(0, n_samples, batch):
X_batch= X[i:i+batch]
y_batch= y[i:i+batch]
batch_size = X_batch.shape[0]
# 计算softmax概率
Z = X_batch @ theta
max_Z = np.max(Z, axis=1, keepdims=True)
exp_Z = np.exp(Z - max_Z) # 形状: (batch_size, num_classes)
sum_exp_Z = np.sum(exp_Z, axis=1, keepdims=True) # 形状: (batch_size, 1)
probs = exp_Z / sum_exp_Z # 形状: (batch_size, num_classes)
# 热编码 y
y_one_hot = np.eye(n_classes)[y_batch] # 形状: (batch_size, num_classes)
# 计算梯度
dZ = (probs - y_one_hot) / batch_size # 形状: (batch_size, num_classes)
dW = X_batch.T @ dZ # 形状: (input_dim, num_classes)
# 更新参数
theta -= lr * dW
再次强调此处使用的是解析梯度，使用数值梯度的版本同步给出，但此时的时间复杂度高达 $O(n^3)$ ，非必要不建议使用
def softmax_regression_epoch_numerical(X, y, theta, lr=0.1, batch=100, epsilon=1e-5):
"""
使用数值梯度方法实现的softmax回归epoch（仅用于教学和理解）
"""
n_samples = X.shape[0]
input_dim, n_classes = theta.shape
for i in range(0, n_samples, batch):
X_batch = X[i:i+batch]
y_batch = y[i:i+batch]
batch_size = X_batch.shape[0]
# 初始化数值梯度矩阵
numerical_grad = np.zeros_like(theta)
# 对每个参数进行扰动，计算数值梯度
for i in range(input_dim):
for j in range(n_classes):
# 保存原始参数值
original_value = theta[i, j]
# 计算 f(θ + ε)
theta[i, j] = original_value + epsilon
Z_plus = X_batch @ theta
loss_plus = softmax_loss(Z_plus, y_batch)
# 计算 f(θ - ε)
theta[i, j] = original_value - epsilon
Z_minus = X_batch @ theta
loss_minus = softmax_loss(Z_minus, y_batch)
# 恢复原始参数值
theta[i, j] = original_value
# 计算中心差分近似梯度
numerical_grad[i, j] = (loss_plus - loss_minus) / (2 * epsilon)
# 使用数值梯度更新参数
theta -= lr * numerical_grad
问题五：双层神经网络随机梯度下降
在完成问题四后，我们得到了一个 SGD， 现在我们考虑一个双层神经网络：
$$Z = W_2^T\cdot\text{ReLU}(W_1^T\cdot X)$$
其中
- $X$: 输入矩阵，形状为 $(N, D)$，其中 $N$ 是批处理大小（batch size），$D$ 是输入维度。
- $W1$: 第一层权重矩阵，形状为 $(D, H)$，其中 $H$ 是隐藏层维度。
- $W2$: 第二层权重矩阵，形状为 $(H, C)$，其中 $C$ 是类别数（输出维度）。
- $y$: 真实标签向量，形状为 $(N,)$。
- $\text{ReLU}(x) = \max(0, x)$: ReLU激活函数。
- 损失函数：交叉熵损失 $L = -\frac{1}{N} \sum_{i=1}^{N} \sum_{j=1}^{C} y_{ij} \log(p_{ij})$，其中 $p_{ij}$ 是softmax概率，$y_{ij}$ 是one-hot编码的真实标签。
考虑继续应用 Softmax loss ，则此时目标为：
$$\text{minimize}_{W_1, W_2} \;\; \frac{1}{m} \sum_{i=1}^m \ell_{\mathrm{ Softmax }}(W_2^T \mathrm{ReLU}(W_1^T x^{(i)}), y^{(i)})$$
也就意味着我们需要每次对 $W_1,W_2$ 进行更新，具体过程为：
- 隐藏层输入计算:
$$H_1 = X \cdot W_1$$
- 输出层(logits)计算:
$$Z = A_1 \cdot W_2$$
反向传播:
- 计算logits的梯度(dZ)
$$\frac{\partial L}{\partial Z} = \frac{1}{N} (P – Y)$$
- 计算W2的梯度(dW2)
$$\frac{\partial L}{\partial W2} = A_1^T \cdot dZ$$
- 计算A1的梯度(dA1)
$$\frac{\partial L}{\partial A_1} = \frac{\partial L}{\partial Z} \cdot W_2^T$$
- 计算W1的梯度(dH1)
$$\frac{\partial L}{\partial H_1} = \frac{\partial L}{\partial A_1} \odot \text{ReLU}'(H_1)$$
- 计算W1的梯度(dW1)
$$\frac{\partial L}{\partial W_1} = X^T \cdot \frac{\partial L}{\partial H_1}$$
更新参数
$$W_1 = W_1\ -\ \text{lr} \cdot dW_1 \\ W_2 = W_2\ -\ \text{lr} \cdot dW_2$$
其中 $\text{lr}$ 是学习率。
那么基于上述推导，我们给出 pass 代码如下：
def nn_epoch(X, y, W1, W2, lr = 0.1, batch=100):
"""
对由权重矩阵W1和W2（无偏置项）定义的双层神经网络执行单轮SGD训练：
logits = ReLU(X * W1) * W2
该函数应使用指定的学习率lr和批处理大小batch（且不随机打乱X的顺序），并原地修改W1和W2矩阵。
参数说明：
X (np.ndarray[np.float32]): 二维输入数组，尺寸为
(样本数 x 输入维度)
y (np.ndarray[np.uint8]): 一维类别标签数组，尺寸为 (样本数,)
W1 (np.ndarray[np.float32]): 第一层权重二维数组，形状为
(输入维度, 隐藏层维度)
W2 (np.ndarray[np.float32]): 第二层权重二维数组，形状为
(隐藏层维度, 类别数)
lr (float): SGD的步长（学习率）
batch (int): SGD小批量样本数量
返回值：
无
"""
def ReLU(x):
return np.maximum(x, 0)
n_samples = X.shape[0]
n_classes = W2.shape[1]
# 批次处理
for i in range(0, n_samples, batch):
X_batch= X[i:i+batch]
y_batch= y[i:i+batch]
# 前向传播
H1 = X_batch @ W1 # 隐藏层输入
A1 = ReLU(H1) # ReLU激活
logits = A1 @ W2 # 输出层
# 计算softmax概率
max_logits = np.max(logits, axis=1, keepdims=True) # 取最大
exp_logits = np.exp(logits - max_logits) # 取指
sum_exp_logits = np.sum(exp_logits, axis=1, keepdims=True) # 求和
probs = exp_logits / sum_exp_logits # 概率
# y 转换为 one-hot 编码
y_onehot = np.eye(n_classes)[y_batch]
# 反向传播
dZ = (probs - y_onehot) / batch # 计算 logits 梯度
dW2 = A1.T @ dZ # 计算 W2 梯度
dA1 = dZ @ W2.T # 计算 A1 梯度
dH1 = dA1 * (H1 > 0) # 计算 H1 梯度
dW1 = X_batch.T @ dH1 # 计算 W1 梯度
# 更新参数
W1 -= lr * dW1
W2 -= lr * dW2
问题6：C++实现 Softmax 回归
这个问题的思路和上一问其实是一样的，只是换成了 C++ 进行实现
void softmax_regression_epoch_cpp(const float *X, const unsigned char *y,
float *theta, size_t m, size_t n, size_t k,
float lr, size_t batch) {
/**
* 对数据执行单轮SGD训练实现softmax回归（C++版本）。
* 按批次顺序处理数据，并原地更新theta矩阵。
*
* 参数：
* X: 输入数据指针，形状为(m, n)
* y: 标签指针，形状为(m,)
* theta: 参数矩阵指针，形状为(n, k)，会被原地修改
* m: 样本数量
* n: 特征维度
* k: 类别数量
* lr: 学习率
* batch: 批处理大小
*/
// 循环处理每个批次
for (size_t start = 0; start < m; start += batch) {
size_t end = std::min(start + batch, m);
size_t current_batch_size = end - start;
// 临时存储当前批次的梯度，初始化为0
std::vector<float> grad(n * k, 0.0f);
// 处理当前批次中的每个样本
for (size_t i = start; i < end; ++i) {
// 计算当前样本的logits Z: 1xk向量
std::vector<float> Z_row(k, 0.0f);
float max_Z = -std::numeric_limits<float>::infinity();
for (size_t j = 0; j < k; ++j) {
for (size_t l = 0; l < n; ++l) {
Z_row[j] += X[i * n + l] * theta[l * k + j];
}
if (Z_row[j] > max_Z) {
max_Z = Z_row[j];
}
}
// 计算exp(Z - max_Z)用于数值稳定性
std::vector<float> exp_Z_row(k);
float sum_exp = 0.0f;
for (size_t j = 0; j < k; ++j) {
exp_Z_row[j] = std::exp(Z_row[j] - max_Z);
sum_exp += exp_Z_row[j];
}
// 计算softmax概率
std::vector<float> probs(k);
for (size_t j = 0; j < k; ++j) {
probs[j] = exp_Z_row[j] / sum_exp;
}
// 构建当前样本的one-hot编码标签
std::vector<unsigned char> one_hot(k, 0);
one_hot[y[i]] = 1;
// 计算当前样本的梯度贡献并累加到grad中
for (size_t j = 0; j < k; ++j) {
float diff = probs[j] - one_hot[j];
for (size_t l = 0; l < n; ++l) {
grad[l * k + j] += diff * X[i * n + l];
}
}
}
// 更新theta：使用平均梯度（除以当前批次大小）和学习率
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < k; ++j) {
theta[i * k + j] -= lr * grad[i * k + j] / current_batch_size;
}
}
}
}
可以看到当前的实现是纯暴力循环，因此其时间复杂度来到了抽象的 $O(n^4)$ 数量级，当然考虑到作为初实现也没有什么问题，但我们也可以对其做出一些优化：
#include <vector>
#include <algorithm>
#include <cmath>
#include <limits>
void softmax_regression_epoch_cpp_optimized(const float *X, const unsigned char *y,
float *theta, size_t m, size_t n, size_t k,
float lr, size_t batch) {
// 为整个批次分配一次性的临时存储，避免在样本循环中反复分配
std::vector<float> Z_batch(batch * k); // 批次的logits
std::vector<float> grad(n * k, 0.0f); // 梯度累加器
std::vector<float> exp_Z(batch * k); // exp(Z)存储
std::vector<float> probs(batch * k); // 概率存储
for (size_t start = 0; start < m; start += batch) {
size_t end = std::min(start + batch, m);
size_t current_batch_size = end - start;
// 重置梯度
std::fill(grad.begin(), grad.end(), 0.0f);
// 批量计算所有样本的logits (矩阵乘法: X_batch @ theta)
for (size_t i = 0; i < current_batch_size; ++i) {
size_t sample_idx = start + i;
for (size_t j = 0; j < k; ++j) {
Z_batch[i * k + j] = 0.0f;
for (size_t l = 0; l < n; ++l) {
// 优化内存访问：连续访问X和theta
Z_batch[i * k + j] += X[sample_idx * n + l] * theta[l * k + j];
}
}
}
// 批量计算softmax
for (size_t i = 0; i < current_batch_size; ++i) {
// 找到当前样本的最大logit值
float max_Z = -std::numeric_limits<float>::infinity();
for (size_t j = 0; j < k; ++j) {
if (Z_batch[i * k + j] > max_Z) {
max_Z = Z_batch[i * k + j];
}
}
// 计算exp(Z - max_Z)和求和
float sum_exp = 0.0f;
for (size_t j = 0; j < k; ++j) {
exp_Z[i * k + j] = std::exp(Z_batch[i * k + j] - max_Z);
sum_exp += exp_Z[i * k + j];
}
// 计算概率
for (size_t j = 0; j < k; ++j) {
probs[i * k + j] = exp_Z[i * k + j] / sum_exp;
}
}
// 批量计算梯度贡献
for (size_t i = 0; i < current_batch_size; ++i) {
size_t sample_idx = start + i;
unsigned char true_class = y[sample_idx];
for (size_t j = 0; j < k; ++j) {
float diff = probs[i * k + j] - (j == true_class ? 1.0f : 0.0f);
for (size_t l = 0; l < n; ++l) {
grad[l * k + j] += diff * X[sample_idx * n + l];
}
}
}
// 更新参数
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < k; ++j) {
theta[i * k + j] -= lr * grad[i * k + j] / current_batch_size;
}
}
}
}
或者更彻底的，使用 Eigen 库：
#include <Eigen/Dense>
#include <vector>
void softmax_regression_epoch_eigen(const float *X, const unsigned char *y,
float *theta, size_t m, size_t n, size_t k,
float lr, size_t batch) {
// 将原始指针转换为Eigen矩阵视图
Eigen::Map<const Eigen::MatrixXf> X_mat(X, m, n);
Eigen::Map<Eigen::MatrixXf> theta_mat(theta, n, k);
for (size_t start = 0; start < m; start += batch) {
size_t end = std::min(start + batch, m);
size_t current_batch_size = end - start;
// 获取当前批次
auto X_batch = X_mat.middleRows(start, current_batch_size);
// 前向传播：计算logits
Eigen::MatrixXf Z = X_batch * theta_mat;
// 计算softmax概率
Eigen::MatrixXf exp_Z = (Z.rowwise() - Z.colwise().maxCoeff()).array().exp();
Eigen::MatrixXf probs = exp_Z.array().colwise() / exp_Z.rowwise().sum().array();
// 创建one-hot标签矩阵
Eigen::MatrixXf y_onehot = Eigen::MatrixXf::Zero(current_batch_size, k);
for (size_t i = 0; i < current_batch_size; ++i) {
y_onehot(i, y[start + i]) = 1.0f;
}
// 计算梯度
Eigen::MatrixXf dZ = (probs - y_onehot) / current_batch_size;
Eigen::MatrixXf grad = X_batch.transpose() * dZ;
// 更新参数
theta_mat -= lr * grad;
}
}
最后，在编译时，针对Windows架构，我们需要在Makefile中添加以下内容：
PYTHON_INCLUDE := $(shell python -m pybind11 --includes)
PYTHON_LIB_DIR := E:/Miniconda3/envs/CMU-10-714/libs # 这里需要修改为自己的Python安装路径
PYTHON_LIB := -lpython39 # 这里需要修改为自己的Python版本
default:
c++ -O3 -Wall -shared -std=c++11 -fPIC \
$(PYTHON_INCLUDE) \
src/simple_ml_ext.cpp \
-o src/simple_ml_ext.pyd \ # 注意后缀名(Copy的时候删掉本条注释)
-L$(PYTHON_LIB_DIR) $(PYTHON_LIB) \
-static-libgcc -static-libstdc++
这是因为：
- Windows的 `pybind11` 编译出的模块文件扩展名是 `.pyd`，而不是 `.so`
- .pyd 文件本质上是一个 DLL，它不仅需要 Python 本身，还可能需要编译器运行时库（MinGW 的 libstdc++-6.dll、libgcc_s_seh-1.dll、libwinpthread-1.dll 等）。
- 如果用 MinGW 编译，运行时库默认不会静态链接，而是需要在运行时找到 .dll 文件。
- Python 导入 .pyd 时，会让 Windows 动态加载器去找依赖 DLL：
- 先在 .pyd 所在目录搜索
- 再在系统路径 (PATH) 搜索
- 如果找不到就报 “找不到指定的模块”。
- 为了避免这些麻烦，我们直接选择静态链接，这样 libgcc 和 libstdc++ 会直接被打进 .pyd，运行时只需要找到 python39.dll就行
至此我们完成了 hw0 的所有内容，这一部分的内容还是很基础的，但写起来真的能把人写累死ㅍ_ㅍ
发表回复