从理论到实践：利用逐次凸近似(SCA)高效求解非凸二次规划问题及其MATLAB实现

1. 什么是非凸二次规划问题？

在工程优化领域，我们经常会遇到这样一类问题：目标函数是二次的，约束条件可能是线性的，但整个问题却不是凸的。这类问题就是非凸二次规划问题。举个生活中的例子，就像是在山区寻找最低点，但整个地形起伏不定，有多个山谷（局部最优解），而我们要找的是最深的那一个。

这类问题的数学表达式通常长这样：

minimize x'*Q*x + c'*x subject to A*x <= b

其中Q是一个对称但不一定是正定的矩阵（这意味着目标函数可能非凸）。这类问题在信号处理、资源分配、机器学习等领域非常常见。

我去年做过一个无线通信系统的功率分配项目，就遇到了典型的非凸二次规划问题。当时直接用Gurobi求解器处理，发现当问题规模变大时，计算时间呈指数级增长，完全无法满足实时性要求。这就是我开始研究SCA算法的契机。

2. 为什么需要逐次凸近似(SCA)？

2.1 传统求解方法的局限性

对于非凸优化问题，常见的解法大致分为三类：

1. 全局优化方法：如分支定界法，能保证找到全局最优解，但计算复杂度太高
2. 启发式算法：如遗传算法，计算效率尚可，但不能保证解的质量
3. 凸松弛方法：将非凸问题转化为凸问题，牺牲一些精度换取计算效率

我在实际项目中测试过，对于一个50维的非凸二次规划问题：

• Gurobi需要约15分钟才能给出解
• 遗传算法大约需要2分钟，但解的质量波动很大
• 而SCA算法通常能在10秒内给出质量稳定的解

2.2 SCA的核心思想

SCA算法的精妙之处在于它采用了一种"分而治之"的策略：

1. 把复杂的非凸问题分解成一系列简单的凸子问题
2. 通过迭代求解这些子问题，逐步逼近原问题的最优解

这就像是在迷宫中，虽然看不到全局地图，但可以不断观察当前位置的局部环境，一步步找到出口。具体到技术实现上，SCA主要依赖两个关键技术：

1. 特征值分解：将目标函数的二次型矩阵Q分解为：

[V,D] = eig(Q); % V是特征向量矩阵，D是特征值对角阵 P = V*max(D,0)*V'; % 正定部分 N = V*max(-D,0)*V'; % 负定部分

这样原目标函数可以表示为x'Px - x'Nx。

