碰撞约束

在碰撞分离阶段除了使用碰撞法线直接分离物体之外，还有一种广泛用于各大刚体模拟引擎中的方法——碰撞约束。使用碰撞约束可以更加优雅的分离物体。

通过之前《物理引擎中的约束》这篇文章我们了解到我们只需要写出我们需要的约束方程，然后找出Jacobian矩阵就可以通过之前那篇文章最后导出的公式求解冲量了。接下来我们将介绍解释碰撞约束的约束方程和推导Jacobian矩阵。

碰撞约束的约束方程

假设物体A与物体B发生碰撞

想一想我们需要什么效果，我们希望物体A与物体B满足分离状态，我们需要使用符号表示出这个状态。这里再次给出约束方程的含义：约束方程是用于描述一个系统中各部件之间相对运动限制的数学表达式。通过这些方程，我们可以定义哪些运动是允许的，哪些是被禁止的。即当方程中的参数满足约束方程我们就不必调整物体状态，否则我们需要调整物体状态。

所以我们可以这么定义碰撞约束的位置约束方程

$$ C = (\mathbf{P}_a-\mathbf{P}_b)\cdot\mathbf{n}\ge0 $$

这个方程表示两点在法向方向上不能穿透（距离 ≥ 0），但可以分离（距离 > 0）或者接触（距离=0）。还可以这么理解当两个物体穿透的时候$$\mathbf{P}_a-\mathbf{P}_b$$与法线的点积小于0（不满足约束方程）剩下的情况（ ≥ 0）就是满足约束的情况了。

找到约束方程后我们还需要找到Jacobian矩阵

Jacobian矩阵

还记得我们《物理引擎中的约束》中提到的我们找到位置约束方程后，我们还需要对时间求导，得到速度约束方程。

对上面得到的碰撞约束方程对时间求一阶导数得到速度约束方程

$$ \dot{C} = \mathbf{n}\cdot\frac{d}{dt}\mathbf{P}_a- \mathbf{n}\cdot\frac{d}{dt}\mathbf{P}_b $$

• $ \frac{d}{dt}\mathbf{P}_a$就是碰撞点a的速度
• $ \frac{d}{dt}\mathbf{P}_b$就是碰撞点b的速度

容易得到点a的速度为

$$ \frac{d}{dt}\mathbf{P}_a=\mathbf{v}_a+(\omega_a\times\mathbf{r}_a) $$

容易得到点b的速度为

$$ \frac{d}{dt}\mathbf{P}_b=\mathbf{v}_b+(\omega_b\times\mathbf{r}_b) $$

带入速度约束方程化简有

$$ \dot{C} = \mathbf{n}\cdot\mathbf{v}_a+\omega_a\cdot(\mathbf{r}_a\times\mathbf{n})-\mathbf{n}\cdot\mathbf{v}_b-\omega_b\cdot(\mathbf{r}_b\times\mathbf{n}) $$

我们将其写成矩阵的形式

$$ \dot{C} = \begin{bmatrix} \mathbf{n} & \mathbf{r}_a\times\mathbf{n} & -n &-\mathbf{r}_b\times\mathbf{n} \end{bmatrix}\begin{bmatrix} \mathbf{v}_a \\ \omega_a\\ \mathbf{v}_b \\ \omega_b \end{bmatrix} $$

其中Jacobian就是

$$ J=\begin{bmatrix} \mathbf{n} & \mathbf{r}_a\times\mathbf{n} & -n &-\mathbf{r}_b\times\mathbf{n} \end{bmatrix} $$

注：Jacobian是我们最后应用不同物体冲量的方向例如

a->ApplyImpulseLinear(Vec2(Jacobian[0]*lambda, Jacobian[1]*lambda)); 
a->ApplyImpulseAngular(Jacobian[2]*lambda);                   
b->ApplyImpulseLinear(Vec2(Jacobian[3]*lambda, Jacobian[4]*lambda)); b->ApplyImpulseAngular(Jacobian[5]*lambda);           

求解冲量大小

根据我们的《物理引擎中的约束》最后的得到的方程我们可以解出冲量

$$ \lambda = -(JM^{-1}J^{T})^{-1}(JV) $$

带bias系数的形式

$$ \lambda = -(JM^{-1}J^{T})^{-1}(JV+b) $$

