通信通道的容量

作者：Robert Gowers、Roger Hill、Sami Al-Izzi、Timothy Pollington和Keith Briggs

来自Boyd和Vandenberghe的《凸优化》，习题4.57，第207-8页

凸优化可用于找到离散无记忆通道的通道容量 \(C\)。考虑一个输入 \(X(t) \in \{1,2,...,n\}\) 和输出 \(Y(t) \in \{1,2,...m\}\) 的通信通道。这意味着随机变量 \(X\) 和 \(Y\) 分别可以取 \(n\) 和 \(m\) 种不同的值。

在离散无记忆通道中，输入和输出之间的关系由转移概率给出：

\(p_{ij} = \mathbb{P}(Y(t)=i | X(t)=j)\)

这些转移概率构成了通道转移矩阵 \(P\)，其中 \(P \in \mathbb{R}^{m\times n}\)。

假设 \(X\) 具有概率分布，用 \(x \in \mathbb{R}^n\) 表示，即：

\(x_j = \mathbb{P}(X(t) = j) \quad j \in \{1,...,n\}\)。

根据香农的理论，通道容量由 \(X\) 和 \(Y\) 之间的最大可能互信息 \(I\) 给出：.. code:: ipython3

#!/usr/bin/env python3 # 作者：R. Gowers, S. Al-Izzi, T. Pollington, R. Hill & K. Briggs

import cvxpy as cp import numpy as np import math from scipy.special import xlogy

def channel_capacity(n, m, P, sum_x=1):
    '''
    通信通道的容量。
    Boyd和Vandenberghe的《凸优化》，练习4.57，第207页

    我们考虑一个通信通道，输入X(t)∈{1,..,n}，输出Y(t)∈{1,...,m}，t=1,2,...。
    输入和输出之间的关系是统计给定的：
    p_(i,j) = ℙ(Y(t)=i|X(t)=j)，i=1,..,m  j=1,...,n

    P矩阵 ∈ ℝ^(m*n) 被称为信道过渡矩阵，并且该信道被称为离散无记忆信道。假设X具有概率分布，表示为x ∈ ℝ^n，即，
    x_j = ℙ(X=j)，j=1,...,n。

    X和Y之间的互信息由下式给出：
    ∑(∑(x_j p_(i,j)log_2(p_(i,j)/∑(x_k p_(i,k)))))

    然后，信道容量C由下式给出：
    C = sup I(X;Y)。
    使用y = Px进行变量变换，得到：
    I(X;Y)=  c^T x - ∑(y_i log_2 y_i)
    其中 c_j = ∑(p_(i,j)log_2(p_(i,j)))
    '''

    # n是不同输入值的数量
    # m是不同输出值的数量
    if n*m == 0:
        print('输入和输出值的范围必须大于零')
        return 'failed', np.nan, np.nan

    # x是输入信号X(t)的概率分布
    x = cp.Variable(shape=n)

    # y是输出信号Y(t)的概率分布
    # P是信道过渡矩阵
    y = P@x

    # I是x和y之间的互信息
    c = np.sum(np.array((xlogy(P, P) / math.log(2))), axis=0)
    I = c@x + cp.sum(cp.entr(y) / math.log(2))

    # 通过最大化互信息进行最大化信道容量
    obj = cp.Maximize(I)
    constraints = [cp.sum(x) == sum_x,x >= 0]

    # 构造并解决问题
    prob = cp.Problem(obj,constraints)
    prob.solve()
    if prob.status=='optimal':
        return prob.status, prob.value, x.value
    else:
        return prob.status, np.nan, np.nan

示例

在这个示例中，我们考虑一个具有两个可能的输入和输出的通信通道，所以:math:n = m = 2。在这种情况下，我们使用以下信道过渡矩阵:

\(P = \pmatrix{0.75,0.25\\0.25,0.75}\)

请注意，:math:`P`的列必须总和为1，且:math:`P`的所有元素必须为正。

np.set_printoptions(precision=3)
n = 2
m = 2
P = np.array([[0.75,0.25],
             [0.25,0.75]])
stat, C, x = channel_capacity(n, m, P)
print('问题状态: ',stat)
print('最优值 C = {:.4g}'.format(C))
print('最优变量 x = \n', x)

问题状态:  optimal
最优值 C = 0.1887
最优变量 x =
 [0.5 0.5]