简介

这篇文章主要描述基于得分匹配（Score matching model）的想法，以及之后主要的修改思路。这种思路是生成模型的一种，与GAN、normal-flow等模型具备同样的功能。

本篇文章大量借鉴棒棒生博客，推荐阅读原文博客。本文章在其基础上加入一些作者本人的思考，并且统一符号，增加阅读流畅性。

逆Ising问题及其求解方法

引言

逆Ising问题涉及从观测数据中推断Ising模型的参数。Ising模型是一种简单的数学模型，用于描述磁性材料中的自旋系统。解决逆Ising问题的方法包括极大似然估计（MLE）、蒙特卡罗方法（MCMC）、平均场近似、伪似然估计（PLE）以及对比散度（CD）。

极大似然估计（MLE）

步骤

1. 定义Hamiltonian: H(s) = −∑_i < jJ_ijs_is_j

2. 定义配分函数: Z = ∑_se^−H(s)

3. 定义似然函数: $$ L(J) = \prod_{k=1}^{N} P(\mathbf{s}^{(k)}|J) = \prod_{k=1}^{N} \frac{e^{-H(\mathbf{s}^{(k)})}}{Z} $$

4. 取对数似然函数: $$ \log L(J) = -\sum_{k=1}^{N} H(\mathbf{s}^{(k)}) - N \log Z $$

5. 优化对数似然函数: 使用数值优化方法（如梯度下降）最大化对数似然函数。

示例代码

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

import numpy as np
from scipy.optimize import minimize

# 观测数据
data = np.array([[1, 1], [1, -1], [-1, 1], [-1, -1]])

# 定义Hamiltonian
def hamiltonian(J, s):
    return -J * s[0] * s[1]

# 定义配分函数
def partition_function(J):
    return 2 * (np.exp(J) + np.exp(-J))

# 定义对数似然函数
def log_likelihood(J, data):
    logL = 0
    for s in data:
        logL += -hamiltonian(J, s)
    logL -= len(data) * np.log(partition_function(J))
    return -logL  # 最小化负对数似然

# 初始猜测
initial_J = 0.1

# 优化参数
result = minimize(log_likelihood, initial_J, args=(data,))
optimal_J = result.x[0]

print(f'MLE Estimated J: {optimal_J}')

蒙特卡罗方法（MCMC）

步骤

1. 初始化模型参数。
2. 生成样本： 使用MCMC方法（如Metropolis-Hastings或Gibbs采样）生成自旋配置样本。
3. 计算统计量： 基于生成的样本，计算系统的期望统计量。
4. 更新参数： 使用梯度方法或其他优化技术调整模型参数。
5. 重复步骤2-4，直到参数收敛。

示例代码

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

import numpy as np

# 观测数据
data = np.array([[1, 1], [1, -1], [-1, 1], [-1, -1]])

# 初始化参数
J = 0.1

# 吉布斯采样函数
def gibbs_sample(J, num_samples, burn_in):
    samples = []
    s = np.random.choice([-1, 1], size=2)
    for _ in range(num_samples + burn_in):
        for i in range(2):
            prob = 1 / (1 + np.exp(-2 * J * s[1 - i]))
            s[i] = 1 if np.random.rand() < prob else -1
        if _ >= burn_in:
            samples.append(s.copy())
    return np.array(samples)

# 计算期望统计量
def calculate_statistics(samples):
    correlations = np.mean(samples[:, 0] * samples[:, 1])
    return correlations

# 更新参数
def update_parameters(J, samples, data):
    data_corr = np.mean(data[:, 0] * data[:, 1])
    sample_corr = calculate_statistics(samples)
    learning_rate = 0.1
    J += learning_rate * (data_corr - sample_corr)
    return J

# 蒙特卡罗迭代
num_iterations = 100
num_samples = 1000
burn_in = 100

for _ in range(num_iterations):
    samples = gibbs_sample(J, num_samples, burn_in)
    J = update_parameters(J, samples, data)
    print(f'Iteration {_}: J = {J}')

print(f'Estimated J: {J}')

平均场近似（Mean Field Approximation）

步骤

1. 定义Hamiltonian。
2. 引入平均场近似： 每个自旋的影响被近似为一个平均场。
3. 计算自旋期望值： 根据平均场计算自旋的期望值。
4. 迭代求解： 通过迭代更新每个自旋的期望值和相互作用参数。
5. 更新相互作用参数： 调整相互作用参数使得模型生成的期望值与观测数据的期望值匹配。

示例代码

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

import numpy as np

# 观测数据
data = np.array([[1, 1], [1, -1], [-1, 1], [-1, -1]])

# 初始化参数
J = 0.1

# 计算观测数据的期望值
mean_s1 = np.mean(data[:, 0])
mean_s2 = np.mean(data[:, 1])

# 定义迭代函数
def mean_field_iteration(J, mean_s1, mean_s2, max_iter=100, tol=1e-5):
    for _ in range(max_iter):
        h1 = J * mean_s2
        h2 = J * mean_s1
        new_mean_s1 = np.tanh(h1)
        new_mean_s2 = np.tanh(h2)
        if np.abs(new_mean_s1 - mean_s1) < tol and np.abs(new_mean_s2 - mean_s2) < tol:
            break
        mean_s1, mean_s2 = new_mean_s1, new_mean_s2
    return mean_s1, mean_s2

# 迭代求解期望值
mean_s1, mean_s2 = mean_field_iteration(J, mean_s1, mean_s2)

# 更新相互作用参数
def update_parameters(J, mean_s1, mean_s2, data):
    observed_corr = np.mean(data[:, 0] * data[:, 1])
    model_corr = mean_s1 * mean_s2
    learning_rate = 0.1
    J += learning_rate * (observed_corr - model_corr)
    return J

# 更新参数
J = update_parameters(J, mean_s1, mean_s2, data)

print(f'Estimated J: {J}')

伪似然估计（PLE）

步骤

1. 定义条件概率： $$\begin{align} P(s_i | \mathbf{s}_{\setminus i}) = \frac{e^{s_i h_i}}{e^{h_i} + e^{-h_i}} \end{align}$$
2. 定义伪似然函数： $$\begin{align} PL(J) = \prod_{i=1}^n \prod_{k=1}^N P(s_i^{(k)} | \mathbf{s}_{\setminus i}^{(k)}) \end{align}$$
3. 取对数伪似然函数： $$\begin{align} \log PL(J) = \sum_{i=1}^n \sum_{k=1}^N \log P(s_i^{(k)} | \mathbf{s}_{\setminus i}^{(k)}) \end{align}$$
4. 优化对数伪似然函数： 使用数值优化方法（如梯度下降）最大化对数伪似然函数。

示例代码

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

import numpy as np
from scipy.optimize import minimize

# 观测数据
data = np.array([[1, 1], [1, -1], [-1, 1], [-1, -1]])

# 定义条件概率
def conditional_probability(s_i, s_j, J):
    return 1 / (1 + np.exp(-2 * J * s_i * s_j))

# 定义对数伪似然函数
def log_pseudolikelihood(J, data):
    logPL = 0
    for s in data:
        logPL += np.log(conditional_probability(s[0], s[1], J)) + np.log(conditional_probability(s[1], s[0], J))
    return -logPL  # 最小化负对数伪似然

# 初始猜测
initial_J = 0.1

# 优化参数
result = minimize(log_pseudolikelihood, initial_J, args=(data,))
optimal_J = result.x[0]

print(f'Optimal J: {optimal_J}')

对比散度（CD）

步骤

1. 定义能量函数： $$\begin{align} E(\mathbf{s}) = -\sum_{i<j} J_{ij} s_i s_j \end{align}$$
2. 初始化模型参数。
3. 正向采样（Positive Phase）： 从观测数据中采样，计算数据分布下的期望值。
4. 负向采样（Negative Phase）： 使用Gibbs采样从模型分布中生成样本，计算模型分布下的期望值。
5. 更新参数： 根据正向采样和负向采样的期望值差异，更新相互作用参数。
6. 迭代，直到参数收敛。

示例代码

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

import numpy as np

# 观测数据
data = np.array([[1, 1], [1, -1], [-1, 1], [-1, -1]])

