第8章率失真理论（Rate-Distortion Theory）

一、导读

8.1 本章学习目标

通过本章的学习，你将能够：

理解率失真问题的基本设定：掌握有损压缩与无损压缩的本质区别

掌握率失真函数 $R(D)$ 的定义：理解其作为优化问题的数学形式

掌握率失真函数的性质：单调性、凸性、连续性

能够计算简单信源的率失真函数：伯努利信源、高斯信源

理解率失真定理的证明思路：典型序列编码和联合典型性译码的应用

掌握平方误差失真下的高斯率失真函数：理解"反向水填充"的概念

理解量化与熵编码的结合：标量量化、矢量量化的区别

能够分析实际压缩系统的率失真性能：理解PSNR、SSIM等质量度量

8.2 为什么本章如此重要

率失真理论是有损压缩的理论基础。从JPEG图片到MP3音频，从H.264视频到深度学习模型压缩，所有有损压缩技术都在率失真权衡的框架内运作。

理论重要性：

率失真函数 $R(D)$ 给出了在给定失真约束下的最小压缩率。这与信道容量 $C$ 给出了可靠通信的最大速率形成对偶。事实上，率失真理论和信道编码理论在数学结构上有很多相似之处。

实践重要性：

现代压缩标准（如JPEG、MPEG、AAC）都隐含地利用了率失真权衡。编码器需要在比特率和重建质量之间做出选择，这正是率失真理论的核心问题。理解 $R(D)$ 有助于设计更高效的压缩算法，以及评估现有算法的性能极限。

方法论重要性：

率失真问题的求解涉及凸优化、变分法、信息投影等高级数学工具。这些工具在机器学习、信号处理、经济学等领域也有广泛应用。

8.3 核心问题引导

在学习本章之前，请思考以下问题：

什么是有损压缩？ 为什么我们愿意接受一定的失真？失真如何量化？

率失真函数 $R(D)$ 的物理意义是什么？ 为什么它是单调递减的？为什么它是凸函数？

对于高斯信源，率失真函数为什么是 $R(D) = \frac{1}{2}\log_2(\sigma^2/D)$？ 这个结果的直观解释是什么？

率失真定理如何证明？ 与信道编码定理的证明有什么相似之处？

量化和熵编码有什么关系？ 为什么先量化再熵编码可以接近率失真极限？

矢量量化为什么比标量量化更高效？ 维度增加如何改善率失真性能？

感知质量与客观失真度量有什么关系？ 为什么PSNR高的图片可能看起来不好？

率失真理论与信道容量有什么对偶关系？

二、核心概念详解

8.1 率失真问题形式化

8.1.1 基本设置

信源：离散或连续随机变量 $X$，服从分布 $p(x)$。

重建字母表：$\hat{\mathcal{X}}$，可能与 $\mathcal{X}$ 相同或不同。

失真度量：函数 $d: \mathcal{X} \times \hat{\mathcal{X}} \to [0, \infty)$，常用的有：

汉明失真：$d(x, \hat{x}) = \mathbb{1}\{x \neq \hat{x}\}$
平方误差：$d(x, \hat{x}) = (x - \hat{x})^2$
绝对误差：$d(x, \hat{x}) = |x - \hat{x}|$

$(n, M, D)$ 率失真码：

编码器：$f: \mathcal{X}^n \to \{1, 2, \ldots, M\}$
译码器：$g: \{1, 2, \ldots, M\} \to \hat{\mathcal{X}}^n$
平均失真：$\frac{1}{n}\sum_{i=1}^n \mathbb{E}[d(X_i, \hat{X}_i)] \leq D$

可达速率：如果存在一系列 $(n, 2^{nR}, D)$ 码，使得 $n \to \infty$ 时失真趋于 $D$，则速率 $R$ 是可达的。

率失真函数：

$$R(D) = \inf\{R : R \text{ is achievable at distortion } D\}$$

8.1.2 率失真函数的优化形式

定理：率失真函数可以表示为：

