***
代码部分可以直接通过Jupyter Notebook来查看

这几天在上Andrew Ng教授开的Coursera系列课程Deep Learning，总觉得光是看视频和做作业还不够，还是得自己动手写写代码，亲自实现课程里提到的算法内容，于是便有了这篇博客，作为自己入门深度学习的里程碑吧。

前馈神经网络

机器学习有两个基本问题，一是回归，二是分类，神经网络大多用于解决分类问题，前馈神经网络（feedforward neural network）是整个神经网络家族中较为常见和较为基础的一种，如下图右上角的DFF所示。图片来源是Cheat Sheets for AI, Neural Networks, Machine Learning, Deep Learning & Big Data。

神经网络中的基本元素是神经元，每层都有一定数量的神经元，神经元组合的多样性决定了神经网络的丰富性。下面是一个简单的前馈神经网络，总共有三层，从左到右分别是输入层、隐层和输出层，输入层的x1和x2表示这个样本只有两个特征（自变量），因为输入层通常不计入内，所以这是一个两层的神经网络，第一层有4个神经元，第二层只有1个。注意，隐层可以不止一层，隐层设置得越多，整个神经网络越庞大。

这个神经网络的工作原理是，给定一个样本的数据，将数据传输到第一层，进行线性变换和激活变换，得到加工过后的数据，这份新数据传到第二层，作为第二层的输入，接着进行线性变换和激活变换，又得到一份新的数据，因为第二层是最后一层了（如果不止两层就一直进行这样的操作直到抵达最后一层为止），所以最终的输出作为我们对该样本的预测值y_hat。

每一个神经元如同工厂的流水车间的机器，它重复做着四件事情：【接受上一层数据作为输入>>线性变换>>激活变换>>输出数据到下一层】，每个神经元中有三个组成部分：权重(weight)矩阵W，偏置(bias)向量b，以及激活函数(activation function) g，用公式表达为下图，其中上标(i)表示这是第i个样本数据，上标[1]和[2]分别表示这是神经网络的第一层与第二层：

$z^{[1](i)} = W^{[1]} x^{(i)} + b^{[1]} $
$a^{[1](i)} = \tanh(z^{[1](i)})$
$z^{[2](i)} = W^{[2]} a^{[1](i)} + b^{[2]}$
$\hat{y}^{(i)} = a^{[2](i)} = \sigma(z^{ [2](i)})$
$y^{(i)}_{prediction} = \begin{cases} 1 & \mbox{if } a^{[2](i)} > 0.5 \\ 0 & \mbox{otherwise } \end{cases}$

公式(1)中，通过简单的线性变换得到了z[1]，z称为prev_activation，接着通过激活变换（这里用的激活函数是tanh函数，下面会讲到），得到了a[1]，称为activation，公式(3)(4)表达的是第二层的线性变换和激活变换，和第一层大同小异，只不过第二层的激活函数用的不是tanh函数而是sigmoid函数。

激活函数有下面三种，都是执行了非线性变换，实现的效果都是将prev_activation转化为activation。

$ sigmoid(z) = \frac{1}{1+e^{-z}} $
$ tanh(z) = \frac{e^{z} - e^{-z}}{e^{z} + e^{-z}} = 2 sigmoid(2z) - 1 $
$
relu(z) =
max(0, z) =
begin{cases}

z, & z>0 \\

0, & z leq 0
end{cases}
$

每一个神经元都可以从上面三种激活函数中选取一个作为自己的激活函数g，经验表明，使用tanh函数的效果总是碾压使用sigmoid函数，所以人们大多使用tanh作为激活函数，近年来人们发现了relu函数，发现它的性能比tanh更好，relu成为了广受欢迎的激活函数。既然sigmoid性能最差，为什么还要介绍它？在一开始的时候说到，神经网络通常用于分类，比方说，给定一张图片，去识别预测它是不是一只猫：

我们的返回值应该是范围在0~1之间的概率值，sigmoid的函数范围是(0, 1), tanh的范围是(-1, 1), relu的是[0, +∞)，使用sigmoid显然更合适些。所以通常一个神经网络的配置是，中间的隐层的所有神经元使用tanh或者relu作为激活函数，输出层的神经元使用sigmoid。

之所以要在线性变换之后进行非线性变化，是因为，如果没有非线性变换，纯粹使用线性变化的话，不管使用了多少层的线性变换，最终的结果通过合并同类项之后仍然是线性变换，100层的神经网络和1层的没有任何差别。神经网络从本质上来说就是一系列的非线性变换。

讲完了基本概念，我们现在的问题是，如何训练这个神经网络。