# 初始化参数
J = 0.1
learning_rate = 0.1
num_iterations = 100
num_samples = 1000
burn_in = 100

# 计算观测数据的期望值
mean_s1s2_data = np.mean(data[:, 0] * data[:, 1])

# 吉布斯采样函数
def gibbs_sample(J, num_samples, burn_in):
    samples = []
    s = np.random.choice([-1, 1], size=2)  # 初始自旋配置
    for _ in range(num_samples + burn_in):
        for i in range(2):
            prob = 1 / (1 + np.exp(-2 * J * s[1 - i]))
            s[i] = 1 if np.random.rand() < prob else -1
        if _ >= burn_in:
            samples.append(s.copy())
    return np.array(samples)

# 计算模型分布的期望值
def calculate_model_expectation(J, num_samples, burn_in):
    samples = gibbs_sample(J, num_samples, burn_in)
    mean_s1s2_model = np.mean(samples[:, 0] * samples[:, 1])
    return mean_s1s2_model

# 对比散度更新
for _ in range(num_iterations):
    mean_s1s2_model = calculate_model_expectation(J, num_samples, burn_in)
    J += learning_rate * (mean_s1s2_data - mean_s1s2_model)
    print(f'Iteration {_}: J = {J}')

print(f'CD Estimated J: {J}')

总结

逆Ising问题涉及从观测数据中推断Ising模型的相互作用参数。本文介绍了五种求解方法：极大似然估计（MLE）、蒙特卡罗方法（MCMC）、平均场近似、伪似然估计（PLE）以及对比散度（CD）。每种方法都有其优缺点，选择合适的方法取决于具体问题的规模、数据特性和计算资源。

波尔兹曼机是物理模型，同时也是早期的机器学习模型。

Reference: * 神经网络的统计力学 – 受限制玻尔兹曼机的统计力学

利用复本方法计算 Sherrington-Kirkpatrick（SK） 模型。该章节主要介绍 Parisi 解。

Reference: * Statistical Physics of Spin Glasses and Information Processing. Nishimori * Replica calculations for the SK model笔误有点多 * RS and RSB solutions for SK model with spin-S

Link: * Hopfield Model * Replica * Sherrington-Kirkpatrick Model

利用复本方法计算 Sherrington-Kirkpatrick（SK） 模型。目的是掌握和熟悉复本方法的使用。

Reference: * Statistical Physics of Spin Glasses and Information Processing. Nishimori * Replica calculations for the SK model笔误有点多

Link: * Hopfield Model * Replica * The Parisi solution of Sherrington-Kirkpatrick Model

Phase Diagram

根据“Hopfield 自由能”中的讨论，可知：

$$\begin{align} & q=\int D z \tanh ^2 \beta(\sqrt{\alpha r} z+m) \label{860} \\ & m=\int D z\langle\xi \tanh \beta(\sqrt{\alpha r} z+m \xi)\rangle=\int D z \tanh \beta(\sqrt{\alpha r} z+m) \label{861} \\ & r=\frac{q}{(1-\beta+\beta q)^2} \label{862}\\ \end{align}$$

通过数值求解 $\eqref{860},\eqref{861},\eqref{862}$ 可以得到Hopfield的相图：

在高温的情况下，热噪音阻碍了模式恢复的过程，因此m = 0, q = 0, r = 0；当温度下降，参数开始在临界线T_g不稳定，可以通过$\eqref{860}$分析得到；随着温度继续下降，Spin Glass相变成亚稳定态，相变线T_M，恢复的模式是局域稳定的，这个相变是一阶相变。当记忆率下降，转变为稳定的态，对应相变线为T_c。

在温度T = 0的时候，每个自旋的平均熵值为$S=-\left.\frac{\partial f}{\partial T}\right|_{T \rightarrow 0}=-\frac{1}{2} \alpha[\ln (1-C)+C /(1-C)]$，其中C = β(1 − q)，可知自旋平均值是负数，这是非物理的。

Hopfield Model with Arbitrary Hebbian Length

Hebbian learning 是一种学习算法和理论，用于解释神经网络中的突触可塑性。其核心思想可以用一句话概括，即“同步触发的神经元会连线在一起”（“Cells that fire together, wire together”）。这个概念最早由加拿大心理学家Donald Hebb在1949年提出，故称为Hebbian学习。

具体来说，Hebbian学习的基本原则如下：

1. 联结权重的更新：如果一个神经元A经常且重复地激活神经元B，那么神经元A和神经元B之间的突触连接会变得更强。这意味着两个神经元之间的联结权重会增加。

2. 时间一致性：为了使联结权重增加，神经元A的发火必须在神经元B的发火之前或同时发生。这种时间上的一致性被认为是形成记忆和学习的基础。

3. 局部性：Hebbian学习是一种局部学习规则，即每个突触的权重更新只依赖于连接的两个神经元的活动，而不依赖于整个网络的状态。

在数学上，Hebbian学习规则可以表示为：

$$\begin{align} \Delta w_{ij} = \eta x_i y_j \end{align}$$

其中，Δw_ij 表示神经元 i 和神经元 j 之间的突触权重变化，η 是学习率，x_i 和 y_j 分别是神经元 i 和神经元 j 的激活状态。

Hebbian学习在神经科学和人工神经网络领域都有重要影响，特别是在理解神经网络如何通过经验和环境进行学习和调整方面。

Hebbian strength 指的是神经元之间突触连接的强度，定义如下的神经元权重：

$$\begin{align} J_{i j}=\frac{1}{N} \sum_{\mu=1}^P\left[c \xi_i^\mu \xi_j^\mu+\gamma \sum_{r=1}^d\left(\xi_i^\mu \xi_j^{\mu+r}+\xi_i^{\mu+r} \xi_j^\mu\right)\right] \end{align}$$

其中c是标准Hebbian strength；γ是不同记忆模式之间的强度；d是模型的Hebbian 长度，之前讨论的是d = 1的情形，只记忆其中一个模式，d = 0是标准的Hopfield Model。ξ_i^μ是二值变量，例如$p\left(\xi_i^\mu= \pm 1\right)=\frac{1}{2} \delta\left(\xi_i^\mu+1\right)+\frac{1}{2} \delta\left(\xi_i^\mu-1\right)$。讨论的是P与N的极限情形，$\alpha=\frac{P}{N}$是memory load。

Computation of the Disorder-Averaged Free Energy

将J重新写为：

$$\begin{align} \mathbf{J}&=\frac{1}{N} \xi^{\mathrm{T}} \mathbf{X} \xi \\ X_{\mu \eta} & =c \delta_{\mu \eta}+\gamma \sum_{r=1}^d\left(\delta_{\mu,(\eta+r) \bmod P}+\delta_{\mu,(\eta-r) \bmod P}\right) \\ & =(c-\gamma) \delta_{\mu \eta}+\gamma \sum_{r=-d}^d \delta_{\mu,(\eta+r) \bmod P} \end{align}$$

X的第m个本征值为，解的思路： $$\begin{align} \lambda_m&=\sum_{k=0}^{P-1} X_{1(k+1)} e^{-2 \pi i m k / P} \\ & =\sum_{k=0}^{P-1} X_{1(k+1)} \cos \left(2 \pi \frac{m k}{P}\right) \\ & =\sum_{k=0}^{P-1}\left[c \delta_{0 k}+\gamma \sum_{r=1}^d\left(\delta_{0,(k+r) \bmod P}+\delta_{0,(k-r) \bmod P)}\right] \cos \left(2 \pi \frac{m k}{P}\right)\right. \\ & =c+\gamma \sum_{r=1}^d\left[\cos \left(-2 \pi \frac{m r}{P}\right)+\cos \left(2 \pi \frac{m r}{P}\right)\right] \\ & =c+2 \gamma \sum_{r=1}^d \cos \left(2 \pi \frac{m r}{P}\right) \end{align}$$ 其中m = 0, 1, …, P − 1。

Hamiltonian和配分函数为： $$\begin{align} \mathcal{H}(\mathbf{s})&=-\frac{1}{2} \sum_{i \neq j} J_{i j} s_i s_j \\ Z&=\operatorname{Tr} \exp \left[\frac{\beta}{2 N} \mathbf{s}^{\mathrm{T}} \boldsymbol{\xi}^{\mathrm{T}} \mathbf{X} \boldsymbol{\xi} \mathbf{s}\right] \end{align}$$