$$R(D) = \min_{p(\hat{x}|x): \mathbb{E}[d(X,\hat{X})] \leq D} I(X; \hat{X})$$

解释：在所有满足失真约束的条件分布中，找到使互信息最小的那个。这个最小互信息就是率失真函数。

直观理解：

$I(X;\hat{X})$ 度量了原始信号和重建信号之间的信息共享
失真约束限制了重建信号不能偏离原始信号太多
最小化互信息意味着在满足失真要求的前提下，尽可能少地传输信息

8.1.3 与信源编码的关系

无损压缩是 $D=0$ 的特例：

$$R(0) = H(X)$$

当 $D$ 增加时，$R(D)$ 减小，因为允许更大的失真意味着可以更激进地压缩。

8.2 率失真函数的性质

8.2.1 单调性

性质：$R(D)$ 是非增函数，即 $D_1 < D_2 \Rightarrow R(D_1) \geq R(D_2)$。

证明：如果 $D_1 < D_2$，则满足 $D_1$ 约束的条件分布集合是满足 $D_2$ 约束的子集。因此最小值在更大的集合上不会更大。

8.2.2 凸性

性质：$R(D)$ 是凸函数。

证明思路：使用互信息的凸性和失真约束的线性性质。对于两个满足失真 $D_1$ 和 $D_2$ 的条件分布，它们的凸组合满足失真 $\lambda D_1 + (1-\lambda)D_2$，且互信息不超过凸组合。

8.2.3 边界值

$R(0) = H(X)$（无损压缩）
$R(D) = 0$ 当 $D \geq D_{\max}$，其中 $D_{\max} = \min_{\hat{x}} \mathbb{E}[d(X, \hat{x})]$

代码示例8-1：率失真函数性质可视化

python

import numpy as np
import matplotlib.pyplot as plt

def binary_entropy(p):
    if p == 0 or p == 1:
        return 0
    return -p * np.log2(p) - (1-p) * np.log2(1-p)

def rate_distortion_bernoulli(p, D_values):
    """伯努利信源的率失真函数（汉明失真）"""
    D_max = min(p, 1-p)
    R_values = []
    
    for D in D_values:
        if D <= D_max:
            R = binary_entropy(p) - binary_entropy(D)
        else:
            R = 0
        R_values.append(max(0, R))
    
    return R_values, D_max

def rate_distortion_gaussian(sigma2, D_values):
    """高斯信源的率失真函数（平方误差失真）"""
    R_values = []
    for D in D_values:
        if D < sigma2:
            R = 0.5 * np.log2(sigma2 / D)
        else:
            R = 0
        R_values.append(R)
    
    return R_values

# 伯努利信源
p = 0.3
D_values_b = np.linspace(0, 0.5, 100)
R_values_b, D_max = rate_distortion_bernoulli(p, D_values_b)

plt.figure(figsize=(10, 6))
plt.plot(D_values_b, R_values_b, 'b-', linewidth=2.5)
plt.axvline(x=D_max, color='r', linestyle='--', alpha=0.5, label=f'D_max={D_max:.4f}')
plt.axhline(y=binary_entropy(p), color='g', linestyle='--', alpha=0.5, 
            label=f'H(X)={binary_entropy(p):.4f}')
plt.xlabel('失真 D', fontsize=12)
plt.ylabel('率 R(D) (bits/symbol)', fontsize=12)
plt.title(f'伯努利({p})信源的率失真函数（汉明失真）', fontsize=14)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('rd_bernoulli.png', dpi=150)
plt.show()

# 高斯信源
sigma2 = 1.0
D_values_g = np.linspace(0.01, 2, 100)
R_values_g = rate_distortion_gaussian(sigma2, D_values_g)

plt.figure(figsize=(10, 6))
plt.plot(D_values_g, R_values_g, 'r-', linewidth=2.5)
plt.axhline(y=0.5*np.log2(sigma2), color='g', linestyle='--', alpha=0.5, 
            label=f'h(X)/2={0.5*np.log2(sigma2*2*np.pi*np.e):.4f}')
