Homec4science
Diffusion Nonlinear (master) Browse

R7172/codes/functions/lbfgs_augmented_lagrangian.m
4bd89d382eedmaster

Last Change
Raw File

lbfgs_augmented_lagrangian.m
View Options

function [Ucur, output, yc] = lbfgs_augmented_lagrangian(f_u, f_grad_u, proj, update_params, ...
gamma_rule, get_feas, params, myFunc)
% LBFGS solver, inner problems are solved exactly and then dual update.
% Augmented Lagrangian:
% L(u, beta, yc) = f(u) + beta/2*||A(u)-b||^2 + <yc, A(u) - b>
%
% 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
% update_params: Function handle determining the update of
% beta(penalty parameter) and sigma(dual step-size)
% gamma_rule: Function handle determining the step size
% 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 : 03.01.2019
tol = 1e-10;
init_hist_track();
init_iters();
tic;
if (params.verbose)
% fprintf('Lbfgs Algorithm running ... \n');
fprintf('Iter. |Feas |Objective |Grad Mapping \n');
fprintf('------------------------------------------------------------------\n');
end;
current_obj = inf;
prev_obj = inf;
beta_old = beta_;
sigma_old = sigma_;
k_t = 1;
feas_old = inf;
gamma_in = 1e-5;
x_cvx = params.x_cvx;
maxIter = 1e5;
memsize = 15;
tot_iter = 0;
m_s = 10;
omega_ = 1.1;
m_0 = 2;
for ii = 2:params.iters
Uold = Ucur;
m_s = m_s*omega_;
% m_s = m_0*ii;
gradFunc = @(u) myFunc(u, yc, beta_);
[Ucur,ff,nfeval, gamma_] = optLBFGS(gradFunc,Uold,maxIter,memsize, 1e3/beta_);
% [Ucur,ff,nfeval, gamma_] = optLBFGS(gradFunc,Uold,m_s,memsize);
% get the current feasibility(A(u) - b), update sigma and yc
feas_cur = get_feas(Ucur);
[beta_, sigma_, k_t] = update_params(norm(feas_cur),feas_old, k_t, ii, beta_old, sigma_old, output.sigma, output.beta);
% beta_
beta_old = beta_;
sigma_old = sigma_;
feas_old = norm(feas_cur);
% update dual
yc = yc + sigma_ * feas_cur;
tot_iter = tot_iter + nfeval;
% 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
% Stopping criteria
if (abs(current_obj-prev_obj) < tol && feas_old < 1e-6)% || ...
%(abs(f_lagrangian(Ucur) - f_lagrangian(Uold)) < 1e-15)
fprintf('Converged: Breaking...\n')
break;
end
prev_obj = current_obj;
track_hist();
verbosity();
if tot_iter > params.iters
fprintf('Max iters reached: Breaking...\n\n')
break;
end
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.beta = nan(params.iters, 1);
output.gamma = nan(params.iters, 1);
output.norm_dual = nan(params.iters, 1);
output.lagr = nan(params.iters, 1);
output.dist_ = nan(params.iters, 1);
output.iter = 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);
output.beta(ii) = beta_;
[output.lagr(ii), ~] = gradFunc(Uold);
output.dist_(ii) = norm(Uav(1:size(Uav, 1)/2,1) - x_cvx(1:size(Uav, 1)/2,1), 'fro')^2;
output.iter(ii) = tot_iter;
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 |Nfeval \n');
fprintf('------------------------------------------------------------------------\n');
end
str = sprintf('%7.0d %3.3e %3.5e %3.3e %d'...
, ii, output.feas(ii), output.obj(ii), output.grad_mapping(ii), nfeval);
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;
gamma_ = params.gamma;
theta = params.theta;
[beta_, sigma_] = update_params(0,0,0,0,0,0,0,0); % set initials
beta_0 = beta_;
end
end
c4science · Help