扣扣技术交流群：460189483
 
参考文章如下，计算过程小白可能看不懂，于是做进一步补充，写出该文
 
https://zhuanlan.zhihu.com/p/20282234?utm_source=qq&utm_medium=social
 

 什么是麦克纳姆轮
 在竞赛机器人和特殊工种机器人中，全向移动经常是一个必需的功能。「全向移动」意味着可以在平面内做出任意方向平移同时自转的动作。为了实现全向移动，一般机器人会使用「全向轮」（Omni Wheel）或「麦克纳姆轮」（Mecanum Wheel）这两种特殊轮子。
 
全向轮：
 
 
 
 
麦克纳姆轮：
 
 
 
 
全向轮与麦克纳姆轮的共同点在于他们都由两大部分组成：轮毂和辊子（roller）。轮毂是整个轮子的主体支架，辊子则是安装在轮毂上的鼓状物。全向轮的轮毂轴与辊子转轴相互垂直，而麦克纳姆轮的轮毂轴与辊子转轴呈 45° 角。理论上，这个夹角可以是任意值，根据不同的夹角可以制作出不同的轮子，但最常用的还是这两种。
 
全向轮与麦克纳姆轮（以下简称「麦轮」）在结构、力学特性、运动学特性上都有差异，其本质原因是轮毂轴与辊子转轴的角度不同。经过分析，二者的运动学和力学特性区别可以通过以下表格来体现。
 
 
 
 
计算过程如下，供参考，学霸可点开大图验算：
 
 
 
 
近年来，麦轮的应用逐渐增多，特别是在 Robocon、FRC 等机器人赛事上。这是因为麦克纳姆轮可以像传统轮子一样，安装在相互平行的轴上。而若想使用全向轮完成类似的功能，几个轮毂轴之间的角度就必须是 60°，90° 或 120° 等角度，这样的角度生产和制造起来比较麻烦。所以许多工业全向移动平台都是使用麦克纳姆轮而不是全向轮，比如这个国产的叉车： 全向移动平台 麦克纳姆轮叉车 美科斯叉车
 
另外一个原因，可能是麦轮的造型比全向轮要酷炫得多，看起来有一种不明觉厉的感觉……
 
 
的确，第一次看到麦轮运转起来，不少人都会惊叹。以下视频直观地说明了麦轮底盘在平移和旋转时的轮子旋转方向。
 

 
 
麦克纳姆轮工作原理
 

 
 
【物理篇-力学专题】E03 S1小车为什么横着走？~滚动摩擦与麦克纳姆轮
 
 
 

 
 
【初中-物理】E08 速度解算？绕圆运动？麦克纳姆轮-进阶
 
麦轮的安装方法
 
麦轮一般是四个一组使用，两个左旋轮，两个右旋轮。左旋轮和右旋轮呈手性对称，区别如下图。
 
 
 
 
安装方式有多种，主要分为：X-正方形（X-square）、X-长方形（X-rectangle）、O-正方形（O-square）、O-长方形（O-rectangle）。其中 X 和 O 表示的是与四个轮子地面接触的辊子所形成的图形；正方形与长方形指的是四个轮子与地面接触点所围成的形状。
 
 
X-正方形：轮子转动产生的力矩会经过同一个点，所以 yaw 轴无法主动旋转，也无法主动保持 yaw 轴的角度。一般几乎不会使用这种安装方式。
X-长方形：轮子转动可以产生 yaw 轴转动力矩，但转动力矩的力臂一般会比较短。这种安装方式也不多见。
O-正方形：四个轮子位于正方形的四个顶点，平移和旋转都没有任何问题。受限于机器人底盘的形状、尺寸等因素，这种安装方式虽然理想，但可遇而不可求。
O-长方形：轮子转动可以产生 yaw 轴转动力矩，而且转动力矩的力臂也比较长。是最常见的安装方式。 
麦轮底盘的正逆运动学模型
 
以O-长方形的安装方式为例，四个轮子的着地点形成一个矩形。正运动学模型（forward kinematic model）将得到一系列公式，让我们可以通过四个轮子的速度，计算出底盘的运动状态；而逆运动学模型（inverse kinematic model）得到的公式则是可以根据底盘的运动状态解算出四个轮子的速度。需要注意的是，底盘的运动可以用三个独立变量来描述：X轴平动、Y轴平动、yaw 轴自转；而四个麦轮的速度也是由四个独立的电机提供的。所以四个麦轮的合理速度是存在某种约束关系的，逆运动学可以得到唯一解，而正运动学中不符合这个约束关系的方程将无解。
 
先试图构建逆运动学模型，由于麦轮底盘的数学模型比较复杂，我们在此分四步进行：
 
①将底盘的运动分解为三个独立变量来描述；
 
②根据第一步的结果，计算出每个轮子轴心位置的速度；
 
③根据第二步的结果，计算出每个轮子与地面接触的辊子的速度；
 
④根据第三部的结果，计算出轮子的真实转速。
   
一、底盘运动的分解
 
我们知道，刚体在平面内的运动可以分解为三个独立分量：X轴平动、Y轴平动、yaw 轴自转。如下图所示，底盘的运动也可以分解为三个量：
 
        表示 X 轴运动的速度，即左右方向，定义向右为正，这是人为设置的底盘水平方向的速度，是已知量；
 
        表示 Y 轴运动的速度，即前后方向，定义向前为正，这是人为设置的底盘垂直方向的速度，是已知量；
 
        表示 Z轴自转的角速度，定义逆时针为正，这是人为设置的底盘原地旋转的速度，是已知量；
 
       以上三个量一般都视为四个轮子的几何中心（矩形的对角线交点）的速度,也是底盘控制时设置的已知量
 
 
 
 
                                                         
 
二、计算出轮子轴心位置的速度（以右上角轮子为例）
 
 
定义：
 
       表示几何中心的合速度（不包括W）向量，由已知量 和 组成的向量合速度，不一定与在一条线上；
 
       为从几何中心指向轮子轴心的矢量；
 
       为轮子轴心的运动速度矢量；
 
      为轮子轴心沿垂直于   的方向（即切线方向)的速度分量；
 
    如下图所示: （手动绘制部分不喜勿喷，怕个别小白不理解，专门增加该部分讲解）
 
                   
 