其中 Tr 表示对分立 s 的求和。使用复本技巧： $$\begin{align} \langle\ln Z\rangle=\lim _{n \rightarrow 0} \frac{\ln \left\langle Z^n\right\rangle}{n}, \end{align}$$ 其中⟨⋅⟩表示对ξ的求和。因此有： $$\begin{align} Z^n=\operatorname{Tr} \exp \left[\frac{\beta}{2 N} \sum_{a=1}^n\left(\mathbf{s}^a\right)^{\mathrm{T}} \boldsymbol{\xi}^{\mathrm{T}} \mathbf{X} \boldsymbol{\xi} \mathbf{s}^a\right] . \end{align}$$

考虑S个凝聚（foreground）模式P − S个非凝聚（background）模式。这里这两个模式的类别，表示什么含义？每次恢复的模式不应该只有一个么？根据以上定义，可以将X分为： $$\begin{align} \mathbf{X}=\left[\begin{array}{ll} \mathbf{X}_{F F} & \mathbf{X}_{F B} \\ \mathbf{X}_{B F} & \mathbf{X}_{B B} \end{array}\right] \end{align}$$ 其中 X_FF ∈ ℝ^S × S, X_BF^T = X_FB ∈ ℝ^{S × (P − S)} and X_BB ∈ ℝ^{(P − S) × (P − S)}。

将X_BB 对角化 X_BB^μν = ∑_σλ_ση_μ^ση_ν^σ，其中 λ_σ 与 η_μ^σ 是本征值与本征态。

使用Hubbard-Stratonovich transformation：$\exp \left[\frac{1}{2} b^2\right]=$ ∫Dxexp [±bx], 其中 $D x=\frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{x^2}{2}\right) d x$，因此有： $$ \begin{align} Z^n=&\operatorname{Tr} \exp \left[ \frac{\beta}{2 N} \sum_{a, i, j, \mu \in B, v \in B} s_i^a \xi_i^\mu X_{\mu \nu} \xi_j^v s_j^a+\frac{\beta}{N} \sum_{a, i, j, \mu \in B, v \in F} s_i^a \xi_i^\mu X_{\mu \nu} \xi_j^v s_j^a +\frac{\beta}{2 N} \sum_{a, i, j, \mu \in F, v \in F} s_i^a \xi_i^\mu X_{\mu \nu} \xi_j^v s_j^a\right] \\ =&\operatorname{Tr} \exp {\left[\frac{\beta}{2 N} \sum_{a, \sigma} \lambda_\sigma\left(\sum_{i, \mu \in B} s_i^a \xi_i^\mu \eta_\mu^\sigma\right)^2+\frac{\beta}{N} \sum_{a, i, j, \mu \in B, v \in F} s_i^a \xi_i^\mu X_{\mu \nu} \xi_j^v s_j^a +\frac{\beta}{2 N} \sum_{a, i, j, \mu \in F, v \in F} s_i^a \xi_i^\mu X_{\mu \nu} \xi_j^\nu s_j^a\right]} \\ = & \operatorname{Tr} \prod_{a, \sigma} \int D x_\sigma^a \exp \left[\sum_{i, \mu \in B} \frac{\xi_i^\mu}{\sqrt{N}}\left(\sum_{a, \sigma} s_i^a \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\frac{\beta}{\sqrt{N}} \sum_{a, j, v \in F} s_i^a X_{\mu \nu} \xi_j^v s_j^a\right)+\frac{\beta}{2 N} \sum_{a, i, j, \mu \in F, v \in F} s_i^a \xi_i^\mu X_{\mu \nu} \xi_j^v s_j^a\right] \end{align} $$

定义： $$\begin{align} \Phi_B=&\exp \left[\sum_{i, \mu \in B} \frac{\xi_i^\mu}{\sqrt{N}}\left(\sum_{a, \sigma} s_i^a \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\frac{\beta}{\sqrt{N}} \sum_{a, j, v \in F} s_i^a X_{\mu \nu} \xi_j^v s_j^a\right)\right] \\ \Phi_F=&\exp \left[\frac{\beta}{2 N} \sum_{a, i, j, \mu \in F, v \in F} s_i^a \xi_i^\mu X_{\mu \nu} \xi_j^v s_j^a\right] \end{align}$$

计算无序平均{ξ_i^μ}，可以有： $$\begin{align} \left\langle Z^n\right\rangle=\left\langle\operatorname{Tr} \prod_{a, \sigma} \int D x_\sigma^a \Phi_B \Phi_F\right\rangle \end{align}$$

首先分析⟨Φ_B⟩，结合两个序参量$q_{a b}=\frac{1}{N} \sum_i^N s_i^a s_i^b \quad\text{for}\quad a \neq b$ 和 $m_\mu^a=\frac{1}{N} \sum_i \xi_i^\mu s_i^a$ 有： $$\begin{aligned} \left\langle\Phi_B\right\rangle=&\exp \left\{\frac{1}{2 N} \sum_{i, \mu \in B}\left[\sum_a s_i^a\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\frac{\beta}{\sqrt{N}} \sum_{j, v \in F} X_{\mu \nu} \xi_j^v s_j^a\right)\right]^2\right\} \\ \left\langle\Phi_B\right\rangle= & \int \prod_{a \neq b} \frac{d q_{a b} d \hat{q}_{a b}}{2 \pi / N} \prod_{a, \mu \in F} \frac{d m_\mu^a d \hat{m}_\mu^a}{2 \pi / N} \\ & \times \exp \left[-\frac{1}{2} N \sum_{a \neq b} \hat{q}_{a b} q_{a b}+\frac{1}{2} \sum_{a \neq b} \hat{q}_{a b} \sum_i s_i^a s_i^b-N \sum_{a, \mu \in F} m_\mu^a \hat{m}_\mu^a+\sum_{a, \mu \in F} \hat{m}_\mu^a \sum_i \xi_i^\mu s_i^a\right] \\ & \times \exp \left[\frac{1}{2} \sum_{\mu \in B} \sum_a\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\beta \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_v^a\right)^2\right] \\ & \times \exp \left[\frac{1}{2} \sum_{\mu \in B} \sum_{a \neq b} q_{a b}\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\beta \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_v^a\right)\right. \\ & \left.\times\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^b+\beta \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_v^b\right)\right] \end{aligned}$$

再结合δ函数傅里叶变换与复本对称的应用：

