时变模型下率失真理论计算速率为负的问题排查求助

问题背景 #

率失真理论是信息论中处理信源压缩与解压失真权衡的核心框架，率失真函数量化了达到指定失真水平所需的最小非负速率（比特/信源符号）。

在时变场景中，信源概率分布和失真约束随时间动态变化，率失真函数需依赖时变的概率与失真参数进行计算。当前实现的时变模型引入了多时间片的信源概率集合，但计算中出现了速率为负的异常结果，违反了速率非负的基本性质。

核心代码实现 #

import numpy as np

def findOptimalChannel(s, time_periods=2):
    
    # Define source
    N = 5
    p_list = [
        [0.45, 0.4, 0.05, 0.05, 0.05],
        [0.35, 0.5, 0.05, 0.05, 0.05]
    ]

    # Initialize matrix
    A_list = []
    for t in range(time_periods):
        A = np.ones((N, N))
        for i in range(0, N):
            for j in range(i):
                A[i, j] = np.exp(s + t * 0.1)
                A[j, i] = np.exp(s + t * 0.1)
        A_list.append(A)

    # initial guess
    q = np.ones(N) / N

    Tu = 1
    Tl = 0
    max_iterations = 1000
    iterations = 0
    while (Tu - Tl > 0.001) and (iterations < max_iterations):
        iterations += 1
        for t in range(time_periods):
            A = A_list[t]
            p = p_list[t]
            b = A.dot(q)
            c = np.zeros(N)
            for k in range(N):
                for j in range(N):
                    c[k] += A[j, k] * p[j] / b[j]

            q = q * c

            Tu = -np.sum(q * np.log(c))
            Tl = -np.max(np.log(c))
    print('out')
    # Compute channel
    Q = np.zeros((N, N))
    for k in range(N):
        for j in range(N):
            Q[k, j] = A[j, k] * q[k] / b[j]

    # Compute distortion
    D = 0
    for j in range(N):
        dummy = 0
        for k in range(N):
            dummy += Q[k, j] * (1 - (j == k))
        D += p[j] * dummy

    # Compute rate
    R = s * D - np.sum(p * np.log(b)) - np.sum(q * np.log(c))

    return R, D, A_list, Q

raw = []
for k in range(len(s)):
    raw.append(findOptimalChannel(s[k]))

潜在问题排查 #

1. Blahut-Arimoto算法的时变场景适配错误 #

标准Blahut-Arimoto算法针对单时间片的信源分布设计，当前代码在多时间片循环中直接覆盖更新q，未对不同时间片的贡献进行加权，导致迭代过程中q的更新逻辑不符合时变场景的联合优化目标，收敛结果异常。

2. 速率计算公式推导错误 #

当前速率公式R = s * D - np.sum(p * np.log(b)) - np.sum(q * np.log(c))存在符号或结构错误：

• 标准拉格朗日形式下，率失真函数的对偶问题中，速率应为互信息，公式需严格对应时变场景的加权互信息计算；
• 拉格朗日乘子s的符号与矩阵A的定义不匹配，若s为负，会导致公式项的符号反转，直接引发负速率。

3. 矩阵A的定义完全偏离理论要求 #

Blahut-Arimoto算法中，矩阵A的标准定义为A[x,y] = exp(-s * d(x,y))，其中d(x,y)是失真函数（此处为指示函数1-(x==y)）。当前代码中A[i,j] = exp(s + t*0.1)，既未引入失真函数，符号也完全错误，会导致后续b、c的计算值偏离合理范围，进而引发速率异常。

4. 失真计算未覆盖所有时间片 #

代码中仅使用最后一个时间片的信源概率p计算失真D，未对所有时间片的失真进行加权平均，导致D的计算值完全错误，直接影响速率R的结果。

5. 迭代收敛条件逻辑错误 #

每次处理单个时间片后就更新Tu和Tl，未在所有时间片处理完成后统一计算收敛指标，导致收敛判断不准确，迭代可能提前终止或无法收敛到最优解。

6. 概率分布未归一化 #

更新q时仅执行q = q * c，未对q进行归一化操作，导致q不再是合法的概率分布（和不为1），后续所有依赖q的计算（如b、c、速率）都会出现偏差。

解决思路 #

1. 修正时变场景的算法框架 #

• 若为独立时间片优化：对每个时间片单独运行Blahut-Arimoto算法，最后按时间权重加权得到总速率；
• 若为联合优化：重新推导多时间片的率失真目标函数，将各时间片的信源概率作为权重，构建联合优化的迭代更新规则。

2. 校准速率与失真的计算公式 #

• 回到率失真理论的拉格朗日对偶推导，确认时变场景下的速率公式：总速率应为各时间片互信息的加权和，即R = sum( w_t * I_t )，其中w_t是时间片t的权重；
• 校验拉格朗日乘子s的符号，确保与矩阵A的定义匹配（A[x,y] = exp(-s*d(x,y))）。

3. 修正矩阵A的定义 #

将矩阵A改为符合理论要求的形式：

for t in range(time_periods):
    A = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            d = 1 - (i == j)  # 失真函数
            A[i, j] = np.exp(-s * d)
    A_list.append(A)

4. 修正失真计算逻辑 #

遍历所有时间片，计算每个时间片的失真并加权：

D = 0
for t in range(time_periods):
    A = A_list[t]
    p = p_list[t]
    # 计算当前时间片的Q
    Q_t = np.zeros((N, N))
    b_t = A.dot(q)
    for k in range(N):
        for j in range(N):
            Q_t[k, j] = A[j, k] * q[k] / b_t[j]
    # 计算当前时间片的失真D_t
    D_t = 0
    for j in range(N):
        dummy = 0
        for k in range(N):
            dummy += Q_t[k, j] * (1 - (j == k))
        D_t += p[j] * dummy
    # 按时间权重累加（此处假设等权重）
    D += D_t / time_periods

5. 修复迭代收敛与概率归一化 #

• 在所有时间片处理完成后，再计算Tu和Tl判断收敛；
• 更新q后立即归一化：q = q / np.sum(q)，确保q是合法的概率分布。

6. 添加数值稳定性处理 #

在计算log操作时，添加极小值避免对数为负无穷：

np.log(b + 1e-12)  # 替换原np.log(b)
np.log(c + 1e-12)  # 替换原np.log(c)

内容的提问来源于stack exchange，提问作者user21412360