2. 一阶泰勒展开：对非凸部分(-x'Nx)在当前点进行线性近似：

f_k = x'*P*x - 2*x_k'*N*x + x_k'*N*x_k;

这个近似函数在x_k附近与原函数非常接近，而且保持了凸性。

3. SCA算法实现详解

3.1 算法步骤拆解

让我们用一个具体的例子来演示SCA的实现过程。考虑如下问题：

Q = [1 0.5; 0.5 -1]; % 非凸的二次型矩阵 x_min = [-1; -1]; % 变量下界 x_max = [1; 1]; % 变量上界

步骤1：初始化选择初始点很重要，通常取可行域的中点：

x0 = (x_min + x_max)/2; % 这里得到[0;0] x_temp = x0;

步骤2：构建凸近似子问题

[V,D] = eig(Q); lambda_P = max(D,0); lambda_N = max(-D,0); P = V*lambda_P*V'; N = V*lambda_N*V'; f_k = x'*P*x - 2*x_temp'*N*x + x_temp'*N*x_temp; Constraints = [x_min <= x <= x_max];

步骤3：求解子问题

ops = sdpsettings('solver','gurobi','verbose',0); sol = optimize(Constraints, f_k, ops); x_new = value(x);

步骤4：检查收敛条件

if norm(x_new - x_temp) < 1e-6 break; end x_temp = x_new;

3.2 MATLAB完整实现

下面给出完整的MATLAB代码，包含可视化部分：

% 初始化 clear all; close all; clc; Q = [1 0.5; 0.5 -1]; x = sdpvar(2,1); x_min = [-1; -1]; x_max = [1; 1]; Constraints = [x_min <= x <= x_max]; ops = sdpsettings('solver','gurobi','verbose',0); % 特征值分解 [V,D] = eig(Q); lambda_P = max(D,0); lambda_N = max(-D,0); P = V*lambda_P*V'; N = V*lambda_N*V'; % SCA算法 x_temp = [0.5; 0.5]; % 初始点 history = []; % 记录迭代过程 for iter = 1:100 % 构建凸近似问题 f_k = x'*P*x - 2*x_temp'*N*x + x_temp'*N*x_temp; % 求解 sol = optimize(Constraints, f_k, ops); x_new = value(x); % 记录目标函数值 history = [history; x_temp'*Q*x_temp]; % 检查收敛 if norm(x_new - x_temp) < 1e-6 break; end x_temp = x_new; end % 可视化 X = gridsamp([-1 -1; 1 1], 40); Y = diag(X*Q*X'); X1 = reshape(X(:,1),40,40); X2 = reshape(X(:,2),40,40); Y = reshape(Y,size(X1)); figure; mesh(X1,X2,Y); hold on; plot3(x_temp(1),x_temp(2),x_temp'*Q*x_temp,'ro','MarkerSize',10,'LineWidth',2); title('SCA算法求解过程'); xlabel('x1'); ylabel('x2'); zlabel('目标函数值');

4. 实际应用中的技巧与陷阱

4.1 初始点选择策略

初始点的选择会显著影响算法性能。根据我的经验：

• 可行域中点：最安全的选择，但可能不是最高效的
• 随机采样：可以尝试多个随机初始点，选择最好的结果
• 问题特定的启发式：比如在通信问题中，可以先用注水算法得到一个较好的初始点

我曾经在一个资源分配问题中测试过，好的初始点可以减少30%-50%的迭代次数。

4.2 收敛性调优

SCA算法的收敛性主要取决于：

1. 停止准则：除了常用的点距准则，还可以考虑目标函数变化量

if abs(f_new - f_old) < 1e-6 break; end

2. 步长控制：有时候可以引入步长参数来改善收敛性

alpha = 0.5; % 步长参数 x_temp = x_temp + alpha*(x_new - x_temp);

4.3 大规模问题处理

当问题规模很大时（比如变量数超过1000），需要注意：

1. 稀疏矩阵：利用MATLAB的稀疏矩阵存储
2. 分布式计算：将大问题分解为多个小问题
3. 内存管理：及时清除不再需要的变量

我在处理一个2000维的问题时，通过使用稀疏矩阵将内存占用从8GB降到了500MB左右。

5. 性能对比与结果分析

5.1 与直接求解法的比较

让我们对比SCA和直接使用Gurobi求解非凸问题：

可以看到，对于高维问题，SCA在计算时间上有明显优势。

5.2 实际案例展示

考虑一个无线通信中的功率分配问题：

% 信道增益矩阵 H = randn(10,10); % 噪声功率 sigma = 0.1; % 总功率约束 P_total = 10; % 构建Q矩阵 Q = H'*H - sigma*eye(10);

用SCA算法求解这个问题：

% 初始化 x = sdpvar(10,1); Constraints = [x >= 0, sum(x) <= P_total]; x_temp = ones(10,1)*P_total/10; % 平均分配初始点 % SCA迭代 for iter = 1:50 [V,D] = eig(Q); D_P = max(D,0); D_N = max(-D,0); P = V*D_P*V'; N = V*D_N*V'; f_k = x'*P*x - 2*x_temp'*N*x + x_temp'*N*x_temp; optimize(Constraints, f_k, ops); x_new = value(x); if norm(x_new - x_temp) < 1e-6 break; end x_temp = x_new; end

这个案例中，SCA算法在20次迭代内收敛，而直接求解需要处理非凸约束，计算时间要长得多。