计算物理期末论文

发表于 2018-12-20 更新于 2020-11-30 分类于物理 Waline：阅读次数：

一、摘要

本文使用二维伊辛模型进行蒙特卡洛模拟，观察了不同温度下系统平均磁矩随时间的演化关系，并且观察到了温度变化时产生的二级相变现象，证明使用如此简单的模型也可以模拟部分物理学现象。

之后以伊辛模型为基础，简要实现了一种神经网络，即 Hopefield Neural Network(HNN),可以根据给定的训练集来生成关系矩阵，并可以用来还原加入噪声的信号。

二、原理

1、二维伊辛模型

伊辛模型是一个以物理学家恩斯特·伊辛为名的数学模型，用于描述物质的铁磁性。该模型中包含了可以用来描述单个原子磁矩的参数，其值只能为+1或-1，分别代表自旋向上或向下，这些磁矩通常会按照某种规则排列，形成晶格，并且在模型中会引入特定交互作用的参数，使得相邻的自旋互相影响。

我们定义总能量

其中为外部磁场。和现实情况一样，我们先生成一个随机分布，然后计算翻转所需能量的变化

当，反转自旋；

当，利用玻尔兹曼分布计算概率，并生成一个随机数来决定是否反转自旋。

这样随着循环次数的增加，系统总会趋向于一个能量最低的地方，至于是局部最低还是全局最低则取决于温度和外磁场，可以使用平均场理论进行粗略的计算。

二维伊辛模型有一个相变点，临界温度满足以下方程：

为模拟方便，这里我们取，则有

2、HNN

大脑中的神经网络可以看作是无数神经元节点相互连接，相互影响所形成的。因此我们可以利用与之类似的伊辛模型来进行模拟，只是这里的相互作用参数不再是一个定值，而是变为一个的矩阵。因此整个系统的总能量可以表示为

对于给定的输入，要使系统能量最低，相互作用参数为

一旦生成了该矩阵，我们就可以对新的输入进行计算，从而使用与 1 中类似的方法得到系统能量最低的状态，也就是经过 HNN 处理后的输出信号。

由于该模型比较简陋，因此在训练时对于相近的输入区分度不是很好，而且会有存储上限（）。

三、模拟

1、二维伊辛模型

下面的代码默认采用的二维格子进行模拟。首先展示了外部磁场为0时，在不同温度下平均磁矩的变化情况。


import numpy as np
import matplotlib.pyplot as plt

class ising(object):
    def __init__(self, lenth = 20,temperature = 1, total_time = 1000):
        self.l = lenth
        self.T = temperature
        self.t = total_time
        self.lattice = np.random.choice((1,-1), (self.l, self.l))  # 随机生成初始格点
        self.magnet = [np.sum(self.lattice) / (self.l ** 2)]

    def energy(self, i, j):
        # 使用周期边界条件，超出边界时即为0
        neighbors = self.lattice[i][(j+1) % (self.l)] + self.lattice[i][(j-1) % (self.l)] + self.lattice[(i+1) % self.l][j] + self.lattice[(i-1) % self.l][j]
        return 2 * self.lattice[i][j] * neighbors

    def run(self):
        for t in range(1, self.t):
            for i in range(self.l):
                for j in range(self.l):
                    energy = self.energy(i, j)
                    if energy <= 0:
                        self.lattice[i][j] *= -1
                    elif np.random.random() < np.exp(-energy / self.T):
                        self.lattice[i][j] *= -1
            mag = np.sum(self.lattice) / (self.l ** 2)
            self.magnet.append(mag)

    def show(self):
        plt.title('Magnetization versus Time  T = %s'% self.T)
        plt.xlabel('Time')
        plt.ylabel('Magnetization')
        plt.ylim((-1.2, 1.2))
        plt.plot(range(self.t), self.magnet)

if __name__ == '__main__':
    k = 221
    plt.figure(figsize = (12, 12), dpi = 80)
    for i in (1, 2, 2.27, 4):
        a = ising(temperature = i)
        a.run()
        plt.subplot(k)
        a.show()
        k += 1
    plt.show()

可以看到，当时，会产生自发的极化现象，也即整个系统会变为有序相。而这个变化是随机的，整个系统的自旋会随机趋向于某一个状态，这是系统表现为铁磁性。

当时，系统会在两个状态之间随机变化。

而当时，系统变为无序的状态，即磁性消失，平均磁场在零附近涨落。

这种有磁性、无磁性两相之间的转变，是一种连续相变，即二级相变。

下面我们绘制系统的平均磁矩随温度变化的图像，可以更直观的反应这种二级相变。

