一、本文概述：

本文给出对偶上升法（dual ascent）求解凸优化问题最优解的代码实例。如果您觉得对您有帮助，请点个赞，加个收藏，谢谢！

二、简单实例

本文以下述实例为例，撰写对偶上升法的迭代步骤，并给出最终可运行的MATLAB代码，以便大家上手学习。
1）优化问题为：
m i n x 1 , x 2   f ( x 1 , x 2 ) = 2 ( x 1 − 1 ) 2 + ( x 2 + 2 ) 2 s . t . { x 1 ≥ 2 x 2 ≥ 0 \begin{aligned} &\mathop{min}\limits_{x_1,x_2} \ f(x_1,x_2)=2(x_1-1)^2+(x_2+2)^2\\ &s.t. \begin{cases} x_1\geq2 \\ x_2\geq0 \end{cases} \end{aligned} ​x1​,x2​min​ f(x1​,x2​)=2(x1​−1)2+(x2​+2)2s.t.{x1​≥2x2​≥0​​

2）拉格朗日函数为：
L ( x 1 , x 2 , λ 1 , λ 2 ) = 2 ( x 1 − 1 ) 2 + ( x 2 + 2 ) 2 + λ 1 ( 2 − x 1 ) + λ 2 ( − x 2 ) \begin{aligned} &L(x_1,x_2,\lambda_1,\lambda_2) \\ &=2(x_1-1)^2+(x_2+2)^2+\lambda_1(2-x_1)+\lambda_2(-x_2)\\ \end{aligned} ​L(x1​,x2​,λ1​,λ2​)=2(x1​−1)2+(x2​+2)2+λ1​(2−x1​)+λ2​(−x2​)​

3）决策变量 x 1 , x 2 x_1,x_2 x1​,x2​解析解为：
∂ L ( x 1 , x 2 , λ 1 , λ 2 ) ∂ x 1 = 0 ∂ L ( x 1 , x 2 , λ 1 , λ 2 ) ∂ x 2 = 0 \begin{aligned} \frac{\partial L(x_1,x_2,\lambda_1,\lambda_2) }{\partial x_1}=0 \\ \frac{\partial L(x_1,x_2,\lambda_1,\lambda_2) }{\partial x_2}=0 \end{aligned} ∂x1​∂L(x1​,x2​,λ1​,λ2​)​=0∂x2​∂L(x1​,x2​,λ1​,λ2​)​=0​
求解上述方程，解得：
x 1 ∗ = λ 1 4 + 1 x 2 ∗ = λ 2 2 − 2 \begin{aligned} x_1^*=\frac{\lambda_1}{4}+1 \\ x_2^*=\frac{\lambda_2}{2}-2 \end{aligned} x1∗​=4λ1​​+1x2∗​=2λ2​​−2​

4）拉格朗日函数更新公式：
第k+1轮迭代更新公式：
λ 1 k + 1 = λ 1 k + t 1 ( 2 − x 1 ) λ 2 k + 1 = λ 2 k + t 1 ( − x 2 ) \begin{aligned} \lambda_1^{k+1}&=\lambda_1^{k}+t_1(2-x_1) \\ \lambda_2^{k+1}&=\lambda_2^{k}+t_1(-x_2) \end{aligned} λ1k+1​λ2k+1​​=λ1k​+t1​(2−x1​)=λ2k​+t1​(−x2​)​

至此，对偶上升方法完毕求解完毕，下面为MATLAB代码！

三、MATLAB代码

clear all;
clc;

x = [100 -120]';   %决策变量初始化   %注：本函数的决策变量为向量！ 

size_x_row  = size(x,1);
size_x_line = size(x,2);
lambda = zeros( size_x_row , size_x_line );   %拉格朗日乘子初始化 
k=1;
G_v = [];  lambda_v = [];  x_v = [];  L_v = [];  f_v = [];
tol = 10^(-5);
t_lambda=0.1;  % 拉格朗日乘子的更新步长
t_x=0.1;      % 决策变量 x 的更新步长
max_iteration = 10000;

for i = 1 : max_iteration
    %% 迭代更新决策变量 x
    x(1) = x(1) - t_x * (  4 * ( x(1) - 1 ) - lambda(1)  );  %% 注意这里可以写解析解，也可以写梯度迭代更新解，也可混合写！
    x(2) = x(2) - t_x * (  2 * ( x(2) + 2 ) - lambda(2)  ); 
    %x(1) = lambda(1)/4+1; 
    %x(2) = lambda(2)/2-2; 
    
    %% 迭代更新拉格朗日乘子 lambda 
    Grad_f = [ (2-x(1)) ; (-x(2)) ] ; 
    lambda = max(   lambda + t_lambda * Grad_f  ,  zeros( size_x_row , size_x_line )    );  
    f =  2*( x(1)-1 )^2 + ( x(2)+2 )^2 ; 
    L =  2*( x(1)-1 )^2 + ( x(2)+2 )^2 + lambda(1)*(2-x(1)) + lambda(2)*(-x(2));
    
    %% 记录历史更新结果
    G_v = [G_v Grad_f];
    x_v = [x_v x];
    lambda_v = [lambda_v lambda];
    f_v = [f_v f]; 
    L_v = [L_v L];
    
    %% 终止条件判断
    % if (k>2) &&( abs(L - L_v(k-1)) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
    if (k>2) &&( abs( f - L ) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
        break;
    end
    k=k+1;
end
if i == max_iteration
    disp(['超过最大迭代次数！']);
end
disp(['迭代次数为：',num2str(k),'次']);
disp(['最优值点为：']);
x 
figure(1)
plot( 2:length(L_v),L_v(2:length(L_v)), '-ro', 'linewidth',2)
disp(['最优值为：',num2str(L),'']);

figure(2)
plot( 2:length(L_v), f_v(2:length(f_v)) - L_v(2:length(L_v)) , '-b*', 'linewidth',2)

四、注意事项及扩展

注意：上述第3)步中，若无法得到解析解，则求梯度亦可。换言之，上述第3)步可用下述代码替代：（计算方法为：拉格朗日函数分别对 x 1 x_1 x1​与 x 2 x_2 x2​求偏导即可）
x 1 k + 1 = x 1 k − t 2 [ 4 ( x 1 − 1 ) − λ 1 ] x 2 k + 1 = x 2 k − t 2 [ 2 ( x 2 + 2 ) − λ 2 ] \begin{aligned} x_1^{k+1}&=x_1^{k} - t_2[4(x_1-1)-\lambda_1] \\ x_2^{k+1}&=x_2^{k} - t_2[2(x_2+2)-\lambda_2] \end{aligned} x1k+1​x2k+1​​=x1k​−t2​[4(x1​−1)−λ1​]=x2k​−t2​[2(x2​+2)−λ2​]​
MATLAB代码中将解析解的情况注释掉了，用数值解求解的。换言之代码的“第18行第19行”与代码中的“第20行第21行”是等价替换的，运行哪一组都可以！！！当然，也可以混合运行！！！

五、运行结果

六、基础回顾

这里需要注意，不要把梯度下降前的更新正负号写反了，对拉格朗日乘子更新时，其梯度前是加号（即：加梯度），而对决策变量更新时，其梯度前为减号（即：减梯度）。其原因是，前者为梯度上升，后者为梯度下降。我们回顾一下，拉更朗日函数的目标是求解下述优化问题：

m a x λ 1 , λ 2 [ m i n x 1 , x 2   L ( x 1 , x 2 , λ 1 , λ 2 ) ] \begin{aligned} \mathop{max}\limits_{\lambda_1,\lambda_2} \left[ \mathop{min}\limits_{x_1,x_2} \ L(x_1,x_2,\lambda_1,\lambda_2)\right] \end{aligned} λ1​,λ2​max​[x1​,x2​min​ L(x1​,x2​,λ1​,λ2​)]​
因此，选取拉格朗日乘子 λ 1 , λ 2 \lambda_1,\lambda_2 λ1​,λ2​以最大化 L ( x 1 , x 2 , λ 1 , λ 2 ) L(x_1,x_2,\lambda_1,\lambda_2) L(x1​,x2​,λ1​,λ2​)，因此对拉格朗日乘子做梯度上升，所以是加号，而选取决策变量 x 1 , x 2 x_1,x_2 x1​,x2​以最小化 L ( x 1 , x 2 , λ 1 , λ 2 ) L(x_1,x_2,\lambda_1,\lambda_2) L(x1​,x2​,λ1​,λ2​)，因此对拉格朗日乘子做梯度下降，所以是减号。

另：若有时间，会陆续更新 “精确直线搜索代码” 及 “ADMM的consensus代码”

更多优化内容，欢迎关注本人微信公众号：

七、后期更正

终止条件应是 原函数f 与 拉格朗日函数L 只差足够小。而并非 拉格朗日函数L 自身收敛。

八、变换函数

甚至更换函数，并在0附近，该方法依然可行

m i n x 1 , x 2   f ( x 1 , x 2 ) = 1 x 1 + x 1 2 + 1 x 2 + x 2 2 s . t . { 0 < x 1 ≤ 1 0 < x 2 ≤ 1 \begin{aligned} &\mathop{min}\limits_{x_1,x_2} \ f(x_1,x_2) = \frac{1}{x_1}+x_1^2 + \frac{1}{x_2}+x_2^2\\ &s.t. \begin{cases} 0<x_1\leq1 \\ 0<x_2\leq1 \end{cases} \end{aligned} ​x1​,x2​min​ f(x1​,x2​)=x1​1​+x12​+x2​1​+x22​s.t.{0<x1​≤10<x2​≤1​​

clear all;
clc;

x = [-10 10]';   %决策变量初始化   %注：本函数的决策变量为向量！ 

lambda = ones( 4 , 1 );   %拉格朗日乘子初始化 
k=1;
G_v = [];  lambda_v = [];  x_v = [];  L_v = [];  f_v = [];
tol = 10^(-5);
t_lambda=0.01;  % 拉格朗日乘子的更新步长
epsilon = 0.1;
t_x_1=0.1;      % 决策变量 x 的更新步长
t_x_2=0.1;      % 决策变量 x 的更新步长
max_iteration = 10000;

for i = 1 : max_iteration
    %% 迭代更新决策变量 x
    x(1) = x(1) - t_x_1 * (  - x(1)^(-2) + 2 * x(1) - lambda(1) + lambda(2)   );  %% 注意这里可以写解析解，也可以写梯度迭代更新解，也可混合写！
    x(2) = x(2) - t_x_2 * (  - x(2)^(-2) + 2 * x(2) - lambda(3) + lambda(4) ); 
        
    %% 迭代更新拉格朗日乘子 lambda 
    Grad_f = [ ( -x(1) ) ; ( x(1)-1 ) ;  ( -x(2) ) ; ( x(2)-1 )  ] ; 
    lambda = max(   lambda + t_lambda * Grad_f  ,  zeros( 4 , 1 )    );  
    f =  ( x(1)^(-1) ) + ( x(1)^(2) ) + ( x(2)^(-1) ) + ( x(2)^(2) ) ; 
    L =  ( x(1)^(-1) ) + ( x(1)^(2) ) + ( x(2)^(-1) ) + ( x(2)^(2) ) + lambda(1)*( -x(1) ) + lambda(2)*( x(1)-1 ) + lambda(3)*( -x(2) ) + lambda(4)*( x(2)-1 ) ;
    
    %% 记录历史更新结果
    G_v = [G_v Grad_f];
    x_v = [x_v x];
    lambda_v = [lambda_v lambda];
    f_v = [f_v f]; 
    L_v = [L_v L];
    
    %% 终止条件判断
    % if (k>2) &&( abs(L - L_v(k-1)) < tol )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
    if ( k > 2 ) && ( abs( (f-L)  ) < tol ) &&  ( Grad_f(1)<=epsilon )  &&  ( Grad_f(2)<=epsilon )  &&  ( Grad_f(3)<=epsilon )  % k == 5000 %k>3&&(f_v(k)-f_v(k-1)<tol)
        break;
    end
    k=k+1;
end
if i == max_iteration
    disp(['超过最大迭代次数！']);
end
disp(['迭代次数为：',num2str(k),'次']);
disp(['最优值点为：']);
x
figure(1)
plot( 2:length(L_v),L_v(2:length(L_v)), '-ro', 'linewidth',2)
disp(['最优值为：',num2str(L),'']);
figure(2)
plot( 2:length(L_v), f_v(2:length(f_v)) - L_v(2:length(L_v)) , '-ro', 'linewidth',2)
figure(3)
plot( 2:length(x_v(1,:)),x_v(1,2:length(x_v)), '-bo', 'linewidth',2)
grid on
xlabel('Iterate')
ylabel('x_1 ')
title(['x_1' ]) 
figure(4)
plot( 2:length(x_v(2,:)),x_v(2,2:length(x_v)), '-go', 'linewidth',2)
grid on
xlabel('Iterate')
ylabel('x_2 ')
title(['x_2' ])