和其他的监督学习方法一样，神经网络同样是需要定义一个损失函数(cost function)，通过对它进行最小化，得到最优的参数W和b，即每个神经元的权重矩阵和偏置向量。神经网络用的损失函数是交叉熵：

$ J = -\frac{1}{m} \sum_{i = 1}^{m}{ [ y^{(i)}log(\hat{y}^{(i)}) + (1-y^{(i)}) log(1-\hat{y}^{(i)}) ] } $

其中m为样本总数，y_hat为最终预测值，是通过每一层的神经元进行层层加工处理得到的，L表示整个神经网络的层数，即[L]表示最后一层：

$ Z^{[1]} = W^{[1]}X + b^{[1]} $
$ A^{[1]} = g^{[1]}(Z^{[1]}) $
$ Z^{[2]} = W^{[2]}X + b^{[2]} $
$ A^{[2]} = g^{[2]}(Z^{[2]}) $
$ ... $
$ A^{[L]} = g^{[L]}(Z^{[L]}) = \hat{Y}$

优化损失函数的方法是梯度下降法，梯度下降的原理这里就不说了，就贴一下参数的更新公式：

$ W^{[l]} = W^{[l]} - \alpha \text{ } dW^{[l]}$
$ b^{[l]} = b^{[l]} - \alpha \text{ } db^{[l]}$
$ \alpha:learning \text{ }rate $
$ dW \text{ | } db: gradients$

以前在看书的时候，看到神经网络部分，总是在说“神经网络训练使用的是反向传播(back propagation)算法”，一直搞不明白这四个字究竟是什么意思，现在总算是有些理解了。

既然有反向传播，就有正向传播，所谓正向，就是沿着神经网络从输入层到输出层，就是沿着每个神经元执行【接受上一层数据作为输入>>线性变换>>激活变换>>输出数据到下一层】的方向，就是下图中从左到右的方向，而反向传播，指的是从输出层到输入层，即从右到左。

所谓的传播，其实就是把整个神经网络运算的过程抽象成了流程图的形式。正向传播时，数据通过第一层神经元，被处理成了新的数据，传输到下一层，再加工成新数据，再传到下一层，这便是数据的传播，而在反向传播时，传播的是梯度，从最后一层（输出层）出发，首先根据正向传播得到的预测值，计算得出最后一个神经元处的A(activation)的梯度，进而得到Z(prev_activation)的梯度，根据公式再算出W, b以及倒数第二层的A的梯度，就像多米诺骨牌，不断进行下去，直到将梯度传播到第一层的神经网络处。

在正向传播时，每个神经元执行了两步操作，线性变化和激活变换，线性变换是由W[l], b[l], A[l-1]得到了z[l]，激活变化是由z[l]得到了A[l]；而在反向传播时，在每个神经元中，根据相应的激活函数由dA[l]得到dz[l]，再由dz[l]得到dW[l], db[l], dA[l-1]，这可以看做是线性变化和激活变换的反向操作。

由此得到了每一层的神经元的W和b的梯度，参数W和b也由此得以更新，反向传播时传播的是梯度，也可以理解为传播的是误差，是当前模型的不足，参数W和b借助这个信息来进行修正，从而优化模型。于是整个神经网络就按照下图的流程进行迭代，直到达到理想的模型效果：

所以整个神经网络的操作流程是这样的：

• 确定神经网络的层数和每层的神经元数
• 初始化参数W和b

• 进行循环迭代

  • 前向传播，得到预测值
  • 计算损失函数
  • 反向传播，得到参数W和b的梯度
  • 对参数W和b进行更新
• 对新数据做预测

代码实现

课程里有给出前馈神经网络的python代码，但是我觉得写得过于繁琐（简单的几步操作居然写了300多行？），太注重形式，太刻意地去遵循代码规范，反倒增加了阅读负担，不太合我心意。

于是我按照自己对算法的理解写了一份代码，进行了精简，只保留最重要的部分，并且加入了L2正则化的部分。为了测试代码的正确性，这里使用的神经网络的层数和神经元数和Coursera上的是一样的，共有四层，每层的神经元数量分别是20/ 7/ 5/ 1，用的数据也是课程的数据。

下面进入代码部分。

首先要调用numpy。为什么要使用numpy？因为不想使用一层又一层的for循环，使用向量化的计算能够大幅度地降低运算时间。

准备好激活函数，隐层统一使用relu，输出层使用sigmoid：

def sigmoid(z):
    '''
    z为prev_activation, size为 nl * m
    '''
    return 1 / (1 + np.exp(-z))

