CVXPY 半定规划（SDP）实战：从标准形式到协方差矩阵补全

科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载半定规划Semidefinite Program, SDP是凸优化中最重要的模型类别之一它在线性规划LP的基础上引入了“矩阵半正定”这一锥约束被广泛应用于协方差估计、组合优化松弛、控制理论与鲁棒优化。本指南以 CVXPY 官方入门示例 doc/source/examples/basic/sdp.rst 为骨架逐步讲解 SDP 的数学标准形式、CVXPY 建模语法矩阵不等式、symmetric变量、trace原子并以随机生成的 SDP 和协方差矩阵补全为例给出可直接复制运行且带源码级原理注解的完整代码。读完本文你将掌握用 CVXPY 声明对称/半正定矩阵变量的正确姿势、X 0与X 0背后的 PSD 约束实现机制、如何用cp.trace高效构建线性矩阵目标与等式约束以及如何把“缺失数据补全”这类实际问题转写成 SDP 并求出数值解。SDP 的标准形式目标、等式约束与半正定锥SDP 是如下形式的优化问题minimize tr(CX) subject to tr(A_i X) b_i, i 1, ..., p X ≽ 0其中tr表示矩阵的迹对角线元素之和X ∈ Sⁿ是优化变量Sⁿ表示所有 n×n 对称矩阵构成的集合C, A_1, ..., A_p ∈ Sⁿ与b_1, ..., b_p ∈ R是给定的问题数据X ≽ 0表示 X 半正定positive semidefinite即对任意向量 z 有 zᵀXz ≥ 0这是一个矩阵不等式matrix inequality。可以看到SDP 与线性规划结构同构目标函数与约束都是变量的线性函数唯一的区别在于可行域从“非负象限”换成了“半正定锥”Sⁿ₊。正因为半正定锥是自对偶的凸锥SDP 拥有完整的对偶理论与多项式时间的内点算法这是它可以被现代求解器高效处理的理论基础。从 CVXPY 源码看半正定约束正是作为一个独立的锥约束类存在的cvxpy/constraints/psd.py 中的PSD(Cone)类实现了X ≽ 0约束其文档字符串明确说明对二维表达式 X 施加 PSD 约束等价于约束其对称部分(X Xᵀ)/2 ≽ 0。该类的is_dcp方法要求被约束的表达式必须是仿射的cvxpy/constraints/psd.py这也是 DCP 规则对锥约束的基本要求。一个经典应用协方差矩阵补全SDP 的一个典型应用是协方差矩阵补全covariance matrix completion。假设已知一个带缺失项的协方差矩阵Σ̃ ∈ Sⁿ₊缺失位置构成集合M ⊂ {1,...,n} × {1,...,n}我们希望补全出一个完整且半正定的协方差矩阵Σ同时保持已知位置的值不变minimize 0 subject to Σ_ij Σ̃_ij, (i,j) ∉ M Σ ≽ 0目标函数取常数 0即我们只关心“是否存在一个半正定矩阵与已知观测一致”这是一个典型的可行性问题feasibility problem。约束分两部分已知位置的元素必须等于观测值等式约束整个矩阵必须半正定锥约束。由于半正定本身就是“协方差矩阵”这一统计概念的数学刻画这个建模方式是金融风险分析、传感器网络、推荐系统等领域中矩阵补全问题的标准起点。求解得到Σ后缺失位置(i,j) ∈ M的补全值直接从Σ.value中读出即可。在 CVXPY 中建模并求解一个随机 SDP生成问题数据下面这段代码构造了一个 n3、p3 的随机 SDP随机生成对称问题数据 C 与 A_i以及随机右端项 b_i# Import packages. import cvxpy as cp import numpy as np # Generate a random SDP. n 3 p 3 np.random.seed(1) C np.random.randn(n, n) A [] b [] for i in range(p): A.append(np.random.randn(n, n)) b.append(np.random.randn())固定np.random.seed(1)是为了保证结果可复现——官方文档给出的最优值 2.654348003008652 正是基于该种子。注意这里 C 和 A_i 虽然由randn生成、本身不一定对称但迹函数tr(AX)只与 A 的对称部分有关因此用于构造tr(A_i X)时并不要求 A_i 预先对称。声明对称矩阵变量与 PSD 约束接下来是 SDP 建模的核心三行# Define and solve the CVXPY problem. # Create a symmetric matrix variable. X cp.Variable((n, n), symmetricTrue) # The operator denotes matrix inequality. constraints [X 0] constraints [ cp.trace(A[i] X) b[i] for i in range(p) ] prob cp.Problem(cp.Minimize(cp.trace(C X)), constraints) prob.solve()这里的关键点在于cp.Variable((n, n), symmetricTrue)symmetricTrue告诉 CVXPY 该变量限定在对称矩阵空间内。从源码看Leaf.__init__会校验只要声明了PSD、NSD、symmetric、diag、hermitian中的任一属性形状必须是方阵否则抛出ValueErrorcvxpy/expressions/leaf.py更关键的是这些“降维属性”dim_reducing_attr彼此互斥——一个变量不能同时声明symmetric和PSDcvxpy/expressions/leaf.py。所以正确做法是变量本身只声明对称性半正定性通过X 0约束施加Leaf.is_symmetric()的实现显示标量或带diag/symmetric/PSD/NSD属性的变量都会被视为对称cvxpy/expressions/leaf.py声明symmetricTrue能显著减少决策变量的自由度从 n² 降到 n(n1)/2降低求解规模。X 0这个语法不是数学符号而是 Python 运算符重载。在 cvxpy/expressions/expression.py 中Expression定义了四个矩阵不等式运算符运算符含义源码实现X 0X 半正定PSDPSD(X - 0)0 XX 半正定PSDPSD(X - 0)X 0X 半负定NSDPSD(0 - X)0 XX 半负定NSDPSD(0 - X)也就是说X 0会被转换成PSD(X - 0)这个约束对象。注意CVXPY不提供严格正定X ≻ 0约束——cvxpy/constraints/psd.py 的文档明确说明在数值求解语境下严格定号没有意义因为浮点求解器只能逼近锥边界。另外PSD 约束要求表达式必须是方阵非方阵会直接报ValueErrorcvxpy/constraints/psd.py。用 trace 原子构造目标与等式约束cp.trace(C X)与cp.trace(A[i] X)使用了迹原子。trace的完整定义在 cvxpy/atoms/affine/trace.py其中包含一个值得注意的数值优化细节当参数是矩阵乘法表达式MulExpression时trace(A B)会被改写为sum(A * B.T)逐元素相乘后求和利用恒等式tr(A B) Σ_ij A_ij·B_ji把 O(n³) 的完整矩阵乘法降到 O(n²) 的对角线计算否则回退到通用的Trace原子其数值实现等价于np.linalg.trace并校验参数必须是方阵最后两维相等且维度 ≥ 2cvxpy/atoms/affine/trace.py。因此文档示例中的C X、A[i] X会自动走这条高效路径。由于 X 是对称变量、trace(C X)与trace(A_i X)均为仿射表达式整个问题满足 DCP 规则。完整运行与输出把上面三段代码合并就是官方示例的完整程序。求解后打印结果# Print result. print(The optimal value is, prob.value) print(A solution X is) print(X.value)在np.random.seed(1)下输出为The optimal value is 2.654348003008652 A solution X is [[ 1.6080571 -0.59770202 -0.69575904] [-0.59770202 0.22228637 0.24689205] [-0.69575904 0.24689205 1.39679396]]两点观察值得留意解矩阵严格对称X.value满足X[i, j] X[j, i]这是symmetricTrue与 PSD 锥共同作用的结果。CVXPY 测试套件中也专门用assert (M.value M.T.value).all()来验证这一性质cvxpy/tests/test_semidefinite_vars.py默认求解器prob.solve()未指定求解器时CVXPY 会按内部优先级自动挑选可用的锥求解器通常优先 SCS/CLARABEL 等。如需显式指定可传入prob.solve(solvercp.SCS)或prob.solve(solvercp.CLARABEL)。源码视角PSD 约束在求解链中的位置从建模到求解X 0约束会经历一条完整的变换链。以 cvxpy/reductions/dcp2cone 为例DCP 问题先被转换为标准的锥规划cone program形式在dcp2cone的 canonicalizers 中凡是需要半正定锥的原子——例如lambda_max使用PSD(t_eye - A)cvxpy/reductions/dcp2cone/canonicalizers/lambda_max_canon.py、log_det使用PSD(X)cvxpy/reductions/dcp2cone/canonicalizers/log_det_canon.py、matrix_frac与sigma_max同样如此——都会显式构造PSD约束对象。在更下游的cone2cone精确变换中cvxpy/reductions/cone2cone/exact.py 会根据目标求解器的接口要求把 PSD 约束原样传递或通过SvecPSDcvxpy/constraints/psd.py转换为长度n(n1)/2的缩放向量化svec表示——这是 SCS 等求解器期望的半正定锥输入格式。PSD约束还实现了residual、dual_residual、dual_violation等属性cvxpy/constraints/psd.py用于评估约束违反程度与对偶可行性这也是调试求解结果时常用的工具。把示例改造成协方差补全问题官方示例展示的是“随机数据 非零目标”的通用 SDP。把它改写成前述协方差补全问题只需三步固定已知子矩阵构造一个带缺失项的对称矩阵Sigma_tilde缺失位置用np.nan占位替换约束对每个(i, j) ∉ M添加Sigma[i, j] Sigma_tilde[i, j]替换目标将目标改为cp.Minimize(0)可行性问题或改为最小化补全矩阵的某种正则项如cp.Minimize(cp.trace(Sigma))鼓励低迹解。一个可直接运行的骨架如下import cvxpy as cp import numpy as np n 3 # 带缺失项的观测矩阵nan 表示缺失位置 Sigma_tilde np.array([ [1.0, np.nan, 0.5], [np.nan, 2.0, 0.1], [0.5, 0.1, 3.0], ]) Sigma cp.Variable((n, n), symmetricTrue) constraints [Sigma 0] for i in range(n): for j in range(n): if not np.isnan(Sigma_tilde[i, j]): constraints.append(Sigma[i, j] Sigma_tilde[i, j]) prob cp.Problem(cp.Minimize(0), constraints) prob.solve() print(Completed covariance matrix:\n, Sigma.value)注意这里同样遵循“变量只声明symmetricTrue半正定用约束”的规则。运行后Sigma.value即为一个与已知观测一致、且半正定的完整协方差矩阵缺失位置(i, j) ∈ M的值就是补全结果。小结SDP 建模要点速查数学形式SDP 线性目标 线性矩阵迹等式约束 X ≽ 0锥约束变量位于对称矩阵空间Sⁿ变量声明对称矩阵用cp.Variable((n, n), symmetricTrue)symmetric、PSD、NSD、diag等属性互斥不可叠加锥约束语法X 0表示半正定X 0表示半负定二者在 cvxpy/expressions/expression.py 中统一转换为PSD约束CVXPY 不提供严格定号约束迹运算cp.trace(A X)会自动改写为 O(n²) 的sum(A * B.T)计算见 cvxpy/atoms/affine/trace.py典型应用协方差矩阵补全、lambda_max/log_det/matrix_frac/sigma_max等原子的锥表示以及 SDP 松弛类问题求解器默认自动选择可用锥求解器也可通过prob.solve(solver...)显式指定。若想继续深入可以阅读同一目录下的姊妹示例 线性规划、二次规划 与 二阶锥规划它们与 SDP 一起构成了 CVXPY 锥规划家族的完整入门路线。赞分享科学计算【免费下载链接】cvxpyA Python-embedded modeling language for convex optimization problems.项目地址https://gitcode.com/gh_mirrors/cv/cvxpy点击查看免费下载相关推荐time_varying_optimization 详解基于 CVXPY 的时变半定规划TV-SDP框架与多项式矩阵 PSD 约束time_varying_optimization 详解基于 CVXPY 的时变半定规划TV SDP框架与多项式矩阵 PSD 约束 导读 time_var人工智能深度学习NLP计算机视觉强化学习使用CVXPY求解半定规划(SDP)问题详解使用CVXPY求解半定规划 SDP 问题详解 什么是半定规划 SDP 半定规划 Semidefinite Programming, SDP 是凸优化领域中的一个科学计算YimMenu Lua 脚本开发gui 表格 API 全解析与实战指南YimMenu Lua 脚本开发gui 表格 API 全解析与实战指南 gui 是 YimMenu Lua 脚本系统中负责 操控菜单 GUI 本体 的核心表格科学计算上一篇ImageGlass图片浏览器终极指南免费开源支持90格式的高效工具下一篇Adobe-GenP终极指南三步解锁Adobe全家桶的完整教程创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考