$$ \begin{aligned} \left\langle\Phi_B\right\rangle=&\int \frac{d q d \hat{q}}{(2 \pi / N)^{n(n-1)}} \frac{d m d \hat{m}}{(2 \pi / N)^{n S}}-N n \sum_{\mu \in F} m_\mu \hat{m}_\mu \exp \left[-\frac{1}{2} N n(n-1) \hat{q} q \quad+\frac{1}{2} \hat{q} \sum_{a \neq b} \sum_i s_i^a s_i^b+\sum_{a, \mu \in F} \hat{m}_\mu \sum_i \xi_i^\mu s_i^a\right] \exp \left[\frac { 1 } { 2 } \sum _ { \mu \in B } \sum _ { a } \left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a\right.\left.\quad+\beta \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_v\right)^2\right] \\ &\times \exp \left[\frac{q}{2} \sum_{\mu \in B} \sum_{a \neq b}\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\beta \sqrt{N} \sum_{v \in F} X_{\mu \nu} m_\nu\right) \times\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^b+\beta \sqrt{N} \sum_{v \in F} X_{\mu \nu} m_v\right)\right] \\ =&\int \frac{d q d \hat{q}}{(2 \pi / N)^{n(n-1)}} \frac{d m d \hat{m}}{(2 \pi / N)^{n S}} \exp \left[-\frac{1}{2} N n(n-1) \hat{q} q+\frac{1}{2} \hat{q} \sum_{a \neq b} \sum_i s_i^a s_i^b-N n \sum_{\mu \in F} m_\mu \hat{m}_\mu +\sum_{a, \mu \in F} \hat{m}_\mu \sum_i \xi_i^\mu s_i^a\right] \\ &\times \exp \left[\frac{1-q}{2} \sum_{\mu \in B} \sum_a\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\beta \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_v\right)^2\right] \\ & \times \exp \left[\frac{q}{2} \sum_{\mu \in B}\left(\sum_{a, \sigma} \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\beta n \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_\nu\right)^2\right] \\ = & \int \frac{d q d \hat{q}}{(2 \pi / N)^{n(n-1)}} \frac{d m d \hat{m}}{(2 \pi / N)^{n S}} \prod_{\mu, a} D y_\mu^a \prod_\mu D z_\mu \\ & \times \exp \left[-\frac{1}{2} N n(n-1) \hat{q} q+\frac{1}{2} \hat{q} \sum_{a \neq b} \sum_i s_i^a s_i^b-N n \sum_{\mu \in F} m_\mu \hat{m}_\mu+\sum_{a, \mu \in F} \hat{m}_\mu \sum_i \xi_i^\mu s_i^a\right] \\ & \times \exp \left[\sqrt{1-q} \sum_{\mu \in B} \sum_a\left(\sum_\sigma \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\beta \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_\nu\right) y_\mu^a\right] \\ & \times \exp \left[\sqrt{q} \sum_{\mu \in B}\left(\sum_{a, \sigma} \eta_\mu^\sigma \sqrt{\beta \lambda_\sigma} x_\sigma^a+\beta n \sqrt{N} \sum_{\nu \in F} X_{\mu \nu} m_\nu\right) z_\mu\right] \\ = & \int \frac{d q d \hat{q}}{(2 \pi / N)^{n(n-1)}} \frac{d m d \hat{m}}{(2 \pi / N)^{n S}} \prod_{\mu, a} D y_\mu^a \prod_\mu D z_\mu \\ & \times \exp \left[-\frac{1}{2} N n(n-1) \hat{q} q+\frac{1}{2} \hat{q} \sum_{a \neq b} \sum_i s_i^a s_i^b-N n \sum_{\mu \in F} m_\mu \hat{m}_\mu+\sum_{a, \mu \in F} \hat{m}_\mu \sum_i \xi_i^\mu s_i^a\right] \\ & \times \exp \left[\sum_{a, \sigma} x_\sigma^a \sqrt{\beta \lambda_\sigma} \sum_{\mu \in B} \eta_\mu^\sigma\left(\sqrt{1-q} y_\mu^a+\sqrt{q} z_\mu\right)\right] \\ & \times \exp \left[\beta \sqrt{N} \sum_{a, \mu \in B} \sum_{v \in F} X_{\mu \nu} m_\nu\left(\sqrt{1-q} y_\mu^a+\sqrt{q} z_\mu\right)\right] \end{aligned} $$

参考原书

$$\begin{align} F\left(r_{\rho \sigma}, q_{\rho \sigma}, m_\rho^1\right)= & -\frac{\alpha \beta^2}{2} \sum_{\rho, \sigma} r_{\rho \sigma} q_{\rho \sigma}-\frac{\alpha}{2} \operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}] -\frac{\beta}{2} \sum_\rho\left(m_\rho^1\right)^2+\left\langle\ln \operatorname{Tr} e^{\beta H_{\xi^1}}\right\rangle_{\xi^1} \label{freeenergy1} \end{align}$$

$$\begin{align} r_{\rho \sigma}&=\frac{1}{\alpha} \sum_{\mu \geq 2} m_\rho^\mu m_\sigma^\mu \\ m_\rho^\mu&=\frac{1}{N} \sum_i \xi_i^\mu S_i^\rho \\ q_{\rho \sigma}&=\frac{1}{N} \sum_i S_i^\rho S_i^\sigma \end{align}$$

$$\begin{align} \langle\ln Z\rangle=\lim _{n \rightarrow 0} \frac{\ln \left\langle Z^n\right\rangle}{n}=\lim _{n \rightarrow 0} \frac{\ln e^{N F\left(\theta^*\right)}}{n}=N \lim _{n \rightarrow 0} \frac{F\left(\theta^*\right)}{n} \label{847} \end{align}$$

Replica-Symmetric Ansätz

讨论到以上的情形，为了继续分析，需要对重叠矩阵做一个近似（考虑最简单的形式）：任意两个纯态应该是对称的。这被称为副本对称（RS）假设。

$$\begin{align} \left\{\begin{array}{l} r_{\rho \sigma}=r, \forall \rho, \sigma \\ m_\rho^1=m, \forall \rho \\ q_{\rho \sigma}=q, \forall \rho \neq \sigma \end{array} \right. \end{align}$$

将$\eqref{freeenergy1}$改写为： $$ \begin{aligned} F(r, q, m)= & -\frac{\alpha \beta^2}{2} r q\left(n^2-n\right)-\frac{\alpha \beta^2}{2} n r-\frac{\alpha}{2} \operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}] \\ & -\frac{\beta}{2} n m^2+\left\langle\ln \operatorname{Tr} e^{\beta H_{\xi}}\right\rangle \end{aligned} $$

结合$\eqref{847}$: $$\begin{align} \langle\ln Z\rangle= & \frac{N \alpha \beta^2 r q}{2}-\frac{N \alpha \beta^2 r}{2}-\frac{\alpha N}{2} \lim _{n \rightarrow 0} \frac{\operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}]}{n}-\frac{\beta N m^2}{2}+N \lim _{n \rightarrow 0} \frac{\left\langle\ln \operatorname{Tr} e^{\beta H_{\xi^1}}\right\rangle}{n} \label{zn3}\\ \beta H_{\xi^1}=&\beta m \xi^1 \sum_\rho S^\rho+\frac{1}{2} \alpha \beta^2 r \sum_{\rho, \sigma} S^\rho S^\sigma \end{align}$$

首先计算$\eqref{zn3}$中最后一项： $$\begin{align} \operatorname{Tr} e^{\beta H_{\xi} 1} & =\operatorname{Tr} e^{\beta m \xi^1 \Sigma_\rho S^\rho+\frac{1}{2} \alpha \beta^2 r\left(\sum_\rho S^\rho\right)^2} \\ & :=\operatorname{Tr} e^{A\left(\Sigma_\rho S^\rho\right)^2+B \Sigma_\rho S^\rho} \\ & =\operatorname{Tr} \sqrt{\frac{A}{\pi}} \int d z e^{-A z^2+2 A z \sum_\rho S^\rho+B \Sigma_\rho S^\rho} \\ & =\sqrt{\frac{A}{\pi}} \int d z e^{-A z^2} \operatorname{Tr} \prod_\rho e^{(2 A z+B) S^\rho} \\ & =\sqrt{\frac{\alpha \beta^2 r}{2 \pi}} \int d z e^{-\frac{1}{2} \alpha \beta^2 r z^2}\left[2 \cosh \left(\alpha \beta^2 r z+\beta m \xi^1\right)\right]^n \\ & =\sqrt{\frac{\alpha \beta^2 r}{2 \pi}} \int d z e^{-\frac{1}{2} \alpha \beta^2 r z^2+n \ln \left[2 \cosh \left(\alpha \beta^2 r z+\beta m \xi^1\right)\right]} \\ & =\sqrt{\frac{1}{2 \pi}} \int d z e^{-\frac{1}{2} z^2+n \ln \left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]} . \end{align}$$

并且可知其极限为： $$\begin{align} \lim _{n \rightarrow 0} \operatorname{Tr} e^{\beta H_{\xi} 1}=\sqrt{\frac{1}{2 \pi}} \int d z e^{-\frac{1}{2} z^2}=1 \end{align}$$

