前言
具体的代码参见我的Github库。
第一步：数据抽象
用惯PyTorch等库的小伙伴应该清楚，机器学习的基本单位通常是“张量”(tensor)。所谓张量，感性来理解，张量就是一个多维的数据结构。举例来说，下面的几种数据结构都算作张量：
- 标量，比如1、2、0等。
- 矩阵。
- 高维矩阵，2*3*4形状的数据。
那么，我们要如何用基本的数据结构来表示张量呢？有这样几种组织方式。
首先就是参考Python，将张量组织为数字、数组或者数组的数组，但是，这就要求我们要手搓一个类似于Python中List的数据结构，此外，也要对标量、矩阵的情况进行特殊优化。因为高维的张量可以用List的List来表示，而数字和数组用这种方式表达，不够直观，不够方便。
其次就是表示为一个成员为Tensor的数组，Tensor可以是一个数字，也可以是一个数组，也可以是一个数组的数组。这种组织结构没有任何问题，但是，一旦涉及到维度的操作（比如三维张量转为矩阵），就很复杂了。
其实我们还有更简单的方法——用一个数组来表示张量，此外，额外设置几个字段来表示形状，如下所示：
1
2
3
4
5
| typedef struct Tensor {
double* data;
int* shape;
int n_dim;
} Tensor;
|
三个字段分别表示存储的具体数据、每个维度的形状以及维度的个数。张量的首要目的是存储数据，其次是存储维度信息，我们完全没有必要按照Python的逻辑来设置我们的机器学习库。
依靠这种表示方式，不仅免去了特殊情况的考量，不论是标量、矩阵还是高维张量，不同之处只在于shape和n_dim两个字段而已，这也方便了我们做涉及维度的操作。当然，这种方法也有一些缺点，拿矩阵来举例子，获取某行某列的元素时，如果依靠数组的数组表示法，编码很简单，而英语这种方法表示，则需要编写额外的运算：
1
2
| double s = Tensor.data[1][2];
double s = Tensor.data[shape[1] + 2];
|
当然，由于张量的维度通常都小于四维，这点开销并不算什么问题。
第二步：反向传播
一个神经网络的优化，依靠于反向传播，首先计算损失函数相对于最终节点的梯度，然后依靠链式法则逐层逐节点计算梯度，最后依靠梯度下降法优化模型的参数。如果您不是很理解，请参考相应的文章。
反向传播有两个要点：
- 要记录该节点的前序节点，便于反方向的运算。
- 要记录该节点的运算方法，即该节点在前序节点的基础上，通过什么运算得到的，从而选择合适的求导公式。
那么，我们的张量表示需要修改修改：
1
2
3
4
5
6
7
8
9
10
11
12
| typedef struct Tensor {
double* data;
int* shape;
int n_dim;
double* grad;
int n_prev;
struct Tensor** prev;
int fn_id;
} Tensor;
|
最后四个字段分别表示：计算出来的梯度、前序节点个数、前序节点表、前序节点的运算ID（即当前节点是通过前序节点的哪种运算得到的）。
对于每一个运算ID，都要编写两个函数：前向传播函数以及反向传播的求导函数，对于张量加法运算而言：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
|
Tensor* AddTensor(Tensor* tensor1, Tensor* tensor2) {
if(!IsSameShape(tensor1, tensor2)) {
perror("Tensor shape not same");
exit(EXIT_OP_SHAPE_FAILURE);
}
Tensor* result = CreateZeroTensor(tensor1->n_dim, tensor1->shape);
int data_num = GetDataNum(result);
for(int i = 0; i < data_num; i++) {
result->data[i] = tensor1->data[i] + tensor2->data[i];
}
result->n_prev = 2;
result->prev = (Tensor**)malloc(sizeof(Tensor*) * 2);
if(result->prev == NULL) {
perror("Prev malloc");
exit(EXIT_MALLOC_FAILURE);
}
result->prev[0] = tensor1;
result->prev[1] = tensor2;
result->fn_id = OP_ADD;
}
void AddBackprop(Tensor* tensor) {
for(int i = 0; i < GetDataNum(tensor); i++) {
tensor->prev[0]->grad[i] = tensor->grad[i];
tensor->prev[1]->grad[i] = tensor->grad[i];
}
}
|
反向传播函数中，我们更新的是前序节点，因为根据链式法则有下面的公式：
$$
\frac{\partial{L}}{\partial{n_{-1}}}=\frac{\partial{L}}{\partial{n}}\frac{\partial{n}}{\partial{n_{-1}}}
$$
其中，$L$代表损失函数，$n_{-1}$代表当前节点的前序节点，$n$代表当前节点。上式的含义是：损失函数相对于当前节点的前序节点的偏导，等于损失函数相对于当前节点的偏导(存储在tensor->grad中)与当前节点数据对前序节点的偏导(运算带来的偏导，加法的偏导为1，故忽略了)。
第三步：梯度下降
通过反向传播，我们计算了每一个节点的梯度，在这之后，我们需要根据梯度优化每个节点的参数，这个过程被称作梯度下降，所依靠的公式为：
$$
P_{+1}=P-\eta G
$$
其中，$P$是未优化的参数，$P_{+1}$是优化后的参数，$\eta$是学习率，代表向优化方向移动的步长，$G$是梯度。这个公式很好理解，计算梯度之后，就确定了函数的参数需要向哪边移动，移动的距离用另外的参数设定。如果不清楚，请参阅相关资料。
实现出来就是：
1
2
3
4
5
6
7
8
9
10
| bool SGDOptimize(Tensor* tensor, double learning_rate) {
int data_num = GetDataNum(tensor);
for(int i = 0; i < data_num; i++) {
tensor->data[i] -= learning_rate * tensor->grad[i];
}
return CheckNAN(tensor);
}
|
第四步：简单的迭代
一个简单的例子：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
| int main(int argc, char* argv[]) {
Tensor* in_tensor = CreateBaseOnArray(2, (int[]){6, 1}, (double[]){1, 2, 3, 4, 5, 6});
Tensor* layer = CreateBaseOnArray(2, (int[]){1, 1}, (double[]){1});
Tensor* gt_tensor = CreateBaseOnArray(2, (int[]){6, 1}, (double[]){6, 12, 18, 24, 30, 36});
double lr = 0.001;
for(int iter = 0; iter < 200; ++iter) {
Tensor* out_tensor = MatMul(in_tensor, layer);
Tensor* loss = MSE(gt_tensor, out_tensor);
Backprop(loss);
SGDOptimize(layer, lr);
printf("%lf\n", layer->grad[0]);
ZeroGradTensor(layer);
free(loss);
free(out_tensor);
}
PrintTensor(layer);
FreeTensor(in_tensor);
FreeTensor(layer);
FreeTensor(gt_tensor);
return 0;
}
|