mag = []
for i in np.arange(0.1, 5, 0.1):
    a = ising(temperature = i, total_time = 1000, lenth = 10)
    a.run()
    mag.append(np.mean(np.abs(a.magnet)))

plt.figure(figsize = (10, 10), dpi = 80)
plt.plot(np.arange(0.1, 5, 0.1), mag, 'o', markersize = 3)
plt.grid(True)
plt.title('Magnetization versus Temperature')
plt.xlabel('Temperature')
plt.ylabel('Magnetization')
plt.ylim((0, 1))
plt.show()

可以很明显的发现在区间内，系统的平均磁矩会由逐渐趋近于，但是变化过程是连续平滑的，因此称之为二级相变。

2、HNN

这部分的模拟中使用了MNIST手写数字数据为基础，其格式为的数值矩阵，总计784个像素值。

我从数据集中提取了0-9十个数字并转化为 npy 文件，为了满足输入的要求，将其进行了二值化处理。由于模型限制，只采用了三个数字进行输入得到关系矩阵，简单测试了生成的神经网络对于信号的还原功能。

import numpy as np
import matplotlib.pyplot as plt

def train(train_data):
    J = np.zeros((784, 784))
    for s in train_data:
        J += np.outer(s, s)
    return J / 3

def recall(J, test_in, steps = 5):
    sgn = np.vectorize(lambda x: -1 if x<0 else +1)
    for i in range(steps):
        test_out = sgn(np.dot(test_in, J))
    return test_out

def destroy(J, possible):

    for i in range(784):
        for j in range(784):
            if np.random.rand() < possible:
                J[i][j] = 0
    return J

def test(J, test_data, steps = 5):
    plt.figure(figsize = (8,12), dpi = 80)
    k = 321
    for n in test_data:
        test_noise = n + np.random.binomial(2, 0.4, 784) # 随机引入噪声信号
        new = recall(J, test_noise, steps) # 经过 HNN 处理后
        plt.subplot(k)
        plt.title('Image with Noise')
        plt.imshow(test_noise.reshape(28, 28), cmap = 'gray')
        k += 1
        plt.subplot(k)
        plt.title('Image after HNN')
        plt.imshow(new.reshape(28, 28), cmap = 'gray')
        k += 1
    plt.show()

if __name__ == '__main__':
    train_data = np.load('./asset/nor_train.npy')[[0, 2, 7]] # 仅使用 0，2，7三个数字进行测试
    J = train(train_data)
    plt.figure(figsize = (8,8), dpi = 80)
    plt.title('Image of J')
    plt.imshow(J.reshape(784, 784), cmap = 'gray')
    plt.show()
    test(J, train_data)

下面验证一下 HNN 的抗干扰能力，首先随机将矩阵中的内容变为0，测试其性能。

J = train(train_data)
J = destroy(J, 0.8)
plt.figure(figsize = (8,8), dpi = 80)
plt.title('Image of J with 80% 0')
plt.imshow(J.reshape(784, 784), cmap = 'gray')
plt.show()
test(J, train_data, 30)

可以看到即使关系矩阵绝大部分信息消失，剩余部分仍然可以较好的还原输入的信号，说明 HNN 抗干扰能力还是很不错的。

下面分别测试和损坏比例下的性能。

J = train(train_data)
J = destroy(J, 0.90)
plt.figure(figsize = (8,8), dpi = 80)
plt.title('Image of J with 90% 0')
plt.imshow(J.reshape(784, 784), cmap = 'gray')
plt.show()
test(J, train_data, 20)

J = train(train_data)
J = destroy(J, 0.95)
plt.figure(figsize = (8,8), dpi = 80)
plt.title('Image of J with 95% 0')
plt.imshow(J.reshape(784, 784), cmap = 'gray')
plt.show()
test(J, train_data, 20)

可以看出其保留性能的极限基本在左右，之后的输出基本就只能靠脑补了

四、总结

本文使用蒙特卡洛方法对伊辛模型进行模拟，得到了无外磁场作用时系统二级相变的直观图像，表明伊辛模型虽然简单，但是物理内涵和应用范围都十分广泛。

之后以伊辛模型为理论基础的 HNN 也验证了这一点。但是我们可以看到这种简单的神经网络模型应用范围其实十分有限，更多的是采用其他更加复杂的神经网络模型。但就是这么一个简单的神经网络，也展现出了其抗干扰能力强、训练简单的优势，可以在此基础上引入改进增强性能，如玻尔兹曼机等。

一、摘要

二、原理

1、二维伊辛模型

2、HNN

三、模拟

1、二维伊辛模型

2、HNN

四、总结

五、致谢