由此可获得最后一项为： $$\begin{align} & \lim _{n \rightarrow 0} \frac{\left\langle\ln \operatorname{Tr} e^{\beta H_{\xi}}\right\rangle}{n} \\ & =\left\langle\lim _{n \rightarrow 0} \frac{\frac{d}{d n} \operatorname{Tr} e^{\beta H_{\xi} 1}}{\operatorname{Tr} e^{\beta H_{\xi} 1}}\right\rangle \\ & =\left\langle\sqrt{\frac{1}{2 \pi}} \lim _{n \rightarrow 0} \frac{d}{d n} \int d z e^{-\frac{1}{2} z^2+n \ln \left[2 \cosh \left(\beta \sqrt{a r} z+\beta m \xi^1\right)\right]}\right\rangle \\ & =\left\langle\sqrt{\frac{1}{2 \pi}} \lim _{n \rightarrow 0} \int d z e^{-\frac{1}{2} z^2} \frac{d}{d n}\left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]^n\right\rangle \\ & =\left\langle\sqrt{\frac{1}{2 \pi}} \lim _{n \rightarrow 0} \int d z e^{-\frac{1}{2} z^2}\left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]^n \ln \left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]\right\rangle \\ & =\left\langle\sqrt{\frac{1}{2 \pi}} \int d z e^{-\frac{1}{2} z^2} \lim _{n \rightarrow 0}\left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]^n \ln \left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]\right\rangle \\ & =\left\langle\sqrt{\frac{1}{2 \pi}} \int d z e^{-\frac{1}{2} z^2} \ln \left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]\right\rangle \\ & =\int D z\left\langle\ln \left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]\right\rangle \end{align}$$

此时用Dz表示对z的高斯积分。

然后计算$\eqref{zn3}$中第三项。由于Q是对称矩阵，将其对角化： $$\begin{align} \mathbf{A Q A}^{-1}=\Lambda=\operatorname{diag}\left(\lambda_1, \lambda_2, \ldots, \lambda_n\right) \end{align}$$

将ln [I − βQ]进行指数展开$\ln (1-x)=-\sum_{n=1}^{\infty} \frac{x^n}{n}$，得到： $$\begin{align} \operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}] & =\operatorname{Tr}\left\{\mathbf{A} \cdot \ln [\mathbf{I}-\beta \mathbf{Q}] \cdot \mathbf{A}^{-1}\right\} \\ & =-\operatorname{Tr}\left\{\sum_{l=1}^{\infty} \frac{\beta^l\left(\mathbf{A Q A}^{-1}\right)^l}{l}\right\} \\ & =-\operatorname{Tr}\left\{\sum_{l=1}^{\infty} \frac{\beta^l(\Lambda)^l}{l}\right\} \\ & =-\sum_{l=1}^{\infty} \frac{\beta^l}{l} \sum_{i=1}^n \lambda_i^l=\sum_{i=1}^n \ln \left[1-\beta \lambda_i\right] \end{align}$$

再结合矩阵恒等式Tr ln K = ln det K，可以计算Q的本征值： $$\begin{align} & \left|\begin{array}{cccc} 1-\lambda & q & \cdots & q \\ q & 1-\lambda & \cdots & q \\ \vdots & \vdots & & \vdots \\ q & q & \cdots & 1-\lambda \end{array}\right| \\ & =\left|\begin{array}{cccc} 1-\lambda+(n-1) q & 1-\lambda+(n-1) q & \cdots & 1-\lambda+(n-1) q \\ q & 1-\lambda & \cdots & q \\ \vdots & \vdots & & \vdots \\ q & q & \cdots & 1-\lambda \end{array}\right| \\ & =[1-\lambda+(n-1) q]\left|\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ q & 1-\lambda & \cdots & q \\ \vdots & \vdots & & \vdots \\ q & q & \cdots & 1-\lambda \end{array}\right| \\ & =[1-\lambda+(n-1) q]\left|\begin{array}{cccc} 1 & 1 & \cdots & 1 \\ 0 & 1-\lambda-q & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots 1-\lambda-q \end{array}\right| \\ & =[1-\lambda+(n-1) q](1-q-\lambda)^{n-1}=0 \end{align}$$

可以得到1个本征值(1 + (n − 1)q)与n − 1个本征值(1 − q)，可以将迹写为： $$\begin{align} \operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}]=\ln (1-\beta+\beta q-n \beta q)+(n-1) \ln (1-\beta+\beta q) \end{align}$$ 从而;

$$\begin{align} \lim _{n \rightarrow 0} \frac{\operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}]}{n} & =\lim _{n \rightarrow 0}\left[\frac{\ln \left(\frac{1-\beta+\beta q-n \beta q}{1-\beta+\beta q}\right)}{n}+\ln (1-\beta+\beta q)\right] \\ & =-\frac{\beta q}{1-\beta+\beta q}+\ln (1-\beta+\beta q) \end{align}$$

得到自由能为：

$$\begin{align} -\beta f & =\frac{1}{N}\langle\ln Z\rangle \\ & =\frac{\alpha \beta^2}{2} r(q-1)-\frac{\alpha}{2}\left[\ln (1-\beta+\beta q)-\frac{\beta q}{1-\beta+\beta q}\right]-\frac{\beta}{2} m^2 +\int D z\left\langle\ln \left[2 \cosh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right]\right\rangle \end{align}$$

求解其极值： $$\begin{align} \left\{\begin{array}{l} \frac{\partial(-\beta f)}{\partial r}=0 \\ \frac{\partial(-\beta f)}{\partial m}=0 \\ \frac{\partial(-\beta f)}{\partial q}=0 \end{array}\right. \end{align}$$

$$\begin{align} q= & -\frac{1}{\beta \sqrt{2 \pi \alpha r}} \int d z e^{-\frac{1}{2} z^2} z\left\langle\tanh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right\rangle+1 \\ = & \frac{1}{\beta \sqrt{2 \pi \alpha r}} \int d z \frac{d e^{-\frac{1}{2} z^2}}{d z}\left\langle\tanh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right\rangle+1 \\ = & \left.\frac{1}{\beta \sqrt{2 \pi \alpha r}} e^{-\frac{1}{2} z^2}\left\langle\tanh \left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right\rangle\right|_{-\infty} ^{+\infty} -\int D z\left\langle 1-\tanh ^2\left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right\rangle+1 \\ = & \int D z\left\langle\tanh ^2\left(\beta \sqrt{\alpha r} z+\beta m \xi^1\right)\right\rangle \\ = & \int D z \tanh ^2 \beta(\sqrt{\alpha r} z+m) \end{align}$$

可以得到联想记忆模型的鞍点方程： $$\begin{align} & q=\int D z \tanh ^2 \beta(\sqrt{\alpha r} z+m) \label{860} \\ & m=\int D z\langle\xi \tanh \beta(\sqrt{\alpha r} z+m \xi)\rangle=\int D z \tanh \beta(\sqrt{\alpha r} z+m) \label{861} \\ & r=\frac{q}{(1-\beta+\beta q)^2} \label{862}\\ \end{align}$$

可以从以上的内容分析相变点。

Zero-Temperature Limit

当T → 0(β → ∞)时候，有：

$$\tanh (\beta x) \rightarrow \operatorname{sign}(x)=\left\{ \begin{array}{ll} 1 & x>0 \\ 0 & x=0 \\ -1 & x<0 \end{array} \right.$$

$\eqref{861}$为：

$$\begin{align} m&=\int D z\operatorname{sign}(\sqrt{\alpha r} z+m)+O(T) \\ & =\operatorname{erf}\left( \frac{m}{\sqrt{2 \alpha r}} \right)+O(T) \end{align}$$

上面这个变换利用正态分布与误差函数之间的关系完成。 另一方面，当β → ∞：

$$\begin{align} 1-q & =\int \frac{d z}{\sqrt{2 \pi}} e^{-\frac{z^2}{2}}\left(1-\tanh ^2 \beta(\sqrt{\alpha r} z+m)\right) \\ & \left.\simeq \frac{1}{\sqrt{2 \pi}} e^{-\frac{z^2}{2}}\right|_{\tanh ^2 \beta(\sqrt{\alpha r} z+m)=0} \int d z\left(1-\tanh ^2 \beta(\sqrt{\alpha r} z+m)\right) \\ & =\frac{1}{\sqrt{2 \pi}} e^{-\frac{m^2}{2 \alpha r}} \frac{1}{\beta \sqrt{\alpha r}} \int d z \frac{\partial}{\partial z} \tanh \beta(\sqrt{\alpha r} z+m) \\ & =\frac{2}{\sqrt{2 \pi}} \frac{1}{\beta \sqrt{\alpha r}} e^{-\frac{m^2}{2 \alpha r}} \end{align}$$

$\eqref{860}$产生q = 1 − CT，其中

$$\begin{align} C \stackrel{\text { def }}{=} \sqrt{\frac{2}{\pi r \alpha}} e^{-\frac{m^2}{2 \alpha r}} \end{align}$$ 将$\eqref{862}$变为r = (1 − C)⁻²。

通过定义辅助变量$y=m / \sqrt{2 \alpha r}$将m 和 r 减少为一个方程： $$\begin{align} \operatorname{erf}(y)=y\left(\sqrt{2 \alpha}+\frac{2}{\sqrt{\pi}} e^{-y^2}\right) . \end{align}$$

其中的一个恒定解为y = m = 0。对于α ≥ α_c = 0.138只有为0的唯一解；当a < α_c时,m ≠ 0的解出现；当α = α_c时，m = 0.967。

通过求解m = erf (y)可以得到m的值。图中纵坐标表示误差$P_{\text{error}}=\frac{1-m}{2}$，横坐标表述存储的比例。可以发现，在图中存在一个跃变的点，当α——c = 0.138的时候，误差跳至0.5，这是一个不连续的转变，此时代表跳转至玻璃相；当α < α_c时，误差很小，表示此时网络可以从之前学习的模式中“恢复”；当α > α_c时，误差为0.5，为瞎猜的几率，可以认为此时网络不能恢复之前学习的模式。

Append

高斯积分和误差函数之间的关系

高斯函数（Gaussian function）与误差函数（erf）之间的关系主要通过正态分布的累积分布函数（CDF）来体现。为了详细说明这一点，我们从标准正态分布及其累积分布函数出发。

标准正态分布

标准正态分布的概率密度函数（PDF）定义为：

$$f(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2 / 2}$$

标准正态分布的累积分布函数

累积分布函数（CDF）是从负无穷到某个值 x 的概率密度函数的积分。对于标准正态分布，CDF 通常记作 Φ(x)：

$$\Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} e^{-t^2 / 2} \, dt$$

误差函数的定义

误差函数（erf）定义为：

$$\text{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} \, dt$$

高斯函数与误差函数的关系

我们可以通过一些变换将标准正态分布的 CDF 表示为误差函数。首先，考虑累积分布函数从负无穷大积分到某个值 x：

$$\Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}} e^{-t^2 / 2} \, dt$$