plt.xlabel('失真 D', fontsize=12)
plt.ylabel('率 R(D) (bits/symbol)', fontsize=12)
plt.title(f'高斯N(0,{sigma2})信源的率失真函数（平方误差）', fontsize=14)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('rd_gaussian.png', dpi=150)
plt.show()

print("率失真函数分析")
print("=" * 60)
print(f"伯努利({p})信源:")
print(f"  H(X) = {binary_entropy(p):.4f}")
print(f"  D_max = {D_max:.4f}")
print(f"\n高斯(0,{sigma2})信源:")
print(f"  R(D) = 0.5*log2({sigma2}/D)")

8.3 伯努利信源的率失真函数

8.3.1 问题设定

信源：$X \sim \text{Bernoulli}(p)$

失真：汉明失真 $d(x, \hat{x}) = \mathbb{1}\{x \neq \hat{x}\}$

8.3.2 推导

我们需要最小化 $I(X;\hat{X})$，约束为 $\mathbb{E}[d(X,\hat{X})] \leq D$。

设 $P(\hat{X} \neq X) = q \leq D$。

互信息：

$$I(X;\hat{X}) = H(X) - H(X|\hat{X})$$

由于对称性，最优条件分布是二元对称信道：

$$P(\hat{X}=1|X=0) = P(\hat{X}=0|X=1) = q$$

此时：

$$H(X|\hat{X}) = H_b(q)$$

因此：

$$R(D) = H_b(p) - H_b(D)$$

当 $D > \min(p, 1-p)$ 时，$R(D) = 0$（可以直接输出常数）。

8.3.3 最优测试信道

达到率失真函数的条件分布 $p^*(\hat{x}|x)$ 对应于一个BSC($D$)信道。

直观解释：为了以失真 $D$ 重建 $X$，我们可以让 $X$ 通过一个交叉概率为 $D$ 的二元对称信道，输出即为 $\hat{X}$。这样平均失真恰好为 $D$，而互信息为 $H_b(p) - H_b(D)$。

8.4 高斯信源的率失真函数

8.4.1 问题设定

信源：$X \sim \mathcal{N}(0, \sigma^2)$

失真：平方误差 $d(x, \hat{x}) = (x - \hat{x})^2$

8.4.2 推导

我们需要最小化 $I(X;\hat{X})$，约束为 $\mathbb{E}[(X-\hat{X})^2] \leq D$。

关键观察：对于给定的方差，高斯分布具有最大微分熵。因此，最优的条件分布应该使 $X-\hat{X}$ 是高斯的。

设 $X = \hat{X} + Z$，其中 $Z \sim \mathcal{N}(0, D)$ 且与 $\hat{X}$ 独立。

则 $\hat{X} \sim \mathcal{N}(0, \sigma^2 - D)$（假设 $D < \sigma^2$）。

互信息：

$$\begin{align*}

I(X;\hat{X}) &= h(X) - h(X|\hat{X}) \\

&= h(X) - h(Z) \\

&= \frac{1}{2}\log_2(2\pi e \sigma^2) - \frac{1}{2}\log_2(2\pi e D) \\

&= \frac{1}{2}\log_2\left(\frac{\sigma^2}{D}\right)

\end{align*}$$

结果：

$$R(D) = \begin{cases} \frac{1}{2}\log_2\left(\frac{\sigma^2}{D}\right) & 0 \leq D \leq \sigma^2 \\ 0 & D > \sigma^2 \end{cases}$$

8.4.3 直观解释

当 $D = \sigma^2$ 时，$R = 0$：可以不发送任何信息，重建为0（均值），失真为方差
当 $D \to 0$ 时，$R \to \infty$：无损压缩高斯信源需要无限比特
每减少一半的失真，需要额外1比特/符号

代码示例8-2：高斯率失真函数应用

python

import numpy as np
import matplotlib.pyplot as plt

