Homec4science
Diffusion Nonlinear (master) Browse

R7172/codes/functions/linearized_augmented_lagrangian.m
3225d2ffcba7master

Last Change
Raw File

linearized_augmented_lagrangian.m
View Options

function [Ucur, output] = linearized_augmented_lagrangian(f_u, f_grad_u, proj, beta_rule, ...
gamma_rule, sigma_rule, get_feas, params)
% Non-convex solver for sdp's. This function is used to solve a nonlinear
% constrained nonconvex problem with an inexact augmented lagrangian method.
% min f(u), subject to A(u) = b, u \in \mathcal{C}
% Augmented Lagrangian:
% L(u, beta, yc) = f(u) + 1/(2*beta)||A(u)-b||^2 + <yc, b-A(u)>
%
% INPUT: f_u: function handle to evaluate the function value
% f_grad_u: function handle to evaluate the gradient of augmented
% lagrangian. Accepts 3 inputs: u, beta, yc
% proj: Projection operator onto set C
% beta_rule: function handle determining the beta update
% gamma_rule: function handle determining the step size
% sigma_rule: function handle determining dual step: ds = 1/sigma
% get_feas: function handle which returns (A(u) - b)
% params: structure containing parameters for the
% algorithm( maxiters, Uinit, tol, savehist etc.)
% OUTPUT: Ucur: Last iterate of the algorithm
% output: structure containing the history of the
% iterates(times, feas, obj, grad_mapping etc.)
%
% Author: Mehmet Fatih Sahin
% Date : 02.11.2018
tol = 1e-15;
init_hist_track();
init_iters();
tic;
if (params.verbose)
fprintf('Nonconvex Algorithm running ... \n');
fprintf('Iter. |Feas |Objective |Grad Mapping \n');
fprintf('------------------------------------------------------------------\n');
end;
current_obj = inf;
prev_obj = inf;
for ii = 1:params.iters
beta_ = beta_rule(ii);
% TODO(mfsahin): add line search here
gamma = gamma_rule(beta_, yc);
Uold = Ucur;
gradU = f_grad_u(Ucur, beta_, yc);
Ucur = Ucur - gamma * gradU;
% projection onto set C
Ucur = proj(Ucur);
% get the current feasibility(A(u) - b), update sigma and yc
feas_cur = get_feas(Ucur);
sigma = sigma_rule(ii, norm(feas_cur));
% update dual
yc = yc + 1/sigma * feas_cur;
% average of the iterates
if params.averaging
Usum = Ucur + Usum;
ysum = yc + ysum;
Uav = 1/(ii+1)*Usum;
ycav = 1/(ii+1)*ysum;
else
Uav = Ucur;
ycav = yc;
end
if abs(current_obj-prev_obj) < tol
fprintf('Converged: Breaking...\n')
break;
end
prev_obj = current_obj;
track_hist();
verbosity();
end
if ii == params.iters
fprintf('Max iters reached: Breaking...\n\n')
end
function [] = init_hist_track()
% Initialization
output.times = nan(params.iters, 1);
output.feas = nan(params.iters, 1);
output.obj = nan(params.iters, 1);
output.grad_mapping = nan(params.iters, 1);
output.sigma = nan(params.iters, 1);
output.gamma = nan(params.iters, 1);
output.norm_dual = nan(params.iters, 1);
end
function [] = track_hist()
% save the history
if (params.savehist == 1)
output.times(ii) = toc; tic;
output.feas(ii) = norm(get_feas(Uav));
output.obj(ii) = f_u(Uav);
output.grad_mapping(ii) = 1/gamma^2 * norm(Ucur-Uold, 2)^2;
current_obj = output.obj(ii);
output.sigma(ii) = sigma;
output.gamma(ii) = gamma;
output.norm_dual(ii) = norm(yc);
end;
end
function [] = verbosity()
if (params.verbose) && (mod(ii, params.verbose)==0 || ii == 1)
if ((mod(ii, 10*params.verbose)==0))
fprintf('Iter. |Feas |Objective |Grad Mapping \n');
fprintf('------------------------------------------------------------------\n');
end
str = sprintf('%d %f %f %f '... %f'...
, ii, output.feas(ii), output.obj(ii), output.grad_mapping(ii));%, norm_gf(ii));
disp(str);
end
end
function [] = init_iters()
% initialize the iterations
Ucur = params.Uinit;
Ucur = Ucur./norm(Ucur, 'fro');
if params.startfeasible
% Start feasible for maxcut
rowN = sqrt(sum(Ucur.^2,2));
Ucur = Ucur./(rowN*ones(1,size(Ucur,2)));
end
yc = params.ycinit; % center point to be updated
ysum = yc;
Usum = Ucur;
end
end
c4science · Help