通过代换 $u = \frac{t}{\sqrt{2}}$，我们得到：

$$du = \frac{dt}{\sqrt{2}} \quad \Rightarrow \quad dt = \sqrt{2} \, du$$

因此，积分变为：

$$\Phi(x) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-t^2 / 2} \, dt = \int_{-\infty}^{x/\sqrt{2}} \frac{1}{\sqrt{2\pi}} e^{-2u^2 / 2} \cdot \sqrt{2} \, du$$

简化后，我们得到：

$$\Phi(x) = \frac{1}{2} \left(1 + \text{erf}\left(\frac{x}{\sqrt{2}}\right)\right)$$

因此，标准正态分布的累积分布函数与误差函数之间的关系为：

$$\Phi(x) = \frac{1}{2} \left(1 + \text{erf}\left(\frac{x}{\sqrt{2}}\right)\right)$$

逆误差函数

误差函数的反函数（逆误差函数）是一个重要工具，用于将累积分布函数反转：

$$x = \sqrt{2} \, \text{erf}^{-1}(2\Phi(x) - 1)$$

示例

假设我们想计算标准正态分布 𝒩(0, 1) 的某个值 x 的累积分布函数。使用 Python，我们可以如下计算：

1
2
3
4
5

import scipy.special as sp

x = 1.0
Phi_x = 0.5 * (1 + sp.erf(x / np.sqrt(2)))
print(Phi_x)

这个示例计算了 x = 1 时的累积分布函数值。

总结

高斯函数（正态分布）与误差函数之间的关系通过正态分布的累积分布函数（CDF）体现。正态分布的累积分布函数可以表示为误差函数的形式：

$$\Phi(x) = \frac{1}{2} \left(1 + \text{erf}\left(\frac{x}{\sqrt{2}}\right)\right)$$

这使得误差函数成为研究正态分布及其相关问题的一个重要工具。

Link: * Hopfield 相图 * Replica

根据前文已知： $$\begin{equation} \begin{aligned} \omega_{ij}&=\frac{1}{N}\xi_i\xi_j \quad i\neq j\\ \omega_{ii}&=0 \end{aligned} \label{omega} \end{equation} $$ $$\begin{equation} H=-\frac{1}{2} \sum_{i, j}^N w_{i j} S_i S_j \label{hamiltonian} \end{equation}$$

Free Energy

假定能够存储P = αN个随机模式。结合$\eqref{omega}$与$\eqref{hamiltonian}$有：

$$\begin{aligned} \left\langle Z^n\right\rangle & =\left\langle\operatorname{Tr} \exp \left[\frac{\beta}{2 N} \sum_{i, j}^N \sum_{\mu=1}^P \sum_{\rho=1}^n \xi_i^\mu \xi_j^\mu S_i^\rho S_j^\rho\right]\right\rangle \\ & =\left\langle\operatorname{Tr} \exp \left[\frac{\beta}{2 N} \sum_{\rho, \mu}\left(\sum_i \xi_i^\mu S_i^\rho\right)\left(\sum_j \xi_j^\mu S_j^\rho\right)\right]\right\rangle \\ & =\left\langle\operatorname{Tr} \exp \left[\frac{\beta N}{2} \sum_{\rho, \mu}\left(\frac{1}{N} \sum_i \xi_i^\mu S_i^\rho\right)^2\right]\right\rangle \\ & =\left\langle\operatorname{Tr} \prod_{\rho, \mu} \exp \left[\frac{\beta N}{2}\left(\frac{1}{N} \sum_i \xi_i^\mu S_i^\rho\right)^2\right]\right\rangle, \end{aligned}$$

其中Tr 表示对所有构型S求和，⟨⋅⟩表示淬火无序平均。使用Hubbard-Stratonovich transformation技巧，将其中指数中的二次项转化为高斯积分：

$$\begin{align} e^{a b^2}=\sqrt{\frac{a}{\pi}} \int e^{-a x^2+2 a b x} d x \end{align}$$

设定这些参数为：

