深度学习学习教程，从入门到精通，直面配分函数 — 完整知识点与代码实现(18）
·
直面配分函数 — 完整知识点与代码实现
一、核心背景：为什么需要直面配分函数？
很多概率模型（如玻尔兹曼机、基于能量的模型）的概率分布形式为：
pmodel(x)=p~model(x)Zp_{\text{model}}(\mathbf{x}) = \frac{\tilde{p}_{\text{model}}(\mathbf{x})}{Z}pmodel(x)=Zp~model(x)
其中 p~\tilde{p}p~ 是未归一化的概率（易计算），Z=∑xp~(x)Z = \sum_{\mathbf{x}} \tilde{p}(\mathbf{x})Z=∑xp~(x) 是配分函数（partition function），求和/积分遍历所有可能的状态，在高维空间中几乎不可计算。
本章的所有方法都在回答同一个问题：如何在无法计算 ZZZ 的情况下训练和使用这些模型？
二、知识点 1：对数似然梯度
2.1 数学推导
对数似然为：
logpmodel(x)=logp~model(x)−logZ\log p_{\text{model}}(\mathbf{x}) = \log \tilde{p}_{\text{model}}(\mathbf{x}) - \log Zlogpmodel(x)=logp~model(x)−logZ
对参数 θ\thetaθ 求梯度：
∇θlogpmodel(x)=∇θlogp~model(x)−∇θlogZ\nabla_\theta \log p_{\text{model}}(\mathbf{x}) = \nabla_\theta \log \tilde{p}_{\text{model}}(\mathbf{x}) - \nabla_\theta \log Z∇θlogpmodel(x)=∇θlogp~model(x)−∇θlogZ
关键推导——第二项可以转化为期望：
∇θlogZ=1Z∇θZ=1Z∑x∇θp~(x)=∑xpmodel(x)∇θlogp~(x)\nabla_\theta \log Z = \frac{1}{Z} \nabla_\theta Z = \frac{1}{Z} \sum_{\mathbf{x}} \nabla_\theta \tilde{p}(\mathbf{x}) = \sum_{\mathbf{x}} p_{\text{model}}(\mathbf{x}) \nabla_\theta \log \tilde{p}(\mathbf{x})∇θlogZ=Z1∇θZ=Z1x∑∇θp~(x)=x∑pmodel(x)∇θlogp~(x)
即：
∇θlogZ=Ex∼pmodel[∇θlogp~model(x)]\nabla_\theta \log Z = \mathbb{E}_{\mathbf{x} \sim p_{\text{model}}} \left[ \nabla_\theta \log \tilde{p}_{\text{model}}(\mathbf{x}) \right]∇θlogZ=Ex∼pmodel[∇θlogp~model(x)]
最终梯度公式（两力平衡）：
∇θlogpmodel(x)=∇θlogp~model(x)⏟正相（增大训练数据概率）−Ex∼pmodel[∇θlogp~model(x)]⏟负相（降低模型生成样本概率）\nabla_\theta \log p_{\text{model}}(\mathbf{x}) = \underbrace{\nabla_\theta \log \tilde{p}_{\text{model}}(\mathbf{x})}_{\text{正相（增大训练数据概率）}} - \underbrace{\mathbb{E}_{\mathbf{x} \sim p_{\text{model}}} \left[ \nabla_\theta \log \tilde{p}_{\text{model}}(\mathbf{x}) \right]}_{\text{负相（降低模型生成样本概率）}}∇θlogpmodel(x)=正相（增大训练数据概率）∇θlogp~model(x)−负相（降低模型生成样本概率）Ex∼pmodel[∇θlogp~model(x)]
- 正相（Positive Phase）：推高训练数据的概率
- 负相（Negative Phase）：压低模型自生成样本的概率
- 难点在于：负相需要从 pmodelp_{\text{model}}pmodel 中采样，而这本身需要 ZZZ
2.2 案例代码：基于能量模型的对数似然梯度计算
import numpy as np
import matplotlib.pyplot as plt
# ============================================
# 知识点1：对数似然梯度的数值验证
# 以一个简单的基于能量的模型为例
# p_model(x) = exp(-E(x;theta)) / Z(theta)
# E(x; theta) = theta * x^2 （简单的二次能量函数）
# ============================================
def energy(x, theta):
"""
计算能量函数 E(x; theta) = theta * x^2
参数:
x: 数据点，numpy数组
theta: 模型参数（标量，控制能量曲面的"陡峭"程度）
返回:
能量值，形状与 x 相同
"""
return theta * x ** 2
def unnormalized_log_prob(x, theta):
"""
计算未归一化的对数概率 log p̃(x) = -E(x; theta)
在基于能量的模型中，p̃(x) = exp(-E(x; theta))
因此 log p̃(x) = -E(x; theta) = -theta * x^2
参数:
x: 数据点
theta: 模型参数
返回:
log p̃(x) 的值
"""
return -energy(x, theta)
def unnormalized_prob(x, theta):
"""
计算未归一化概率 p̃(x) = exp(-E(x; theta))
用于数值积分求配分函数 Z
参数:
x: 数据点
theta: 模型参数
返回:
p̃(x) 的值
"""
return np.exp(unnormalized_log_prob(x, theta))
def partition_function(theta, x_range=(-10, 10), n_points=10000):
"""
通过数值积分计算配分函数 Z = ∫ p̃(x) dx
注意：仅在一维连续空间可行，高维空间此方法不可行
参数:
theta: 模型参数
x_range: 积分范围
n_points: 积分点数（越多越精确）
返回:
配分函数 Z 的近似值
"""
# 在 x_range 范围内均匀生成采样点
x_points = np.linspace(x_range[0], x_range[1], n_points)
# 计算每个点的未归一化概率
y_points = unnormalized_prob(x_points, theta)
# 使用梯形法则进行数值积分
# np.trapz 使用梯形规则近似定积分
Z = np.trapz(y_points, x_points)
return Z
def log_partition_function(theta):
"""
计算 log Z，直接用 log 包裹数值积分结果
参数:
theta: 模型参数
返回:
log Z 的近似值
"""
return np.log(partition_function(theta))
def log_likelihood(data, theta):
"""
计算整个数据集的平均对数似然
L(theta) = (1/N) * sum_i [log p̃(x_i; theta) - log Z(theta)]
参数:
data: 训练数据，一维numpy数组
theta: 模型参数
返回:
平均对数似然值
"""
# 对每个数据点计算 log p̃(x_i)
log_p_tilde = unnormalized_log_prob(data, theta)
# 减去 log Z 得到真正的对数似然
log_Z = log_partition_function(theta)
# 返回平均对数似然
return np.mean(log_p_tilde - log_Z)
def gradient_log_unnormalized_prob(x, theta):
"""
计算正相梯度：∇_θ log p̃(x; theta)
由于 log p̃(x) = -theta * x^2
所以 ∂/∂θ log p̃(x) = -x^2
参数:
x: 数据点
theta: 模型参数（此处不需要，因为梯度不依赖theta）
返回:
正相梯度值
"""
return -x ** 2
def true_log_likelihood_gradient(data, theta):
"""
计算真实的对数似然梯度（通过数值微分验证）
使用有限差分法：∂L/∂θ ≈ (L(θ+ε) - L(θ-ε)) / (2ε)
参数:
data: 训练数据
theta: 模型参数
返回:
梯度的数值近似
"""
epsilon = 1e-5 # 微小扰动量
# 中心差分法，比前向差分更精确
grad = (log_likelihood(data, theta + epsilon) -
log_likelihood(data, theta - epsilon)) / (2 * epsilon)
return grad
def analytical_log_likelihood_gradient(data, theta, x_range=(-10, 10), n_samples=50000):
"""
通过解析公式计算对数似然梯度
∇_θ L = 正相 - 负相
正相 = (1/N) * sum_i ∇_θ log p̃(x_i)
负相 = E_{x~p_model} [∇_θ log p̃(x)]
负相通过从 p_model 采样来近似
采样方法：使用拒绝采样（rejection sampling）
参数:
data: 训练数据
theta: 模型参数
x_range: 采样范围
n_samples: 用于近似负相的采样数
返回:
梯度的估计值
"""
# === 正相计算 ===
# 对每个训练数据点计算 ∇_θ log p̃(x_i)
positive_phase = np.mean(gradient_log_unnormalized_prob(data, theta))
# === 负相计算（核心难点）===
# 需要从 p_model(x) 采样，使用拒绝采样
# 提议分布使用均匀分布 U(-5, 5)
# 接受条件：u < p̃(x) / (M * q(x))，其中 M 是上界常数
max_x = 5.0 # 采样范围边界
# 在网格上找到 p̃(x) 的最大值，作为拒绝采样的上界常数
x_grid = np.linspace(-max_x, max_x, 10000)
p_tilde_max = np.max(unnormalized_prob(x_grid, theta))
# 拒绝采样过程
samples = [] # 存储接受的样本
attempts = 0 # 记录尝试次数
while len(samples) < n_samples:
# 从均匀分布 U(-max_x, max_x) 采样候选点
x_candidate = np.random.uniform(-max_x, max_x)
# 计算接受概率：p̃(x) / p̃_max
# 因为 q(x) 是均匀分布，所以归一化常数被约掉
accept_prob = unnormalized_prob(x_candidate, theta) / p_tilde_max
# 以概率 accept_prob 接受该候选点
if np.random.uniform() < accept_prob:
samples.append(x_candidate)
attempts += 1
# 防止无限循环
if attempts > n_samples * 100:
print("警告：拒绝采样效率过低")
break
samples = np.array(samples) # 转为numpy数组
# 对采样样本计算梯度的平均值，即负相
negative_phase = np.mean(gradient_log_unnormalized_prob(samples, theta))
# 最终梯度 = 正相 - 负相
gradient = positive_phase - negative_phase
return gradient
# ============ 实验运行 ============
np.random.seed(42) # 固定随机种子以保证结果可复现
# 生成训练数据：从 N(0, 1) 中采样
# 真实数据分布是一个标准正态分布
train_data = np.random.normal(0, 1, size=500)
# 测试不同 theta 值下的对数似然
theta_values = np.linspace(0.1, 3.0, 50) # 从0.1到3.0均匀取50个值
ll_values = [] # 存储对应的对数似然值
for theta in theta_values:
ll = log_likelihood(train_data, theta)
ll_values.append(ll)
# 打印最优 theta
best_theta = theta_values[np.argmax(ll_values)]
print(f"数值方法找到的最优 theta: {best_theta:.4f}")
print(f"理论最优 theta: 0.5000 (因为真实分布是N(0,1)，能量函数theta*x^2对应高斯精度=theta)")
# 验证梯度计算的正确性
test_theta = 1.0
analytical_grad = analytical_log_likelihood_gradient(train_data, test_theta)
numerical_grad = true_log_likelihood_gradient(train_data, test_theta)
print(f"\n在 theta={test_theta} 处:")
print(f" 解析梯度（采样近似）: {analytical_grad:.6f}")
print(f" 数值梯度（有限差分）: {numerical_grad:.6f}")
print(f" 差异: {abs(analytical_grad - numerical_grad):.6f}")
三、知识点 2：随机最大似然与对比散度
3.1 核心思想
负相需要从 pmodelp_{\text{model}}pmodel 采样，直接采样不可行。两种主要近似方法：
马尔可夫链蒙特卡罗 (MCMC)：构造一个以 pmodelp_{\text{model}}pmodel 为平稳分布的马尔可夫链，运行足够步后样本近似服从 pmodelp_{\text{model}}pmodel。
吉布斯采样（Gibbs Sampling）：对于玻尔兹曼机，逐分量采样：
xi∼pmodel(xi∣x−i)x_i \sim p_{\text{model}}(x_i \mid x_{-i})xi∼pmodel(xi∣x−i)
| 方法 | 思想 | 优点 | 缺点 |
|---|---|---|---|
| 随机最大似然 (SML) | 每步用MCMC生成样本，且持久化马尔可夫链的状态（从上一步继续） | 理论上渐近一致 | 需要足够多的MCMC步数 |
| 对比散度 (CD-k) | 从训练数据初始化MCMC链，只运行 kkk 步（通常 k=1k=1k=1） | 速度快，实践效果好 | 梯度有偏，kkk 越小偏差越大 |
CD-k 的关键步骤：
- 从训练样本 x\mathbf{x}x 出发（初始化MCMC链）
- 运行 kkk 步吉布斯采样得到 x′\mathbf{x}'x′
- 近似梯度：∇θlogp≈∇θlogp~(x)−∇θlogp~(x′)\nabla_\theta \log p \approx \nabla_\theta \log \tilde{p}(\mathbf{x}) - \nabla_\theta \log \tilde{p}(\mathbf{x}')∇θlogp≈∇θlogp~(x)−∇θlogp~(x′)
3.2 案例代码：受限玻尔兹曼机（RBM）上的 CD-1 与 SML
import numpy as np
# ============================================
# 知识点2：对比散度(CD-1) 和 随机最大似然(SML)
# 在受限玻尔兹曼机(RBM)上的实现
# ============================================
class RBM:
"""
受限玻尔兹曼机 (Restricted Boltzmann Machine)
结构：可见层 v (n_visible个二值单元) + 隐藏层 h (n_hidden个二值单元)
能量函数：E(v, h) = -b^T v - c^T h - v^T W h
联合概率：p(v, h) = exp(-E(v, h)) / Z
RBM的特殊结构使得：给定一层，另一层的各单元条件独立
因此可以逐单元进行吉布斯采样
"""
def __init__(self, n_visible, n_hidden, learning_rate=0.01):
"""
初始化RBM
参数:
n_visible: 可见层单元数
n_hidden: 隐藏层单元数
learning_rate: 学习率
"""
self.n_visible = n_visible # 可见层维度
self.n_hidden = n_hidden # 隐藏层维度
self.lr = learning_rate # 学习率
# 初始化权重和偏置
# 权重使用小随机数初始化（打破对称性）
self.W = np.random.randn(n_visible, n_hidden) * 0.01 # 权重矩阵
self.b = np.zeros(n_visible) # 可见层偏置
self.c = np.zeros(n_hidden) # 隐藏层偏置
# SML (Persistent Contrastive Divergence) 用的持久化马尔可夫链状态
self.persistent_chain = None # 初始为空，训练时初始化
def sigmoid(self, x):
"""
sigmoid激活函数 σ(x) = 1 / (1 + exp(-x))
用于计算条件概率 p(h_j=1|v) 和 p(v_i=1|h)
参数:
x: 输入值，numpy数组
返回:
σ(x) 的值，范围在 (0, 1) 之间
"""
# np.clip 防止数值溢出（exp过大导致inf）
x = np.clip(x, -500, 500)
return 1.0 / (1.0 + np.exp(-x))
def sample_bernoulli(self, probs):
"""
从伯努利分布中采样
对于每个元素 p_i，以概率 p_i 生成 1，以概率 1-p_i 生成 0
参数:
probs: 概率值，numpy数组，每个元素在 [0,1] 之间
返回:
与 probs 同形状的 0/1 二值数组
"""
# 生成均匀随机数，如果小于对应概率则为1
return (np.random.random(probs.shape) < probs).astype(np.float64)
def prob_hidden_given_visible(self, v):
"""
计算隐藏层的条件概率 p(h|v)
对于RBM：p(h_j=1|v) = σ(c_j + sum_i W_{ij} * v_i) = σ(c + v^T W)
这是RBM的关键性质：隐藏单元之间条件独立
参数:
v: 可见层状态，形状 (batch_size, n_visible)
返回:
p(h=1|v)，形状 (batch_size, n_hidden)
"""
# v @ W: 矩阵乘法，形状 (batch, n_visible) × (n_visible, n_hidden) → (batch, n_hidden)
# + self.c: 广播加上隐藏层偏置
return self.sigmoid(v @ self.W + self.c)
def prob_visible_given_hidden(self, h):
"""
计算可见层的条件概率 p(v|h)
对于RBM：p(v_i=1|h) = σ(b_i + sum_j W_{ij} * h_j) = σ(b + W h)
参数:
h: 隐藏层状态，形状 (batch_size, n_hidden)
返回:
p(v=1|h)，形状 (batch_size, n_visible)
"""
# h @ W.T: 矩阵乘法，形状 (batch, n_hidden) × (n_hidden, n_visible) → (batch, n_visible)
return self.sigmoid(h @ self.W.T + self.b)
def gibbs_step(self, v):
"""
执行一步吉布斯采样：v → h → v'
1. 给定v采样h：从p(h|v)采样
2. 给定h采样v'：从p(v|h)采样
这是对比散度和SML的核心操作
参数:
v: 当前可见层状态
返回:
v_new: 新的可见层状态
h_prob: 隐藏层概率（用于梯度计算）
h_sample: 隐藏层采样结果
"""
# 步骤1：从可见层计算隐藏层概率并采样
h_prob = self.prob_hidden_given_visible(v) # 计算 p(h|v)
h_sample = self.sample_bernoulli(h_prob) # 从 p(h|v) 采样
# 步骤2：从隐藏层重建可见层概率并采样
v_prob = self.prob_visible_given_hidden(h_sample) # 计算 p(v|h)
v_new = self.sample_bernoulli(v_prob) # 从 p(v|h) 采样
return v_new, h_prob, h_sample
def contrastive_divergence(self, v_data, k=1):
"""
对比散度算法 (CD-k)
核心思想：从训练数据初始化马尔可夫链，只运行k步
近似梯度公式：
∂L/∂W ≈ (1/N) * sum_i [v_i * p(h|v_i)^T] - (1/N) * sum_i [v'_i * p(h|v'_i)^T]
∂L/∂b ≈ mean(v_data) - mean(v_recon)
∂L/∂c ≈ mean(p(h|v_data)) - mean(p(h|v_recon))
参数:
v_data: 训练数据batch，形状 (batch_size, n_visible)
k: 吉布斯采样步数（CD-1中k=1）
返回:
各参数的梯度
"""
batch_size = v_data.shape[0] # 获取batch大小
# === 正相（Positive Phase）===
# 在数据分布上计算统计量
h_prob_data = self.prob_hidden_given_visible(v_data) # p(h|v_data)
# 正相的梯度贡献：数据期望 <v_i * h_j>_data
# 外积求和：v_data.T @ h_prob_data
positive_W = v_data.T @ h_prob_data / batch_size # 形状 (n_visible, n_hidden)
positive_b = np.mean(v_data, axis=0) # 可见层偏置的正梯度
positive_c = np.mean(h_prob_data, axis=0) # 隐藏层偏置的正梯度
# === 负相（Negative Phase）===
# 从数据开始运行k步吉布斯采样
v_current = v_data.copy() # 复制数据作为MCMC链的起点
for step in range(k):
# 每步执行完整的吉布斯采样 v → h → v'
v_current, h_prob_recon, _ = self.gibbs_step(v_current)
# k步后的隐藏层概率（用于负相梯度计算）
h_prob_recon = self.prob_hidden_given_visible(v_current)
# 负相的梯度贡献：模型分布期望 <v_i * h_j>_model
negative_W = v_current.T @ h_prob_recon / batch_size
negative_b = np.mean(v_current, axis=0)
negative_c = np.mean(h_prob_recon, axis=0)
# 计算各参数梯度 = 正相 - 负相
grad_W = positive_W - negative_W
grad_b = positive_b - negative_b
grad_c = positive_c - negative_c
return grad_W, grad_b, grad_c
def train_cd(self, data, epochs=100, batch_size=32, k=1):
"""
使用CD-k训练RBM
参数:
data: 训练数据，形状 (n_samples, n_visible)
epochs: 训练轮数
batch_size: 批大小
k: CD的步数
"""
n_samples = data.shape[0] # 训练样本总数
print(f"开始CD-{k}训练，共{epochs}轮...")
for epoch in range(epochs):
# 打乱数据顺序（随机梯度下降的要求）
indices = np.random.permutation(n_samples)
data_shuffled = data[indices]
total_loss = 0 # 累计重建误差
n_batches = 0 # batch计数
# 按batch遍历数据
for i in range(0, n_samples, batch_size):
# 取出当前batch
batch = data_shuffled[i:i + batch_size]
# 计算CD梯度
grad_W, grad_b, grad_c = self.contrastive_divergence(batch, k=k)
# 更新参数：θ ← θ + lr * ∇L（最大化对数似然等价于梯度上升）
self.W += self.lr * grad_W
self.b += self.lr * grad_b
self.c += self.lr * grad_c
# 计算重建误差（用于监控训练进度，不是梯度的一部分）
# 重建步骤：data → hidden → recon
h_prob = self.prob_hidden_given_visible(batch)
h_sample = self.sample_bernoulli(h_prob)
recon = self.prob_visible_given_hidden(h_sample)
# 使用均方误差衡量重建质量
loss = np.mean((batch - recon) ** 2)
total_loss += loss
n_batches += 1
# 每10轮打印一次训练信息
if (epoch + 1) % 10 == 0:
avg_loss = total_loss / n_batches
print(f" Epoch {epoch+1}/{epochs}, 平均重建误差: {avg_loss:.4f}")
def train_sml(self, data, epochs=100, batch_size=32, k=1):
"""
使用随机最大似然(SML) / 持续对比散度(PCD) 训练RBM
SML与CD的核心区别：
- CD：每个batch都从训练数据重新初始化MCMC链
- SML：维护一个"持久化"的马尔可夫链，每个batch从上一步结束的位置继续
优势：SML能更好地探索模型分布的全局结构
参数:
data: 训练数据
epochs: 训练轮数
batch_size: 批大小
k: 每步的吉布斯采样步数
"""
n_samples = data.shape[0]
# 初始化持久化链：从训练数据中随机选取一个样本
# 这个链会在整个训练过程中持续演化
self.persistent_chain = data[np.random.randint(0, n_samples, batch_size)].copy()
print(f"开始SML训练，共{epochs}轮...")
for epoch in range(epochs):
indices = np.random.permutation(n_samples)
data_shuffled = data[indices]
total_loss = 0
n_batches = 0
for i in range(0, n_samples, batch_size):
batch = data_shuffled[i:i + batch_size]
# === 正相：在数据上计算 ===
h_prob_data = self.prob_hidden_given_visible(batch)
positive_W = batch.T @ h_prob_data / batch.shape[0]
positive_b = np.mean(batch, axis=0)
positive_c = np.mean(h_prob_data, axis=0)
# === 负相：在持久化链上计算 ===
# 关键区别：从persistent_chain继续，而非从数据开始
v_persistent = self.persistent_chain[:batch.shape[0]].copy()
for step in range(k):
# 在持久化链上执行k步吉布斯采样
v_persistent, _, _ = self.gibbs_step(v_persistent)
# 更新持久化链状态（供下一个batch使用）
self.persistent_chain[:batch.shape[0]] = v_persistent.copy()
# 计算负相梯度
h_prob_persistent = self.prob_hidden_given_visible(v_persistent)
negative_W = v_persistent.T @ h_prob_persistent / batch.shape[0]
negative_b = np.mean(v_persistent, axis=0)
negative_c = np.mean(h_prob_persistent, axis=0)
# 梯度更新
self.W += self.lr * (positive_W - negative_W)
self.b += self.lr * (positive_b - negative_b)
self.c += self.lr * (positive_c - negative_c)
# 重建误差监控
h_prob = self.prob_hidden_given_visible(batch)
h_sample = self.sample_bernoulli(h_prob)
recon = self.prob_visible_given_hidden(h_sample)
loss = np.mean((batch - recon) ** 2)
total_loss += loss
n_batches += 1
if (epoch + 1) % 10 == 0:
avg_loss = total_loss / n_batches
print(f" Epoch {epoch+1}/{epochs}, 平均重建误差: {avg_loss:.4f}")
def reconstruct(self, v):
"""
用训练好的RBM重建输入数据
v → p(h|v) → p(v|h) → 重建的v
参数:
v: 输入数据
返回:
重建后的概率值
"""
h_prob = self.prob_hidden_given_visible(v) # 编码
v_recon = self.prob_visible_given_hidden(h_prob) # 解码
return v_recon
# ============ 生成模拟数据 ============
# 生成简单的二值模式数据（类似手写数字的简化版）
def generate_patterns(n_samples):
"""
生成两种二值模式的数据
模式1：左侧为1，右侧为0（左偏模式）
模式2：左侧为0，右侧为1（右偏模式）
加入噪声使其更真实
参数:
n_samples: 样本数
返回:
二值数据矩阵，形状 (n_samples, 16)
"""
n_visible = 16 # 4x4 像素的简化图像
data = np.zeros((n_samples, n_visible))
for i in range(n_samples):
if np.random.random() < 0.5:
# 模式1：左侧高激活
pattern = np.array([1,1,1,1, 1,1,0,0, 0,0,0,0, 0,0,0,0])
else:
# 模式2：右侧高激活
pattern = np.array([0,0,0,0, 0,0,0,0, 0,0,1,1, 1,1,1,1])
# 加入随机翻转噪声（10%概率翻转每个像素）
noise_mask = np.random.random(n_visible) < 0.1
pattern = pattern.copy()
pattern[noise_mask] = 1 - pattern[noise_mask] # 翻转被噪声选中的位
data[i] = pattern
return data
# 生成训练数据
np.random.seed(42)
train_data = generate_patterns(500)
# ============ 训练与对比 ============
# 用CD-1训练
print("=" * 50)
print("方法1: CD-1 训练")
print("=" * 50)
rbm_cd = RBM(n_visible=16, n_hidden=8, learning_rate=0.05)
rbm_cd.train_cd(train_data, epochs=50, batch_size=32, k=1)
# 用SML训练
print("\n" + "=" * 50)
print("方法2: SML 训练")
print("=" * 50)
rbm_sml = RBM(n_visible=16, n_hidden=8, learning_rate=0.05)
rbm_sml.train_sml(train_data, epochs=50, batch_size=32, k=1)
# 对比重建效果
test_sample = train_data[0:1] # 取第一个样本做测试
print(f"\n原始数据: {test_sample[0].astype(int)}")
print(f"CD-1重建: {(rbm_cd.reconstruct(test_sample)[0] > 0.5).astype(int)}")
print(f"SML重建: {(rbm_sml.reconstruct(test_sample)[0] > 0.5).astype(int)}")
四、知识点 3：伪似然（Pseudolikelihood）
3.1 核心思想
直接计算 p(x)p(\mathbf{x})p(x) 需要 ZZZ，但条件概率 p(xi∣x−i)p(x_i \mid \mathbf{x}_{-i})p(xi∣x−i) 中配分函数会被约掉：
p(xi∣x−i)=p(x)∑xi′p(xi′,x−i)=p~(x)∑xi′p~(xi′,x−i)p(x_i \mid \mathbf{x}_{-i}) = \frac{p(\mathbf{x})}{\sum_{x_i'} p(x_i', \mathbf{x}_{-i})} = \frac{\tilde{p}(\mathbf{x})}{\sum_{x_i'} \tilde{p}(x_i', \mathbf{x}_{-i})}p(xi∣x−i)=∑xi′p(xi′,x−i)p(x)=∑xi′p~(xi′,x−i)p~(x)
分母只需对 xix_ixi 的可能取值求和（如二值情况下只有两项），不需要全局配分函数。
伪似然目标函数：
PL(θ)=∑i=1nlogpθ(xi∣x−i)\text{PL}(\theta) = \sum_{i=1}^{n} \log p_\theta(x_i \mid \mathbf{x}_{-i})PL(θ)=i=1∑nlogpθ(xi∣x−i)
- 优点：完全避免配分函数的计算，梯度无偏
- 缺点：假设各维度条件独立（不适用于需要捕获长程依赖的场景），训练和评估不一致（优化的是伪似然而非真正的似然）
3.2 案例代码
import numpy as np
# ============================================
# 知识点3：伪似然 (Pseudolikelihood)
# 在基于能量的模型中实现伪似然训练
# ============================================
class EnergyBasedModel:
"""
基于能量的模型（Ising模型的一种推广）
能量函数：E(x; W, b) = -x^T W x - b^T x
其中 W 是对称权重矩阵（无自连接），b 是偏置
未归一化概率：p̃(x) = exp(-E(x; W, b))
"""
def __init__(self, n_units, learning_rate=0.01):
"""
初始化基于能量的模型
参数:
n_units: 变量维度（单元数）
learning_rate: 学习率
"""
self.n = n_units
self.lr = learning_rate
# 初始化权重为小随机数，对称矩阵
# 乘以0.01防止初始能量差异过大
self.W = np.random.randn(n_units, n_units) * 0.01
# 使权重矩阵对称：W = (W + W^T) / 2
self.W = (self.W + self.W.T) / 2
# 对角线置零（无自连接）
np.fill_diagonal(self.W, 0)
# 偏置初始化为零
self.b = np.zeros(n_units)
def energy(self, x):
"""
计算能量函数 E(x) = -x^T W x - b^T x
参数:
x: 状态向量，形状 (batch_size, n_units)
返回:
每个样本的能量值，形状 (batch_size,)
"""
# x @ self.W @ x.T 的对角线元素等于每个样本的 x^T W x
# 更高效的方式：逐元素乘法再求和
interaction = -0.5 * np.sum((x @ self.W) * x, axis=1) # 交互项
bias = -x @ self.b # 偏置项
return interaction + bias
def unnormalized_log_prob(self, x):
"""
计算 log p̃(x) = -E(x)
参数:
x: 状态向量
返回:
log p̃(x) 的值
"""
return -self.energy(x)
def conditional_prob(self, x, i):
"""
计算条件概率 p(x_i = 1 | x_{-i})
关键推导：
p(x_i=1|x_{-i}) / p(x_i=0|x_{-i}) = exp(-E(x_i=1,x_{-i}) + E(x_i=0,x_{-i}))
令 ΔE = E(x_i=1, x_{-i}) - E(x_i=0, x_{-i})
则 p(x_i=1|x_{-i}) = σ(-ΔE) = 1/(1+exp(ΔE))
对于Ising模型：
ΔE = -(sum_j W_{ij} * x_j + b_i)
所以 p(x_i=1|x_{-i}) = σ(sum_j W_{ij} * x_j + b_i)
参数:
x: 当前状态，形状 (batch_size, n_units)
i: 要计算条件概率的变量索引
返回:
p(x_i=1|x_{-i})，形状 (batch_size,)
"""
# 计算 x_i 的局部场（所有邻居对i的影响之和）
# h_i = sum_j W_{ij} * x_j + b_i
# 注意：W_{ii}=0，所以不需要特殊处理
local_field = x @ self.W[:, i] + self.b[i]
# 条件概率 = sigmoid(local_field)
return 1.0 / (1.0 + np.exp(-np.clip(local_field, -500, 500)))
def pseudolikelihood(self, x):
"""
计算伪似然：PL(x) = sum_i log p(x_i | x_{-i})
对于二值变量，x_i ∈ {0, 1}：
log p(x_i | x_{-i}) = x_i * log(p_i) + (1-x_i) * log(1-p_i)
其中 p_i = p(x_i=1 | x_{-i})
这就是对每个维度分别计算二元交叉熵
参数:
x: 数据，形状 (batch_size, n_units)
返回:
每个样本的伪似然值，形状 (batch_size,)
"""
batch_size = x.shape[0]
pl = np.zeros(batch_size) # 初始化伪似然累加器
for i in range(self.n):
# 计算第i个变量的条件概率
p_i = self.conditional_prob(x, i)
# 二元交叉熵：x_i * log(p_i) + (1-x_i) * log(1-p_i)
# 添加小常数1e-10防止log(0)
log_p = x[:, i] * np.log(p_i + 1e-10) + \
(1 - x[:, i]) * np.log(1 - p_i + 1e-10)
pl += log_p # 累加到伪似然
return pl
def pseudolikelihood_gradient(self, x):
"""
计算伪似然关于参数 W 和 b 的梯度
∂PL/∂W_{ij} = sum_k [∂/∂W_{ij} log p(x_k | x_{-k})]
对于 k = i 或 k = j，梯度不为零（因为 W_{ij} 只影响 x_i 和 x_j 的条件概率）
更简洁的实现：对每个维度 i，
∂PL/∂W_{ij} = (x_i - p_i) * x_j （对于所有 j ≠ i）
∂PL/∂b_i = x_i - p_i
参数:
x: 数据，形状 (batch_size, n_units)
返回:
grad_W: W的梯度
grad_b: b的梯度
"""
batch_size = x.shape[0]
grad_W = np.zeros_like(self.W) # 权重梯度，形状 (n, n)
grad_b = np.zeros_like(self.b) # 偏置梯度，形状 (n,)
for i in range(self.n):
# 计算 p_i = p(x_i=1|x_{-i})
p_i = self.conditional_prob(x, i) # 形状 (batch_size,)
# "残差"：实际值 - 预测概率
residual = x[:, i] - p_i # 形状 (batch_size,)
# 对偏置的梯度：∂PL/∂b_i = (1/N) * sum(residual)
grad_b[i] += np.mean(residual)
# 对权重的梯度：∂PL/∂W_{ij} = (1/N) * sum(residual * x_j)
# 这是一个向量外积的期望
# residual[:, None] 将其变为列向量，x 是 (batch, n)
grad_W[i, :] += np.mean(residual[:, None] * x, axis=0)
# 使梯度矩阵对称（因为W是对称的）
grad_W = (grad_W + grad_W.T) / 2
# 对角线置零（与W的约束一致）
np.fill_diagonal(grad_W, 0)
return grad_W, grad_b
def train_pseudolikelihood(self, data, epochs=100, batch_size=32):
"""
使用伪似然训练基于能量的模型
参数:
data: 训练数据，形状 (n_samples, n_units)
epochs: 训练轮数
batch_size: 批大小
"""
n_samples = data.shape[0]
print("开始伪似然训练...")
for epoch in range(epochs):
# 打乱数据
indices = np.random.permutation(n_samples)
data_shuffled = data[indices]
total_pl = 0 # 累计伪似然
n_batches = 0
for i in range(0, n_samples, batch_size):
batch = data_shuffled[i:i + batch_size]
# 计算伪似然梯度
grad_W, grad_b = self.pseudolikelihood_gradient(batch)
# 梯度上升更新（最大化伪似然）
self.W += self.lr * grad_W
self.b += self.lr * grad_b
# 记录伪似然值（用于监控）
pl = np.mean(self.pseudolikelihood(batch))
total_pl += pl
n_batches += 1
if (epoch + 1) % 10 == 0:
avg_pl = total_pl / n_batches
print(f" Epoch {epoch+1}/{epochs}, 平均伪似然: {avg_pl:.4f}")
# ============ 实验 ============
np.random.seed(42)
# 生成具有空间相关性的二值数据
def generate_correlated_data(n_samples, n_units):
"""
生成具有相邻单元相关性的二值数据
使用链式依赖结构：x_i 与 x_{i-1} 正相关
参数:
n_samples: 样本数
n_units: 变量维度
返回:
二值数据矩阵
"""
data = np.zeros((n_samples, n_units))
# 第一个单元随机
data[:, 0] = (np.random.random(n_samples) > 0.5).astype(float)
for i in range(1, n_units):
# 后续单元以80%概率与前一个单元相同，20%概率相反
# 这创建了相邻单元之间的正相关
same_as_prev = np.random.random(n_samples) < 0.8
data[:, i] = np.where(same_as_prev, data[:, i-1], 1 - data[:, i-1])
return data
# 生成训练数据
train_data = generate_correlated_data(1000, n_units=10)
# 训练模型
model = EnergyBasedModel(n_units=10, learning_rate=0.05)
model.train_pseudolikelihood(train_data, epochs=50, batch_size=64)
# 查看学到的权重矩阵
print("\n学到的权重矩阵W（应显示相邻单元间的正权重）:")
print(np.round(model.W[:5, :5], 3)) # 打印前5x5的子矩阵
五、知识点 4：得分匹配与比率匹配
5.1 得分匹配（Score Matching）
得分（Score） 定义为对数概率密度的梯度：
s(x)=∇xlogpmodel(x)s(\mathbf{x}) = \nabla_\mathbf{x} \log p_{\text{model}}(\mathbf{x})s(x)=∇xlogpmodel(x)
关键观察：得分的计算不依赖配分函数，因为 ∇xlogZ=0\nabla_\mathbf{x} \log Z = 0∇xlogZ=0：
s(x)=∇xlogp~(x)s(\mathbf{x}) = \nabla_\mathbf{x} \log \tilde{p}(\mathbf{x})s(x)=∇xlogp~(x)
得分匹配目标：令模型得分匹配数据得分。目标函数：
J(θ)=12Epdata[∥smodel(x;θ)−sdata(x)∥2]J(\theta) = \frac{1}{2} \mathbb{E}_{p_{\text{data}}} \left[ \| s_{\text{model}}(\mathbf{x}; \theta) - s_{\text{data}}(\mathbf{x}) \|^2 \right]J(θ)=21Epdata[∥smodel(x;θ)−sdata(x)∥2]
通过分部积分，可以避免直接计算 sdatas_{\text{data}}sdata：
J(θ)=Epdata[tr(∇xsmodel(x))+12∥smodel(x)∥2]J(\theta) = \mathbb{E}_{p_{\text{data}}} \left[ \text{tr}(\nabla_\mathbf{x} s_{\text{model}}(\mathbf{x})) + \frac{1}{2} \| s_{\text{model}}(\mathbf{x}) \|^2 \right]J(θ)=Epdata[tr(∇xsmodel(x))+21∥smodel(x)∥2]
5.2 比率匹配（Ratio Matching）
适用于离散数据。目标是匹配概率比率 p(x)p(¬ix)\frac{p(\mathbf{x})}{p(\neg_i \mathbf{x})}p(¬ix)p(x)，其中 ¬ix\neg_i \mathbf{x}¬ix 是翻转第 iii 位后的向量。
J(θ)=∑i=1nEpdata[(11+exp(si(x)−si(¬ix)))2]J(\theta) = \sum_{i=1}^{n} \mathbb{E}_{p_{\text{data}}} \left[ \left( \frac{1}{1 + \exp(s_i(\mathbf{x}) - s_i(\neg_i \mathbf{x}))} \right)^2 \right]J(θ)=i=1∑nEpdata[(1+exp(si(x)−si(¬ix))1)2]
其中 si(x)=logp~(¬ix)−logp~(x)s_i(\mathbf{x}) = \log \tilde{p}(\neg_i \mathbf{x}) - \log \tilde{p}(\mathbf{x})si(x)=logp~(¬ix)−logp~(x)。
5.3 案例代码
import numpy as np
# ============================================
# 知识点4：得分匹配 (Score Matching) 与 比率匹配 (Ratio Matching)
# ============================================
class GaussianEBM:
"""
基于能量的连续模型（高斯形式）
p̃(x) = exp(-0.5 * x^T A x - b^T x)
其中 A 是正定矩阵，b 是偏置向量
这对应一个多元高斯分布
"""
def __init__(self, n_dim, learning_rate=0.001):
"""
初始化模型
参数:
n_dim: 数据维度
learning_rate: 学习率
"""
self.n_dim = n_dim
self.lr = learning_rate
# A: 正定矩阵参数（初始化为单位阵）
self.A = np.eye(n_dim)
# b: 偏置向量
self.b = np.zeros(n_dim)
def score(self, x):
"""
计算模型得分 s(x) = ∇_x log p̃(x)
推导：
log p̃(x) = -0.5 * x^T A x - b^T x
∇_x log p̃(x) = -A x - b
参数:
x: 数据点，形状 (batch_size, n_dim)
返回:
得分向量，形状 (batch_size, n_dim)
"""
# -(x @ A^T + b)，由于A是对称的，A^T = A
return -(x @ self.A + self.b)
def score_jacobian(self, x):
"""
计算得分的雅可比矩阵 ∇_x s(x) = ∂s_i/∂x_j
对于 s(x) = -Ax - b
雅可比矩阵 ∂s/∂x = -A（与x无关）
但在一般情况下雅可比矩阵依赖x
这里返回 -A 的 batch 广播版本
参数:
x: 数据点，形状 (batch_size, n_dim)
返回:
雅可比张量，形状 (batch_size, n_dim, n_dim)
"""
batch_size = x.shape[0]
# 对于线性能量模型，雅可比是常数 -A
return np.tile(-self.A, (batch_size, 1, 1))
def score_matching_loss(self, x):
"""
计算得分匹配的损失函数
J(θ) = E[tr(∇_x s(x)) + 0.5 * ||s(x)||^2]
这是经过分部积分后的形式，不需要知道数据得分
参数:
x: 数据batch，形状 (batch_size, n_dim)
返回:
损失值（标量）
"""
batch_size = x.shape[0]
# 计算模型得分
s = self.score(x) # 形状 (batch_size, n_dim)
# 计算得分雅可比矩阵
jac = self.score_jacobian(x) # 形状 (batch_size, n_dim, n_dim)
# 项1: tr(∇_x s(x)) = 对每个样本计算雅可比矩阵的迹
# np.trace 沿最后两个轴计算迹
trace_term = np.trace(jac, axis1=1, axis2=2) # 形状 (batch_size,)
# 项2: 0.5 * ||s(x)||^2
norm_term = 0.5 * np.sum(s ** 2, axis=1) # 形状 (batch_size,)
# 总损失 = 两项之和的均值
loss = np.mean(trace_term + norm_term)
return loss
def score_matching_gradient(self, x):
"""
计算得分匹配损失对参数A和b的梯度
对于线性能量模型有解析梯度：
∂J/∂A = -I + E[(Ax+b)(Ax+b)^T] + E[A x x^T]（简化形式）
这里使用自动微分思想手动推导
参数:
x: 数据batch
返回:
grad_A, grad_b
"""
batch_size = x.shape[0]
# 得分 s(x) = -Ax - b
s = self.score(x) # (batch_size, n_dim)
# 雅可比的迹 = -tr(A) = -sum(A的对角元素)
trace_A = np.trace(self.A)
# ∂J/∂A:
# 项1: ∂tr(∇_x s)/∂A = -I（因为tr(-A)对A的梯度是-I）
grad_A_trace = -np.eye(self.n_dim)
# 项2: ∂(0.5||s||^2)/∂A = E[s * (-x^T)] 因为 s = -Ax-b
# = -E[s @ x^T] 但还需要考虑对称性
grad_A_norm = s.T @ x / batch_size # (n_dim, n_dim)
grad_A = grad_A_trace + grad_A_norm
# ∂J/∂b:
# 项1: tr(∇_x s)/∂b = 0（迹不依赖b）
# 项2: ∂(0.5||s||^2)/∂b = E[s * (-1)] = -E[s]
grad_b = -np.mean(s, axis=0)
return grad_A, grad_b
def train_score_matching(self, data, epochs=200, batch_size=64):
"""
使用得分匹配训练模型
参数:
data: 训练数据
epochs: 训练轮数
batch_size: 批大小
"""
n_samples = data.shape[0]
print("开始得分匹配训练...")
for epoch in range(epochs):
indices = np.random.permutation(n_samples)
data_shuffled = data[indices]
total_loss = 0
n_batches = 0
for i in range(0, n_samples, batch_size):
batch = data_shuffled[i:i + batch_size]
# 计算梯度
grad_A, grad_b = self.score_matching_gradient(batch)
# 梯度下降（最小化损失）
self.A -= self.lr * grad_A
self.b -= self.lr * grad_b
# 确保A保持对称
self.A = (self.A + self.A.T) / 2
loss = self.score_matching_loss(batch)
total_loss += loss
n_batches += 1
if (epoch + 1) % 50 == 0:
print(f" Epoch {epoch+1}/{epochs}, 损失: {total_loss/n_batches:.4f}")
class DiscreteRatioMatching:
"""
比率匹配模型，用于离散二值数据
比率匹配的核心思想：
匹配 log p̃(¬_i x) / p̃(x) 的比率，其中 ¬_i x 是翻转第i位的向量
"""
def __init__(self, n_dim, learning_rate=0.01):
"""
初始化比率匹配模型
参数:
n_dim: 数据维度
learning_rate: 学习率
"""
self.n_dim = n_dim
self.lr = learning_rate
# 使用简单的线性模型: log p̃(x) = b^T x + x^T W x
self.W = np.random.randn(n_dim, n_dim) * 0.01
self.W = (self.W + self.W.T) / 2 # 对称
np.fill_diagonal(self.W, 0)
self.b = np.zeros(n_dim)
def log_unnormalized_prob(self, x):
"""
计算 log p̃(x) = x^T W x + b^T x
参数:
x: 二值数据，形状 (batch_size, n_dim)
返回:
log p̃(x) 的值，形状 (batch_size,)
"""
# 交互项 + 偏置项
return np.sum((x @ self.W) * x, axis=1) + x @ self.b
def ratio_matching_loss(self, x):
"""
计算比率匹配损失
对于每个样本和每个维度 i：
1. 构造 ¬_i x（翻转第i位）
2. 计算 s_i = log p̃(¬_i x) - log p̃(x)
3. 损失项 = (1/(1+exp(s_i)))^2
直觉：如果模型正确，则翻转一位后的概率应与原概率有一定比率关系
参数:
x: 数据batch，形状 (batch_size, n_dim)
返回:
损失值
"""
batch_size = x.shape[0]
log_p_x = self.log_unnormalized_prob(x) # log p̃(x)
total_loss = 0 # 累计损失
for i in range(self.n_dim):
# 构造翻转第i位的样本 ¬_i x
x_flipped = x.copy() # 复制原数据
x_flipped[:, i] = 1 - x[:, i] # 翻转第i位（0→1, 1→0）
# 计算 log p̃(¬_i x)
log_p_flipped = self.log_unnormalized_prob(x_flipped)
# 计算 s_i = log p̃(¬_i x) - log p̃(x)
s_i = log_p_flipped - log_p_x
# σ(-s_i)^2 = (1/(1+exp(s_i)))^2
# 使用 sigmoid 函数
sigma_neg = 1.0 / (1.0 + np.exp(np.clip(s_i, -500, 500)))
loss_i = sigma_neg ** 2
total_loss += np.mean(loss_i)
return total_loss / self.n_dim # 对所有维度取平均
def ratio_matching_gradient(self, x):
"""
计算比率匹配损失的梯度（数值近似版，更稳健）
使用有限差分法
参数:
x: 数据batch
返回:
grad_W, grad_b
"""
epsilon = 1e-5
grad_W = np.zeros_like(self.W)
grad_b = np.zeros_like(self.b)
# 当前损失
loss_orig = self.ratio_matching_loss(x)
# 对b的梯度
for i in range(self.n_dim):
self.b[i] += epsilon
loss_plus = self.ratio_matching_loss(x)
self.b[i] -= epsilon
grad_b[i] = (loss_plus - loss_orig) / epsilon
# 对W的梯度（只计算上三角，然后对称化）
for i in range(self.n_dim):
for j in range(i + 1, self.n_dim):
self.W[i, j] += epsilon
self.W[j, i] += epsilon
loss_plus = self.ratio_matching_loss(x)
self.W[i, j] -= epsilon
self.W[j, i] -= epsilon
grad_W[i, j] = (loss_plus - loss_orig) / epsilon
grad_W[j, i] = grad_W[i, j] # 对称
return grad_W, grad_b
def train(self, data, epochs=50, batch_size=32):
"""
使用比率匹配训练模型
参数:
data: 训练数据
epochs: 训练轮数
batch_size: 批大小
"""
n_samples = data.shape[0]
print("开始比率匹配训练...")
for epoch in range(epochs):
indices = np.random.permutation(n_samples)
data_shuffled = data[indices]
total_loss = 0
n_batches = 0
for i in range(0, n_samples, batch_size):
batch = data_shuffled[i:i + batch_size]
grad_W, grad_b = self.ratio_matching_gradient(batch)
# 梯度下降
self.W -= self.lr * grad_W
self.b -= self.lr * grad_b
self.W = (self.W + self.W.T) / 2
np.fill_diagonal(self.W, 0)
loss = self.ratio_matching_loss(batch)
total_loss += loss
n_batches += 1
if (epoch + 1) % 10 == 0:
print(f" Epoch {epoch+1}/{epochs}, 损失: {total_loss/n_batches:.4f}")
# ============ 得分匹配实验 ============
np.random.seed(42)
# 生成二维高斯数据（真实分布）
# 真实协均值矩阵 Sigma = [[2, 0.5], [0.5, 1]]
true_mean = np.array([1.0, -0.5])
true_cov = np.array([[2.0, 0.5], [0.5, 1.0]])
train_data = np.random.multivariate_normal(true_mean, true_cov, size=1000)
# 训练得分匹配模型
model_sm = GaussianEBM(n_dim=2, learning_rate=0.001)
model_sm.train_score_matching(train_data, epochs=200, batch_size=64)
# 比较学到的参数与真实参数
# 真实精度矩阵 A_true = Sigma^{-1}
true_precision = np.linalg.inv(true_cov)
print(f"\n真实精度矩阵 A:\n{np.round(true_precision, 4)}")
print(f"学到的矩阵 A:\n{np.round(model_sm.A, 4)}")
# ============ 比率匹配实验 ============
# 生成简单二值数据
binary_data = np.zeros((500, 6))
for i in range(500):
if np.random.random() < 0.5:
binary_data[i] = [1,1,0,0,0,0]
else:
binary_data[i] = [0,0,1,1,0,0]
model_rm = DiscreteRatioMatching(n_dim=6, learning_rate=0.01)
model_rm.train(binary_data, epochs=30, batch_size=32)
print(f"\n学到的偏置 b:\n{np.round(model_rm.b, 3)}")
六、知识点 5：去噪得分匹配（Denoising Score Matching）
6.1 核心思想
得分匹配需要计算 ∇xs(x)\nabla_\mathbf{x} s(\mathbf{x})∇xs(x)（雅可比矩阵），这在高维空间计算代价很高。
去噪得分匹配通过向数据添加噪声来简化问题：
- 对数据 x\mathbf{x}x 添加噪声：x~=x+ϵ,ϵ∼N(0,σ2I)\tilde{\mathbf{x}} = \mathbf{x} + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma^2 I)x~=x+ϵ,ϵ∼N(0,σ2I)
- 目标变为匹配去噪后的得分：
J(θ)=12EpdataEqσ(x~∣x)[∥sθ(x~)−∇x~logqσ(x~∣x)∥2]J(\theta) = \frac{1}{2} \mathbb{E}_{p_{\text{data}}} \mathbb{E}_{q_\sigma(\tilde{\mathbf{x}} | \mathbf{x})} \left[ \| s_\theta(\tilde{\mathbf{x}}) - \nabla_{\tilde{\mathbf{x}}} \log q_\sigma(\tilde{\mathbf{x}} | \mathbf{x}) \|^2 \right]J(θ)=21EpdataEqσ(x~∣x)[∥sθ(x~)−∇x~logqσ(x~∣x)∥2]
其中 ∇x~logqσ(x~∣x)=−x~−xσ2\nabla_{\tilde{\mathbf{x}}} \log q_\sigma(\tilde{\mathbf{x}} | \mathbf{x}) = -\frac{\tilde{\mathbf{x}} - \mathbf{x}}{\sigma^2}∇x~logqσ(x~∣x)=−σ2x~−x。
优势：
- 不需要计算得分的雅可比矩阵
- 概率密度 q(x~∣x)q(\tilde{x}|x)q(x~∣x) 已知且简单（高斯噪声）
- 在高维空间中更实用
6.2 案例代码
import numpy as np
# ============================================
# 知识点5：去噪得分匹配 (Denoising Score Matching)
# ============================================
class DenoisingScoreMatching:
"""
去噪得分匹配模型
使用神经网络参数化得分函数 s_θ(x) ≈ ∇_x log p_data(x)
核心流程：
1. 从训练数据采样 x
2. 给 x 加噪声得到 x̃ = x + σ*ε
3. 目标得分 = -(x̃ - x)/σ² (已知的高斯噪声得分)
4. 训练网络预测这个目标得分
"""
def __init__(self, n_dim, hidden_dim=64, noise_std=0.1, learning_rate=0.001):
"""
初始化去噪得分匹配模型
参数:
n_dim: 数据维度
hidden_dim: 隐藏层维度
noise_std: 噪声标准差 σ（核心超参数）
learning_rate: 学习率
"""
self.n_dim = n_dim
self.hidden_dim = hidden_dim
self.sigma = noise_std # 噪声标准差
self.lr = learning_rate
# 简单的两层神经网络参数
# 第一层：输入层 → 隐藏层
self.W1 = np.random.randn(n_dim, hidden_dim) * np.sqrt(2.0 / n_dim)
self.b1 = np.zeros(hidden_dim)
# 第二层：隐藏层 → 输出层（输出得分向量）
self.W2 = np.random.randn(hidden_dim, n_dim) * np.sqrt(2.0 / hidden_dim)
self.b2 = np.zeros(n_dim)
def relu(self, x):
"""
ReLU激活函数: f(x) = max(0, x)
参数:
x: 输入
返回:
max(0, x)
"""
return np.maximum(0, x)
def relu_derivative(self, x):
"""
ReLU的导数: f'(x) = 1 if x > 0, else 0
参数:
x: 输入
返回:
导数值
"""
return (x > 0).astype(float)
def forward(self, x):
"""
前向传播，计算得分预测 s_θ(x)
架构: x → [W1, b1] → ReLU → [W2, b2] → s_θ(x)
参数:
x: 输入数据，形状 (batch_size, n_dim)
返回:
score_pred: 预测的得分，形状 (batch_size, n_dim)
cache: 中间结果（用于反向传播）
"""
# 第一层的线性变换
z1 = x @ self.W1 + self.b1 # 形状 (batch, hidden_dim)
# ReLU激活
a1 = self.relu(z1) # 形状 (batch, hidden_dim)
# 第二层的线性变换（输出层无激活函数）
score_pred = a1 @ self.W2 + self.b2 # 形状 (batch, n_dim)
# 缓存中间结果用于反向传播
cache = {'x': x, 'z1': z1, 'a1': a1}
return score_pred, cache
def add_noise(self, x):
"""
给数据添加高斯噪声
x̃ = x + σ * ε, ε ~ N(0, I)
参数:
x: 原始数据，形状 (batch_size, n_dim)
返回:
x_noisy: 加噪后数据
noise: 添加的噪声（用于计算目标得分）
"""
noise = np.random.randn(*x.shape) # 生成标准正态噪声
x_noisy = x + self.sigma * noise # 加噪
return x_noisy, noise
def target_score(self, x, x_noisy):
"""
计算目标得分（噪声条件概率的梯度）
对于 q(x̃|x) = N(x̃; x, σ²I)
∇_{x̃} log q(x̃|x) = -(x̃ - x)/σ²
这是已知的解析表达式，不需要估计
参数:
x: 原始数据
x_noisy: 加噪后数据
返回:
目标得分向量
"""
return -(x_noisy - x) / (self.sigma ** 2)
def dsm_loss(self, x):
"""
计算去噪得分匹配损失
L(θ) = (1/2) * E[||s_θ(x̃) - ∇log q(x̃|x)||²]
直觉：让网络在噪声输入上预测"指向真实数据的方向"
参数:
x: 原始数据batch
返回:
loss: 损失值
cache: 中间结果
"""
# 步骤1: 添加噪声
x_noisy, noise = self.add_noise(x)
# 步骤2: 神经网络预测得分
score_pred, cache = self.forward(x_noisy)
# 步骤3: 计算目标得分
score_target = self.target_score(x, x_noisy)
# 步骤4: MSE损失
diff = score_pred - score_target # 形状 (batch, n_dim)
loss = 0.5 * np.mean(np.sum(diff ** 2, axis=1)) # 对维度求和，对batch求平均
cache['x_noisy'] = x_noisy
cache['score_pred'] = score_pred
cache['score_target'] = score_target
cache['diff'] = diff
return loss, cache
def backward(self, cache):
"""
反向传播计算梯度
参数:
cache: 前向传播和损失计算的中间结果
返回:
各参数的梯度
"""
batch_size = cache['x'].shape[0]
diff = cache['diff'] # score_pred - score_target，形状 (batch, n_dim)
# ∂L/∂score_pred = diff（MSE的梯度）
d_score = diff / batch_size # 形状 (batch, n_dim)
# 第二层梯度
# score = a1 @ W2 + b2
# ∂L/∂W2 = a1^T @ d_score
grad_W2 = cache['a1'].T @ d_score # 形状 (hidden_dim, n_dim)
grad_b2 = np.sum(d_score, axis=0) # 形状 (n_dim,)
# 传播到隐藏层
# ∂L/∂a1 = d_score @ W2^T
d_a1 = d_score @ self.W2.T # 形状 (batch, hidden_dim)
# ReLU的梯度
d_z1 = d_a1 * self.relu_derivative(cache['z1']) # 形状 (batch, hidden_dim)
# 第一层梯度
# ∂L/∂W1 = x^T @ d_z1
# 注意：这里的x实际是x_noisy（网络的输入是加噪数据）
grad_W1 = cache['x_noisy'].T @ d_z1 # 形状 (n_dim, hidden_dim)
grad_b1 = np.sum(d_z1, axis=0) # 形状 (hidden_dim,)
return grad_W1, grad_b1, grad_W2, grad_b2
def train(self, data, epochs=200, batch_size=64):
"""
训练去噪得分匹配模型
参数:
data: 训练数据
epochs: 训练轮数
batch_size: 批大小
"""
n_samples = data.shape[0]
print(f"开始去噪得分匹配训练 (σ={self.sigma})...")
for epoch in range(epochs):
indices = np.random.permutation(n_samples)
data_shuffled = data[indices]
total_loss = 0
n_batches = 0
for i in range(0, n_samples, batch_size):
batch = data_shuffled[i:i + batch_size]
# 前向传播 + 损失计算
loss, cache = self.dsm_loss(batch)
# 反向传播
grad_W1, grad_b1, grad_W2, grad_b2 = self.backward(cache)
# 梯度下降更新
self.W1 -= self.lr * grad_W1
self.b1 -= self.lr * grad_b1
self.W2 -= self.lr * grad_W2
self.b2 -= self.lr * grad_b2
total_loss += loss
n_batches += 1
if (epoch + 1) % 50 == 0:
print(f" Epoch {epoch+1}/{epochs}, 损失: {total_loss/n_batches:.4f}")
def predict_score(self, x):
"""
用训练好的模型预测得分
参数:
x: 数据点
返回:
预测的得分 ∇_x log p(x)
"""
score, _ = self.forward(x)
return score
def langevin_dynamics_sample(self, n_samples, n_steps=100, step_size=0.01):
"""
使用朗之万动力学从模型中采样
x_{t+1} = x_t + (step_size/2) * s_θ(x_t) + sqrt(step_size) * ε
这是基于得分的生成方法
参数:
n_samples: 生成的样本数
n_steps: 朗之万动力学的步数
step_size: 步长
返回:
生成的样本
"""
# 从先验分布（标准正态）初始化
x = np.random.randn(n_samples, self.n_dim) * 3.0
for step in range(n_steps):
# 计算当前得分
score = self.predict_score(x)
# 朗之万更新
noise = np.random.randn(*x.shape)
x = x + (step_size / 2) * score + np.sqrt(step_size) * noise
# 可选：逐步退火步长（改善采样质量）
if step == n_steps // 2:
step_size *= 0.5 # 后半程减半步长
return x
# ============ 实验 ============
np.random.seed(42)
# 生成混合高斯数据
def generate_mog_data(n_samples):
"""
生成双峰高斯混合数据
混合权重各0.5，两个分量的均值不同
"""
data = []
for _ in range(n_samples):
if np.random.random() < 0.5:
# 分量1：均值 [-2, 0]
sample = np.random.randn(2) * 0.5 + np.array([-2, 0])
else:
# 分量2：均值 [2, 0]
sample = np.random.randn(2) * 0.5 + np.array([2, 0])
data.append(sample)
return np.array(data)
train_data = generate_mog_data(2000)
# 训练模型（测试不同噪声标准差）
for sigma in [0.05, 0.1, 0.5]:
print(f"\n{'='*50}")
print(f"噪声标准差 σ = {sigma}")
print(f"{'='*50}")
model = DenoisingScoreMatching(
n_dim=2,
hidden_dim=32,
noise_std=sigma,
learning_rate=0.001
)
model.train(train_data, epochs=100, batch_size=64)
# 采样生成
samples = model.langevin_dynamics_sample(n_samples=500, n_steps=200, step_size=0.01)
print(f" 生成样本均值: {np.mean(samples, axis=0).round(3)}")
print(f" 生成样本标准差: {np.std(samples, axis=0).round(3)}")
七、知识点 6：噪声对比估计（NCE）
7.1 核心思想
NCE 将密度估计问题转化为一个二分类问题：
- 定义带参数的未归一化模型 pmodel(x;θ)=p~model(x;θ)p_\text{model}(\mathbf{x}; \theta) = \tilde{p}_\text{model}(\mathbf{x}; \theta)pmodel(x;θ)=p~model(x;θ)（ZZZ 被吸收进参数）
- 引入一个已知的噪声分布 pnoise(x)p_\text{noise}(\mathbf{x})pnoise(x)
- 生成"混合"数据：真实数据（标签 y=1y=1y=1）和噪声样本（标签 y=0y=0y=0）
- 训练分类器区分两者
关键公式：
p(y=1∣x)=pmodel(x)pmodel(x)+k⋅pnoise(x)p(y=1|\mathbf{x}) = \frac{p_\text{model}(\mathbf{x})}{p_\text{model}(\mathbf{x}) + k \cdot p_\text{noise}(\mathbf{x})}p(y=1∣x)=pmodel(x)+k⋅pnoise(x)pmodel(x)
其中 kkk 是每个真实样本对应的噪声样本数。
NCE 目标函数（每个真实样本的贡献）：
Ji=logp(y=1∣xi)+k⋅Ex′∼pnoise[logp(y=0∣x′)]J_i = \log p(y=1 | \mathbf{x}_i) + k \cdot \mathbb{E}_{\mathbf{x}' \sim p_\text{noise}} [\log p(y=0 | \mathbf{x}')]Ji=logp(y=1∣xi)+k⋅Ex′∼pnoise[logp(y=0∣x′)]
关键优势：如果将模型写为 pmodel(x)=1Z(θ)p~(x;θ)p_\text{model}(\mathbf{x}) = \frac{1}{Z(\theta)} \tilde{p}(\mathbf{x}; \theta)pmodel(x)=Z(θ)1p~(x;θ)，NCE 的梯度会驱使 Z(θ)→1Z(\theta) \to 1Z(θ)→1，因此 ZZZ 也被隐式学习了。
7.2 案例代码
import numpy as np
# ============================================
# 知识点6：噪声对比估计 (Noise Contrastive Estimation, NCE)
# ============================================
class NCEModel:
"""
噪声对比估计模型
目标：估计一个概率分布 p_model(x) = p̃(x; θ) / Z(θ)
方法：通过区分真实数据和噪声样本来间接学习 θ 和 Z
模型假设：
- p̃(x; θ) 使用指数族形式
- 噪声分布 p_noise 选用已知的简单分布（如高斯）
"""
def __init__(self, n_dim, noise_mean=None, noise_cov=None,
learning_rate=0.01, n_noise_samples=10):
"""
初始化NCE模型
参数:
n_dim: 数据维度
noise_mean: 噪声分布的均值（默认零向量）
noise_cov: 噪声分布的协方差（默认单位阵）
learning_rate: 学习率
n_noise_samples: 每个真实样本对应的噪声样本数 k
"""
self.n_dim = n_dim
self.lr = learning_rate
self.k = n_noise_samples # 噪声样本倍数
# 模型参数：使用二次型能量模型
# log p̃(x) = -0.5 * x^T A x - b^T x + c
# A 必须正定以保证概率可积
self.A = np.eye(n_dim) * 0.5 # 初始化为半倍单位阵
self.b = np.zeros(n_dim) # 偏置
self.log_Z = 0.0 # log配分函数（NCE会自动学习它！）
# 噪声分布参数
if noise_mean is None:
self.noise_mean = np.zeros(n_dim)
else:
self.noise_mean = noise_mean
if noise_cov is None:
self.noise_cov = np.eye(n_dim) * 4.0 # 较大方差的噪声
else:
self.noise_cov = noise_cov
# 预计算噪声分布精度矩阵和归一化常数
self.noise_precision = np.linalg.inv(self.noise_cov)
self.noise_log_norm = -0.5 * (n_dim * np.log(2 * np.pi) +
np.log(np.linalg.det(self.noise_cov)))
def log_unnormalized_model(self, x):
"""
计算模型的未归一化对数概率
log p̃(x; θ) = -0.5 * x^T A x - b^T x
参数:
x: 数据，形状 (batch_size, n_dim)
返回:
log p̃(x) 的值，形状 (batch_size,)
"""
return -0.5 * np.sum((x @ self.A) * x, axis=1) - x @ self.b
def log_model_prob(self, x):
"""
计算模型的（近似）对数概率
log p_model(x) ≈ log p̃(x) - log Z
在NCE中，log Z 也被学习
参数:
x: 数据
返回:
log p_model(x)
"""
return self.log_unnormalized_model(x) - self.log_Z
def log_noise_prob(self, x):
"""
计算噪声分布的对数概率
log p_noise(x) = -0.5 * (x-μ)^T Σ^{-1} (x-μ) + log_norm
参数:
x: 数据，形状 (batch_size, n_dim)
返回:
log p_noise(x)，形状 (batch_size,)
"""
diff = x - self.noise_mean # 偏差
# 二次型
quad = np.sum((diff @ self.noise_precision) * diff, axis=1)
return -0.5 * quad + self.noise_log_norm
def nce_probability(self, x, log_p_model, log_p_noise):
"""
计算NCE分类概率 p(y=1|x)
p(y=1|x) = p_model(x) / (p_model(x) + k * p_noise(x))
在对数空间中实现以避免数值溢出：
log_ratio = log p_model(x) - log p_noise(x)
p(y=1|x) = sigmoid(log_ratio - log(k))
参数:
x: 数据
log_p_model: log p_model(x)
log_p_noise: log p_noise(x)
返回:
p(y=1|x) 的值
"""
# log(p_model/p_noise) - log(k) = log_ratio - log(k)
log_ratio = log_p_model - log_p_noise - np.log(self.k)
# sigmoid
return 1.0 / (1.0 + np.exp(-np.clip(log_ratio, -500, 500)))
def nce_loss_and_gradient(self, real_data):
"""
计算NCE损失和梯度
损失函数（对每个真实样本x_i和k个噪声样本x'_j）：
L_i = -log p(y=1|x_i) - sum_j log p(y=0|x'_j)
总损失 L = (1/N) * sum_i L_i
梯度通过链式法则计算
参数:
real_data: 真实数据，形状 (batch_size, n_dim)
返回:
loss: 损失值
grad_A, grad_b, grad_log_Z: 各参数梯度
"""
batch_size = real_data.shape[0]
# === 生成噪声样本 ===
# 从噪声分布中为每个真实样本采样k个噪声样本
noise_samples = np.random.multivariate_normal(
self.noise_mean, self.noise_cov,
size=(batch_size, self.k)
) # 形状 (batch_size, k, n_dim)
# 重塑为 (batch_size * k, n_dim)
noise_flat = noise_samples.reshape(-1, self.n_dim)
# === 计算对数概率 ===
# 真实数据
log_p_real = self.log_model_prob(real_data) # 形状 (batch_size,)
log_q_real = self.log_noise_prob(real_data) # 形状 (batch_size,)
# 噪声样本
log_p_noise = self.log_model_prob(noise_flat) # 形状 (batch_size * k,)
log_q_noise = self.log_noise_prob(noise_flat) # 形状 (batch_size * k,)
# === 计算NCE分类概率 ===
p_real = self.nce_probability(real_data, log_p_real, log_q_real) # p(y=1|真实数据)
p_noise = self.nce_probability(noise_flat, log_p_noise, log_q_noise) # p(y=1|噪声样本)
# === 损失计算 ===
# -log p(y=1|x_real) - sum log p(y=0|x_noise)
# = -log p(y=1|x_real) - sum log(1-p(y=1|x_noise))
loss_real = -np.mean(np.log(p_real + 1e-10)) # 真实数据部分
loss_noise = -np.mean(np.log(1 - p_noise + 1e-10)) # 噪声数据部分
loss = loss_real + loss_noise
# === 梯度计算 ===
# 关键推导：
# 对于真实样本：∂L/∂θ = -(1 - p(y=1|x_real)) * ∂log p̃(x_real)/∂θ
# 对于噪声样本：∂L/∂θ = p(y=1|x_noise) * ∂log p̃(x_noise)/∂θ
# 真实数据的梯度权重
weight_real = (1 - p_real) / batch_size # 形状 (batch_size,)
# 噪声数据的梯度权重
weight_noise = p_noise / (batch_size * self.k) # 形状 (batch_size * k,)
# log p̃(x) = -0.5 x^T A x - b^T x
# ∂log p̃/∂A = -0.5 * x x^T (矩阵形式)
# ∂log p̃/∂b = -x
# A 的梯度
# 真实数据贡献
grad_A_real = np.zeros_like(self.A)
for i in range(batch_size):
grad_A_real += weight_real[i] * (-0.5 * np.outer(real_data[i], real_data[i]))
# 噪声数据贡献
grad_A_noise = np.zeros_like(self.A)
for i in range(len(noise_flat)):
grad_A_noise += weight_noise[i] * (-0.5 * np.outer(noise_flat[i], noise_flat[i]))
grad_A = grad_A_real + grad_A_noise # 总梯度
# b 的梯度
grad_b_real = -real_data.T @ weight_real # 形状 (n_dim,)
grad_b_noise = -noise_flat.T @ weight_noise # 形状 (n_dim,)
grad_b = grad_b_real + grad_b_noise
# log_Z 的梯度
# NCE的独特之处：它也学习Z！
# ∂L/∂log_Z = (1/N) * sum[(1-p(y=1|x_real))] - (k/N) * sum[p(y=1|x_noise)]
# 简化理解：当模型概率过高时，log_Z增大；反之减小
grad_log_Z = np.mean(1 - p_real) - np.mean(p_noise)
return loss, grad_A, grad_b, grad_log_Z
def train(self, data, epochs=300, batch_size=64):
"""
使用NCE训练模型
参数:
data: 真实训练数据
epochs: 训练轮数
batch_size: 批大小
"""
n_samples = data.shape[0]
print(f"开始NCE训练 (k={self.k}, 噪声方差={np.diag(self.noise_cov).mean():.2f})...")
for epoch in range(epochs):
indices = np.random.permutation(n_samples)
data_shuffled = data[indices]
total_loss = 0
n_batches = 0
for i in range(0, n_samples, batch_size):
batch = data_shuffled[i:i + batch_size]
loss, grad_A, grad_b, grad_log_Z = self.nce_loss_and_gradient(batch)
# 梯度下降（最小化损失）
self.A -= self.lr * grad_A
self.b -= self.lr * grad_b
self.log_Z -= self.lr * grad_log_Z
# 确保A对称
self.A = (self.A + self.A.T) / 2
total_loss += loss
n_batches += 1
if (epoch + 1) % 50 == 0:
print(f" Epoch {epoch+1}/{epochs}, "
f"损失: {total_loss/n_batches:.4f}, "
f"log_Z: {self.log_Z:.4f}")
# ============ NCE实验 ============
np.random.seed(42)
# 生成二维高斯数据
# 真实分布 N([1, -1], [[1.5, 0.3], [0.3, 0.8]])
true_mean = np.array([1.0, -1.0])
true_cov = np.array([[1.5, 0.3], [0.3, 0.8]])
train_data = np.random.multivariate_normal(true_mean, true_cov, size=1000)
# 设置NCE
# 噪声分布选一个较宽的高斯，覆盖数据所在区域
noise_cov = np.eye(2) * 5.0 # 大方差
model = NCEModel(
n_dim=2,
noise_mean=np.zeros(2),
noise_cov=noise_cov,
learning_rate=0.005,
n_noise_samples=20 # 每个真实样本对应20个噪声样本
)
model.train(train_data, epochs=300, batch_size=64)
# 结果比较
true_precision = np.linalg.inv(true_cov)
print(f"\n真实精度矩阵:\n{np.round(true_precision, 4)}")
print(f"学到的精度矩阵 A:\n{np.round(model.A, 4)}")
print(f"真实 log Z: {0.5*(np.log(np.linalg.det(true_cov)) + 2*np.log(2*np.pi)):.4f}")
print(f"学到的 log Z: {model.log_Z:.4f}")
print(f"\n真实均值: {true_mean}")
print(f"模型隐含均值: {np.round(-np.linalg.solve(model.A, model.b), 4)}")
八、知识点 7：估计配分函数
8.1 核心思想
有时我们仍然需要直接估计 ZZZ（例如计算似然值或进行模型比较）。
直接方法：对离散空间，遍历所有状态：
Z=∑xp~(x)Z = \sum_\mathbf{x} \tilde{p}(\mathbf{x})Z=x∑p~(x)
但这通常不可行（状态数指数增长）。
主要方法：
| 方法 | 思想 | 特点 |
|---|---|---|
| 桥式采样 | 利用两个分布之间的样本比率 | 需要从两个分布都能采样 |
| 退火重要性采样 (AIS) | 在 p0p_0p0 和 pmodelp_\text{model}pmodel 之间构造一系列中间分布 | 最常用，通过 MCMC 在各中间分布之间转移 |
| 变分上界 | 利用 logZ\log ZlogZ 的上界进行估计 | 给出 logZ\log ZlogZ 的上界估计 |
AIS 的核心思想：
构造温度序列 0=β0<β1<⋯<βT=10 = \beta_0 < \beta_1 < \cdots < \beta_T = 10=β0<β1<⋯<βT=1，定义中间分布：
pβt(x)=p~model(x)βt⋅p0(x)1−βtZβtp_{\beta_t}(\mathbf{x}) = \frac{\tilde{p}_\text{model}(\mathbf{x})^{\beta_t} \cdot p_0(\mathbf{x})^{1-\beta_t}}{Z_{\beta_t}}pβt(x)=Zβtp~model(x)βt⋅p0(x)1−βt
其中 p0p_0p0 是已知配分函数的简单分布。
利用重要性采样的比率关系：
Z1Z0=Ep0[∏t=1Twtwt−1]\frac{Z_1}{Z_0} = \mathbb{E}_{p_0} \left[ \prod_{t=1}^{T} \frac{w_t}{w_{t-1}} \right]Z0Z1=Ep0[t=1∏Twt−1wt]
8.2 案例代码：退火重要性采样（AIS）
import numpy as np
# ============================================
# 知识点7：退火重要性采样 (AIS) 估计配分函数
# ============================================
class AISPartitionEstimator:
"""
使用退火重要性采样(AIS)估计基于能量模型的配分函数
AIS 原理：
1. 选择一个已知Z的简单分布 p_0 (如高斯)
2. 在 p_0 和 p_model 之间插值一系列分布
3. 从 p_0 开始，依次在各插值分布间做MCMC转移
4. 累积重要性权重来估计 Z_model/Z_0
"""
def __init__(self, n_dim, energy_fn, beta_schedule=None):
"""
初始化AIS估计器
参数:
n_dim: 数据维度
energy_fn: 能量函数 E(x)，可调用
beta_schedule: 退火温度调度 [0, β1, β2, ..., 1]
"""
self.n_dim = n_dim
self.energy_fn = energy_fn
# 设置退火调度：如果没有指定，使用均匀间隔的100个温度
if beta_schedule is None:
self.betas = np.linspace(0, 1, 100)
else:
self.betas = beta_schedule
# 初始分布 p_0: 标准正态 N(0, I)
# 其配分函数 Z_0 = (2π)^{d/2}
self.log_Z_0 = 0.5 * n_dim * np.log(2 * np.pi)
def log_p_beta(self, x, beta):
"""
计算退火分布的（未归一化）对数概率
log p_β(x) = -β * E(x) - (1-β) * 0.5 * ||x||²
= -β * E(x) - (1-β) * E_0(x)
其中 E_0(x) = 0.5 * ||x||² 是标准正态的能量
参数:
x: 状态，形状 (n_dim,) 或 (batch, n_dim)
beta: 退火温度参数
返回:
log p_β(x)（未归一化）
"""
log_p_model = -self.energy_fn(x) # log p̃_model(x)
log_p_0 = -0.5 * np.sum(x ** 2) # log p_0(x) (标准正态的未归一化对数概率)
return beta * log_p_model + (1 - beta) * log_p_0
def log_weight_ratio(self, x, beta_from, beta_to):
"""
计算从 beta_from 到 beta_to 的重要性权重比
log(w_new/w_old) = log p_{β_to}(x) - log p_{β_from}(x)
参数:
x: 当前状态
beta_from: 起始温度
beta_to: 目标温度
返回:
log 权重比
"""
return self.log_p_beta(x, beta_to) - self.log_p_beta(x, beta_from)
def mcmc_step(self, x, beta, step_size=0.1):
"""
在退火分布 p_β 上执行一步Metropolis-Hastings MCMC
提议分布：高斯随机游走 q(x'|x) = N(x; x, step_size² I)
接受概率：α = min(1, p_β(x')/p_β(x))
参数:
x: 当前状态，形状 (n_dim,)
beta: 当前退火温度
step_size: 提议分布的标准差
返回:
x_new: 新状态（可能被拒绝，等于x）
accepted: 是否接受
"""
# 提议新状态：在当前状态附近高斯采样
x_proposal = x + np.random.randn(self.n_dim) * step_size
# 计算接受比（对数空间计算防止溢出）
log_accept_ratio = (self.log_p_beta(x_proposal, beta) -
self.log_p_beta(x, beta))
# Metropolis接受/拒绝
if np.log(np.random.uniform()) < log_accept_ratio:
return x_proposal, True # 接受
else:
return x, False # 拒绝
def run_ais(self, n_chains=100, n_mcmc_steps=10, step_size=0.5):
"""
运行AIS算法
流程：
1. 从 p_0 初始化 n_chains 条链
2. 对每个退火阶段 t=1,...,T:
a. 更新温度从 β_{t-1} 到 β_t
b. 累积重要性权重
c. 在 p_{β_t} 上做MCMC步以更新样本
3. 利用样本权重估计 Z_model/Z_0
参数:
n_chains: 并行运行的MCMC链数（越多估计越稳定）
n_mcmc_steps: 每个温度阶段的MCMC步数
step_size: MCMC步长
返回:
log_Z_estimate: log Z 的估计值
log_Z_std: log Z 估计的标准差（不确定性）
"""
# 步骤1：从初始分布 p_0 = N(0, I) 采样
chains = np.random.randn(n_chains, self.n_dim) # (n_chains, n_dim)
# 初始化 log 权重（每条链独立）
log_weights = np.zeros(n_chains)
# 步骤2：逐步退火
for t in range(1, len(self.betas)):
beta_prev = self.betas[t - 1] # 上一个温度
beta_curr = self.betas[t] # 当前温度
beta_diff = beta_curr - beta_prev # 温度增量
for chain_idx in range(n_chains):
# 2a: 更新重要性权重
# w_t = w_{t-1} * p_{β_t}(x) / p_{β_{t-1}}(x)
# 即 log w_t = log w_{t-1} + log_ratio
log_weight_change = self.log_weight_ratio(
chains[chain_idx], beta_prev, beta_curr
)
log_weights[chain_idx] += log_weight_change
# 2b: 在 p_{β_t} 上执行 MCMC 步
for _ in range(n_mcmc_steps):
chains[chain_idx], _ = self.mcmc_step(
chains[chain_idx], beta_curr, step_size
)
# 步骤3：估计 log(Z_model/Z_0)
# 使用 log-sum-exp 技巧数值稳定地计算
# log(mean(exp(log_weights))) = log_sum_exp(log_weights) - log(n_chains)
max_log_w = np.max(log_weights)
log_mean_weights = (max_log_w +
np.log(np.mean(np.exp(log_weights - max_log_w))))
# log Z_model = log Z_0 + log(mean(weights))
log_Z_estimate = self.log_Z_0 + log_mean_weights
# 使用 bootstrap 估计标准差
n_bootstrap = 1000
bootstrap_estimates = []
for _ in range(n_bootstrap):
# 有放回重采样
boot_indices = np.random.choice(n_chains, size=n_chains, replace=True)
boot_weights = log_weights[boot_indices]
boot_max = np.max(boot_weights)
boot_log_mean = boot_max + np.log(np.mean(np.exp(boot_weights - boot_max)))
bootstrap_estimates.append(self.log_Z_0 + boot_log_mean)
log_Z_std = np.std(bootstrap_estimates)
return log_Z_estimate, log_Z_std
class ToyIsingModel:
"""
二维 Ising 模型（玩具版）
在 L x L 的方格上，每个格点有自旋 σ_i ∈ {-1, +1}
能量：E(σ) = -J * sum_{<i,j>} σ_i * σ_j - h * sum_i σ_i
其中 <i,j> 表示最近邻对
"""
def __init__(self, L, J=1.0, h=0.0):
"""
初始化Ising模型
参数:
L: 方格边长（总格点数 L^2）
J: 耦合强度（正值=铁磁，负值=反铁磁）
h: 外磁场强度
"""
self.L = L # 方格边长
self.n = L * L # 总格点数
self.J = J # 耦合常数
self.h = h # 外场
# 预计算最近邻对的列表
self.neighbors = []
for i in range(L):
for j in range(L):
idx = i * L + j # 当前格点的一维索引
# 右邻居（周期性边界）
right = i * L + (j + 1) % L
# 下邻居（周期性边界）
down = ((i + 1) % L) * L + j
self.neighbors.append((idx, right))
self.neighbors.append((idx, down))
def energy(self, spins):
"""
计算Ising模型能量
E = -J * sum_{<i,j>} σ_i σ_j - h * sum_i σ_i
参数:
spins: 自旋配置，形状 (n,) 或 (batch, n)，值为 ±1
返回:
能量值
"""
if spins.ndim == 1:
spins = spins.reshape(1, -1)
batch_size = spins.shape[0]
energy = np.zeros(batch_size)
# 交互项
for (i, j) in self.neighbors:
energy -= self.J * spins[:, i] * spins[:, j]
# 外场项
energy -= self.h * np.sum(spins, axis=1)
return energy
def exact_log_partition(self):
"""
计算精确的 log Z（仅适用于小系统）
遍历所有 2^n 个状态
仅当 n <= 20 时可行
"""
if self.n > 20:
print(f"警告：系统太大 (n={self.n})，精确计算不可行")
return None
log_Z = -np.inf # 初始化为负无穷
total_states = 2 ** self.n
# 使用 log-sum-exp 技巧
log_unnorm_probs = []
for state_idx in range(total_states):
# 将整数转换为二进制自旋配置
# 0 → -1, 1 → +1
bits = np.array([(state_idx >> k) & 1 for k in range(self.n)])
spins = 2 * bits - 1 # 转换为 ±1
# 计算该状态的能量
E = self.energy(spins.reshape(1, -1))[0]
log_unnorm_probs.append(-E) # log p̃ = -E
log_unnorm_probs = np.array(log_unnorm_probs)
# log Z = log sum exp(log_p̃) = max + log sum exp(log_p̃ - max)
max_val = np.max(log_unnorm_probs)
log_Z = max_val + np.log(np.sum(np.exp(log_unnorm_probs - max_val)))
return log_Z
# ============ AIS 实验 ============
np.random.seed(42)
# 创建小规模Ising模型（4x4 = 16个格点）
L = 4
ising = ToyIsingModel(L=L, J=0.5, h=0.1) # 铁磁耦合+弱外场
print(f"Ising 模型: {L}x{L} 格点, J={ising.J}, h={ising.h}")
# 计算精确配分函数
exact_log_Z = ising.exact_log_partition()
print(f"\n精确 log Z: {exact_log_Z:.6f}")
# 使用AIS估计配分函数
print("\n运行AIS估计...")
ais = AISPartitionEstimator(
n_dim=ising.n,
energy_fn=ising.energy,
beta_schedule=np.linspace(0, 1, 200) # 200个退火阶段
)
# 多次运行AIS以评估稳定性
n_trials = 5
ais_estimates = []
for trial in range(n_trials):
log_Z_est, log_Z_std = ais.run_ais(
n_chains=500, # 500条并行链
n_mcmc_steps=5, # 每阶段5步MCMC
step_size=1.0 # MCMC步长
)
ais_estimates.append(log_Z_est)
print(f" Trial {trial+1}: log Z = {log_Z_est:.4f} ± {log_Z_std:.4f}")
ais_mean = np.mean(ais_estimates)
ais_std = np.std(ais_estimates)
print(f"\nAIS 估计: log Z = {ais_mean:.4f} ± {ais_std:.4f}")
print(f"精确值: log Z = {exact_log_Z:.6f}")
print(f"绝对误差: {abs(ais_mean - exact_log_Z):.6f}")
print(f"相对误差: {abs(ais_mean - exact_log_Z) / abs(exact_log_Z) * 100:.4f}%")
九、知识总结对比表
| 方法 | 是否需要 Z | 适用空间 | 计算复杂度 | 梯度是否有偏 |
|---|---|---|---|---|
| 对数似然梯度 | 需要采样近似 | 连续/离散 | 中 | 采样引入方差 |
| CD-1 | 不需要 | 离散为主 | 低 | 有偏 |
| SML/PCD | 不需要 | 离散为主 | 中 | 渐近无偏 |
| 伪似然 | 不需要 | 离散 | 低 (O(n)) | 无偏但目标不同 |
| 得分匹配 | 不需要 | 连续 | 中(需雅可比) | 无偏 |
| 去噪得分匹配 | 不需要 | 连续 | 低 | 无偏 |
| 比率匹配 | 不需要 | 离散 | 高 (O(n)) | 无偏 |
| NCE | 隐式学习 Z | 连续/离散 | 低 | 渐近无偏 |
| AIS | 直接估计 Z | 通用 | 高 | 有偏但可控 |
每个方法都是在计算可行性和估计质量之间做权衡。实践中，CD-1 在 RBM 中最常用，NCE 在语言模型中有重要应用，去噪得分匹配在生成模型中越来越受欢迎（如 score-based diffusion models）。
更多推荐
1
0- 0
已为社区贡献49条内容
所有评论(0)