def gaussian_rd(sigma2, D):
    """高斯信源率失真函数"""
    if D >= sigma2:
        return 0
    return 0.5 * np.log2(sigma2 / D)

# 应用场景：图像压缩
# 假设像素值归一化到[0,1]，方差约为0.1

sigma2_img = 0.1

# 不同质量级别对应的失真
quality_levels = np.linspace(0.01, 0.15, 50)
rates = [gaussian_rd(sigma2_img, D) for D in quality_levels]

plt.figure(figsize=(10, 6))
plt.plot(quality_levels, rates, 'b-', linewidth=2.5)
plt.xlabel('失真 D (MSE)', fontsize=12)
plt.ylabel('率 R(D) (bits/pixel)', fontsize=12)
plt.title('图像压缩的率失真曲线（高斯近似）', fontsize=14)
plt.grid(True, alpha=0.3)
plt.axhline(y=8, color='r', linestyle='--', alpha=0.3, label='未压缩: 8 bits/pixel')
plt.legend()
plt.tight_layout()
plt.savefig('image_rd.png', dpi=150)
plt.show()

print("图像压缩率失真分析")
print("=" * 60)
print(f"信号方差: {sigma2_img}")
for D in [0.01, 0.02, 0.05, 0.1]:
    R = gaussian_rd(sigma2_img, D)
    psnr = 10 * np.log10(1/D)  # 峰值信噪比
    print(f"D={D:.4f}: R={R:.4f} bits/pixel, PSNR≈{psnr:.1f} dB")

8.5 率失真定理

8.5.1 定理陈述

定理（率失真定理）：对于离散无记忆信源和单字母失真度量：

正面（Achievability）：对于任意 $R > R(D)$，存在码率为 $R$、平均失真不超过 $D$ 的编码方案。

反面（Converse）：任何码率 $R < R(D)$ 的编码方案必然有平均失真大于 $D$。

8.5.2 证明思路（正面）

随机码本构造：

选择达到 $R(D)$ 的最优条件分布 $p^*(\hat{x}|x)$

生成 $2^{nR}$ 个重建序列 $\hat{x}^n$，每个符号i.i.d.服从 $p^(\hat{x}) = \sum_x p(x)p^(\hat{x}|x)$

这些重建序列构成"码本"

编码：

给定源序列 $x^n$，寻找码本中与之联合典型的 $\hat{x}^n$
如果找到，输出其索引；否则报错

译码：输出索引对应的重建序列

错误分析：

由联合典型性引理，对于给定的 $x^n$，随机生成的 $\hat{x}^n$ 与之联合典型的概率约 $2^{-nI(X;\hat{X})}$
如果有 $2^{nR}$ 个候选，找到联合典型对的概率高（当 $R > I(X;\hat{X})$）
联合典型性保证平均失真接近 $D$

8.5.3 证明思路（反面）

假设存在码率 $R < R(D)$、失真 $\leq D$ 的编码方案。

则 $I(X^n;\hat{X}^n) \leq nR < nR(D)$。

但由率失真函数的定义，任何失真 $\leq D$ 的重建必须满足 $I(X;\hat{X}) \geq R(D)$。

矛盾。因此这样的编码方案不存在。

8.6 量化与熵编码

8.6.1 标量量化

均匀量化：将输入范围划分为等宽区间，每个区间用一个代表值（通常是区间中心）表示。

量化步长 $\Delta$ 与失真的关系（高斯信源）：

$$D \approx \frac{\Delta^2}{12}$$

熵编码：量化后的索引可以用熵编码（如Huffman）进一步压缩。

8.6.2 矢量量化

思想：一次量化多个样本（一个向量），而不是逐个量化。

优势：

可以利用样本间的相关性
随着维度增加，性能趋近率失真极限

挑战：

码本设计复杂度高
搜索最优码字计算量大

8.6.3 率失真优化的量化

Lloyd-Max算法：迭代优化量化器的判决边界和代表值，最小化给定比特率下的失真。

步骤：