其中V为

$$ V=\begin{bmatrix} v_{ax} \\ v_{ay}\\ \omega_a\\ v_{bx} \\ v_{by} \\ \omega_b \end{bmatrix} $$

M的逆为

$$ M^{-1} = \begin{bmatrix} \frac{1}{m_a} & & & & & \\ & \frac{1}{m_a} & & & & \\ & & \frac{1}{I_a} & & & \\ & & & \frac{1}{m_b} & & \\ & & & & \frac{1}{m_b} & \\ & & & & &\frac{1}{I_b} \end{bmatrix} $$

JV是一个数，$JM^{-1}J^{T}$也是一个数，所以求出来$\lambda$也是一个数我们的公式是正确的

注:这里可以使用高斯-赛德尔（Gauss-Seidel）方法求解线性方程组

Baumgarte 稳定性方法（Baumgarte Stabilization Method） 是一种常用于刚体动力学中处理约束系统稳定性的数值方法，特别是在使用Lagrange 乘子法求解带有约束的系统时。它的主要目标是防止数值积分过程中因误差累积导致的约束漂移（constraint drift）问题。

如果不使用 Baumgarte 稳定性方法 直接应用碰撞约束求解出来的冲量可能会导致在数值积分中由于舍入误差或逼近误差，约束无法完全维持。最后导致系统状态随时间逐渐“漂移”出满足约束的区域。

Baumgarte 稳定性方法核心思想

Baumgarte 方法通过引入一个微分项来“惩罚”偏离约束的行为，从而稳定系统。

在碰撞约束中我们一般使用

$$ \dot C = JV = 0 $$

表示一个速度约束方程

利用 Baumgarte 修正，构造一个新的方程：

$$ \dot C = JV + b = 0 $$

其中

$$ b = \frac{\beta}{\Delta t}C $$

• $\Delta t为时间步长$
• $\beta$是一个系数，可以跟据你需要的修正速度调整，大小在0到1之间，越大修正的越快

这一项我们通常叫做bias系数

引入bias系数后我们约束求解的冲量结果为

$$ \vec {\lambda}_{imp} = (JM^{-1}J^T)^{-1}-(JV_0+b) $$

优点

• 简单易实现，只需对加速度约束项做改写。
• 适用于显式或半隐式积分方法。
• 可用于任意类型的约束系统（位置、速度等）。

缺点

• 参数选择敏感。若 $\beta$ 选得不合适，可能导致系统刚性太强或振荡。
• 并不是完全消除漂移，而是“压制”它；如果时间步过大仍可能不稳定。

碰撞响应

​ 在游戏物理引擎中两个物体在碰撞检测阶段结束后一般下一阶段就是要处理碰撞了，也就是碰撞处理阶段，碰撞处理阶段可以细分成两个阶段，先是经过分离阶段，再就是碰撞响应了。

​ 碰撞响应要处理的就是如何实现类似现实世界中的弹性碰撞效果，对此我们需要更新物体碰撞后的运动状态，使得物体的运动状态与现实世界的相似。

传统冲量法

​ 还记得我们高中学习物理是怎么解决这类碰撞问题的吗？我们通常使用动量定理那一套方法，所以在这里我们还是使用这套方法来处理。

通过碰撞检测，我们可以得到以下信息：

• 碰撞点位置：$\mathbf{p}_a,\mathbf{p}_b$
• 碰撞法线：$\mathbf{n}$
• 穿透深度：$depth$

并且我们还具有物体此时物体的动力学状态。

• 质心位置：$\mathbf{p}_A,\mathbf{p}_B$
• 线速度：$\mathbf{v}_{Alinear},\mathbf{v}_{Blinear}$
• 角速度：$\mathbf{\omega}_A,\mathbf{\omega}_B$
• 质量：$\mathbf{m}_A,\mathbf{m}_B$
• 转动惯量：$\mathbf{I}_A,\mathbf{I}_B$
• 恢复系数：$e$

如果我们能计算出碰撞时的冲量，那么我们就可以通过冲量来改变速度而达成目的。

令A，B的广义速度为

$$ \mathbf v_A =\mathbf{v}_{Alinear}+\mathbf{\omega}_A\times{\mathbf{r_a}}\\ \mathbf v_B =\mathbf{v}_{Blinear}+\mathbf{\omega}_B\times{\mathbf{r_b}} $$