那么可以计算出轮子中心速度：  =   +  =   +   x 
 
分别计算 几何中心X、Y 轴的分量为：
 
  =   -       
 
 =   + 
 
其中 、表示为在X轴、Y轴的速度数值分量，表示的值大小，不是向量
 
 与 表示 r在 X轴、Y轴的投影长度，也就是两个轮子之间的水平距离（宽度）一半、垂直距离（长度）的一半
 
 =  X cosθ  =  X  = W x r X = 
 
 =  x sinθ =  x = w x r x= 
 
所以可求得几何中心X、Y轴速度数值（非向量）分别为：
 
  =   -      =   - 
 
 =   +  =  + 
 
同理可以算出其他三个轮子轴心的速度。
 
                               
 
三、计算辊子的速度
 
 
根据轮子轴心的速度，可以分解出沿辊子方向的速度 和垂直于辊子方向的速度。其中 是可以无视的（思考题：为什么垂直方向的速度可以无视？），而
 
 =  *    此处表示向量点乘向量，其中  是沿辊子方向的单位向量。这一步表示向量在向量上的投影，是数值大小，不是向量噢，这样就可以求得沿辊子方向的速度数值大小。例如向量a*b表示向量b在向量a上的投影，是一个长度单位，不是向量
 
向量几何知识了解参考 https://www.jianshu.com/p/6b37baa326ec
 
 
 =   +    ，其中、表示、的向量，上面我们求出来的是、的数值，那向量如何表示呢？
 
可以通过膜*单位方向向量实现，如下图所示，设X轴单位方向向量，坐标就是(1,0)，Y轴单位方向向量，坐标就是(0,1)，那么的方向向量就可以表示为*，的方向向量就可以表示为*，且*=1，*=1，原因是两个向量点乘结果是数值，不是向量噢
 
是单位向量，可以表示为（-，），或者 = -* + *
 
 =  *  =（* + *）*（-* + *）
 
 = -*** +  * ** =  -*  +  * 
 
所以最终求得沿辊子方向的速度数值大小 =  -*  +  * 
 
                            
 
四、计算轮子的速度
 
从与地面接触的辊子速度到轮子线转速的计算比较简单：原因是辊子与轮子夹角是45°，
 
 =   =（ -*  +  * ）* = - + 
 
根据上面求出的           =   - ，        =   + ，由此可得
 
   =-（  - ）+ （ + ） =   -  +  (    + )
 