固定代表值，优化判决边界（最近邻规则）

固定判决边界，优化代表值（质心规则）

重复直到收敛

代码示例8-3：Lloyd-Max算法实现

python

import numpy as np
import matplotlib.pyplot as plt

def lloyd_max_algorithm(data, n_levels, max_iter=100):
    """
    Lloyd-Max量化器设计
    
    参数:
        data: 训练数据
        n_levels: 量化级数
        max_iter: 最大迭代次数
    
    返回:
        boundaries: 判决边界
        representatives: 代表值
    """
    
    # 初始化：均匀分布的代表值
    min_val, max_val = data.min(), data.max()
    representatives = np.linspace(min_val, max_val, n_levels)
    
    for iteration in range(max_iter):
        # 步骤1：分配（最近邻）
        assignments = np.argmin(np.abs(data.reshape(-1, 1) - representatives.reshape(1, -1)), axis=1)
        
        # 步骤2：更新代表值（质心）
        new_representatives = np.zeros(n_levels)
        for i in range(n_levels):
            mask = assignments == i
            if np.sum(mask) > 0:
                new_representatives[i] = np.mean(data[mask])
            else:
                new_representatives[i] = representatives[i]
        
        # 检查收敛
        if np.max(np.abs(new_representatives - representatives)) < 1e-6:
            break
        
        representatives = new_representatives
    
    # 计算判决边界
    boundaries = (representatives[:-1] + representatives[1:]) / 2
    boundaries = np.concatenate([[-np.inf], boundaries, [np.inf]])
    
    return boundaries, representatives

def quantize(data, boundaries, representatives):
    """量化数据"""
    quantized = np.zeros_like(data)
    for i in range(len(representatives)):
        mask = (data >= boundaries[i]) & (data < boundaries[i+1])
        quantized[mask] = representatives[i]
    return quantized

# 测试
np.random.seed(42)
data = np.random.normal(0, 1, 10000)

n_levels = 8
boundaries, reps = lloyd_max_algorithm(data, n_levels)

# 量化
quantized = quantize(data, boundaries, reps)

# 计算失真
distortion = np.mean((data - quantized)**2)

# 计算熵
from collections import Counter
counts = Counter(np.round(quantized, 6))
probs = np.array(list(counts.values())) / len(quantized)
entropy = -np.sum(probs * np.log2(probs))

print("Lloyd-Max量化器")
print("=" * 60)
print(f"量化级数: {n_levels}")
print(f"失真 (MSE): {distortion:.6f}")
print(f"熵: {entropy:.4f} bits")
print(f"代表值: {reps.round(4)}")

# 可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

ax1.hist(data, bins=50, density=True, alpha=0.5, label='原始数据')
for r in reps:
    ax1.axvline(x=r, color='r', linestyle='--', alpha=0.3)
ax1.set_xlabel('值', fontsize=12)
ax1.set_ylabel('密度', fontsize=12)
ax1.set_title('Lloyd-Max量化器代表值', fontsize=14)
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.scatter(data[:500], quantized[:500], alpha=0.5, s=10)
ax2.plot([-3, 3], [-3, 3], 'r--', alpha=0.5, label='理想')
ax2.set_xlabel('原始值', fontsize=12)
ax2.set_ylabel('量化值', fontsize=12)
ax2.set_title('量化散点图', fontsize=14)
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lloyd_max.png', dpi=150)
plt.show()

8.7 并行高斯信源与水填充对偶

8.7.1 问题设定

考虑 $k$ 个独立的高斯信源，第 $i$ 个的方差为 $\sigma_i^2$。总比特率约束为 $R$。

优化问题：分配比特率 $R_i$ 给每个信源，最小化总失真：

$$\min_{R_1, \ldots, R_k} \sum_{i=1}^k D_i(R_i) \quad \text{s.t.} \quad \sum_{i=1}^k R_i \leq R$$

其中 $D_i(R_i) = \sigma_i^2 2^{-2R_i}$ 是高斯信源的失真率函数。