$$\begin{align} \left\{\begin{array}{l} b \rightarrow \frac{1}{N} \sum_i \xi_i^\mu S_i^\rho \\ x \rightarrow m_\rho^\mu \\ a \rightarrow \frac{\beta N}{2} \end{array} \right. \end{align}$$

通过对m_ρ^μ积分得到： $$\begin{align} \left\langle Z^n\right\rangle & =\left\langle\operatorname{Tr} \int \prod_{\rho, \mu} \sqrt{\frac{\beta N}{2 \pi}} d m_\rho^\mu \exp \left[-\frac{\beta N}{2}\left(m_\rho^\mu\right)^2+\beta m_\rho^\mu \sum_i \xi_i^\mu S_i^\rho\right]\right\rangle \\ & =\left\langle\operatorname{Tr} \int \prod_{\rho', \mu'} \sqrt{\frac{\beta N}{2 \pi}} d m_{\rho'}^{\mu'} \exp \left[-\frac{\beta N}{2} \sum_{\rho, \mu}\left(m_\rho^\mu\right)^2+\beta \sum_{\rho, \mu} m_\rho^\mu \sum_i \xi_i^\mu S_i^\rho\right]\right\rangle \\ & =\left\langle\operatorname{Tr} \int \prod_{\rho', \mu'} \sqrt{\frac{\beta N}{2 \pi}} d m_{\rho'}^{\mu'} \exp \left[-\frac{\beta N}{2} \sum_{\mu \geq 2} \sum_\rho\left(m_\rho^\mu\right)^2+\right.\left.\beta \sum_{\mu \geq 2} \sum_\rho m_\rho^\mu \sum_i \xi_i^\mu S_i^\rho-\frac{\beta N}{2} \sum_\rho\left(m_\rho^1\right)^2+\beta \sum_\rho m_\rho^1 \sum_i \xi_i^1 S_i^\rho\right]\right\rangle . \end{align}$$

最后假设只有第一个模式被恢复μ = 1，在表达式中将其单独分出来，因此交错项m_ρ^μ ∼ O(1)。

接下来考虑没有被恢复的模式(μ ≥ 2)，由于没有被恢复，因此⟨∑_iξ_i^μS_i^ρ⟩_ξ = 0并且⟨(∑_iξ_i^μS_i^ρ)²⟩_ξ = N + ⟨∑_i ≠ jξ_i^μξ_j^μS_i^ρS_j^ρ⟩_ξ = N。m_ρ^μ(μ ≥ 2)的大小数量级为（可以将这个过程认为是一维链上的随机游走过程，计算的是绝对值的期望）： $$\begin{align} m_\rho^\mu=\frac{1}{N} \sum_i \xi_i^\mu S_i^\rho \approx O\left(\frac{1}{\sqrt{N}}\right) \end{align}$$

为了将整个数量级变为O(1)，因此将m_ρ^μ进行变换$m_\rho^\mu \rightarrow \frac{m_\rho^\mu}{\sqrt{\beta N}}$。变化之后为：

$$ \begin{align} \left\langle Z^n\right\rangle= & \left(\frac{1}{\sqrt{2 \pi}}\right)^{n(P-1)}\left\langle\operatorname { T r } \int \prod _ { \rho' , \mu' > 1 } d m _ { \rho' } ^ { \mu' } \prod _ { \rho } \sqrt { \beta N } d m _ { \rho } ^ { 1 } \operatorname { e x p } \left[-\frac{1}{2} \sum_{\mu \geq 2} \sum_\rho\left(m_\rho^\mu\right)^2+ \sqrt{\frac{\beta}{N}} \sum_{\mu \geq 2} \sum_\rho m_\rho^\mu \sum_i \xi_i^\mu S_i^\rho-\frac{\beta N}{2} \sum_\rho\left(m_\rho^1\right)^2+\beta \sum_\rho m_\rho^1 \sum_i \xi_i^1 S_i^\rho\right]\right\rangle \label{znstart} \end{align} $$

接下来将其中的每一部分进行分解化简。

对于μ ≥ 2的部分，由于模式ξ是随机选取的，在模式的数量较多的情况下可以认为存在⟨exp (Aξ)⟩_{{ξ = ±1}} = exp (−A) + exp (A) = 2cosh (A) ∝ exp (ln cosh (A))，再结合近似$\ln \cosh x=\frac{x^2}{2}+\cdots$ 在 x → 0：

$$\begin{aligned} & \left\langle\exp \left[\sqrt{\frac{\beta}{N}} \sum_{\mu \geq 2} \sum_\rho m_\rho^\mu \sum_i \xi_i^\mu S_i^\rho\right]\right\rangle_{\xi_i^\mu: \mu>1} \\ & \propto \exp \left[\sum_{\mu \geq 2, i} \ln \cosh \left(\sqrt{\frac{\beta}{N}} \sum_\rho m_\rho^\mu S_i^\rho\right)\right] \\ & \cong \exp \left[\sum_{\mu \geq 2} \sum_i \frac{\beta}{2 N}\left(\sum_\rho m_\rho^\mu S_i^\rho\right)^2\right] \end{aligned}$$

然后将平方项(m_ρ^μ)²改写： $$\begin{align} \sum_\rho\left(m_\rho^\mu\right)^2&=\sum_{\rho, \sigma} m_\rho^\mu \delta_{\rho \sigma} m_\sigma^\mu \\ \frac{1}{N} \sum_i\left(\sum_\rho m_\rho^\mu S_i^\rho\right)^2 & =\frac{1}{N} \sum_i \sum_\rho m_\rho^\mu S_i^\rho \sum_\sigma m_\sigma^\mu S_i^\sigma \\ & =\sum_{\rho, \sigma} m_\rho^\mu \frac{1}{N} \sum_i S_i^\rho S_i^\sigma m_\sigma^\mu \\ & :=\sum_{\rho, \sigma} m_\rho^\mu q_{\rho \sigma} m_\sigma^\mu \end{align}$$

为了进一步化简该表达式：

$$\begin{align} \kappa_{\rho \sigma}=\delta_{\rho \sigma}-\frac{\beta}{N} \sum_i S_i^\rho S_i^\sigma:=\delta_{\rho \sigma}-\beta q_{\rho \sigma} \end{align}$$

写为矩阵形式：

$$\begin{align} \mathbf{K}&=\mathbf{I}-\beta \mathbf{Q} \\ q_{\rho \sigma}&=\left\{\begin{array}{ll} \frac{1}{N} \sum_i S_i^\rho S_i^\sigma & \rho \neq \sigma \\ 1 & \rho=\sigma \end{array}\right. \label{q} \end{align}$$

K, Q 是对称的 n × n 矩阵。因此有：

$$\begin{equation} \begin{aligned} \left\langle Z^n\right\rangle & \propto \operatorname{Tr} \int \prod_{\rho, \sigma} d q_{\rho \sigma} \delta\left(q_{\rho \sigma}-\frac{1}{N} \sum_i S_i^\rho S_i^\sigma\right) \\ & \times \prod_{\mu \geq 2, \rho} d m_\rho^\mu \exp \left[-\frac{1}{2} \sum_{\mu \geq 2} \sum_{\rho, \sigma} m_\rho^\mu \kappa_{\rho \sigma} m_\sigma^\mu\right] \\ & \times\left\langle\int \prod_\rho d m_\rho^1 \exp \left[-\frac{\beta N}{2} \sum_\rho\left(m_\rho^1\right)^2+\beta \sum_\rho m_\rho^1 \sum_i \xi_i^1 S_i^\rho\right]\right\rangle_{\xi^1} \end{aligned} \label{zn} \end{equation}$$

式$\eqref{zn}$中第一行的作用，将$\eqref{q}$关系明确的写进计算公式中，忽略不相关因子。利用多变量高斯积分：

$$ \int_{R^n} d \mathbf{m} e^{-\mathbf{M}^{\mathrm{T}} \mathbf{K M}}=\sqrt{\frac{\pi^n}{\operatorname{det}(\mathbf{K})}}, $$

可以得到：

$$\begin{align} \int \prod_{\mu \geq 2, \rho} d m_\rho^\mu \exp \left[-\frac{1}{2} \sum_{\mu \geq 2} \sum_{\rho, \sigma} m_\rho^\mu \kappa_{\rho \sigma} m_\sigma^\mu\right]=\frac{C}{(\operatorname{det} \mathbf{K})^{\frac{P-1}{2}}} \end{align}$$ 其中C是常量。

由于det (e^K) = e^Tr K，可得det K = e^Tr ln K，因此有：

$$\begin{align} (\operatorname{det} \mathbf{K})^{-\frac{P-1}{2}}=e^{-\frac{P-1}{2} \operatorname{Tr} \ln \mathbf{K}}=e^{-\frac{P-1}{2} \operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}]} \end{align}$$

再结合Dirac函数的傅里叶分解：

$$\begin{align} \delta(x)=\frac{1}{2 \pi} \int_{-\infty}^{+\infty} e^{-\mathrm{i} k x} d k, \end{align}$$

将$\eqref{zn}$中第一行Dirac函数进行变换：

$$ \begin{aligned} \delta\left(q_{\rho \sigma}-\frac{1}{N} \sum_i S_i^\rho S_i^\sigma\right) \propto \int_{-\infty}^{+\infty} d r_{\rho \sigma} \exp \left[-i r_{\rho \sigma} q_{\rho \sigma}+ i r_{\rho \sigma} \frac{1}{N} \sum_i S_i^\rho S_i^\sigma\right] \end{aligned} $$

将r进行重标度$r_{\rho \sigma} \rightarrow -\frac{\mathrm{i} N \alpha \beta^2}{2} r_{\rho \sigma}$，使得r_ρσ ∼ O(1)，这种操作对计算结果没有影响。上式变换忽略了因子项。 这样操作使得最后的结果变得简洁。其中α = P/N。

$$\begin{align} \delta\left(q_{\rho \sigma}-\frac{1}{N} \sum_i S_i^\rho S_i^\sigma\right) \propto \int_{-\infty}^{+\infty} d r_{\rho \sigma} \exp \left[-\frac{N \alpha \beta^2}{2} r_{\rho \sigma} q_{\rho \sigma}+\frac{\alpha \beta^2}{2} \sum_{i} r_{\rho \sigma} S_i^\rho S_i^\sigma\right] \label{fuliye} \end{align}$$

