如何向量化CVXPY代码加速编译？解决大规模K值运行问题

核心问题是代码用循环逐个构建子表达式，导致CVXPY编译时需要处理大量零散计算节点，当K超过40时复杂度急剧上升。按照警告提示，向量化代码是唯一有效的解决途径，具体优化步骤如下：

1. 将字典转换为矩阵格式 #

把e_gp、f_gp从字典转为numpy二维数组，r_gp转为一维数组，利用矩阵乘法完成批量计算，避免循环生成子表达式：

# 转换e_gp为K×K矩阵
E = np.array([[e_gp[(k, kp)] for kp in range(K)] for k in range(K)])
# 转换f_gp为K×K矩阵
F = np.array([[f_gp[(k, kp)] for kp in range(K)] for k in range(K)])
# 转换r_gp为K维数组
R = np.array([r_gp[k] for k in range(K)])

2. 向量化计算t_e、t_f #

用矩阵与向量的乘法替换循环求和，直接得到整体向量结果，大幅减少子表达式数量：

# t_e是排除自身的干扰和：等价于(E @ p)减去自身项
t_e = E @ p - np.diag(E) * p
# t_f是全部项的和，直接用矩阵乘法
t_f = F @ p
# t_r直接用转换后的数组R
t_r = R

3. 简化约束表达式 #

保持约束的向量化写法，避免拆分单个元素处理：

S = (t_e + t_f + t_r) / p
constraints = [p >= p_min, p <= p_max, S <= 1/sinr_t]

完整优化后代码 #

import numpy as np
import cvxpy as cp

K = 50  # 可测试大K值

# 原始字典数据（实际使用时可按K扩展）
e_gp = {(0, 0): 0.01, (0, 1): 0.01, (0, 2): 0.01, (0, 3): 0.01, (0, 4): 0.01, (1, 0): 0.01, (1, 1): 0.01, (1, 2): 0.01, (1, 3): 0.01, (1, 4): 0.01, (2, 0): 0.01, (2, 1): 0.01, (2, 2): 0.01, (2, 3): 0.01, (2, 4): 0.01, (3, 0): 0.01, (3, 1): 0.01, (3, 2): 0.01, (3, 3): 0.01, (3, 4): 0.01, (4, 0): 0.01, (4, 1): 0.01, (4, 2): 0.01, (4, 3): 0.01, (4, 4): 0.01}
f_gp = {(0, 0): 0.0022083736662195513, (0, 1): 1.6361300004264737e-07, (0, 2): 6.673065375868695e-07, (0, 3): 1.8216969961107378e-05, (0, 4): 7.747569454751729e-06, (1, 0): 0.0007991844189121604, (1, 1): 0.0958477227950784, (1, 2): 0.0006366120217946077, (1, 3): 0.005370938447086471, (1, 4): 0.0007190951574825933, (2, 0): 0.0001461922258284424, (2, 1): 8.781670998829123e-06, (2, 2): 0.013922894884250001, (2, 3): 0.0007882566113873586, (2, 4): 5.609072687532707e-05, (3, 0): 1.6647180495599578e-05, (3, 1): 9.701212848634301e-06, (3, 2): 8.415947875397424e-05, (3, 3): 0.002340059870467834, (3, 4): 6.260839608528006e-06, (4, 0): 0.002396438491055845, (4, 1): 2.777751959882801e-05, (4, 2): 0.00011551904109760868, (4, 3): 0.0005849270607107239, (4, 4): 0.013415428308115004}
r_gp = {0: 3.9230406757174874e-05, 1: 0.04336223814165744, 2: 0.0009650237242044663, 3: 5.605495010604486e-05, 4: 0.0008930969228883886}

# 转换为矩阵/数组格式（兼容K扩展的情况）
E = np.array([[e_gp[(k, kp)] if (k, kp) in e_gp else 0.0 for kp in range(K)] for k in range(K)])
F = np.array([[f_gp[(k, kp)] if (k, kp) in f_gp else 0.0 for kp in range(K)] for k in range(K)])
R = np.array([r_gp[k] if k in r_gp else 0.0 for k in range(K)])

sinr_t = cp.Variable(pos=True)
p = cp.Variable(shape=(K,), pos=True)
p_min = 1e-10 * np.ones(K)
p_max = 1.0 * np.ones(K)

# 向量化计算
t_e = E @ p - np.diag(E) * p
t_f = F @ p
t_r = R

objective = cp.Maximize(sinr_t)
S = (t_e + t_f + t_r) / p
constraints = [p >= p_min, p <= p_max, S <= 1/sinr_t]

problem = cp.Problem(objective, constraints)
sol = problem.solve(solver=cp.MOSEK, gp=True, verbose=False)

额外优化建议 #

• 若K非常大（如超过100），可预先计算排除自身的干扰矩阵，避免重复操作：

  E_no_self = E - np.diag(np.diag(E))
  t_e = E_no_self @ p
  

• 确保CVXPY为最新版本，新版对向量化表达式的编译优化更高效。

内容的提问来源于stack exchange，提问作者José Rezende