其中

$$ \mathbf{r_a}=\mathbf{p_a}-\mathbf{p_A}\\ \mathbf{r_b}=\mathbf{p_b}-\mathbf{p_B} $$

根据动量定理以及角动量定理有

$$ J_A \mathbf n= m_A\mathbf{v}^{'}_{Alinear}-m_A\mathbf{v}_{Alinear} \\ $$

$$ J_B\mathbf n = m_B\mathbf{v}^{'}_{Blinear}-m_B\mathbf{v}_{Blinear} \\ $$

$$ J_A\mathbf n\times{\mathbf r_a} = I_A\mathbf{\omega}^{'}_{A}-I_A\mathbf{\omega}_{A} \\ $$

$$ J_B\mathbf n\times{\mathbf r_b} = I_B\mathbf{\omega}^{'}_{B}-I_B\mathbf{\omega}_{B} $$

由牛顿第三定律可知

$$ J_A = - J_B $$

令

$$ J = J_A = -J_B $$

则碰撞后的广义速度为

$$ \mathbf v_A^{'} =\mathbf{v}^{'}_{Alinear}+\mathbf{\omega}_A^{'}\times{\mathbf{r_a}} \\ $$

$$ \mathbf v_B ^{'}=\mathbf{v}^{'}_{Blinear}+\mathbf{\omega}_B^{'}\times{\mathbf{r_b}} $$

联立动量定理、牛顿第三定律、碰撞后的广义速度可得

$$ \mathbf v_A^{'} =\mathbf{v}_{A}+\frac{J\mathbf{n}}{m_A}+(\mathbf I_A^{-1}(J\mathbf{n\times{r_a}}))\times \mathbf{r_a} \\ $$

$$ \mathbf v_B^{'} =\mathbf{v}_{B}-\frac{J\mathbf{n}}{m_B}-(\mathbf I_B^{-1}(J\mathbf{n\times{r_b}}))\times \mathbf{r_b} $$

由于

$$ e = \frac{(\mathbf v_A^{'} -\mathbf v_B^{'})\cdot\mathbf n}{(\mathbf v_B -\mathbf v_A)\cdot\mathbf n} \\ $$

$$ (\mathbf v_A^{'} -\mathbf v_B^{'})\cdot\mathbf n = -e({\mathbf v_A -\mathbf v_B})\cdot\mathbf n $$

将碰撞后的广义速度两式子相减带入恢复系数的式子中有

$$ (\frac{J\mathbf n}{m_A}+\frac{J\mathbf n}{m_B}+(\mathbf{I}_A^{-1}(J\mathbf n\times{\mathbf r_a}))\times{\mathbf r_a}+(\mathbf{I}_B^{-1}(J\mathbf n\times{\mathbf r_b}))\times{ \mathbf r_b})\cdot\mathbf n=-(1+e)(\mathbf{v}_a-\mathbf{v}_b)\cdot\mathbf n $$

由于我们只改变碰撞法线方向上的冲量所以我们在两边点乘上碰撞法线即可，然后化简得到

$$ J = \frac{-(1+e)(\mathbf{v}_A-\mathbf v_B)\cdot{\mathbf{n}}}{\frac{1}{m_A}+\frac{1}{m_B}+\mathbf{n}\cdot((\mathbf{I}_A^{-1}(\mathbf {r_a}\times {\mathbf{n}}))\times{\mathbf r_a}+(\mathbf{I}_B^{-1}(\mathbf {\mathbf r_b}\times {\mathbf{n}}))\times{\mathbf r_b})} $$

注：需要用到混合积（标量三重积）

特别的，在二维情况下

$$ J = \frac{-(1+e)(\mathbf{v}_A-\mathbf v_B)\cdot{\mathbf{n}}}{\frac{1}{m_A}+\frac{1}{m_B}+\mathbf{n}\cdot(\frac{(\mathbf {r_a}\times {\mathbf{n}})^2}{I_A}+\frac{(\mathbf {r_b}\times {\mathbf{n}})^2}{I_B})} $$

摩擦力效果

​ 要实现摩擦力的效果需要将上公式中的 $\mathbf n$ 全部替换为$\mathbf n$的正交单位向量即可

碰撞约束的响应(弹性碰撞)

​ 在碰撞约束完成后物体不会穿透其他物体，但是也不会发生弹性碰撞，所以我们在碰撞约束中的Baumgarte稳定性方法系数bias上加上 $e(\mathbf v_A-\mathbf v_B)\cdot \mathbf n$ 即可

即bias为

$$ b = \frac{\beta}{\Delta t}C+e(\mathbf v_A-\mathbf v_B)\cdot \mathbf n $$

摩擦力效果

​ 要实现摩擦力的效果需要将上约束方程中的 $\mathbf n$ 全部替换为$\mathbf n$的正交单位向量即可

引用

https://math.stackexchange.com/questions/4501575/full-derivation-of-impulse-formula-for-collision-response

折叠表达式（Fold Expression）是 C++17 引入的特性，用于简化对可变参数模板（Variadic Template）的操作。它通过将参数包（Parameter Pack）与运算符结合，避免了递归展开模板的繁琐代码。

基本语法

形式 语法 展开方式 示例

一元左折叠 (... op args) 左结合，按 ((a op b) op c) 展开
一元右折叠 (args op ...) 右结合，按 (a op (b op c)) 展开
二元左折叠 (init op ... op args) 左结合，带初始值 (0 + ... + args) → (((0 + a) + b) + c)
二元右折叠 (args op ... op init) 右结合，带初始值 (args + ... + 0) → (a + (b + (c + 0)))

示例代码

1. 求和函数

   template<typename... Args>
   auto sum(Args... args) {
    return (args + ...); //右折叠 展开为：args1 + (args2 + (...))
   }
   
   int main() {
    std::cout << sum(1, 2, 3, 4); // 输出 10
   }

2. 检查所有参数是否为真

   template<typename... Args>
   bool all_true(Args... args) {
    return (args && ...); //右折叠 展开为：args1 && args2 && ... && argsN
   }
   
   int main() {
    std::cout << std::boolalpha << all_true(true, true, false); // 输出 false
   }

3. 拼接字符串

   template<typename... Args>
   std::string concat(Args... args) {
    return (std::string{} + ... + args); //左折叠 展开为：((str1 + str2) + str3)...
   }
   
   int main() {
    std::cout << concat("Hello", ", ", "C++", "!"); // 输出 "Hello, C++!"
   }

注意事项

• 空参数包处理：空参数包需指定默认值，例如 (args + ... + 0)。
• 结合性差异：左折叠（从左侧开始计算）和右折叠（从右侧开始计算）可能因运算符特性导致结果不同。
• 运算符限制：仅支持二元运算符，且需注意运算符优先级。

在C++11中引入了可变参数模板，顾名思义参数的数量是可以改变的，我们先来看一段代码

#include <iostream>
using namespace std;

void myPrint()
{
    cout << "the last one";
}

template <typename T, typename... Args>
void myPrint(T first, Args... rest)
{
    cout << "process:" << first << endl;
    myPrint(rest...);
}

int main()
{
    myPrint(123, "hello", 3.14);
    return 0;
}

输出

process:123
process:hello
process:3.14
the last one

我们慢慢分析这段代码，这段代码先是定义了一个myPrint函数，这个函数没有参数，当参数包为空时会调用该函数。
再是定义了一个函数模板，函数模板，模板参数中有一个typename... Args，这个参数是可变参数模板的关键。
在C++17之前可变参数模板是通过递归展开的。我们可以看到myPrint函数中调用了myPrint函数，这表明该可变参数函数模板通过递归的方式展开参数，而上述提到的无参数的函数模板就是该递归的终止条件。
对于递归展开的可变参数函数模板，必须确保每次递归参数包严格减少，必须定义终止函数即无参数同名函数。

上述代码的调用过程分析。

在main函数中我们给myPrint函数传入参数包{123,"hello",3.14}

1. 第一次调用
   T=int，Args = {const char*,double}，rest... = {"hello",3.14}
   cout<<123<<endl;
   myPrint("hello",3.14);
2. 第二次调用
   T = const char*，Args = {double}，rest... ={3.14}
   cout<<"hello"<<endl;
   myPrint(3.14);
3. 第三次调用
   T = double，Args = {}，rest... ={}
   cout<<3.14<<endl;
   myPrint();

4. 第四次调用
   cout << "the last one";

   底层机制

   递归展开在编译期间生成多个函数实现：

5. myPrint(int,const char*,double)
6. myPrint(const char*,double)
7. myPrint(double)

8. myPrint()
   每个实例处理不同数量的类型。

   缺点

   参数包过大会导致编译时间加长，生成代码体积变大。

   其他方案

   在C++17以后可以使用折叠表达式简化代码，避免了现实递归。

   #include <iostream>
   using namespace std;
   
   template <typename... Args>
   void myPrint(Args... rest)
   {
    (cout <<...<<rest);
   }
   
   int main()
   {
    myPrint(123, "hello", 3.14);//输出123hello3.14
    return 0;
   }

   在下一篇文章中介绍折叠表达式

前言

在我们介绍人工神经元的时候提到了人工神经元模型是模仿生物体的神经元进行设计的模型。在那一篇文章中我们了解到生物神经元中的电信号在达到阈值时才会输出电信号给下一个神经元。人工神经元为了模仿生物神经元这一特征在人工神经元中加入了激活函数（Activation functions）

这篇文章你将会了解：

• 什么是激活函数
• 常用的激活函数
• 使用python编写带有激活函数的神经元完成非线性回归

什么是激活函数

由于我们之前的神经元只有一段线性函数y=wx+b，没有激活函数的神经网络实质上是一个线性回归模型，只能解决线性可分的问题，而对于线性不可分的任务我们之前的神经元就无能为力了，这时候我们就引入了激活函数，使得我们的函数变得不线性也就是变成非线性函数，这可以增加模型泛化能力。

例如我们不可能用一条直线去拟合这种分布的数据：

而曲线可以

所以激活函数能够帮助你的神经元完成更加复杂的工作。

并且我们人类思考问题的过程并不是完美的拟合，而是离散的分类，我们会给食物能否吃饱分成多类：能吃饱，不能吃饱，能吃几分饱(这也是我们根据经验得出的新类)，而我们并不会去使用一个食物的大小和饱腹程度的函数去判断是否能吃饱。

我们把之前线性函数的输出放入激活函数中计算后的值做出最后的输出

常见的激活函数

sigmoid

$$ sigmoid(x) = \frac{1}{1+e^{-x} } $$

sigmoid激活函数特点：

• 把(−∞,+∞)的数映射到(0,1)之间，因此它对每个神经元的输出进行了归一化；
• Sigmoid 函数非常合适用于将预测概率作为输出的模型。因为概率的取值范围是 0 到 1。
• 函数处处可导。
• 函数很快逼近y=0和y=1，导数无限小，这会导致权重更新的非常慢
• 一般用于二分类的输出层

sigmoid的导函数

$$ \frac{\partial sigmoid(x)}{\partial x} = sigmoid(x)(1-sigmoid(x)) $$

Tanh

$$ Tanh(x) = \frac{e^{x}-e^{-x}}{e^{x}+e^{-x} } $$

Tanh激活函数特点：

• 函数是将取值为 (−∞,+∞) 的数映射到(−1,1) 之间；
• 在0附近的导数相比sigmoid更大，收敛速度比sigmoid快
• 函数处处可导。
• 函数很快逼近y=0和y=1，导数无限小，这会导致权重更新的非常慢

Tanh的导函数

$$ \frac{\partial Tanh(x)}{\partial x} = 1-Tanh(x)^{2} $$

ReLU

$$ ReLU(x) = \begin{cases} 0 & x\leqslant 0 \\ x & x>0 \end{cases} $$

ReLU激活函数特点：

• 计算简单高效，相比sigmoid、tanh没有指数运算
• 相比sigmoid、tanh更符合生物学神经激活机制
• 收敛速度较快，大约是 sigmoid、tanh 的 6 倍
• 在x<0时导数为0，这会导致权重参数无法得到更新，这被称为神经元死亡，或者梯度消失

ReLU的导函数

$$ \frac{\partial ReLU(x)}{\partial x} = \begin{cases} 0 & x\leqslant 0 \\ 1 & x>0 \end{cases} $$

使用python编写带有激活函数的神经元完成非线性回归

只使用Sigmoid完成

# 导入numpy库
import numpy as np
from matplotlib import pyplot as plt


def sigmoid(x):
    return 1/(1+np.exp(-x))


def sigmoid_deriv(a):
    return a*(1-a)


def getdata(count):
    """获取指定数量的数据"""
    xs = np.random.rand(count)
    xs = np.sort(xs)
    ys = np.zeros(count)
    for i in range(count):
        x = xs[i]
        yi = 0.7*x+(0.5-np.random.rand())/50+0.5
        if yi > 0.8:
            ys[i] = 1
    return xs, ys


data = getdata(100)
xs = data[0]
ys = data[1]

w = np.random.rand(1)
b = np.random.rand(1)
z = w*xs+b
a = sigmoid(z)

plt.xlim(-0.1, 1)
plt.ylim(-0.1, 3)
plt.scatter(xs, ys)
plt.plot(xs, a)
plt.show()

for j in range(1000):
    for i in range(len(xs)):
        x = xs[i]
        y = ys[i]
        # 前向传播
        z = w*x+b
        a = sigmoid(z)#加入激活函数
        e = (y - a)**2

        alpha = 0.05

        # 反向传播
        deda = -2*(y-a)
        dadz = sigmoid_deriv(a)#加入激活函数后本质还是复合函数求导
        dzdw = x
        dzdb = 1

        dedw = deda*dadz*dzdw
        w = w - alpha*dedw
        dedb = deda*dadz*dzdb
        b = b - alpha*dedb

z = w*xs+b
a = sigmoid(z)
plt.xlim(0, 1)
plt.ylim(0, 3)
plt.scatter(xs, ys)
plt.plot(xs, a)
plt.show()

在上一篇文章中我们使用了化简后的Rosenblatt感知器进行回归预测，在这篇文章中我们将了解损失函数、梯度下降算法以及梯度下降算法的的实现，且加入偏置项系数b。

损失函数

损失函数(Loss Function)：定义在单个样本上，算的是一个样本的误差.

代价函数(Cost Function)：定义在整个训练集上，算的是所有样本的误差，也就是损失函数的平均。

在机器学习中几乎所有的模型都会有损失函数，损失函数描述的是：

预测值和样本数据的实际目标值之间的误差与预测函数的系数之间的关系

均方误差

均方误差代价函数其表达式如下：

$$ e(w,b) = \frac{1}{n}\sum_{i=0}^{n}(y_{i}-(wx_{i}+b))^2 $$

其中y_i表示样本数据的实际目标值表示wx_i+b预测值根据样本数据x_i计算出的预测值。

我们可以发现这个图像是一个开口向上的二元二次函数图像，既然是一个开口向上的二元二次函数，那么我们是不是可以找到一个w和b使得整体误差最小(全局最优解)呢？是的在机器学习中正是这样。

寻找最低点的方法有

• 正规方程

  $$ w=\frac{\sum_{i=0}^{n}(x_{i}y_{i})}{\sum_{i=0}^{n}x_{i}^2} $$

  推导过程是使用了二元一次函数的顶点坐标公式

• 梯度下降

梯度下降

梯度下降(Gradient descent)的可以类比为一个下山的过程，如图所示函数看似为一座山，山中全是雾，有一个人在山上任意一个位置想要快速找到水，由于山上全是雾，人只能看见脚下的路，那么人可以判断脚下土地的坡度判断自己是否为下坡路，每走一步调整自己的方向，每次向坡度最大的方向走，就走到最低点的位置喝到水。

数学中实现

在数学中具有一个概念叫做梯度，梯度描述的就是一个面上某个点的瞬时变化率。和斜率类似，其实梯度就是斜率更广的说法。

既然描述的是面上某个点的瞬时变化率，那么我们就可以对他进行求导。我们要得到w最小同时b最小的一个点，我们只需要分别对w和b求偏导，当某个点的梯度大于零，对应的偏导数大于零我们只需要让对应的值减去alpha(alpha>0)倍的偏导数，偏导数为正，对应值往小的调整；当某个点的梯度小于零，对应的偏导数小于零我们只需要让对应的值减去alpha(alpha>0)倍的偏导数，偏导数为负，对应值往大的调整。

• 为什么要在偏导数上乘以一个alpha

  alpha叫做学习率，在上一篇文章中遇到过。在调整w和b的过程中如果没有学习率alpha它将可能无法或者很慢的回到最低点，乘以一个alpha可以加快或者调慢速度

$$ \xi =\xi -\alpha \times \frac{\partial f(\xi )}{\partial \xi } $$

随机梯度下降算法

随机梯度下降算法(SGD)取出其中一个样本的数据依次进行梯度下降

反向传播

我们已经知道了梯度下降算法，但是我们如何获得求得损失函数中w和b的偏导呢。其实高中数学足以

观察公式：

$$ e(w,b) = (y_{i}-(wx_{i}+b))^2 $$

使用复合函数的求导法则：

$$ \frac{\partial e}{\partial w}=2(y_{i}-b-wx_{i})\cdot (-x_{i}) $$

$$ \frac{\partial e}{\partial b}=2(y_{i}-b-wx_{i})\cdot (-1) $$

调整w和b

$$ w =w -\alpha \times \frac{\partial e}{\partial w } $$

$$ b =b -\alpha \times \frac{\partial e}{\partial b } $$

代码实现

将使用SGD实现

首先先引入需要用到的库

import numpy
from matplotlib import pyplot as plt

引入numpy数学库以方便我们的数学计算

引入matplotlib中的pyplot模块进行画图操作

生成数据以便进行训练

def getdata(count):
    """获取指定数量的数据"""
    x = np.random.rand(count)  # 生成区间[0,1)的随机数
    x = np.sort(x)  # 进行从小到大排序
    y = [-1.5*xx+np.random.rand()/3+1 for xx in x]  # 使用列表推导式生成对应的值
    return x, y

创建一个名为getdata的函数随机生成有共同特征的数据，有一个参数count来接收获取数据的数量

获取100个数据

data = getdata(100)
xs = data[0]
ys = data[1]

编写预测函数

w = np.random.rand(1)
b = np.random.rand(1)
z = w*xs+b

绘制图像

plt.xlim(0, 1)
plt.ylim(0, 3)
plt.scatter(xs, ys)
plt.plot(xs, z)
plt.show()

随机梯度下降算法和反向传播

for j in range(500):
    for i in range(len(xs)):
        x = xs[i]
        y = ys[i]

        z = w*x+b
        e = (y - z)**2#损失函数
        alpha = 0.05
        
        #反向传播和随机梯度下降
        dedz = -2*(y-z)
        dzdw = x
        dzdb = 1
        dedw = dedz*dzdw
        w = w - alpha*dedw
        dedb = dedz*dzdb
        b = b - alpha*dedb
        
        """
        你也可以这么写
        dedw = -2*(y-z)*x
        w = w - alpha*dedw
        dedb = -2*(y-z)*1
        b = b - alpha*dedb
        """

再次绘制拟合后的图像

z = w*xs+b
plt.xlim(0, 1)
plt.ylim(0, 3)
plt.scatter(xs, ys)
plt.plot(xs, z)
plt.show()

可以看到已经完美拟合

完整代码

# 导入numpy库
import numpy as np
from matplotlib import pyplot as plt


def getdata(count):
    """获取指定数量的数据"""
    x = np.random.rand(count)  # 生成区间[0,1)的随机数
    x = np.sort(x)  # 进行从小到大排序
    y = [-1.5*xx+np.random.rand()/3+1 for xx in x]  # 使用列表推导式生成对应的值
    return x, y


data = getdata(100)
xs = data[0]
ys = data[1]

w = np.random.rand(1)
b = np.random.rand(1)
z = w*xs+b

plt.xlim(0, 1)
plt.ylim(0, 3)
plt.scatter(xs, ys)
plt.plot(xs, z)
plt.show()

for j in range(500):
    for i in range(len(xs)):
        x = xs[i]
        y = ys[i]

        z = w*x+b
        e = (y - z)**2
        alpha = 0.05

        dedz = -2*(y-z)
        dzdw = x
        dzdb = 1

        dedw = dedz*dzdw
        w = w - alpha*dedw
        dedb = dedz*dzdb
        b = b - alpha*dedb
        """
        你也可以这么写
        dedw = -2*(y-z)*x
        w = w - alpha*dedw
        dedb = -2*(y-z)*1
        b = b - alpha*dedb
        """

z = w*xs+b
plt.xlim(0, 1)
plt.ylim(0, 3)
plt.scatter(xs, ys)
plt.plot(xs, z)
plt.show()

Rosenblatt感知器，感知器也叫预测机

那么感知器有什么用呢？感知器可以拟合任何的线性函数，任何线性分类或线性回归问题都可以用感知器来解决。

使用python实现Rosenblatt感知器

pass：此次python实验的预测函数不含偏置项系数，所以函数图像过原点，只能拟合因果是成正比的数据或只能完成可被过原点的直线分类的数据

算法图

线性回归问题

1. 首先先引入需要用到的库

   import numpy
   from matplotlib import pypolt 

   引入numpy数学库以方便我们的数学计算

   引入matplotlib中的pyplot模块进行画图操作

   numpy使用手册：https://www.numpy.org.cn/

   matplotlib使用手册：https://www.matplotlib.org.cn/

2. 获取数据以便进行训练

   def getdata(count):
       """获取指定数量的数据"""
       x = np.random.rand(count)#生成区间[0,1)的随机数
       x = np.sort(x)#进行从小到大排序
       y = [1.5*xx+np.random.rand()/3 for xx in x]#使用列表推导式生成对应的值
       return x,y

   创建一个名为getdata的函数随机生成有共同特征的数据，有一个参数count来接收获取数据的数量

3. 获取100个数据

   xs,ys = getdata(100)#获取100个数据

4. 编写预测函数

   w = float("%.2f"%np.random.uniform(0,1))#随机生成权重参数
   alpha = 0.05#设置学习率
   y_predict =  w*xs#简单的预测函数

5. 绘制图像以便我们观察

   plt.scatter(xs,ys)#绘制散点图
   plt.plot(xs,y_predict)#绘制预测函数的图像
   plt.show()#显示图像

   

6. 算法的实现

   for i in range(100):#对每个数据进行一次权重调整
           x = xs[i]
           y = ys[i]
           #取出一个数据
           y_predict =  w*x#计算相应数据的预测值
           e = y - y_predict#计算相应数据的预测值与实际值的误差（使用差值评估误差）
           w = w + alpha*e*xs#根据误差调整权重参数
           y_predict = w * xs#重新配置预测函数

7. 多次训练

   why：因为学习率的存在，我们的预测函数只进行了一次在总体上的调整，无法完美的拟合。但是如果学习率设置太大会造成预测曲线不能收敛

   for _ in range(100):#进行100次训练
       for i in range(100):#对每个数据进行一次权重调整
           x = xs[i]
           y = ys[i]
           #取出一个数据
           y_predict =  w*x#计算相应数据的预测值
           e = y - y_predict#计算相应数据的预测值与实际值的误差（使用差值评估误差）
           w = w + alpha*e*xs#根据误差调整权重参数
           y_predict = w * xs#重新配置预测函数

8. 再次画出预测曲线

   plt.plot(xs,y_predict)#绘制预测函数的图像
   plt.show()#显示图像

   

   可以看到完美的拟合了数据

线性回归问题代码

import numpy as np#引入numpy数学计算库
from matplotlib import pyplot as plt#引入matplotlib画图库

def getdata(count):
    """获取指定数量的数据"""
    x = np.random.rand(count)
    x = np.sort(x)
    y = [1.5*xx+np.random.rand()/3 for xx in x]
    return x,y

xs,ys = getdata(100)#获取100个数据

w = float("%.2f"%np.random.uniform(0,1))#随机生成权重参数
alpha = 0.05#设置学习率
y_predict =  w*xs#简单的预测函数

plt.scatter(xs,ys)#绘制散点图
plt.plot(xs,y_predict)#绘制预测函数的图像
plt.show()

for _ in range(100):#进行100次训练
    for i in range(100):#对每个数据进行一次权重调整
        x = xs[i]
        y = ys[i]
        #取出一个数据
        y_predict =  w*x#计算相应数据的预测值
        e = y - y_predict#计算相应数据的预测值与实际值的误差（使用差值评估误差）
        w = w + alpha*e*xs#根据误差调整权重参数
        y_predict = w * xs#重新配置预测函数
    #画图
    plt.clf()
    plt.xlim(0,1.1)
    plt.ylim(0,1.85)
    plt.scatter(xs,ys)
    plt.plot(xs,y_predict)
    plt.pause(0.01)
   
plt.show()