
科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载本文围绕 CVXPY 官方文档中的经典应用示例 Channel_capacity_BV4.57.rst 展开讲解如何把香农信道容量问题转化为凸优化问题并用 CVXPY 的 DCP 接口与entr熵原子直接求解。读完本文你将掌握离散无记忆信道DMC的建模方法、互信息与信道容量之间的凸优化关系、cp.entr原子的底层原理含其指数锥规范形变换并能独立实现任意规模n×m信道的容量计算。问题背景离散无记忆信道的信道容量考虑一个离散无记忆信道discrete memoryless channel其输入随机变量$X(t) \in {1,2,\dots,n}$输出随机变量 $Y(t) \in {1,2,\dots,m}$即输入可取 $n$ 种符号、输出可取 $m$ 种符号。信道输入与输出之间的统计关系由转移概率唯一刻画$$p_{ij} \mathbb{P}\big(Y(t)i \mid X(t)j\big)$$这些转移概率构成信道转移矩阵$P \in \mathbb{R}^{m \times n}$。假设输入符号服从概率分布 $x \in \mathbb{R}^n$即$$x_j \mathbb{P}(X(t)j), \quad j1,\dots,n$$由香农信息论信道容量 $C$ 定义为输入分布 $x$ 与输出分布 $y$ 之间互信息的最大值$$C \sup_{x} I(X;Y)$$其中互信息可写为输出熵与条件熵之差。由于 $y Px$互信息的表达式为$$I(X;Y) -\sum_{i1}^{m} y_i \log_2 y_i \sum_{j1}^{n}\sum_{i1}^{m} x_j p_{ij}\log_2 p_{ij}$$凸优化建模把容量问题写成 DCP关键观察是$x\log x$ 在 $x \geq 0$ 上是凸函数因此 $-y\log y$即熵项是凹函数互信息关于 $x$ 是凹的。最大化凹函数等价于最小化凸函数故信道容量问题可以写成如下凸优化问题$$\begin{array}{ll} \text{minimize} -I(X;Y) \ \text{subject to} \sum_{i1}^{n} x_i 1, \quad x \succeq 0 \end{array}$$约束 $\sum x_i 1$ 与 $x \succeq 0$ 来自 $x$ 是一个概率分布这一事实。为了让 CVXPY 的表达更简洁示例做了一次变量代换$y Px$$y$ 即输出符号的概率分布并预计算常数向量 $c$$$I(X;Y) c^T x - \sum_{i} y_i \log_2 y_i, \qquad c_j \sum_{i1}^{m} p_{ij}\log_2 p_{ij}$$这样目标函数中与输入分布 $x$ 相关的部分退化为线性项 $c^T x$非线性部分只剩 $-\sum y_i \log_2 y_i$恰好可以用 CVXPY 内置的凹原子entr(y) -y\log y自然对数直接表达再除以 $\log 2$ 完成以 2 为底的对数换算。CVXPY 完整实现channel_capacity 函数示例将整个求解过程封装为一个可复用的函数channel_capacity(n, m, P, sum_x1)完整代码如下#!/usr/bin/env python3 # author: 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_x1): Boyd and Vandenberghe, Convex Optimization, exercise 4.57 page 207 Capacity of a communication channel. We consider a communication channel, with input X(t)∈{1,..,n} and output Y(t)∈{1,...,m}, for t1,2,... .The relation between the input and output is given statistically: p_(i,j) ℙ(Y(t)i|X(t)j), i1,..,m j1,...,n The matrix P ∈ ℝ^(m*n) is called the channel transition matrix, and the channel is called a discrete memoryless channel. Assuming X has a probability distribution denoted x ∈ ℝ^n, i.e., x_j ℙ(Xj), j1,...,n The mutual information between X and Y is given by ∑(∑(x_j p_(i,j)log_2(p_(i,j)/∑(x_k p_(i,k))))) Then channel capacity C is given by C sup I(X;Y). With a variable change of y Px this becomes I(X;Y) c^T x - ∑(y_i log_2 y_i) where c_j ∑(p_(i,j)log_2(p_(i,j))) # n is the number of different input values # m is the number of different output values if n*m 0: print(The range of both input and output values must be greater than zero) return failed, np.nan, np.nan # x is probability distribution of the input signal X(t) x cp.Variable(shapen) # y is the probability distribution of the output signal Y(t) # P is the channel transition matrix y Px # I is the mutual information between x and y c np.sum(np.array((xlogy(P, P) / math.log(2))), axis0) I cx cp.sum(cp.entr(y) / math.log(2)) # Channel capacity maximised by maximising the mutual information obj cp.Maximize(I) constraints [cp.sum(x) sum_x, x 0] # Form and solve problem 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各关键步骤说明输入校验当n*m 0时输入或输出符号集为空函数直接返回失败状态与nan避免后续构造空问题。决策变量x cp.Variable(shapen)是 $n$ 维输入符号概率分布y Px是 $m$ 维输出符号概率分布由转移矩阵的线性作用给出为 CVXPY 的矩阵乘法运算符。常数向量 c利用scipy.special.xlogy(P, P)逐元素计算 $p_{ij}\log p_{ij}$在 $p_{ij}0$ 处由xlogy正确返回 0除以 $\log 2$ 转成以 2 为底再沿axis0即对每个输入符号 $j$ 累加所有输出 $i$求和得到 $c \in \mathbb{R}^n$。目标函数I cx cp.sum(cp.entr(y) / math.log(2))。这里cp.entr(y)逐元素计算 $-y_i\log y_i$自然对数cp.sum沿所有输出累加除以 $\log 2$ 换算为 $\log_2$。由于entr是凹原子且cx是仿射函数最大化I满足 DCP 规则。约束与求解约束[cp.sum(x) sum_x, x 0]保证 $x$ 为概率分布。sum_x参数默认取 1允许灵活调整例如归一化到其他常数。构造cp.Problem(obj, constraints)后调用prob.solve()最后根据prob.status返回求解状态、最优值即信道容量 $C$与最优输入分布 $x$。深入源码cp.entr 原子是如何实现的示例的核心是熵原子cp.entr。该原子定义在 cvxpy/atoms/elementwise/entr.py并在 cvxpy/atoms/init.py 中被导出为顶层 API。从源码可以看出它的关键性质数值计算numeric方法调用-xlogy(x, x)实现逐元素 $-x\log x$并在定义域外结果为nan处返回-np.inf从而在 DCP 求解时自动把不可行点推向可行域。曲率判定is_atom_convex()返回False、is_atom_concave()返回True即entr是凹原子sign_from_args返回符号未知。这决定了它只能出现在最大化目标或凹约束左侧等允许凹函数的 DCP 位置。定义域_domain返回[self.args[0] 0]即 $x \geq 0$ 被自动作为隐式约束加入问题。梯度_grad给出 $-(\log x 1)$ 的逐元素梯度供需要导数信息的求解器如 NLP 求解器使用。仓库中的单元测试 cvxpy/tests/test_atoms.py#L1045-L1055 验证了entr对稠密与稀疏常数输入都给出正确结果例如对x0.5有0.5*log(2)对x2有-2*log(2)且 $x0$ 处值为 0。底层原理entr 如何被规范化为指数锥entr之所以能被现代锥求解器高效处理是因为 CVXPY 的 DCP 到锥dcp2cone规约把它改写为指数锥约束。见 cvxpy/reductions/dcp2cone/canonicalizers/entr_canon.pydef entr_canon(expr, args, solver_contextNone): x args[0] shape expr.shape t Variable(shape) # -x\log(x) t x\exp(t/x) 1 ones Constant(np.ones(shape)) constraints [ExpCone(t, x, ones)] return t, constraints即约束 $-x\log x \geq t$ 等价于 $x e^{t/x} \leq 1$这正是指数锥ExpCone(t, x, ones)的定义。指数锥本身实现在 cvxpy/constraints/exponential.py类ExpCone要求三个参数均为仿射、实值且形状一致。因此用entr建模后任何支持指数锥的求解器如 SCS、CLARABEL、ECOS 等都可以直接求解该问题。运行示例2×2 对称信道示例考虑输入、输出符号数均为 2 的信道$n m 2$信道转移矩阵取$$P \begin{pmatrix} 0.75 0.25 \ 0.25 0.75 \end{pmatrix}$$注意两个前提条件P 的每一列必须和为 1列随机矩阵保证对任意输入分布输出仍是概率分布且P 的所有元素必须为正。调用与输出如下np.set_printoptions(precision3) n 2 m 2 P np.array([[0.75, 0.25], [0.25, 0.75]]) stat, C, x channel_capacity(n, m, P) print(Problem status: , stat) print(Optimal value of C {:.4g}.format(C)) print(Optimal variable x \n, x)预期输出Problem status: optimal Optimal value of C 0.1887 Optimal variable x [0.5 0.5]结果解读该二元对称信道BSC的容量约为0.1887 bit/符号最优输入分布为均匀分布 $x [0.5, 0.5]$。这与信息论结论一致——对称信道的容量在等概率输入时取得。扩展与验证更大规模信道与求解器选择任意 n×m 信道channel_capacity的参数化设计使其天然支持任意规模的转移矩阵只需传入正确的n、m与形状为(m, n)的P。求解器验证仓库中的熵相关测试 cvxpy/tests/nlp_tests/test_entropy_related.py 展示了两种求解路径的交叉验证同一样本既可用 NLP 求解器 IPOPTsolvercp.IPOPT, nlpTrue求解也可用锥求解器 CLARABELsolvercp.CLARABEL求解两者最优解偏差不超过 $10^{-4}$。这为信道容量问题提供了可靠的求解器冗余验证思路。参数化复用sum_x参数允许在不改动代码的情况下把分布约束从 $\sum x 1$ 改为任意常数便于研究非归一化情形或与其他问题组合。该示例还收录于文档索引 doc/source/examples/index.rst与熵最大化max_entropy、水填充water_filling等应用同属 CVXPY 的高级应用示例集可作为学习凸优化建模与 DCP 规则的标准参考案例。赞分享科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载相关推荐CVXPY 共识优化Consensus Optimization实战用 ADMM 分布式求解可分凸问题CVXPY 共识优化Consensus Optimization实战用 ADMM 分布式求解可分凸问题 本文以 CVXPY 官方示例文档 consensu科学计算CVXPY 实战用乘子法Method of Multipliers求解带等式约束的凸优化问题CVXPY 实战用乘子法Method of Multipliers求解带等式约束的凸优化问题 本文以 CVXPY 官方文档示例 doc/source/ex科学计算Holo1.5-3B性能优化终极指南10个技巧提升模型推理速度和准确率 Holo1.5 3B性能优化终极指南10个技巧提升模型推理速度和准确率 Holo1.5 3B作为一款先进的视觉语言模型VLM专为计算机使用代理设计上一篇炉石传说终极插件指南50功能完整配置与实战教程下一篇如何在5分钟内掌握Mermaid Live Editor免费在线图表制作终极指南创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考