计算物理期末论文

一、摘要

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

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

二、原理

1、二维伊辛模型

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

我们定义总能量

$$E=-J\sum_{<ij>}{s_i s_{j}}-\mu H\sum_{i} s_{i}$$

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

$$E_{flip}=-2E$$

当 $E_{flip}<0$,反转自旋;

当$E_{flip}>0$,利用玻尔兹曼分布计算概率 $p=exp(-E_{flip}/k_B T)$ ,并生成一个随机数来决定是否反转自旋。

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

二维伊辛模型有一个相变点,临界温度 $T_c$ 满足以下方程:

$$\sinh{\frac{2J}{kT_c}}= 1$$

为模拟方便,这里我们取 $J/k=1$,则有

$$T_c=\frac{2}{\ln(1+{\sqrt 2)}}\approx 2.269$$

2、HNN

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

$$E=-\sum_{i,j}{J_{ij}s_i s_j}$$

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

$$J_{i,j}=\frac{1}{M}\sum_m s_i(m)s_j(m)$$

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

由于该模型比较简陋,因此在训练时对于相近的输入区分度不是很好,而且会有存储上限($\sim 0.13N$)。

三、模拟

1、二维伊辛模型

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
'''
@Author: Editing
@Date: 2018-12-16 15:32:49
@LastEditors: Editing
@LastEditTime: 2018-12-20 15:49:27
@Description: Final Homework
'''

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()

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

当 $T\approx T_c$ 时,系统会在两个状态之间随机变化。

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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<T<3$ 区间内,系统的平均磁矩会由 $1$ 逐渐趋近于 $0$,但是变化过程是连续平滑的,因此称之为二级相变。

2、HNN

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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 的抗干扰能力,首先随机将 $J$ 矩阵中 $80\%$ 的内容变为0,测试其性能。

1
2
3
4
5
6
7
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 抗干扰能力还是很不错的。

下面分别测试 $90\%$ 和 $95\%$ 损坏比例下的性能。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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)

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

四、总结

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

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

五、致谢