CVXPY编译缓慢且报‘子表达式过多’警告的技术优化咨询

问题背景 #

使用CVXPY求解简单无约束二次优化问题时，矩阵规模较小时可正常求解，但规模增大后编译耗时剧增。例如：

• N_t=1时，编译+求解总耗时约12秒（求解器仅耗时30ms）
• N_t=2时，总耗时增至200秒

简化版Python代码如下：

import numpy as np
import cvxpy as cp
from scipy.linalg import dft
import time as time


start_time = time.time()

N_r = 8
N_t = 1
N_s = 10

x = np.random.randint(2, size = N_s*N_r*N_t)

H = dft(N_s*N_r*N_t)*1j
H_H = H.conj().T

err = np.random.random(N_r)
Err = np.kron(np.diag(np.ones(N_s)), np.kron(np.diag(np.ones(N_t)), np.diag(err)))

y = H_H@Err@H@x

var = cp.Variable(shape = (N_r), complex = True)
Err_est = cp.kron(np.diag(np.ones(N_s)), cp.kron(np.diag(np.ones(N_t)), cp.diag(var)))

res = cp.sum_squares(H_H@Err_est@H@x - y)

prob = cp.Problem(cp.Minimize(res))
prob.solve(solver = cp.CVXOPT, verbose = True)

print(f'Result: {np.abs(err-var.value)}')
print(f'\nProcessing duration: {time.time()-start_time:.2f} seconds')

当调整N_t、N_r或N_s增大矩阵规模时，CVXPY会弹出以下警告：

UserWarning: Objective contains too many subexpressions. Consider vectorizing your CVXPY code to speed up compilation.
warnings.warn("Objective contains too many subexpressions. ")

UserWarning: Constraint #0 contains too many subexpressions. Consider vectorizing your CVXPY code to speed up compilation.
warnings.warn(f"Constraint #{i} contains too many subexpressions. ")

查看CVXPY源码发现debug_tools.py中设置了MAX_NODES=10000的阈值，但不清楚节点计数规则及如何减少节点以支持N_r=20、N_t=10的场景。

同时用Matlab CVX实现了相同问题，使用SDPT3求解器，耗时远低于Python版本：

tic
N_r = 8;
N_t = 1;
N_s = 10;

x = randi([0 1], N_s*N_r*N_t,1);

H = dftmtx(N_s*N_r*N_t)*1i;
H_H = H';

err = randn(1,N_r);
Err = kron(diag(ones(1,N_s)), kron(diag(ones(N_t,1)), diag(err)));

y = H_H*Err*H*x;

cvx_begin
    variable varx(N_r)
    Err_est = kron(diag(ones(1,N_s*N_t)),diag(varx));
    
    minimize(sum_square_abs(H_H*Err_est*H*x - y,1)) 
cvx_end

disp('Result: '); 
disp(err.'-varx)

toc

核心问题 #

• 为何Matlab实现速度快这么多？
• CVXPY的节点是如何计算的？
• 如何在不改变问题的前提下减少节点数量？
• 未来如何避免此类问题？
• 如何降低编译耗时，达到Matlab CVX的编译速度？

用户尝试将cp.sum_squares()替换为cp.norm()，编译耗时有所降低，但大场景下（如N_r=16、N_t=8）仍需约4000秒，Matlab仅需11秒，且Python内核易崩溃。已在CVXPY GitHub仓库提交Issue#2205，但当前及未来版本无法解决，寻求替代方案。

解决方案 #

1. Matlab CVX速度更快的原因 #

• 符号处理引擎差异：Matlab CVX基于成熟的Maple符号计算内核，在表达式化简、节点合并上的优化更彻底，能高效处理大规模kronecker乘积和二次型展开。
• 编译机制不同：Matlab CVX编译阶段直接生成紧凑的问题描述，避免CVXPY中大量子表达式的冗余存储与遍历；CVXPY依赖Python动态特性，子表达式的生成与管理开销更高。
• 底层实现优化：Matlab矩阵运算本身经过高度优化，CVX生成求解器输入时的矩阵操作效率优于CVXPY的纯Python实现。

2. CVXPY节点计数规则 #

CVXPY的节点指表达式树中的每个操作单元：

• 每个变量、常数是一个节点
• 每个运算（如@、kron、sum_squares）会生成新节点，且每个子表达式会被拆分为多个节点
• MAX_NODES=10000是编译性能预警阈值，节点数超过该值时，编译时间会呈指数增长

例如，多层kron操作会生成大量中间子表达式节点，二次型展开也会产生大量元素级运算节点。

3. 不改变问题前提下减少节点数量 #

（1）预计算常数矩阵，避免动态生成子表达式

将原代码中H_H@Err_est@H@x拆解为常数矩阵与变量的线性组合，预计算合并后的矩阵，避免CVXPY动态生成多层kron和矩阵乘法的子节点：

# 预计算常数矩阵，将目标转化为A@var的形式
Hx = H @ x
# 构造kron(ones(Ns*Nt), var)的展开矩阵
kron_var_vec = np.kron(np.ones(N_s*N_t), np.eye(N_r)).T  # 形状(N_r, N_s*N_r*N_t)
A = H_H @ (Hx * kron_var_vec)  # 形状(N_s*N_r*N_t, N_r)
res = cp.sum_squares(A @ var - y)

（2）直接构造对角向量替代多层kron

原代码中Err_est是块对角矩阵，直接构造对角向量而非多次调用kron：

err_est_diag = np.kron(np.ones(N_s*N_t), var)
res = cp.sum_squares(H_H @ (err_est_diag * Hx) - y)

这种方式减少了kron操作生成的子节点，同时利用元素级乘法替代矩阵乘法，降低节点数量。

（3）使用cp.norm替代cp.sum_squares

cp.norm(..., 2)的内部实现更紧凑，生成的节点数更少，已验证能降低编译耗时。

4. 未来避免此类问题的方法 #

• 优先预计算常数部分：将所有与变量无关的矩阵运算提前用NumPy完成，只保留变量的线性/二次组合部分给CVXPY处理。
• 避免嵌套kron操作：块对角矩阵尽量直接构造对角向量，而非通过多次kron生成矩阵。
• 选择低开销的表达式形式：优先使用向量运算代替矩阵运算，用cp.norm代替cp.sum_squares，减少子表达式数量。
• 显式构造二次型：对于二次优化问题，直接推导目标函数的二次型矩阵（P）和线性项（q），构造0.5*cp.quad_form(var, P) + q.T @ var + c的形式，完全避免复杂表达式生成。

5. 达到Matlab CVX编译速度的替代方案 #

（1）使用CVXOPT直接构造QP问题

跳过CVXPY编译层，手动推导目标函数的二次型矩阵，用CVXOPT原生接口求解：

import cvxopt as cvx
# 目标函数sum_squares(A@var - y) = var'*A'*A*var - 2y'*A*var + y'y
P = cvx.matrix(A.conj().T @ A)
q = cvx.matrix(-2 * y.conj().T @ A)
sol = cvx.solvers.qp(P, q)
var_value = np.array(sol['x']).flatten()

（2）利用闭式解求解

该问题是无约束二次优化，可直接用闭式解计算，完全不需要调用优化求解器：

var_value = np.linalg.inv(A.conj().T @ A) @ A.conj().T @ y

（3）使用其他Python优化库

如scipy.optimize.minimize（针对无约束问题）、jaxopt（基于JAX的自动微分优化库，编译速度更快）。

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