def relu(z):
    '''
    z为prev_activation, size为 nl * m
    '''
    return np.maximum(0, z)

设置好神经网络的结构：

n0, m = X.shape
n1 = 20
n2 = 7
n3 = 5
n4 = 1
layers_dims = [n0, n1, n2, n3, n4] # [12288, 20, 7, 5, 1] 
L = len(layers_dims) - 1  # 4层神经网络，不计输入层

构建神经网络模型：

##### neural network model
def neural_network(X, Y, learning_rate=0.01, num_iterations=2000, lambd=0):
    m = X.shape[1]
    ### initialize forward propagation
    param_w = [i for i in range(L+1)]
    param_b = [i for i in range(L+1)]
    np.random.seed(10)
    for l in range(1, L+1):
        if l < L:
            param_w[l] = np.random.randn(layers_dims[l], layers_dims[l - 1]) * np.sqrt(2 / layers_dims[l - 1])
        if l == L:
            param_w[l] = np.random.randn(layers_dims[l], layers_dims[l - 1]) * 0.01
        param_b[l] = np.zeros((layers_dims[l], 1))

    activations = [X, ] + [i for i in range(L)]
    prev_activations = [i for i in range(L+1)]

    dA = [i for i in range(L+1)]
    dz = [i for i in range(L+1)]
    dw = [i for i in range(L+1)]
    db = [i for i in range(L+1)]

    for i in range(num_iterations):
        ### forward propagation
        for l in range(1, L+1):
            prev_activations[l] = np.dot(param_w[l], activations[l-1]) + param_b[l]
            if l < L:
                activations[l] = relu(prev_activations[l])
            else:
                activations[l] = sigmoid(prev_activations[l])

        cross_entropy_cost = -1/m * (np.dot(np.log(activations[L]), Y.T) \
                                     + np.dot(np.log(1-activations[L]), 1-Y.T))
        regularization_cost = 0
        for l in range(1, L+1):
            regularization_cost += np.sum(np.square(param_w[l])) * lambd/(2*m)
        cost = cross_entropy_cost + regularization_cost

        ### initialize backward propagation
        dA[L] =  np.divide(1-Y, 1-activations[L]) - np.divide(Y, activations[L])
        assert dA[L].shape == (1, m)

        ### backward propagation
        for l in reversed(range(1, L+1)):
            if l == L:
                dz[l] = dA[l] * activations[l] * (1-activations[l])
            else:
                dz[l] = dA[l].copy()
                dz[l][prev_activations[l] <= 0] = 0

            dw[l] = 1/m * np.dot(dz[l], activations[l-1].T) + param_w[l] * lambd/m
            db[l] = 1/m * np.sum(dz[l], axis=1, keepdims=True)
            dA[l-1] = np.dot(param_w[l].T, dz[l])

            assert dz[l].shape == prev_activations[l].shape
            assert dw[l].shape == param_w[l].shape
            assert db[l].shape == param_b[l].shape
            assert dA[l-1].shape == activations[l-1].shape

            param_w[l] = param_w[l] - learning_rate * dw[l]
            param_b[l] = param_b[l] - learning_rate * db[l]

        if i % 100 == 0:
            print("cost after iteration {}: {}".format(i, cost))

Andrew Ng教授是用一个dict来保存每层神经元的参数的，比如说，在调用第三层的参数W3和b3时，他的写法是:parameters['W' + str(3)]和parameters['b' + str(3)]，这样写没有错误，虽然直观，但是很麻烦，我的做法是，分别使用list来保存W和b，根据位置读取，对应上面的就是param_w[3]和param_b[3]。

注意到，python（以及其他编程语言）是从0开始计数的，而不是从1开始，这意味着，当神经网络共有4层的时候，我的param_w和param_b的长度是5，我对activation、prev_activation以及各个梯度dA/dz/dw/db都是用长度为5的list来保存的，而Ng教授记录这些变量的长度都是4。

我认为我的写法是更容易理解的写法，因为整个神经网络包括输入层在内一共有5层，但是因为习惯上我们不计输入层，所以这是个4层网络，如果我们把输入层称为第0层，输出层称为第4层，没有什么问题，并且恰好符合了编程语言从0开始计数的习惯。所以在我的写法下，activations[3]就表示第三层的输出，param_b[2]表示第二层的偏置向量，不会有什么误解，很直观。

而对于Ng教授的做法，获取每个参数时是基于字符串的，比如说grads["db" + str(4)]表示第四层的b的梯度，而在用for循环遍历每一层神经元的时候，又是基于位置的，这么一来，你就会在+1 和 -1之后迷失自我，即便你确保了自己没有出错，你也已经花费了不少精力来判断这个地方究竟是该+1还是-1还是保持原样。

