记一个神经网络中出现的混沌图样

2022 年 10 月 13 日 PaperWeekly

©作者 | 叶耿杰

单位 | 武汉大学

研究方向 | 凝聚态物理

自己写的 numpy BP 神经网络（拟合 [0,2] 上的函数 sin(2 pi x)/4+0.5 ）第一次跑通，兴奋之余胡乱修改了网络结构。看着 loss 曲线时而下降，时而振荡，突然，有一个网络（代码附在文末）啥也没学到，但给出了熟悉的图案：

赶紧把 loss 曲线的数据小心翼翼的导出来，作散点图，果不其然：

原来训练集为间隔 0.1 采样，反复试验均无分岔图案；当训练集改为间隔 0.05 采样时，出现了分岔图案。

反复运行了好几次，都能出现这种图样，说明与参数初始化无关。

抽去了一个 2 节点 /sigmoid 激活函数的层，仍然出现了类似的分岔图（代码是抽去之后的）。

没有照抄大佬的代码，以下代码说不定有错；但经过试验，取间隔 0.1 采样的数据集，很小的网络规模（1,3，3,1，全 sigmoid），还是能看出它在试图拟合一个函数的。

第一次更新：

补充一些实验结果：

@Horizony大佬和 @镇戎大佬所言甚妙：将待拟合的函数换成常函数 y=0.5，仍然出现相同的图样。

分岔的发生与采样点的密度有关；而分岔“合并”的速率与学习率的衰减有关：学习率衰减越快，分岔合并得也越快。不负责任地推测，假设那儿已经有了一张完整的分岔图，而加密采样点起“平移”作用，加快/减慢学习率衰减起伸/缩作用。

注意到上图的散点图出现了一些台阶的特征，这是不必要的设置——每 100 次迭代衰减一次学习率造成的。去除这个设置后，散点图显得更平滑（相应地可以将迭代次数减少到百分之一）。

混沌很可能是 sigmoid 激活函数导致的，而 relu 不起作用：单层，甚至单个 sigmoid 神经元足以产生分岔图样，虽然并不像 Logistic。

解析计算似乎有希望了？

目前为止，最令我惊讶的是，我本以为最无关的参数——采样间隔，反而是控制着分岔图样出现的最关键参数，实在匪夷所思。

会不会有这么一种可能：每个神经网络的 loss 曲线都是某个混沌映射的反向的图样，只是平时因为位置或伸缩不合适而看不出来？

第二次更新：

非常感谢各位大佬推荐的论文和从动力系统的角度做出的解答…… 同时深刻地意识到自己实在太菜，看懂它们可能还得几年。(´；ω；｀)

再补充一些现象层面的东西：

正如 @年轻诗人所说的，学习率增长就能得到正向的分岔图。

学习率和采样间隔似乎满足某种尺度不变性。比如下面四张图：

分别是固定学习率每次迭代增长 0.01%，而样本间距与初始学习率的组合分别为

(0.05,0.3),

(0.005,0.03),

(0.0005,0.003),

(0.00005,0.0003)

的 loss 曲线（只使用单个神经元）。虽然 3 -周期的位置有细微的差别，但是比起其他的参数组合，它们确实非常相似。

如果这样的假设是正确的，那么一方面，分岔可能在学习率很低的情况下发生，只要训练集密度够高；另一方面，训练集的元素甚至可以只有一个，只要学习率高到离谱，同样可以触发分岔。

对于后者，实验表明确实如此。由于不用迭代上百次就能求得周期点（中的一个）的位置，这可能可以用于编写快速生成分岔图的算法。鉴于其只有一个神经元、过高学习率以及不再要求网络学到训练集的特征的特点，可能已经偏离了神经网络的问题，因此，请详见： https://zhuanlan.zhihu.com/p/567118016

第三次更新：

上面那个算法并不准确，大家还是忘了它吧

按照 @年轻诗人的意见，训练集大小与学习率之间的关系是因为我在对 MSE 求导的时候没有正确地归一化（忘记除以样本数目 n），导致训练集扩大等效于学习率增大相同的倍数。

同时他敏锐地发现起作用的是最后一层，并给出了权重迭代的表达式（见评论区）。这个迭代之所以能生成（类似）Logistic 映射的分岔图，是因为它本身，在取绝对值的意义下，与 Logistic 映射的迭代式的图象非常接近，而学习率恰好占据了 Logistic 映射中参数的位置。

@Peter Griffin 的回答对我大有启发。对于一个优化问题，大的步长确实会带来混沌的可能。以一维梯度下降为例，在极值点的附近，我们可以对势函数做简谐近似，即展开到二阶项。由于常数项、一次项系数可以通过平移消去，二次项系数可以通过重新选取长度量纲而归一，我们实际上只要考虑

代入梯度下降的公式

可得

这个数列的命运随取值不同，只有单调收敛、振荡收敛、周期为 2 的振荡（），和振荡发散四种。

而在极值点附近稍远一点的地方，简谐近似失效，我们可以将势函数展到三次项，考虑

其极小值点仍为，迭代关系为