这样求出来的是轮子与地面接触的一点的线速度，如果要求轮子转速(单位rpm），需要再除以轮子半径
 
    =    ，单位是弧度/秒，转化到rpm（转/分钟），需要30*/PI  ,到此，轮子速度已经求出来啦
 
          
 
根据上图定义，可知 = a ,  = b，由此简化可得轮子线速度为：
 
    =   -  +  
 
由此计算获得四个轮子线速度如下
 
    =   -  +  
 
    =   +  -  
 
    =   -  -  
 
    =   +  +  
 
以上方程组就是O-长方形麦轮底盘的逆运动学模型，而正运动学模型可以直接根据逆运动学模型中的三个方程解出来，此处不再赘述。
 
另一种计算方式
 
「传统」的推导过程虽然严谨，但还是比较繁琐的。这里介绍一种简单的逆运动学计算方式。
 
我们知道，全向移动底盘是一个纯线性系统，而刚体运动又可以线性分解为三个分量。那么只需要计算出麦轮底盘在「沿X轴平移」、「沿Y轴平移」、「绕几何中心自转」时，四个轮子的速度，就可以通过简单的加法，计算出这三种简单运动所合成的「平动+旋转」运动时所需要的四个轮子的转速。而这三种简单运动时，四个轮子的速度可以通过简单的测试，或是推动底盘观察现象得出。
 
当底盘沿着 X 轴平移时：
 
 =  -  
 
 =  +  
 
 =  -  
 
 =  +  
 
当底盘沿着 Y 轴平移时：
 
 =  +  
 
 =  +  
 
 =  +  
 
 =  +  
 
当底盘绕几何中心自转时：
 
    = +  
 
    =  -  
 
    =  -  
 
    =  +  
 
将以上三个方程组相加，得到的恰好是根据「传统」方法计算出的结果。这种计算方式不仅适用于O-长方形的麦轮底盘，也适用于任何一种全向移动的机器人底盘。
 
Makeblock 麦轮底盘的组装
 
理论分析完成，可以开始尝试将其付诸实践了。
 
第一步，组装矩形框架。
 
 
第二步，组装电机模块。
 
由于麦轮底盘的四个轮子速度有约束关系，必须精确地控制每个轮子的速度，否则将会导致辊子与地面发生滑动摩擦，不仅会让底盘运动异常，还会让麦轮的寿命减少。所以必须使用编码电机。
 
 
 
 
第三步，将电机模块安装到框架上。
 
 
 
 
 
 
第四步，将麦轮安装到框架上。
 
 
 
 
第五步，安装电路板并接线。
 
编码电机必须配上相应的驱动板才能正常工作。这里使用的 Makeblock 编码电机驱动板，每一块板可以驱动两个电机。接线顺序在下文中会提及，也可以随意接上，在代码中定义好对应的顺序即可。
 
 
 
 
第六步，装上电池。
 
 
 
 
至此，一个能独立运行的麦轮底盘就完成了。
 
 
 
控制程序
 
根据麦轮的底盘的运动学模型，要完全控制它的运动，需要有三个控制量：X轴速度、Y轴速度、自转角速度。要产生这三个控制量，有很多种方法，本文将使用一个 USB 游戏手柄，左边的摇杆产生平移速度，右边的摇杆产生角速度。
 
首先将一个 USB Host 模块连接到 Orion 主板的 3 口。
 
。
 
 
 
 
然后插上一个无线 USB 游戏手柄。
 
 
 
 
 
 
然后再添加其他细节，就大功告成啦！
 
 
 
 
 
 
其他细节：
 
 
 
#include <Wire.h>
#include <SoftwareSerial.h>
#include "MeOrion.h"

MeUSBHost joypad(PORT_3);
//	手柄代码（红灯亮模式）
//	默认：128-127-128-127-15-0-0-128
//	左一：128-127-128-127-15-1-0-128
//	右一：128-127-128-127-15-2-0-128
//	左二：128-127-128-127-15-4-0-128
//	右二：128-127-128-127-15-8-0-128
//	三角：128-127-128-127-31-0-0-128 （0001 1111）
//	方形：128-127-128-127-143-0-0-128 （1000 1111）
//	叉号：128-127-128-127-79-0-0-128 （0100 1111）
//	圆圈：128-127-128-127-47-0-0-128 （0010 1111）
//	向上：128-127-128-127-0-0-0-128 （0000 0000）
//	向下：128-127-128-127-4-0-0-128 （0000 0100）
//	向左：128-127-128-127-6-0-0-128 （0000 0110）
//	向右：128-127-128-127-2-0-0-128 （0000 0010）
//	左上：128-127-128-127-7-0-0-128 （0000 0111）
//	左下：128-127-128-127-5-0-0-128 （0000 0101）
//	右上：128-127-128-127-1-0-0-128 （0000 0001）
//	右下：128-127-128-127-3-0-0-128 （0000 0011）
//	选择：128-127-128-127-15-16-0-128
//	开始：128-127-128-127-15-32-0-128
//	摇杆：右X-右Y-左X-左Y-15-0-0-128


MeEncoderMotor motor1(0x02, SLOT2);
MeEncoderMotor motor2(0x02, SLOT1);
MeEncoderMotor motor3(0x0A, SLOT2);
MeEncoderMotor motor4(0x0A, SLOT1);

//  底盘：a = 130mm, b = 120mm

float linearSpeed = 100;
float angularSpeed = 100;
float maxLinearSpeed = 200;
float maxAngularSpeed = 200;
float minLinearSpeed = 30;
float minAngularSpeed = 30;

void setup()
{
    //  要上电才能工作，不能只是插上 USB 线来调试。
	motor1.begin();
	motor2.begin();
	motor3.begin();
	motor4.begin();
	
	Serial.begin(57600);
	joypad.init(USB1_0);
}

void loop()
{
    Serial.println("loop:");
	//setEachMotorSpeed(100, 50, 50, 100);
	if(!joypad.device_online)
    {
        //  若一直输出离线状态，重新拔插 USB Host 的 RJ25 线试一下。
        Serial.println("Device offline.");
        joypad.probeDevice();
        delay(1000);
    }
    else
    {
        int len = joypad.host_recv();
        parseJoystick(joypad.RECV_BUFFER);
        delay(5);
    }
	//delay(500);
}


void setEachMotorSpeed(float speed1, float speed2, float speed3, float speed4)
{
	motor1.runSpeed(speed1);
	motor2.runSpeed(-speed2);
	motor3.runSpeed(-speed3);
	motor4.runSpeed(-speed4);
}

void parseJoystick(unsigned char *buf)   //Analytic function, print 8 bytes from USB Host
{
    //  输出手柄的数据，调试用
    // int i = 0;
    // for(i = 0; i < 7; i++)
    // {
    //     Serial.print(buf[i]);  //It won't work if you connect to the Makeblock Orion.
    //     Serial.print('-');
    // }
    // Serial.println(buf[7]);
    // delay(10);

    //  速度增减
    switch (buf[5])
    {
        case 1:
            linearSpeed += 5;
            if (linearSpeed > maxLinearSpeed)
            {
                linearSpeed = maxLinearSpeed;
            }
            break;
        case 2:
            angularSpeed += 5;
            if (angularSpeed > maxAngularSpeed)
            {
                angularSpeed = maxAngularSpeed;
            }
            break;
        case 4:
            linearSpeed -= 5;
            if (linearSpeed < minLinearSpeed)
            {
                linearSpeed = minLinearSpeed;
            }
            break;
        case 8:
            angularSpeed -= 5;
            if (angularSpeed < minAngularSpeed)
            {
                angularSpeed = minAngularSpeed;
            }
            break;
        default:
            break;
    }


    
    if ((128 != buf[0]) || (127 != buf[1]) || (128 != buf[2]) || (127 != buf[3]))
    {
        //  处理摇杆
        float x = ((float)(buf[2]) - 127) / 128;
        float y = (127 - (float)(buf[3])) / 128;
        float a = (127 - (float)(buf[0])) / 128;
        mecanumRun(x * linearSpeed, y * linearSpeed, a * angularSpeed);
    }
    else
    {
        switch (buf[4])
        {
            case 0:
                mecanumRun(0, linearSpeed, 0);
                break;
            case 4:
                mecanumRun(0, -linearSpeed, 0);
                break;
            case 6:
                mecanumRun(-linearSpeed, 0, 0);
                break;
            case 2:
                mecanumRun(linearSpeed, 0, 0);
                break;
            case 7:
                mecanumRun(-linearSpeed/2, linearSpeed/2, 0);
                break;
            case 5:
                mecanumRun(-linearSpeed/2, -linearSpeed/2, 0);
                break;
            case 1:
                mecanumRun(linearSpeed/2, linearSpeed/2, 0);
                break;
            case 3:
                mecanumRun(linearSpeed/2,  -linearSpeed/2, 0);
                break;
            default:
                mecanumRun(0, 0, 0);
                break;
        }
    }
}

void mecanumRun(float xSpeed, float ySpeed, float aSpeed)
{
    float speed1 = ySpeed - xSpeed + aSpeed; 
    float speed2 = ySpeed + xSpeed - aSpeed;
    float speed3 = ySpeed - xSpeed - aSpeed;
    float speed4 = ySpeed + xSpeed + aSpeed;
    
    float max = speed1;
    if (max < speed2)   max = speed2;
    if (max < speed3)   max = speed3;
    if (max < speed4)   max = speed4;
    
    if (max > maxLinearSpeed)
    {
        speed1 = speed1 / max * maxLinearSpeed;
        speed2 = speed2 / max * maxLinearSpeed;
        speed3 = speed3 / max * maxLinearSpeed;
        speed4 = speed4 / max * maxLinearSpeed;
    }
    
    setEachMotorSpeed(speed1, speed2, speed3, speed4);
}
 

 
 
Makeblock 麦克纳姆轮 全向移动机器人
 
 PS：答疑
 
1. 这轮子这么灵活，为啥不都装成麦克纳姆轮呢，多方便啊
 
答：1. 如果上面的算法你掌握了，那么就明白所有的理论计算都是基于小车轮子在一个平面的，那么如果有悬挂的，或者地面不平坦的时候，很明显计算过程中的速度会有偏差，就会导致你的小车底盘跑的不准、不直，如果你因此计算里程计就更差了。2. 效率上不行，为啥呢？记得上面那个问题吗？棍子垂直方向的速度不考虑，为啥？因为这个速度主要用于棍子自转，对小车运行没有影响，那这就势必造成了能量浪费（相比普通轮子），所以效率就降低了 。3 .因为结构特性原因，不利于爬坡、越障，所以比较适合于平地使用，大大限制了应用场景
 
2. 程序中为啥有速度反向
 
因为电机运行方向的原因

*注:本博客参考李宏毅老师2020年机器学习课程. 视频链接
 


 
目录
 
1 反向传播算法数学基础
1.1 前向传播
1.2 反向传播
  
2 矩阵运算
3 模型实现
3.1 模型结构设计
3.2 模型测试
 

 
 
1 反向传播算法数学基础
 
反向传播算法是用于对神经网络中的各个网络参数计算偏导值的一种算法，其核心是链式求导法则。
 注：本节涉及到较复杂的数学计算，了解思想即可。
 

 (图一)
 
在图一中，画出了一个神经网络的前两层结构的部分神经元。图中，$x_i$表示网络的输入。$w_{ij}$表示第i层神经元的参数，j仅用作区分不同的参数。$z_{ij}$表示第i层神经元经过激活函数之前的值，$a_{ij}$表示第i层第j个神经元的输出。$C$表示网络计算得到的损失函数值（Cost）.由此可知：
 $$z_{11}=\sum _{j=1}^n{w}_{1j}{x}_j\text{\hspace{1em}(1)}$$
 $$a=\sigma \left(z\right)\text{\hspace{1em}(2)}$$
 其中$\sigma$表示激活函数。

 如果要计算$\frac{\partial C}{\partial w_{11}}$，根据链式求导法则，有：
 $$\frac{\partial C}{\partial w_{11}}=\frac{\partial z_{11}}{\partial w_{11}}\frac{\partial C}{\partial z_{11}}\text{\hspace{1em}(3)}$$
 在等式3的右边，计算前一项的过程称为前向传播，计算后一项的过程称为反向传播。
 
1.1 前向传播
 
前向传播即计算$\frac{\partial z_{11}}{\partial w_{11}}$，根据式1可知，$\frac{\partial z_{11}}{\partial w_{11}}=x_1$，实际上：
 $$\frac{\partial z_{2i}}{\partial w_{2i}}=a_{11},i=1,2,...m\text{\hspace{1em}(4)}$$
 照此推理，只要知道上一层神经元的输出，就可以求得这一层神经元的输入z对这一层神经元的参数w的偏导数值。在顺序计算损失函数值得过程中，所有的$\frac{\partial z}{\partial w}$都可以计算完毕，因此称为前向传播。
 
1.2 反向传播
 
$\frac{\partial C}{\partial z_{11}}$的计算要稍微复杂一些，首先可以根据式2，得：
 $$\frac{\partial C}{\partial z_{11}}=\frac{\partial C}{\partial a_{11}}{\sigma }^{\prime }(z_{11})\text{\hspace{1em}(5)}$$
 由于$z_{11}$是一个已知数值，因此${\sigma }^{\prime }(z_{11})$是一个常数。对于上式等号右侧得第一项，从图1中可以看出，$a_{11}$得变化首先会改变$z_{2j},j=1,2,...,m$的值变化，最终导致C的值变化，那么：
 $$\frac{\partial C}{\partial a_{11}}=\sum _{i=1}^m\frac{\partial C}{\partial z_{2i}}{w}_{2i}\text{\hspace{1em}(6)}$$
 将式5和式6联立化简得：
 $$\frac{\partial C}{\partial z_{11}}={\sigma }^{\prime }(z_{11})\sum _{i=1}^m\frac{\partial C}{\partial z_{2i}}{w}_{2i}\text{\hspace{1em}(7)}$$
 上式中的${\sigma }^{\prime }(z_{11})$和$w_{2i}$的值均为已知，于是不难发现，式7的定义是递归的，也就是说，要计算损失函数值C对某一层神经元的输入$z_i$的偏导值，必须先知道C对下一层神经元的输入$z_{i+1}$的偏导值。那么自然会想到从最后一层开始，向前计算。
 
 (图2)
 
在图2中给出了神经网络的最后一层的结构。
 $$y_i=\frac{1}{1+e^{-z_i}},i=1,2,...,t\text{\hspace{1em}(8)}$$
 当使用交叉熵作为损失函数时：
 $$C=-\sum _{i=1}^t(\widehat{y_i}\ln y_i+(1-\widehat{y_i})\ln (1-y_i))\text{\hspace{1em}(9)}$$
 则
 $$\frac{\partial C}{\partial z_1}=\frac{\partial C}{\partial y_1}\frac{\partial y_1}{\partial z_1}=-(\frac{\widehat{y_1}}{y_1}-\frac{1-\widehat{y_1}}{1-y_1}){\sigma }^{\prime }(y_1)=y_1-\widehat{y_1}\text{\hspace{1em}(10)}$$
 综上，结合式3、4、7、10，便可以求得网络中任何一个神经元对于其参数w的偏导值。对于式3，我们可以写出它的表达式：
 $$\frac{\partial C}{\partial w_{11}}=x_1{\sigma }^{\prime }(z_{11})\sum _{i=1}^m\frac{\partial C}{\partial z_{2i}}{w}_{2i}\text{\hspace{1em}(11)}$$
 
2 矩阵运算
 
下面我们将上述计算过程写成矩阵运算的形式，这种形式更适合在实际的编码过程中使用各种矩阵运算库进行加速运算（如numpy、tensorflow、PyTorch等）。
 
首先回到图1，我们将神经网络的输入写成$X_1\in {\mathbb{R}}^{n\times 1}$，表示它是一个n行1列的矩阵，对应输入具有n个维度。
对于神经网络的第一层，假设第一层有k个神经元，则每一个神经元都应该具有n个参数w。因此将第一层神经元的所有参数w写成矩阵$W_1\in {\mathbb{R}}^{k\times n}$，表示它是一个k行n列的矩阵，其中每一行表示同一个神经元对输入$X$的不同维度的对应参数，每一列表示不同神经元对输入$X$的相同维度的对应参数，下标1代表它是第一个隐藏层的参数。
同理，第一层的每一个神经元还有一个参数b，将该曾所有的b写成矩阵为$B_1\in {\mathbb{R}}^{k\times 1}$。
第一个隐藏层的激活函数前的输出$Z$也是一个矩阵，表示为$Z_1\in {\mathbb{R}}^{k\times 1}$，且$Z_1=W_1{X}_1+B_1$。该层的输出表示为$A_1\in {\mathbb{R}}^{k\times 1}$，且$A_1=\sigma (Z_1)$，同时该层的输出也是第二层的输入，因此$A_1=X_2$
其余各层的运算与上述过程类似，不再赘述。
对于神经网络的最后一层，如图2，将得到整个网络的输出$Y$，与之前不同的是，将激活函数$\sigma$换成了Softmax函数，假设网络输出包含t维，那么$Y\in {\mathbb{R}}^{t\times 1}$，$Y$的各个维度的值的计算遵循公式8。
 
因此，矩阵运算计算第二层输入的过程可以表示如下：
 

 (图3)
 
其中，dot表示矩阵乘法。在前向传播的过程中，不断重复以上过程，就可以得到所有的层的输入$X$。该过程如下图4，前向传播的任务就是完成矩阵$X$的计算：
 

 (图4)
 
将损失函数$C$对第i层(非最后一层)的参数$W_i$求导的计算式11加以扩展，改写为矩阵运算形式如下：
 $$\frac{\partial C}{\partial W_i}=[{\sigma }^{\prime }(Z_i)\ast (W_{i+1}^T\frac{\partial C}{\partial Z_{i+1}})]X_i^T\text{\hspace{1em}(12)}$$
 上式中，$W_i$为第i层的网络权重，$X_i$为第i层神经元的输入，$Z_i=W_i{X}_i+B_i$为第i层神经元激活前的输出值，$\ast$表示元素乘法（即对应位置元素相乘），$T$表示矩阵转置。
 假设该网络的第i层有m个神经元，第i+1层有n个神经元，第i-1层有p个神经元,那么式11中各变量的维度分别如下表：
 
变量
维度
变量
维度
$X_i$
(p,1)
$Z_i$
(m,1)
$Z_{i+1}$
(n,1)
$W_i$
(m,p)
$\frac{\partial C}{\partial W_i}$
(m,p)
$\frac{\partial C}{\partial Z_{i+1}}$
(n,1)
$W_{i+1}^T\frac{\partial C}{\partial Z_{i+1}}$
(m,1)
$W_{i+1}^T$
(m,n)
上述讨论中我们没有考虑对B的偏导数，实际上，如果将式3中等号右边第一项替换为$\frac{\partial z_{11}}{\partial b_{11}}$即可得出对B的偏导，并且这一项始终为1，因此：
 $$\frac{\partial C}{\partial B_i}={\sigma }^{\prime }(Z_i)\ast (W_{i+1}^T\frac{\partial C}{\partial Z_{i+1}})\text{\hspace{1em}(13)}$$
 
3 模型实现
 
根据上述数学推理，设计并实现一个神经网络已经变得可能。下面我们开始着手设计实现这个数据结构。
 
3.1 模型结构设计
 
首先明确我们要设计的神经网络的结构，假设该模型有m层，m最小为3，因为除了输入层与输出层之外，网络应该包含至少一个隐藏层。
 每个隐藏层都有自己的W和B。若第$i\in [1,m-1]$层有$n_i$个神经元，那么其参数满足：$W_i\in {\mathbb{R}}^{n_i\times n_{i-1}},B_i\in {\mathbb{R}}^{n_i\times 1}$，整个网络的参数一共2(m-2)个矩阵（其中m-2个矩阵为权重W的矩阵，m-2个矩阵为偏置B的矩阵）。网络结构定义如下：
 
import numpy as np
class Network():
    def __init__(self, *layers):
        # 至少包含输入层，输出层，和一个隐藏层，所以至少有三层
        assert(len(layers) > 2)
        # W的维度为(当前层神经元数，上一层神经元数)
        self.W = [np.random.randn(y, x)
                  for x, y in zip(layers[:-1], layers[1:])]
        # B的维度为(当前层神经元数，1)
        self.B = [np.random.randn(y, 1)
                  for x, y in zip(layers[:-1], layers[1:])]
 
根据式12可知，在前向传播的过程中，需要记录$X_i$和$Z_i$的值，以方便计算偏导值。
 
def forward(network, X):
    Zs = []
    Xs = []
    for W, B in zip(network.W, network.B):
        Xs.append(X)
        Z = W.dot(X)+B
        X = sigmoid(Z)
        Zs.append(Z)
    Y = sigmoid(Z)
    return Xs, Zs, Y
# sigmoid函数
def sigmoid(Z):
    A = 1/(1+np.exp(-Z))
    return A
 
反向传播的计算$\frac{\partial C}{\partial Z_i}$的表达式就是式12等式右边第一项，计算过程中需要保存所有的$\frac{\partial C}{\partial Z_i}$。
 
def backward(network, pcpz, Zs):
    P_z = [pcpz]
    for W, Z in zip(network.W[::-1][:-1], Zs[::-1][1:]):
        pcpz = sigma_(Z)*(W.T.dot(pcpz))
        P_z.append(pcpz)
    P_z.reverse()
    return P_z
# sigmoid函数的导数
def sigma_(Z):
    return sigmoid(Z)*(1-sigmoid(Z))
 
在开始反向传播之前，首先需要按照式10计算$\frac{\partial C}{\partial Z_i}$的最后一项：
 
def get_last_pcpz(Y0, Y1):
    return Y1-Y0
 
下面给出训练过程的代码：
 
def start_train(train_x, train_y,epoch,lr,*hidden_layers):
    nn = Network(train_x[0].shape[0], *hidden_layers, train_y[0].shape[0])
    loss = [] #记录损失函数值的变化
    for i in range(epoch):
        # 初始化W和B在当轮迭代中的变化量
        delta_W = [np.zeros_like(w) for w in nn.W]
        delta_B = [np.zeros_like(b) for b in nn.B]
        l = 0
        # 对训练数据中的每一个样本进行计算
        for X, Y in zip(train_x, train_y):
            Xs, Zs, Y1 = forward(nn, X) #前向传播
            pcpz = get_last_pcpz(Y, Y1) #计算最后一项
            P_z = backward(nn, pcpz, Zs) #反向传播
            l += np.abs(Y-Y1).sum() #计算当前样本的损失函数值
            # 更新当前轮次的W和B由于该样本所带来的变化量
            delta_B = [db+pz for db, pz in zip(delta_B, P_z)] 
            delta_W = [dw+pz.dot(x.T) for dw, pz, x in zip(delta_B, P_z, Xs)]
        loss.append(l)
        # 更新模型参数
        nn.W = [w-dw*lr/len(train_x) for w, dw in zip(nn.W, delta_W)]
        nn.B = [b-db*lr/len(train_x) for b, db in zip(nn.B, delta_B)]
    return nn, loss
 
3.2 模型测试
 
使用一组简单的数据对所实现的模型进行测试。该组数据由2个样本构成，其值分别为(1,1)和(0,0)，我们将前者归为类别1，后者归为类别2，对于类别1，我们用(1,0)表示，类别2，我们用(0,1)表示，下面对着一组简单的数据进行训练，并给出训练过程中loss的变化情况。
 
import matplotlib.pyplot as plt
train_x = [np.ones((2, 1)), np.zeros((2, 1))]
train_y = [np.asarray([[1], [0]]), np.asarray([[0], [1]])]
nn, loss = start_train(train_x, train_y,3000,1,3,4)
for X, Y in zip(train_x, train_y):
    Xs, Zs, Y1 = forward(nn, X)
    print("input:{} output: {} class:{}".format(X.tolist(),Y1.tolist(),np.argmax(Y1)+1))
plt.plot(loss)
plt.show()
 
input:[[1.0], [1.0]] output: [[0.9989013915763175], [0.00127671782739844]] class:1
input:[[0.0], [0.0]] output: [[0.0008225495145913559], [0.9990508703257248]] class:2
 

一、梯度下降
 
1. 原理讲解
 
        在讲梯度下降之前，先向大家推荐一门课程《数值分析》，有的书也将它称为《计算方法》。
 
 
 
数值分析(numerical analysis)，为数学的一个分支，是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科。它以数字计算机求解数学问题的理论和方法为研究对象，为计算数学的主体部分。
 
 
        这门课程详细得介绍了如何用计算机求解积分、微分、导数、方程根、微分方程、一元线性方程组、函数极值。那么我们今天要讲解的梯度下降方法就是《数值分析》里面的一个模块。梯度下降算法是用来求一个函数的极小值点。举个例子：
 
求    的函数极小值。函数图像如图：
 
 
 

  
  
 
 
在数学上，求y的极小值点可以通过对y求一次导数，根据导数的零点和函数的单调性能够精确地算出极小值点地值。我们在这里手算以一下：
 
 ,   解出x1 = 3.270083,    x2 = -1.936750。结合函数图像可以知道极小值点是x1。
 
但是，计算机解决的问题是普遍性的，通用性的。这道题的y毕竟是一个普通的多项式函数，如果现在一个函数带上积分、带上三角函数，带上e，带上指数对数函数，此时人工手算的方法就远远不够了。
 
梯度下降方法就是通过迭代的方式，每次进一小步寻找最小值。
 
 
注意到x1的左侧A点导数值小于0，x1的右侧B点导数值大于0，x1处导数值等于0. 随机初始化一个点p，当p在点A时，p = p + δ，当p在B点时，p = p - δ，其中δ = 前进步长，也就是可以写成 p = p - 学习率 * p点导数值。当p点充分接近于x1时，可以停止迭代，那么经过迭代，极小值点就可以找出来。
 
公式化为：  (r是很小的一个数，表示学习率）。当然最后迭代的结果取决于选择的点的位置，当初始点p在A, B之间，则能够迭代出x1，当初始点在c附近，则无法迭代出结果，也就是说算法不会收敛，只会向左一泻千里。当然，数值分析里还有很多算法都取决于初始值的位置。
 
2.代码实现
 
说这么多，我们用代码实现一下这个算法，选用cpp语言。
 
#include<cstdio>
#include<cmath>

using namespace std;

long long N;        //迭代次数
double p;           //初始位置
double r;           //学习率


double ff(double x){        // y的导数
    return 3.0 * pow(x, 2) - 4.0 * x - 19;
}


int main(){
    scanf("%lld%lf%lf",&N, &p, &r);
    for(int i = 0;i<N;i++){
        p = p - r * ff(p);      //p处的导数
    }
    printf("%lf", p);
    return 0;
}
 
运行一下程序，显示执行结果：
 

  程序运行结果
 
N
p
r
极小值点
100
0
0.01
3.270083
100
1
0.05
3.270083
10
0
0.2
-4989812654385.430664
10
-5
0.01
-1447.283969
 
 
从上面的表格可以看到，当我们的参数选择的合适的时候，是可以计算出极小值点x1的，但是如果学习率过大，会出现点p在极值点处左右摇摆的情况，而不能收敛到x1。而当p被“甩”到左侧时，就会出现最后一种情况，函数永远无法收敛。当然，现在有很多方法能够有效避免这种情况，比如学习率随着迭代逐渐改变，而不使用一个固定的学习率，或者随机初始化几个点，然后找到能够收敛的点等等措施。
 
3. 用梯度下降求解线性回归
 
        线性是一个好的特性，线性回归同样也是一个最简单的回归，我们可以用手算或者某种公式就能够算出需要评估的斜率和截距。回顾用手算线性回归方法：已知<x1, y1>, ...... , <xn, yn>，拟合出直线 y = ax + b。
 
<1,6> , <2, 8> , <3, 9.5>, <4, 11.5>, <5, 14>
 
首先确定距离函数L(也称为损失函数）， ，我们要寻找使L取极小值的a和b。把xi和yi看成常数，那么L可以看成关于a和b的函数，也即求L的极小值点。
 
我们有两种方法，一种是手算求导找到函数零点后判断极小值点，另一种方法就是使用我们的梯度下降方法。两种方法都算一遍。计算结果如图：
 
    
 
 
令，解出a和b即可。结果为：
 
  （抄百度百科的）
 
计算出a = 1.95, b = 3.95
 
用梯度下降方法计算：
 
回顾梯度下降公式：
 
 
 
那么每次迭代都小部分更新a和b，当迭代次数充分后，可以计算出a和b的数值解。
 
下面上代码：
 
#include<iostream>
#include<unordered_map>

using namespace std;

/**
    用梯度下降方法求线性回归(不是最小二乘)

    输入：N L, 表示数据结点数和训练轮次
    输出：a, b， 表示y = ax + b

**/

double a, b;          // y = ax + b
double r1, r2;      //a, b 学习率
long long N;        //顶点个数
long long L;        //迭代次数
unordered_map<double, double> data;

void init(){        //初始化
    a = 0.5;        
    b = 0.5;
    r1 = r2 = 0.005;    
}

void train(){
    double na = 0, nb = 0;
    for(pair<double, double> p : data){
        double xi = p.first, y_bar = p.second;
        na += (xi * (a * xi + b - y_bar));
        nb += (a * xi + b - y_bar);
    }
    a = a - (r1 * na);
    b = b - (r2 * nb);
    return;
}

int main(){
    init();
    cin>>N>>L;     

    double x, y;
    for(long long i = 0;i<N;i++){
        cin>>x>>y;
        data[x] = y;
    }

    for(long long i = 0;i<L;i++){
        train();
    }

    cout<<"a = "<<a<<endl<<"b = "<<b<<endl;
    return 0;
}
 
执行结果如图：
 
 
可以看到a = 1.95， b = 3.95，计算结果和解析解一致。 
 
所以我们也可以用梯度下降的方法来求解线性回归（不过都能通过手算出解析解了，为啥还有梯度下降多此一举呢，手动斜眼）。
 
练习1：用梯度下降求解的极值点
 
练习2：用梯度下降求解点集{<0.5,1.455>, <1, 9.777>, <1.5, 29.224>, <2,  403.429>, <2.2, 1422.257>} 拟合中a的结果。
 
二、反向传播
 
反向传播是在神经网络中，求解参数的一个方法，其中用到了梯度下降求解极小值和函数求导的链式法则两个工具。
 
我们先考虑一个极简单的情况，其中神经元内不含激活函数,（也就是激活函数是y = x)：
 
 
 
 
已知输出结果有下面公式：
 
 
损失函数为：   
 
其中 
 
   (1)
 
    (2)
 
我们需要做的是，寻找当L最小的时候, wi,j的值。换句话说，对于函数L，其自变量为wi,j。当然采用梯度下降法。
 
回忆梯度下降法公式（用w1,1举例子）：
 
 
每轮迭代时，都需要计算        ，如何计算呢？这明显是一个多元函数的嵌套，求导符合链式法则。
 
则    ,  等式右边第一项求导后是  ， 等式第二项求导可根据(1)式进行，求导结果是xi,1。
 
那么二者相乘即可算出 
 
从上面的公式中可以看出，每轮迭代时，w1,1的更新公式是
 
同理， 也可以按这种方法，一步一步迭代出来。
 
 
 
考虑激活函数
 
令 
 
激活函数采用sigmoid函数：
 
 
 
还是观察w1,1的更新公式： 
 
对比原来的更新公式，其实就多了一个中间变量zi,1。
 
上式第一项不变，第二项是对sigmoid求导，第三项不变，最后可以整理成一个求解参数式子，即为反向传播。
 
为什么叫反向传播呢？
 
这里只是包含一层的神经网络，没法体现出反向传播的直观，但是如果是深层结构，我们更新完最内层的w后，其求导结果还可以继续使用，然后再通过内层对外层的结果，算出L对外层w的导数，实现了误差函数的结果反向传播。
 
 
 
用C++实现简单的神经网络
 
为了方便，采用随机梯度下降方法，同时隐藏层的激活函数采用sigmoid函数，学习率固定为0.2.
 
下图是首先需要的公式推导，相应的偏导数已经计算出来。
 
 
 
 
代码如下：
 
#include<iostream>
#include<bits/stdc++.h>

/**
    神经网络(2,3,2),2输入，2输出，最后一维是偏置，恒为1，激活函数采用sigmoid，输出不包含softmax

    输入：
    N    L  , N点的个数，L迭代轮次
    然后输入N个(x1,x2,y1,y2)

    然后输入测试样例

    输入如：

20 10000
0 1 1 0
1 2 1 0
1 3 1 0
-1 2 1 0
-1 3 1 0
2 10 1 0
-3 6 1 0
3 7 1 0
0 5 1 0
-1 8 1 0
0 -1 0 1
1 -2 0 1
1 -3 0 1
-1 -2 0 1
-1 -3 0 1
2 -10 0 1
3 -6 0 1
-3 -2 0 1
0 -5 0 1
2 -3 0 1

//测试样例
2 10
2 -5
0 -9
0 4


**/
using namespace std;


int N;          //顶点数
int L;          //学习轮次

double W[3][3];
double U[2][3];
double rate = 0.02;


struct elem{
    double x[3];            //
    double y[2];
};

elem data[10005];

void init(){
    for(int i = 0;i<3;i++){
        for(int j = 0;j<3;j++){
            W[i][j] = 0.5;
        }
    }

    for(int i = 0;i<3;i++){
        U[1][i] = 0.5;
        U[2][i] = 0.5;
    }
}

double h[3];        //记录前向传播的值
double z[3];
double o[3];

void train(int n){          //
    //前向传播
    for(int i = 0;i<3;i++){
        double d = 0;
        for(int j = 0;j<3;j++){
            d += W[i][j] * data[n].x[j];
        }
        h[i] = d;
    }
    for(int i = 0;i<3;i++){
        z[i] = 1.0/(1.0 +pow(exp(1.0),-h[i]));
    }

    for(int i = 0;i<2;i++){
        double d = 0;
        for(int j = 0;j<3;j++){
            d += U[i][j] * z[j];
        }
        o[i] = d;
    }

    //梯度下降更新
    double uu[3][2];        //  dL/du
    fill(uu[0], uu[0]+3 * 2, 0);

    //更新uij
    for(int i = 0;i<3;i++){
        for(int j = 0;j<2;j++){
            uu[i][j] = (o[i]-data[n].y[i]) * z[j];
            U[i][j] = U[i][j] - rate * uu[i][j];
        }
    }

    //更新wij
    double ww[3][3];    //    dL/wij
    fill(ww[0], ww[0]+3*3, 0);


    for(int i = 0;i<3;i++){
        for(int j = 0;j<3;j++){
            for(int k = 0;k<2;k++){
                double tmp = (1.0/(1.0 +pow(exp(1.0),-h[i]))) * (1-((1.0/(1.0 +pow(exp(1.0),-h[i])))));
                ww[i][j] += (o[k]-data[n].y[k])*U[k][i]*tmp*data[n].x[j];
            }
            W[i][j] = W[i][j] - rate * ww[i][j];
        }
    }
    return;
}

void printParam(){
    cout<<"W:"<<endl;
    for(int i = 0;i<3;i++){
        for(int j = 0;j<3;j++){
            cout<<W[i][j]<<" ";
        }
        cout<<endl;
    }
    cout<<endl<<"U:"<<endl;
    for(int i = 0;i<2;i++){
        for(int j = 0;j<3;j++){
            cout<<U[i][j]<<endl;
        }
        cout<<endl;
    }
}

void res(elem e){
    for(int i = 0;i<3;i++){
        double d = 0;
        for(int j = 0;j<3;j++){
            d += W[i][j] * e.x[j];
        }
        h[i] = d;
    }

    for(int i = 0;i<3;i++){
        z[i] = 1.0/(1.0 +pow(exp(1.0),-h[i]));
    }

    for(int i = 0;i<2;i++){
        double d = 0;
        for(int j = 0;j<3;j++){
            d += U[i][j] * z[j];
        }
        o[i] = d;
    }

}

int main(){
    cin>>N>>L;
    for(int i = 0;i<N;i++){
        cin>>data[i].x[0]>>data[i].x[1]>>data[i].y[0]>>data[i].y[1];
        data[i].x[2] = 1;
    }

    init();

    for(int i = 0;i<L;i++){
        for(int j = 0;j<N;j++){
            train(j);
        }
    }

    printParam();
    for(int i = 0 ;i<100;i++){
        elem e;
        cin>>e.x[0]>>e.x[1];
        e.x[2] = 1;
        res(e);
        cout<<o[0]<<"   "<<o[1]<<endl;
    }

    return 0;
}

 
 
 
 
 
 
 
 
 
 
 
 

BP算法是深度学习中最重要也是最基础的算法之一，也是Hinton本人的图灵奖工作之一。更重要的是，BP是深度学习岗位面试必问题之一。本文会以最容易理解的方式来介绍反向传播，避免大量公式。
 

1. 前向传播和反向传播
 
前向还是反向区别在于计算顺序。我们通常把inference称为前向，因为计算顺序是从浅层到深层，一层一层地传递。而反向传播就是从深层到浅层。
 
 
先看一个二层感知机网络，如上图。每条虚线代表一个权重w，每个节点会带有一个b。在一开始会给所有权重初始化赋值，当模型没有被训练的时候也能有前向传播的输出。
 
每一次模型前向传播的输出y，我们会拿y和标准答案进行比较，就会有一个量化差异。
 
我们把这个量化差异称为损失函数，用公式就可以表示。训练模型就是为了最下化输出值和标准答案的差异，即 最小化损失函数：   （公式 1）
 
损失函数可以根据具体任务选取，比如分类任务可以使用交叉熵，超分辨率任务可以使用MSE等等。在本文中，我们不讨论损失函数的选取。咱们只需要明白，损失函数就是这一轮模型输出和标准答案(即标签)差异的量化结果。
 

2. 梯度
 
事实上，方向传播中梯度的定义和咱们以前在高数里学习的梯度是有区别的。不过思想是互通的。
 
上节提到的“最小化损失函数”，咱在这节要将具体方法。只需要把每层参数加以“斧正”，则最后的输出就会越接近标准答案(标注值)，从而损失函数值会越低。
 
那么参数该如何“斧正”呢？每次斧正多少呢？答：使用梯度来斧正，每次斧正的量就是一个负梯度。
 
用公式来表示就是（w表示网络参数）：
 
在本来的参数值上再加一个增量值（该增量可以为负），就变为更新后的参数。所以，我们只用确定了就直到该怎么优化网络了。当然，这也是反向传播算法的精髓。公式如下：
 
（公式 2）
 
用一句话来解释：直接通过偏导寻找参数对于损失函数的负梯度。
 
怎么计算偏导？
 
链式法则走一波，放心，本文是不会介绍链式法则的。
 
下图截自西瓜书：误差反向传播算法
 
 

3. 异形层反向传播
 
掌握第2节就可以明白Hinton老爷子最初的反向传播算法了。可是，随着神经网络的发展，咱们的网络层可不止多层感知机呀。当遇到卷积层，最大池化，平均池化，空洞卷积，反卷积，子像素卷积等等异形网络层的时候，怎么传播呢？
 
尤其是池化层和pixel shuffle这种并不可导的网络层。在神经网络里，我们通常可以用一个函数来建模一个网络层，如果所有的网络层都是可导函数，则我们可以通过链式法则计算任何一层的梯度完成训练。事实上，并不是所有网络层都是可导的。
 
有一篇关于卷积层和池化层的反向传播文章《卷积池化反向传播》，大家可以先看着。
 
我先留个坑，以后时间充足再来详细解析异形层的反向传播。