Home > solver > TEBD_v2_2.m

TEBD_v2_2

PURPOSE ^

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

SYNOPSIS ^

function results = TEBD_v2_2(par,prob)

DESCRIPTION ^

TEBD v2.2: main solver. Only generates dual variables for the last iteration of the algorithm. 
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_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 %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_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_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';
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_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     [x_tilde,Z_tilde] = prob.resolvent1(TEBD_ops(x,'X-alpha*(Y+X2)',gradf_x,theta*lambda,u),Z_tilde,lambda*theta);
0101     
0102     % perform 2nd prox using Moreau's identity
0103     [u_plus_lambda_x_div_lambda, u_plus_lambda_x] = TEBD_ops(u ,'(X+Y*alpha)/alpha',x_tilde,lambda);
0104     [resolv2] = prob.resolvent2(u_plus_lambda_x_div_lambda,lambda);
0105     u_tilde = TEBD_ops(u_plus_lambda_x ,'-', resolv2, lambda);
0106     
0107     % compute v\in T(x) or v\in T^\epsilon(x)
0108     [v_1, v_1_div_theta,dh1] = TEBD_ops(x,'[(X-Y)/alpha +(X2-Y2)*alpha2]/alpha2-X3-Y3',x_tilde,lambda,u_tilde,u,theta,gradf_x,u_tilde);
0109     [v_2,dh2] = TEBD_ops(u,'(X-Y)/alpha+X2',u_tilde,lambda,x_tilde);
0110     
0111     if par.adaptive == 0
0112         lambda_a = lambda;
0113     else
0114         %% find adaptive (aggresive) stepsize lambda
0115         [lambda_a] = GetAdaptiveLambdaScaled(prob,lambda,sigma,x,x_tilde,v_1,u,u_tilde,v_2,gradf_x,theta);
0116     end
0117     
0118     % perform the extragradient step
0119     x = TEBD_ops(x ,'-', v_1,lambda_a);
0120     u = TEBD_ops(u ,'-', v_2,lambda_a);
0121     
0122     %% save the time until this iteration
0123     secs = toc(tstart);
0124     
0125     %% Obtain the residuals relevant to each problem
0126     obj = prob.evalf2(x_tilde,u_tilde,dh2,dh1);
0127     residuals = prob.residuals2(x_tilde,u_tilde,v_1_div_theta,v_2,obj);
0128     
0129     %% Save History: save all the information to evaluate the performance
0130     runhist = TEBD_SaveHistory(obj,runhist,iter,residuals,secs,lambda_a/lambda,theta);
0131     
0132     %% check all stopping conditions
0133     [stopstatus,printstatus,statusstr,runhist] = TEBD_Check_StopCrit(par,residuals,runhist,iter,stopstatus);
0134     if printstatus
0135         TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0136         fprintf(statusstr);
0137     end
0138     
0139     if residuals.resnorm < par.tol || iter >= par.maxit
0140         
0141         %% print last iterate
0142         TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0143         fprintf(statusstr);
0144         
0145         %% Save results
0146         [~,y_tilde] = prob.resolvent2(u_plus_lambda_x_div_lambda,lambda);
0147         obj_d = prob.evalf(x_tilde,y_tilde,Z_tilde);
0148         residuals_d = prob.residuals(x_tilde,u_tilde,v_1_div_theta,v_2,obj_d,y_tilde,Z_tilde,u_plus_lambda_x,lambda);
0149         results = TEBD_GetFinalValues_2(runhist,iter,par,prob,x,u,u_tilde,x_tilde,y_tilde,Z_tilde,residuals_d,secs);
0150         break
0151     end
0152     
0153     %% print current progress
0154     if  iter == 1 || mod(iter,par.printit)==0
0155         
0156         %% print  iterate
0157         TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0158     end
0159     
0160     %% Dynamic scaling of the inner product
0161     if dynamicscaling == 1 && (max(residuals.res1,residuals.res2) > par.tol2||(residuals.res1>residuals.res2*5)||(residuals.res2>residuals.res1*5))
0162         if isdynamicscaling == 0
0163             fprintf(' ****Starting dynamic scaling at iteration = %d ****\n', iter);
0164             isdynamicscaling = 1;
0165         end
0166         if mod(iter,dyn_scale_update_it)==0
0167             wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0168             wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0169             if  wrel_norm2/wrel_norm1 > scaleratio
0170                 theta = theta/(scalecorrection^2);
0171                 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0172             elseif  wrel_norm1/wrel_norm2 > scaleratio
0173                 theta = theta*(scalecorrection^2);
0174                 TEBD_PrintCurrentIterate(par,runhist,iter,theta,residuals);
0175             end
0176             lambda = sigma/sqrt(theta);
0177         elseif mod(iter,dyn_scale_update_it) == 1
0178             % Short memory agregation of the decition residuals
0179             wrel_norm2 =  residuals.res2^(1/dyn_scale_update_it);
0180             wrel_norm1 = residuals.res1^(1/dyn_scale_update_it);
0181         else
0182             % Short memory agregation of the decition residuals
0183             wrel_norm2 = wrel_norm2 * residuals.res2^(1/dyn_scale_update_it);
0184             wrel_norm1 = wrel_norm1 * residuals.res1^(1/dyn_scale_update_it);
0185         end
0186     elseif isdynamicscaling == 1 && dynamicscaling == 1
0187         fprintf(' ****Stopping dynamic scaling at iteration = %d ****\n', iter);
0188         isdynamicscaling = 0;
0189     end
0190     
0191 end
0192 
0193 fprintf('\n\n ****End 2EBD-HPE algorithm ****\n');
0194 %%---------------------------------------------------------------
0195 %% end of main loop
0196 %%---------------------------------------------------------------
0197 %%
0198 %% *****************************************************************************
0199 end
0200 
0201 function [lambda_a] = GetAdaptiveLambdaScaled(prob,lambda,sigma,x,x_tilde,v_1,u,u_tilde,v_2,gradf_x,theta)
0202 %Compute the error
0203 eps = prob.error1(x,x_tilde,gradf_x);
0204 
0205 lambda_sq = lambda^2;
0206 
0207 %Compute the residual v
0208 v_x = v_1;
0209 v_u = v_2;
0210 Norm_v_x_sq = prob.norm(v_x)^2;
0211 Norm_v_u_sq = prob.norm(v_u)^2;
0212 
0213 
0214 Delta_x = TEBD_ops(x_tilde,'-',x);
0215 Delta_u = TEBD_ops(u_tilde,'-',u);
0216 Norm_Delta_x = prob.norm(Delta_x);
0217 Norm_Delta_u = prob.norm(Delta_u);
0218 Norm_Delta_x_sq = Norm_Delta_x^2;
0219 Norm_Delta_u_sq = Norm_Delta_u^2;
0220 
0221 tr_v_Delta_x = prob.trace(v_x,Delta_x);
0222 tr_v_Delta_u = prob.trace(v_u,Delta_u);
0223 
0224 alpha = sigma^2*(Norm_Delta_x_sq/theta+Norm_Delta_u_sq);
0225 
0226 a = (Norm_v_x_sq)/theta + Norm_v_u_sq;
0227 b = (lambda*Norm_v_x_sq+tr_v_Delta_x)/theta + lambda*Norm_v_u_sq + tr_v_Delta_u + eps;
0228 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;
0229 discriminant = sqrt(b^2-a*c);
0230 if b <= 0
0231     beta = (-b+discriminant)/a;
0232 else
0233     beta = -c/(b+discriminant);
0234 end
0235 if beta <=  1e-13
0236     beta = 0;
0237 else
0238     %fprintf('%d beta = %8.3e\n',iter,beta);
0239 end
0240 lambda_a = lambda + beta;
0241 end

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