Home > solver > TEBD_NC_v2_2.m

TEBD_NC_v2_2

PURPOSE ^

TEBD v2.2: main solver. Only generates dual variables for the last iteration of the algorithm.

SYNOPSIS ^

function results = TEBD_NC_v2_2(par,prob)

DESCRIPTION ^

TEBD v2.2: main solver. Only generates dual variables for the last iteration of the algorithm.
Does NOT allows cell defined variables
%
%
% Solves two-easy-block structured convex optimization problems of the form:
%
%                              min f(x)+h_1(x)+h_2(x)
% results = TEBD_NC_v2_2(par,prob)
%
% input:
%        - par                : (struct) parameters of the algorithm
%        - prob               : (struct) problem definition
%
% output:
%        - results            : (struct) contains the results after running 2EBD-HPE using the given parameters (par) on the given problem (prob)
%
% [1] R. D.C Monteiro, C. Ortiz and B. F. Svaiter. A first-order
%     block-decomposition method for solving two-easy-block structured
%     semidefinite programs, 2012.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2EBD-HPE:
% Modified by C.Ortiz
% Last Modified: 5/31/2013

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 %TEBD v2.2: main solver. Only generates dual variables for the last iteration of the algorithm.
0002 %Does NOT allows cell defined variables
0003 %%
0004 %%
0005 %% Solves two-easy-block structured convex optimization problems of the form:
0006 %%
0007 %%                              min f(x)+h_1(x)+h_2(x)
0008 %% results = TEBD_NC_v2_2(par,prob)
0009 %%
0010 %% input:
0011 %%        - par                : (struct) parameters of the algorithm
0012 %%        - prob               : (struct) problem definition
0013 %%
0014 %% output:
0015 %%        - results            : (struct) contains the results after running 2EBD-HPE using the given parameters (par) on the given problem (prob)
0016 %%
0017 %% [1] R. D.C Monteiro, C. Ortiz and B. F. Svaiter. A first-order
0018 %%     block-decomposition method for solving two-easy-block structured
0019 %%     semidefinite programs, 2012.
0020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0021 %% 2EBD-HPE:
0022 %% Modified by C.Ortiz
0023 %% Last Modified: 5/31/2013
0024 function results = TEBD_NC_v2_2(par,prob)
0025 
0026 
0027 sigma   = par.sigma;
0028 
0029 
0030 
0031 x = prob.x0;
0032 u = prob.u0;
0033 
0034 
0035 
0036 par.AlgName = '2EBD-v2.2-NoCell';
0037 fprintf('\n\n ****Begin 2EBD V2.2. iterations = %d ****\n', par.maxit);
0038 fprintf(par.AlgName);
0039 fprintf('\n');
0040 fprintf(TEBD_LegendAppend(par));
0041 fprintf('\n');
0042 
0043 tstart = tic;
0044 %%-----------------------------------------
0045 
0046 
0047 
0048 
0049 if isfield(par,'balancedinit') && par.balancedinit >0
0050     theta = TEBD_GetInitialScaling_NC_v2_2(prob,x,u,sigma);
0051 else
0052     theta = 1;
0053 end
0054 
0055 lambda = sigma/sqrt(theta);
0056 
0057 
0058 %%
0059 %% dynamicscaling = 1: dynamic scaling of the primal inner product (dynamic change of the weight in primal and dual intermediate stepsizes)
0060 %                   0: no dynamic scaling of the primal inner product
0061 dynamicscaling = par.dynamicscaling;
0062 
0063 %dynamic scaling flag
0064 isdynamicscaling = 0;
0065 
0066 if dynamicscaling == 1
0067     %% dyn_scale_update_it: # of iterations to update the primal inner product
0068     %                       scaling
0069     dyn_scale_update_it = par.dyn_scale_updateiteration;
0070     
0071     %% scaleratio: accepted ratio between primal and dual relative residuals
0072     %              for of the dynamic inner product scaling
0073     scaleratio = par.scaleratio;
0074     
0075     %% scalecorrection: change made in the primal inner product if the primal
0076     %                   and dual relative residuals are not within the accepted
0077     %                   ratio for the dynamic inner product scaling
0078     scalecorrection = par.scalecorrection;
0079 end
0080 
0081 %% print headers
0082 TEBD_PrintHeaders;
0083 
0084 %% save time preprocessing
0085 runhist.preprocesstime = toc(tstart);
0086 
0087 
0088 stopstatus = 1;
0089 
0090 tstart = tic;
0091 
0092 Z_tilde.pos = 1;
0093 
0094 %%
0095 for iter = 1:par.maxit;
0096     %% Perform main iteration
0097     
0098     % perform 1st prox
0099     gradf_x = prob.gradf(x);
0100     gradf_x_plus_u = gradf_x+u;
0101     [x_tilde,Z_tilde] = prob.resolvent1(x - gradf_x_plus_u*(theta*lambda),Z_tilde,lambda*theta);
0102     
0103     % perform 2nd prox using Moreau's identity
0104     u_plus_lambda_x = u +x_tilde*lambda;
0105     [resolv2] = prob.resolvent2(u_plus_lambda_x/lambda,lambda);
0106     u_tilde = u_plus_lambda_x - resolv2*lambda;
0107     
0108     % compute v\in T(x) or v\in T^\epsilon(x)
0109     v_1 = ((x-x_tilde)/lambda) + ((u_tilde-u)*theta);
0110     v_2 = (u-u_tilde)/lambda;
0111     
0112     if par.adaptive == 0
0113         lambda_a = lambda;
0114     else
0115         %% find adaptive (aggresive) stepsize lambda
0116         [lambda_a] = GetAdaptiveLambdaScaled(prob,lambda,sigma,x,x_tilde,v_1,u,u_tilde,v_2,gradf_x,theta);
0117     end
0118     
0119     % perform the extragradient step
0120     x = x - v_1*lambda_a;
0121     u = u - v_2*lambda_a;
0122     
0123     %% save the time until this iteration
0124     secs = toc(tstart);
0125     
0126     %% Obtain the residuals relevant to each problem
0127     dh2 = v_2+x_tilde;
0128     dh1 = v_1/theta - gradf_x - u_tilde;
0129     obj = prob.evalf2(x_tilde,u_tilde,dh2,dh1);
0130     residuals = prob.residuals2(x_tilde,u_tilde,v_1/theta,v_2,obj);
0131     
0132     %% Save History: save all the information to evaluate the performance
0133     runhist = TEBD_SaveHistory(obj,runhist,iter,residuals,secs,lambda_a/lambda,theta);
0134     
0135     %% check all stopping conditions
0136     [stopstatus,printstatus,statusstr,runhist] = TEBD_Check_StopCrit(par,residuals,runhist,iter,stopstatus);
0137     if printstatus
0138         TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0139         fprintf(statusstr);
0140     end
0141     
0142     if residuals.resnorm < par.tol || iter >= par.maxit
0143         
0144         %% print last iterate
0145         TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0146         fprintf(statusstr);
0147         
0148         %% Save results
0149         [~,y_tilde] = prob.resolvent2(u_plus_lambda_x/lambda,lambda);
0150         obj_d = prob.evalf(x_tilde,y_tilde,Z_tilde);
0151         residuals_d = prob.residuals(x_tilde,u_tilde,v_1/theta,v_2,obj_d,y_tilde,Z_tilde,u_plus_lambda_x,lambda);
0152         results = TEBD_GetFinalValues_2(runhist,iter,par,prob,x,u,u_tilde,x_tilde,y_tilde,Z_tilde,residuals_d,secs);
0153         break
0154     end
0155     
0156     %% print current progress
0157     if  iter == 1 || mod(iter,par.printit)==0
0158         
0159         %% print  iterate
0160         TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0161     end
0162     
0163     %% Dynamic scaling of the inner product
0164     if dynamicscaling == 1 && (max(residuals.res1,residuals.res2) > par.tol2||(residuals.res1>residuals.res2*5)||(residuals.res2>residuals.res1*5))
0165         if isdynamicscaling == 0
0166             fprintf(' ****Starting dynamic scaling at iteration = %d ****\n', iter);
0167             isdynamicscaling = 1;
0168         end
0169         if mod(iter,dyn_scale_update_it)==0
0170             wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0171             wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0172             if  wrel_norm2/wrel_norm1 > scaleratio
0173                 theta = theta/(scalecorrection^2);
0174                 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0175             elseif  wrel_norm1/wrel_norm2 > scaleratio
0176                 theta = theta*(scalecorrection^2);
0177                 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0178             end
0179             lambda = sigma/sqrt(theta);
0180         elseif mod(iter,dyn_scale_update_it) == 1
0181             % Short memory agregation of the decition residuals
0182             wrel_norm2 =  residuals.res2^(1/dyn_scale_update_it);
0183             wrel_norm1 = residuals.res1^(1/dyn_scale_update_it);
0184         else
0185             % Short memory agregation of the decition residuals
0186             wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0187             wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0188         end
0189     elseif isdynamicscaling == 1 && dynamicscaling == 1
0190         fprintf(' ****Stopping dynamic scaling at iteration = %d ****\n', iter);
0191         isdynamicscaling = 0;
0192     end
0193     
0194 end
0195 
0196 fprintf('\n\n ****End 2EBD-HPE algorithm ****\n');
0197 %%---------------------------------------------------------------
0198 %% end of main loop
0199 %%---------------------------------------------------------------
0200 %%
0201 %% *****************************************************************************
0202 end
0203 
0204 function [lambda_a] = GetAdaptiveLambdaScaled(prob,lambda,sigma,x,x_tilde,v_1,u,u_tilde,v_2,gradf_x,theta)
0205 %Compute the error
0206 eps = prob.error1(x,x_tilde,gradf_x);
0207 
0208 lambda_sq = lambda^2;
0209 
0210 %Compute the residual v
0211 v_x = v_1;
0212 v_u = v_2;
0213 Norm_v_x_sq = norm(v_x)^2;
0214 Norm_v_u_sq = norm(v_u)^2;
0215 
0216 
0217 Delta_x = (x_tilde-x);
0218 Delta_u = (u_tilde-u);
0219 Norm_Delta_x = norm(Delta_x);
0220 Norm_Delta_u = norm(Delta_u);
0221 Norm_Delta_x_sq = Norm_Delta_x^2;
0222 Norm_Delta_u_sq = Norm_Delta_u^2;
0223 
0224 tr_v_Delta_x = v_x'*Delta_x;
0225 tr_v_Delta_u = v_u'*Delta_u;
0226 
0227 alpha = sigma^2*(Norm_Delta_x_sq/theta+Norm_Delta_u_sq);
0228 
0229 a = (Norm_v_x_sq)/theta + Norm_v_u_sq;
0230 b = (lambda*Norm_v_x_sq+tr_v_Delta_x)/theta + lambda*Norm_v_u_sq + tr_v_Delta_u + eps;
0231 c = (lambda_sq*Norm_v_x_sq+2*lambda*tr_v_Delta_x+Norm_Delta_x_sq)/theta+lambda_sq*Norm_v_u_sq+Norm_Delta_u_sq+2*lambda*(tr_v_Delta_u+eps)-alpha;
0232 discriminant = sqrt(b^2-a*c);
0233 if b <= 0
0234     beta = (-b+discriminant)/a;
0235 else
0236     beta = -c/(b+discriminant);
0237 end
0238 if beta <=  1e-13
0239     beta = 0;
0240 else
0241     %fprintf('%d beta = %8.3e\n',iter,beta);
0242 end
0243 lambda_a = lambda + beta;
0244 end

Generated on Tue 06-Aug-2013 17:07:59 by m2html © 2005