下面是我在做Ng教授布置的编程作业时，需要填写的反向传播的部分，要求我在里面填入5行代码来完成，这导致了我在 l 和 l+1和 l-1三者之间纠结犹豫了很久才终于填写正确。我个人认为不是一段user-friendly的代码。

for l in reversed(range(L - 1)):
    # lth layer: (RELU -> LINEAR) gradients.
    # Inputs: "grads["dA" + str(l + 1)], current_cache". 
    # Outputs: "grads["dA" + str(l)] , grads["dW" + str(l + 1)] , grads["db" + str(l + 1)]

    ### START CODE HERE ### (approx. 5 lines)
    current_cache = caches[l]
    dA_prev_temp, dW_temp, db_temp = linear_activation_backward(grads['dA' + str(l + 1)], current_cache, activation="relu")
    grads["dA" + str(l)] = dA_prev_temp
    grads["dW" + str(l + 1)] = dW_temp
    grads["db" + str(l + 1)] = db_temp
    ### END CODE HERE ###

Ng教授的代码中，对于不同的子函数，层数L时而等于5，时而等于4，很是混乱，而我的层数L始终等于4，符合python从0开始计数的习惯，使得整个算法写起来很快，也很好理解，一看就明白当前迭代到了哪一层。

另外，在写算法的时候，要重点关注每个参数的size，注意看它是几乘几的矩阵，很多时候出bug都来自于参数的size弄错了，下面是各参数的size表，里面的12288和209表示我们使用的训练集的数据有209个样本（图片），每个样本有12288个特征（每张图片是64像素*64像素的，所以有64*64*3=12288）。

关于size的规律是，梯度和参数的size相同，同一层的W和dW的size相同，b和db的size相同，其实也很好理解，看参数更新的公式就知道了；另外，同一层的A/z/dA/dz的size都是相同的，这个要记住。可以这么理解，原数据的size是p×m，即m个p维的样本（注意在神经网络中，每一列代表一个样本，和其他机器学习方法中每一行代表一个样本是反过来的），经过每一层的神经元处理之后，得到了新的数据，可以理解为是降维的过程，由一开始的12288×209的数据，通过第一层后变为了20×209的数据，即把高维的12288个特征浓缩为20个新特征，再浓缩为7个特征，最后得到了1×209的预测值。

最后是预测函数，将训练好的神经网络用于一套新的数据上：

def predict(X_new, parameters, threshold=0.5):
    param_w = parameters["param_w"]
    param_b = parameters["param_b"]

    activations = [X_new, ] + [i for i in range(L)]
    prev_activations = [i for i in range(L + 1)]
    m = X_new.shape[1]

    for l in range(1, L + 1):
        prev_activations[l] = np.dot(param_w[l], activations[l - 1]) + param_b[l]
        if l < L:
            activations[l] = relu(prev_activations[l])
        else:
            activations[l] = sigmoid(prev_activations[l])
    prediction = (activations[L] > threshold).astype("int")
    return prediction

经验与收获

• 这是我第一次用python写一个具体的算法，虽然一直在python做数据分析，但是一直没有写过一个逻辑完整的、能应用于实际场合的算法（自己动手写完整算法的经历，在今天之前，是用R语言写了逐步回归和随机森林），今天算是填补了我在python上的这块空白。
• 使用assert语句来确保每个参数的shape正确，能够减少出现bug的几率
• numpy中有一些操作和函数是element-wise的，要留意

• 整个算法没一会儿就写完了，但是训练的时候cost一直没有收敛，检查了老半天，才发现是由两个问题导致的：

  • relu的导函数写错了，返回值是0与1没错，但是我用来判断的参数竟然是A而不是z，真是写得莫名其妙；
  • 初始化不正确（我后来才发现初始化是神经网络的关键），对于权重W我一直使用的写法是param_w[l] = np.random.randn(layers_dims[l], layers_dims[l - 1]) * 0.01，导致cost没能实现收敛，应该是由vanishing gradients导致的，后来我使用了He初始化和Xavier初始化，即把上面的0.01改成np.sqrt(2 / layers_dims[l - 1])或者np.sqrt(1 / layers_dims[l - 1])，效果明显。
• 怎么知道自己的最终算法写对了？我和课程使用的是同一套数据，如果使用我的算法，在相同的神经网络结构下，和Ng教授使用的算法得到的精确度差不多，就说明没问题了。
• 比起写代码，理解代码背后的算法内容更重要。