ADMM算法求解一个简单的例子

Posted 2023-01-19 机器学习的小学生

求解下面的带有等式约束和简单的边框约束(box constraints)的优化问题：

$$\\beginequation \\beginaligned \\min_x,y (x-1)^2+(y-2)^2 \\\\ s.t. 0\\leq x \\leq 3,\\\\ 1 \\leq y \\leq 4, \\\\ 2x+3y=5 \\endaligned \\endequation$$

% 求解下面的最小化问题：
% min_x,y (x-1)^2 + (y-2)^2
% s.t. 0 \\leq x \\leq 3
%      1 \\leq y \\leq 4
%      2x + 3y = 5

% augumented lagrangian function:
% L_rho(x,y,lambda) = (x-1)^2 + (y-2)^2 + lambda*(2x + 3y -5) + rho/2(2x + 3y -5)^2
% solve x  min f(x) = (x-1)^2 +  lambda*(2x + 3y -5) + rho/2(2x + 3y -5)^2,s.t. 0<= x <=3
% solve y  min f(y) = (y-2)^2 + lambda*(2x + 3y -5) + rho/2(2x + 3y -5)^2, s.t. 1<= y <=4
% sovle lambda :update lambda = lambda + rho(2x + 3y -5)
% rho ,we set rho = min(rho_max,beta*rho) beta is a constant ,we set to 1.1,rho0 = 0.5;

% x0,y0 都为1是一个可行解。
param.x0 = 1;
param.y0 = 1;
param.lambda = 1;
param.maxIter = 30;
param.beta = 1.1; % a constant
param.rho = 0.5;
param.rhomax = 2000;

[Hx,Fx] = getHession_F('f1');
[Hy,Fy] = getHession_F('f2');

param.Hx = Hx;
param.Fx = Fx;
param.Hy = Hy;
param.Fy = Fy;

% solve problem using admm algrithm
[x,y] = solve_admm(param);

% disp minimum
disp(['[x,y]:' num2str(x) ',' num2str(y)]);

function [H,F] = getHession_F(fn)
% fn : function name
% H :hessian matrix
% F :一次项系数
syms x y lambda rho;
if strcmp(fn,'f1')
    f = (x-1)^2 + lambda*(2*x + 3*y -5) + rho/2*(2*x + 3*y -5)^2;
    H = hessian(f,x);
    F = (2*lambda + (rho*(12*y - 20))/2 - 2);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    fcol = collect(f,'x');
    disp(fcol);
elseif strcmp(fn,'f2')
    f = (y-2)^2 + lambda*(2*x + 3*y -5) + rho/2*(2*x + 3*y -5)^2;
    H = hessian(f,y);
    F = (3*lambda + (rho*(12*x - 30))/2 - 4);
    fcol = collect(f,'y');
    disp(fcol);
end
% fcol = collect(f,'x');
% fcol = collect(f,'y');
% disp(fcol);
end

function [x,y] = solve_admm(param)

x = param.x0;
y = param.y0;
lambda = param.lambda;
beta = param.beta;
rho = param.rho;
rhomax = param.rhomax;
Hx = param.Hx;
Fx = param.Fx;
Hy = param.Hy;
Fy = param.Fy;
%%
options = optimoptions('quadprog','Algorithm','interior-point-convex');
xlb = 0;
xub = 3;
ylb = 1;
yub = 4;
maxIter = param.maxIter;
i = 1;
% for plot
funval = zeros(maxIter-1,1);
iterNum = zeros(maxIter-1,1);
while 1
    if i == maxIter
        break;
    end
    % solve x
    Hxx = eval(Hx);
    Fxx = eval(Fx);
    x = quadprog(Hxx,Fxx,[],[],[],[],xlb,xub,[],options);% descend
    % solve y
    Hyy = eval(Hy);
    Fyy = eval(Fy);
    y = quadprog(Hyy,Fyy,[],[],[],[],ylb,yub,[],options);%descend
    % update lambda
    lambda = lambda + rho*(2*x + 3*y -5); % ascend
    % rho = min(rhomax,beta*rho);
    funval(i) = compute_fval(x,y);
    iterNum(i) = i;
    i = i + 1;
end

plot(iterNum,funval,'-r');

end

function  fval  = compute_fval(x,y)
fval = (x-1)^2 + (y-2)^2;
end

结果：

[x,y]:0.53846,1.3077

从上图我们可以看到，算法在第5次迭代后就几乎收敛。