8.7.2 反向水填充解

使用拉格朗日乘子法，最优解满足：

$$R_i = \max\left(0, \frac{1}{2}\log_2\frac{\sigma_i^2}{\mu}\right)$$

其中 $\mu$ 由总比特率约束确定。

直观解释：这与信道容量的水填充类似，但是"反向"的——我们将更多比特分配给方差更大的信源（因为它们更难压缩）。

代码示例8-4：反向水填充

python

import numpy as np
import matplotlib.pyplot as plt

def reverse_water_filling(variances, total_rate):
    """
    反向水填充算法
    
    参数:
        variances: 各信源的方差
        total_rate: 总比特率约束
    
    返回:
        比特率分配，阈值
    """
    k = len(variances)
    var = np.array(variances)
    
    # 排序方差
    sorted_indices = np.argsort(var)[::-1]  # 降序
    sorted_var = var[sorted_indices]
    
    # 二分搜索阈值
    def compute_rates(mu):
        rates = np.maximum(0, 0.5 * np.log2(sorted_var / mu))
        return rates
    
    # 找到合适的mu
    mu_low = 1e-10
    mu_high = np.max(sorted_var)
    
    for _ in range(100):
        mu_mid = (mu_low + mu_high) / 2
        rates_mid = compute_rates(mu_mid)
        if np.sum(rates_mid) > total_rate:
            mu_low = mu_mid
        else:
            mu_high = mu_mid
    
    mu = (mu_low + mu_high) / 2
    rates = compute_rates(mu)
    
    # 恢复原始顺序
    final_rates = np.zeros(k)
    for i, idx in enumerate(sorted_indices):
        final_rates[idx] = rates[i]
    
    return final_rates, mu

# 测试
variances = [1, 2, 3, 4, 5]
total_rate = 5

rates, mu = reverse_water_filling(variances, total_rate)
distortions = [v * 2**(-2*r) if r > 0 else v for v, r in zip(variances, rates)]

print("反向水填充比特分配")
print("=" * 60)
print(f"方差: {variances}")
print(f"总比特率: {total_rate}")
print(f"阈值 μ: {mu:.4f}")
print(f"比特分配: {rates.round(4)}")
print(f"总比特: {np.sum(rates):.4f}")
print(f"失真: {[round(d, 4) for d in distortions]}")
print(f"总失真: {sum(distortions):.4f}")

# 可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

channels = np.arange(len(variances))
ax1.bar(channels, variances, color='blue', alpha=0.5, label='方差')
ax1.bar(channels, rates, color='red', alpha=0.7, label='比特分配')
ax1.set_xlabel('信源索引', fontsize=12)
ax1.set_ylabel('值', fontsize=12)
ax1.set_title('反向水填充比特分配', fontsize=14)
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.bar(channels, distortions, color='green', alpha=0.7)
ax2.set_xlabel('信源索引', fontsize=12)
ax2.set_ylabel('失真', fontsize=12)
ax2.set_title('各信源失真', fontsize=14)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('reverse_water_filling.png', dpi=150)
plt.show()

三、程序员的视角

8.1 在实际压缩标准中的应用

JPEG：

DCT变换 → 频域系数

量化（根据量化表，高频系数量化更粗）

Zig-zag扫描 + 游程编码 + Huffman编码

MP3/AAC：

心理声学模型确定掩蔽阈值

MDCT变换

量化以满足掩蔽约束

熵编码

视频编码（H.264/H.265）：

运动估计 + 残差计算

DCT/整数变换

量化 + 熵编码（CAVLC/CABAC）

8.2 质量度量

PSNR（峰值信噪比）：

$$\text{PSNR} = 10 \log_{10}\left(\frac{\text{MAX}^2}{\text{MSE}}\right)$$

SSIM（结构相似性）：考虑亮度、对比度、结构的相似性。

VMAF（视频多方法评估融合）：结合多种指标的感知质量度量。

8.3 常见陷阱

PSNR不等于感知质量：高PSNR的图片可能看起来不好