可以看到，适当选取系数，它完全可以变成 Logistic 映射。用 Mathematica 硬解方程，这是一个 8 次方程，其中两个根是不动点，两个根是周期 2 的点，4 个根是周期 3 的点。决定周期 3 有实根的条件是，根式中的

其解为大于等于有关的某个值。这就说明，当学习率足够大时，我们就得到了产生任意周期和产生混沌的充分条件——实的、连续的、在极值点附近有界的、具有周期 3 的映射。

（用这种方式算出的值似乎大于模拟出来的出现混沌时对应的值）

而值在一段时间中不变，给了充足的时间，以演化到收敛、周期或混沌的轨道上，从而形成清晰的分岔图。

由此可以猜测，在一个优化问题中，目标附近的非谐效应（几乎总是存在）、过大的学习率和过慢的学习率下降是产生分岔和混沌的诱因。

import numpy as np
import copy as cp
import matplotlib.pyplot as plt


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


def dsig(x):
    s = sig(x)
    return np.exp(-x) * s * s


def relu(x):
    if x < 0: return 0.1 * x
    return x


def drelu(x):
    if x >= 0: return 1
    return 0.1


Sig = (sig, dsig)
Relu = (np.vectorize(relu), np.vectorize(drelu))


class ConnectLayer:
    def __init__(self, inp, output, func_tup, lrate_tup):
        self.inp = inp
        self.output = output
        self.func = func_tup[0]
        self.dfunc = func_tup[1]
        self.lrate = lrate_tup[0]
        self.lrate_decay = lrate_tup[1]
        self.decay_cnt = 0
        self.w = np.random.random((output, inp))
        self.b = np.random.random((output, 1))
        self.yjs_cache = None
        self.xjs_cache = None

    def forward(self, xs):
        # n = xs.shape[1] 
        self.xjs_cache = xs
        yjs = self.w @ xs + self.b
        self.yjs_cache = yjs
        return self.func(yjs)

    def backF(self, upstream):  # f对y求导
        return self.dfunc(self.yjs_cache) * upstream

    def backPass(self, upstream):  # y对x求导
        return self.w.T @ upstream

    def refw(self, upstream):  # y对w求导
        return upstream @ self.xjs_cache.T 

    def refb(self, upstream):  # y对b求导
        # 注意：不要直接使用np.sum，否则行列不稳定
        db = np.sum(upstream, axis=1)
        return db.reshape(self.b.shape)

    def backward(self, upstream): 
        M = self.backF(upstream)
        self.w -= self.lrate * self.refw(M)
        self.b -= self.lrate * self.refb(M)

        self.decay_cnt += 1
        if self.decay_cnt == 100:
            self.lrate *= self.lrate_decay
            self.decay_cnt = 0

        return self.backPass(M)


def g(x):
    return np.sin(2 * np.pi * x) / 4 + 0.5


def genData():
    xs = np.arange(0, 2.001, 0.05)
    n = len(xs)
    ys = []
    for x in xs:
        ys.append(g(x))
    ys = np.array(ys)
    return xs.reshape(1, n), ys.reshape(1, n)


def mse(fjs, yjs): 
    m, n = yjs.shape
    delta = fjs - yjs
    return np.sum(delta * delta) / m / n


def dmse(fjs, yjs):
    return 2 * (fjs - yjs)


Mse = (mse, dmse)


class ScalarLayer:
    def __init__(self, inp, answer, loss_tup):
        self.inp = inp
        self.answer = answer
        self.loss = loss_tup[0]
        self.dloss = loss_tup[1]

    def forward(self, fyjs):
        return self.loss(fyjs, self.answer)

    def backward(self, fyjs):  # L对f求导
        return self.dloss(fyjs, self.answer)


class nn:
    def __init__(self, layer_msg, funcs_msg, loss, xs, ys, lrate):
        self.layer_msg = layer_msg
        self.layers = []
        self.xs = xs
        self.ys = ys
        for i in range(len(layers_msg) - 1):
            self.layers.append(ConnectLayer(
                inp=layers_msg[i],
                output=layers_msg[i + 1],
                func_tup=funcs_msg[i],
                lrate_tup=lrate
            ))
        self.outlet = ScalarLayer(1, ys, loss)

    def train(self, TURNS):
        ls = []
        for t in range(TURNS):
            # 单步训练
            data = cp.deepcopy(self.xs)
            for lay in self.layers:
                data = lay.forward(data)
            L = self.outlet.forward(data)
            data = self.outlet.backward(data)
            for lay in self.layers[::-1]:
                data = lay.backward(data)
            ls.append(L)
        return ls

    def test(self):
        newxs = np.hstack((self.xs - 1, self.xs + 1))
        data = cp.deepcopy(newxs)
        for lay in self.layers:
            data = lay.forward(data)
        return newxs, data


if __name__ == '__main__':
    xs, ys = genData()
    layers_msg = [1, 4, 2, 2, 1]
    funcs_msg = [Relu, Sig, Sig, Sig]
    model = nn(layers_msg, funcs_msg, Mse, xs, ys, lrate=(0.3, 0.999))
    ls = model.train(40000)
    plt.plot(ls)
    plt.show()
    print('end')