接下来将$\eqref{zn}$中的最后一项与$\eqref{fuliye}$中最后一项组合一起： $$\begin{align} \left\langle\operatorname{Tr} \exp \left[\beta \sum_\rho m_\rho^1 \sum_i \xi_i^1 S_i^\rho+\frac{\alpha \beta^2}{2} \sum_{i} r_{\rho \sigma} S_i^\rho S_i^\sigma\right]\right\rangle_{\xi^1} = & \left\langle\exp \left\{\sum_i \ln \operatorname{Tr} \exp \left(\beta \sum_\rho m_\rho^1 \xi_i^1 S^\rho+\frac{\alpha \beta^2}{2} r_{\rho \sigma} S^\rho S^\sigma\right)\right\}\right\rangle_{\xi^1} \\ = & \exp \left\{N\left\langle\ln \operatorname{Tr} \exp \left(\beta \sum_\rho m_\rho^1 \xi^1 S^\rho+\frac{\alpha \beta^2}{2} r_{\rho \sigma} S^\rho S^\sigma\right)\right\rangle_{\xi^1}\right\} \\ := & \exp \left\{N\left\langle\ln \operatorname{Tr} \exp \left(\beta H_{\xi^1}\right)\right\rangle_{\xi^1}\right\} \end{align}$$

定义，因为上面分析的是一个q_ρσ的δ函数展开，原来表达式中在积分符号中含有∏_ρσ，在指数中将会增加求和符号∑_ρσ： $$\begin{align} \beta H_{\xi^1}=\frac{1}{2} \alpha \beta^2 \sum_{\rho \sigma}r_{\rho \sigma} S^\rho S^\sigma+\beta \sum_\rho m_\rho^1 \xi^1 S^\rho, \end{align}$$

最终$\eqref{zn}$转化为： $$\begin{equation} \begin{aligned} & \left\langle Z^n\right\rangle \propto \int \prod_\rho d m_\rho^1 \prod_{\rho, \sigma} d q_{\rho \sigma} d r_{\rho \sigma} \exp \left[-\frac{N}{2} \alpha \beta^2 \sum_{\rho, \sigma} r_{\rho \sigma} q_{\rho \sigma}\right] \\ & \times \exp \left[-\frac{P-1}{2} \operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}]\right] \exp \left[-\frac{\beta N}{2} \sum_\rho\left(m_\rho^1\right)^2+N\left\langle\ln \operatorname{Tr} e^{\beta H_{\xi^1}}\right\rangle_{\xi^1}\right] \end{aligned} \label{zn2} \end{equation}$$

因为假定N足够大，使用 Laplace’s Method：

$$ \int_a^b e^{N f(z)} d z \approx \sqrt{\frac{2 \pi}{-N f^{\prime \prime}\left(z_0\right)}} e^{N f\left(z_0\right)} $$

其中Z₀表示f(z)的极值点，这样就把积分转化为在极值点起作用的标量。采用近似 ⟨Zⁿ⟩ ∼ e^NF(θ*)，其中用θ表示序参量。

F(θ^*) = max_θF(θ)

淬火无序平均值变为： $$\begin{align} \langle\ln Z\rangle=\lim _{n \rightarrow 0} \frac{\ln \left\langle Z^n\right\rangle}{n}=\lim _{n \rightarrow 0} \frac{\ln e^{N F\left(\theta^*\right)}}{n}=N \lim _{n \rightarrow 0} \frac{F\left(\theta^*\right)}{n} \label{847} \end{align}$$ 当n足够大时，近似有P − 1 ≃ P = αN：

$$\begin{align} F\left(r_{\rho \sigma}, q_{\rho \sigma}, m_\rho^1\right)= & -\frac{\alpha \beta^2}{2} \sum_{\rho, \sigma} r_{\rho \sigma} q_{\rho \sigma}-\frac{\alpha}{2} \operatorname{Tr} \ln [\mathbf{I}-\beta \mathbf{Q}] -\frac{\beta}{2} \sum_\rho\left(m_\rho^1\right)^2+\left\langle\ln \operatorname{Tr} e^{\beta H_{\xi^1}}\right\rangle_{\xi^1} \label{freeenergy1} \end{align}$$

序参量分析

通过以上讨论已经得到与相变相关的序参量(r_ρσ, q_ρσ, m_ρ¹)，即其满足的方程$\eqref{freeenergy1}$。在继续进行极值分析前，先回顾一下其中每一个变量的含义。α是存储模式的比例；β是逆温度；r是在进行傅里叶变换的时候引入的变量；q是在进行变量替换时候引入的变量；Q是q矩阵形式；m是在进行高斯积分引入的变量；H_ξ是将S重新组合之后的表达形式。

接下来计算F(r_ρσ, q_ρσ, m_ρ¹)的极值，首先处理q_ρσ，结合$\eqref{zn}$：

$$\begin{align} \frac{\partial F}{\partial q_{\rho \sigma}}=0 \Rightarrow \frac{\partial}{\partial q_{\rho \sigma}}\left[-\frac{N \alpha \beta^2}{2} \sum_{\rho, \sigma} r_{\rho \sigma} q_{\rho \sigma}+\frac{\beta}{2}(\sqrt{\beta N})^2 \sum_{\mu \geq 2} \sum_{\rho, \sigma} m_\rho^\mu q_{\rho \sigma} m_\sigma^\mu\right]=0 \text {, } \end{align}$$

得到解为： $$\begin{align} r_{\rho \sigma}=\frac{1}{\alpha} \sum_{\mu \geq 2} m_\rho^\mu m_\sigma^\mu \end{align}$$

从表达式中可以看出，之前进行标度的变换，目的是为了使得最终的共轭关系变得更直观。

然后处理m_ρ^μ，结合最开始的$\eqref{znstart}$：

$$\begin{align} \frac{\partial F}{\partial m_\rho^\mu}=0 \Rightarrow \frac{\partial}{\partial m_\rho^\mu}\left[-\frac{\beta N}{2}\left(m_\rho^\mu\right)^2+\beta m_\rho^\mu \sum_i \xi_i^\mu S_i^\rho\right]=0 \end{align}$$

可以得到： $$\begin{align} m_\rho^\mu=\frac{1}{N} \sum_i \xi_i^\mu S_i^\rho \end{align}$$

可以看出m_ρ^μ在处理不同的模式之间影响的作用，也称作交叠项。

最后处理r_ρσ，结合$\eqref{zn2}$： $$\begin{align} \frac{\partial F}{\partial r_{\rho \sigma}}=0 \Rightarrow \frac{\partial}{\partial r_{\rho \sigma}}\left[-\frac{N \alpha \beta^2}{2} \sum_{\rho, \sigma} r_{\rho \sigma} q_{\rho \sigma}+\frac{\alpha \beta^2}{2} \sum_{i, \rho, \sigma} r_{\rho \sigma} S_i^\rho S_i^\sigma\right]=0, \end{align}$$

得到序参量： $$\begin{align} q_{\rho \sigma}=\frac{1}{N} \sum_i S_i^\rho S_i^\sigma \end{align}$$

q_ρσ通常理解为两个纯态之间的相互重叠。如果一个单一状态主导了相空间，爱德华兹-安德森序参量就用来表征该状态的大小。

Append

行列式与迹的关系证明

1
2

M = {{a, b}, {c, d}};
Sum[Det[MatrixPower[M, k]]/(k!), {k, 0, Infinity}] // FullSimplify

Reference: * Hopfield模型的统计物理视频 * Analytic and Algorithmic Solution of Random Satisfiability Problems * Gibbs states and the set of solutions of random constraint satisfaction problems

讨论Hopfield Model 的相关内容，包含其本身的一些性值以及其在组合优化问题中的应用。

Reference: * Neural networks and physical systems with emergent collective computational abilities * “Neural” Computation of Decisions in Optimization Problems * 神经网络的统计力学 * Neural Networks for Combinatorial Optimization: A Review of More Than a Decade of Research

Link: * Hopfield 自由能求解