忽略 perceptual shaping：人眼对不同频率的敏感度不同

块效应：低比特率下DCT块的边界可见

四、代码实战

8.1 完整的率失真性能分析工具

python

import numpy as np
import matplotlib.pyplot as plt

class RateDistortionAnalyzer:
    """率失真性能分析器"""
    
    def __init__(self, data):
        self.data = np.array(data).flatten()
        self.variance = np.var(self.data)
    
    def theoretical_rd_gaussian(self, D_values):
        """理论高斯率失真函数"""
        return [0.5 * np.log2(self.variance / D) if D < self.variance else 0 
                for D in D_values]
    
    def empirical_rd_uniform_quantization(self, n_bits_values):
        """
        经验率失真曲线（均匀量化）
        
        参数:
            n_bits_values: 不同的比特数列表
        
        返回:
            rates, distortions
        """
        rates = []
        distortions = []
        
        for n_bits in n_bits_values:
            n_levels = 2**n_bits
            
            # 均匀量化
            min_val, max_val = self.data.min(), self.data.max()
            delta = (max_val - min_val) / n_levels
            representatives = min_val + delta/2 + np.arange(n_levels) * delta
            
            # 量化
            indices = np.clip(((self.data - min_val) / delta).astype(int), 0, n_levels-1)
            quantized = representatives[indices]
            
            # 失真
            D = np.mean((self.data - quantized)**2)
            
            # 实际熵
            from collections import Counter
            counts = Counter(indices)
            probs = np.array(list(counts.values())) / len(self.data)
            H = -np.sum(probs * np.log2(probs))
            
            rates.append(H)
            distortions.append(D)
        
        return rates, distortions
    
    def plot_rd_curve(self, D_values=None, n_bits_range=range(1, 9)):
        """绘制率失真曲线"""
        
        if D_values is None:
            D_values = np.linspace(0.01 * self.variance, self.variance, 100)
        
        # 理论曲线
        R_theoretical = self.theoretical_rd_gaussian(D_values)
        
        # 经验曲线
        rates_emp, dist_emp = self.empirical_rd_uniform_quantization(list(n_bits_range))
        
        fig, ax = plt.subplots(figsize=(10, 6))
        ax.plot(D_values, R_theoretical, 'b-', linewidth=2.5, label='理论极限（高斯）')
        ax.plot(dist_emp, rates_emp, 'ro-', linewidth=2, markersize=8, label='均匀量化')
        ax.set_xlabel('失真 D (MSE)', fontsize=12)
        ax.set_ylabel('率 R(D) (bits/sample)', fontsize=12)
        ax.set_title('率失真性能分析', fontsize=14)
        ax.legend(fontsize=11)
        ax.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.savefig('rd_analysis.png', dpi=150)
        plt.show()


# 使用示例
if __name__ == "__main__":
    np.random.seed(42)
    data = np.random.normal(0, 1, 10000)
    
    analyzer = RateDistortionAnalyzer(data)
    analyzer.plot_rd_curve()
    
    print(f"数据方差: {analyzer.variance:.4f}")

五、深入思考

5.1 历史背景

率失真理论由香农在1959年的论文中系统建立。Berger、Gray、Neuhoff等人后续做出了重要贡献。Lloyd-Max算法分别由Lloyd（1957，内部报告）和Max（1960）独立提出。

5.2 前沿方向

神经压缩：使用深度学习进行端到端压缩，学习非线性变换和概率模型。

感知率失真理论：结合感知质量度量，而非简单的MSE。

分布式率失真：多个相关信源的联合压缩（Slepian-Wolf编码）。

5.3 延伸阅读

Cover & Thomas《信息论基础》第10章、第13章

Gray "Source Coding Theory"

Gersho & Gray "Vector Quantization and Signal Compression"

六、实验笔记

实验1：量化器性能比较

目标：比较均匀量化、Lloyd-Max量化和高斯率失真极限的性能。

步骤：

生成高斯数据

设计不同比特数的均匀量化器和Lloyd-Max量化器

测量失真和实际熵

绘制率失真曲线，与理论极限比较

实验2：矢量量化增益

目标：验证矢量量化相对于标量量化的增益。

步骤：

实现2维和4维矢量量化

比较相同比特率下的失真

分析增益随维度的变化

七、本章小结

本章介绍了率失真理论的基本概念和应用。

率失真函数 $R(D)$ 给出了给定失真约束下的最小压缩率。它是单调递减的凸函数，边界值为 $R(0)=H(X)$ 和 $R(D_{\max})=0$。

伯努利信源的率失真函数为 $H_b(p)-H_b(D)$，高斯信源的为 $\frac{1}{2}\log_2(\sigma^2/D)$。这些闭式解提供了重要的直觉。

率失真定理证明了 $R(D)$ 是可达的极限，证明使用了随机码本和联合典型性译码。

量化是实现率失真权衡的实际方法。Lloyd-Max算法优化了标量量化器，而矢量量化通过利用样本间相关性可以获得更好的性能。

反向水填充解决了并行信源的比特分配问题，与信道容量的水填充形成对偶。

从程序员的角度看，理解率失真理论有助于设计高效的压缩算法，选择合适的质量参数，以及评估压缩性能的理论上界。

八、习题

习题8.1（伯努利信源）：计算 $Bernoulli(0.2)$ 信源在汉明失真 $D=0.1$ 下的率失真函数值。

习题8.2（高斯信源）：高斯信源方差为4，要求失真不超过0.5。计算所需的最小比特率。

习题8.3（性质证明）：证明率失真函数 $R(D)$ 是凸函数。

习题8.4（编程题）：实现Lloyd-Max算法，比较它与均匀量化在不同分布下的性能。

习题8.5（反向水填充）：三个高斯信源方差分别为1、4、9。总比特率约束为3 bits。计算最优比特分配和各信源失真。

习题8.6（挑战题）：推导拉普拉斯信源在绝对误差失真下的率失真函数。

习题8.7（应用题）：在设计图像压缩系统时，如何在PSNR和感知质量之间权衡？

习题8.8（分析题）：比较标量量化和矢量量化的优缺点。为什么实际系统中主要使用标量量化？

习题8.9（开放性问题）：研究神经图像压缩方法的最新进展，撰写一篇500字左右的分析报告。

习题8.10（综合题）：结合率失真理论和信道容量，分析一个端到端通信系统的性能极限（源压缩 + 信道传输）。

第8章 率失真理论（Rate-Distortion Theory）

一、导读

8.1 本章学习目标

8.2 为什么本章如此重要

8.3 核心问题引导

二、核心概念详解

8.1 率失真问题形式化

8.1.1 基本设置

8.1.2 率失真函数的优化形式

8.1.3 与信源编码的关系

8.2 率失真函数的性质

8.2.1 单调性

8.2.2 凸性

8.2.3 边界值

8.3 伯努利信源的率失真函数

8.3.1 问题设定

8.3.2 推导

8.3.3 最优测试信道

8.4 高斯信源的率失真函数

8.4.1 问题设定

8.4.2 推导

8.4.3 直观解释

8.5 率失真定理

8.5.1 定理陈述

8.5.2 证明思路（正面）

8.5.3 证明思路（反面）

8.6 量化与熵编码

8.6.1 标量量化

8.6.2 矢量量化

8.6.3 率失真优化的量化

8.7 并行高斯信源与水填充对偶

8.7.1 问题设定

8.7.2 反向水填充解

三、程序员的视角

8.1 在实际压缩标准中的应用

8.2 质量度量

8.3 常见陷阱

四、代码实战

8.1 完整的率失真性能分析工具

五、深入思考

5.1 历史背景

5.2 前沿方向

5.3 延伸阅读

六、实验笔记

实验1：量化器性能比较

实验2：矢量量化增益

七、本章小结

八、习题

第8章率失真理论（Rate-Distortion Theory）