Home > solver > DSA_BD.m

DSA_BD

PURPOSE ^

DSA-BD: main SDP solver simple.

SYNOPSIS ^

function [obj,X,y,Z,runhist] = DSA_BD(blk,At,C,b,par,X0,y0,Z0)

DESCRIPTION ^

 DSA-BD: main SDP solver simple.
%
% Explained in detail in Software_files/DSA_BD_Quick_Guide.pdf
%
% Solves problems in the standard primal and dual SDP of the form:
%
% (P) min{ <C,X> : A(X) = b, X positive semidefinite}
% (D) min{ <b,y> : At(y) + Z = C, Z positive semidefinite}
%
%*************************************************************************
% C.Ortiz
% Last Modified : 3/31/2011

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % DSA-BD: main SDP solver simple.
0002 %%
0003 %% Explained in detail in <a title="Software_files/DSA_BD_Quick_Guide.pdf" href="http://www2.isye.gatech.edu/~cod3/CamiloOrtiz/Software_files/DSA_BD_Quick_Guide.pdf" class="style_1">Quick user guide</a>
0004 %%
0005 %% Solves problems in the standard primal and dual SDP of the form:
0006 %%
0007 %% (P) min{ <C,X> : A(X) = b, X positive semidefinite}
0008 %% (D) min{ <b,y> : At(y) + Z = C, Z positive semidefinite}
0009 %%
0010 %%*************************************************************************
0011 %% C.Ortiz
0012 %% Last Modified : 3/31/2011
0013 function [obj,X,y,Z,runhist] = DSA_BD(blk,At,C,b,par,X0,y0,Z0)
0014 
0015 
0016 if isfield(par,'tol'); tol = par.tol;else tol = 1e-3; end;
0017 sigma   = par.sigma;
0018 sigma_2 = sigma^2;
0019 
0020 
0021 
0022 y = y0;
0023 if iscell(Z0); Z = Z0; else Z = {Z0}; end;
0024 if iscell(X0); X = X0; else X = {X0}; end;
0025 
0026 
0027 
0028 m = length(b);
0029 
0030 fprintf('\n\n ****Begin DSA_BD algorithm. iterations = %d ****\n', par.maxit);
0031 fprintf(DSA_BD_LegendAppend(par));
0032 fprintf([' n = ',num2str(par.n),', m = ',num2str(par.m),'\n']);
0033 
0034 % tstart = cputime;
0035 tstart = clock;
0036 %%-----------------------------------------
0037 
0038 normb = norm(b);
0039 normC = DSA_BD_ops(C,'norm');%norm(C)
0040 
0041 
0042 
0043 %% normalize =  1:divide Ax=b by normb and A*y+z=c by normC
0044 %               2:Precondition A  by L inverse where LtL = AAt
0045 %               3-Precondition A  by SL inverse where SLtLSt = AAt
0046 normalize = par.normalize;
0047 
0048 %% adaptivelambda   0: no agressive choice of extragradient step size
0049 %                   1: agressive choice of extragradient step size
0050 adaptivelambda = par.adaptivelambda;
0051 
0052 %% ad_iteration:    the frequency in which the adaptive agressive choice of
0053 %                   extragradient step size
0054 ad_iteration = par.adaptiveiteration;
0055 
0056 %% dynamicscaling = 1: dynamic scaling of the primal inner product (dynamic change of the weight in primal and dual intermediate stepsizes)
0057 %                   0: no dynamic scaling of the primal inner product
0058 dynamicscaling = par.dynamicscaling;
0059 
0060 %% dyn_scale_update_it: # of iterations to update the primal inner product
0061 %                       scaling
0062 dyn_scale_update_it = par.dyn_scale_updateiteration;
0063 
0064 %% scaleratio: accepted ratio between primal and dual relative residuals
0065 %              for of the dynamic inner product scaling
0066 scaleratio = par.scaleratio;
0067 
0068 %% scalecorrection: change made in the primal inner product if the primal
0069 %                   and dual relative residuals are not within the accepted
0070 %                   ratio for the dynamic inner product scaling
0071 scalecorrection = par.scalecorrection;
0072 
0073 %% scalerelnorm: 0: Checks the ratio of the primal and dual relative residuals
0074 %                   before normalizing for the dynamic inner product scaling
0075 %                1: Checks the ratio of the primal and dual relative residuals
0076 %                   after normalizing for the dynamic inner product scaling
0077 scalerelnorm = par.scalerelnorm;
0078 
0079 %% balancedinit =   0: No balanced initialization
0080 %                   1: Set y0 = argmin norm(Ainv(U)y-c)^2 where U=AAt.
0081 %                   2: Choose theta such that the primal and dual relative
0082 %                      residuals after one iterations are O(1).
0083 %                   3: Perform 1 and 2
0084 balancedinit = par.balancedinit;
0085 
0086 %% Normalize the problem
0087 L.R = sparse(m,m);
0088 L.Rt = L.R';
0089 L.S = sparse(m,m);
0090 L.St = L.S';
0091 b_orig = b;
0092 if normalize  == 1
0093     b = b/(normb+1);
0094     C = DSA_BD_ops(C,'/',(normC+1));%C/(normC+1)
0095     %A = DSA_BD_ops(A, '/',(normb+1));%A /(normb+1);
0096     At = DSA_BD_ops(At,'/',(normb+1));%At/(normb+1);
0097     %     for i = 1:m
0098     %         At{i} = At{i}/(normb+1);
0099     %     end
0100     y=y0*(normC+1)/(normb+1);
0101 end
0102 if normalize  >= 2 || balancedinit > 1
0103     [L] = DSA_BD_cholAAt(blk,At,m);
0104     runhist.nnz_L_A = nnz(L.R);
0105     fprintf('# non-zeros of chol L_A: %d\n',runhist.nnz_L_A);
0106     if L.p > 0
0107         fprintf(' rows of A linearly dependent. p = %4.0d, m = %4.0d\n',L.p,m);
0108         error(' nothing done');
0109     end
0110     if normalize >= 2
0111         b  = L.Rt\(L.St* b);%b  = L_A\ b;
0112     end
0113 end
0114 %% *************************************************************************
0115 % Calculate the Lipschitz Constant and lambda
0116 % -----------------------------------------
0117 %
0118 Lip       = DSA_BD_computeLipschitzConstant(blk,At,par);%computeLipschitzConstant3(A,At,n,par);
0119 
0120 
0121 lambda = sigma/Lip;
0122 
0123 
0124 
0125 %% *************************************************************************
0126 %  Balance the initial iterates and initial primal inner product scaling
0127 % -----------------------------------------
0128 %
0129 if balancedinit == 0
0130     theta = 1;
0131 else
0132     if balancedinit ~= 1
0133         theta = 1;
0134         if normalize < 2
0135             y = (L.S*(L.R\(L.Rt\(L.St*DSA_BD_AXfun(blk,At,C)))));
0136         elseif normalize >= 2
0137             y = (L.Rt\(L.St*DSA_BD_AXfun(blk,At,C)));%-(L_A\(St*AXfun3(A,C)));
0138         end
0139         X = DSA_BD_ops(X,'zeros');
0140         if ~iscell(X); X = {X}; end;
0141     end
0142     if balancedinit ~= 2
0143         theta = 1;
0144         AX_prev = DSA_BD_AXfun(blk,At,X,normalize,L);%AXfun3(A,X,normalize,L_A,St);
0145         y_prev = y;
0146         X_prev = X;
0147         lambda_x = lambda/theta;
0148         lambda_y = lambda*theta;
0149         [y,X,y_tilde,X_tilde,AX_tilde,AX_tilde_ORIG,Aty_tilde] = ...
0150             PerformIteration(blk,y_prev,X_prev,AX_prev,lambda_y,lambda_x,lambda,theta,sigma_2,C,b,At,L, ...
0151             normalize,adaptivelambda,1,ad_iteration,par);
0152         [rel_norm1,rel_norm2] ...
0153             = ComputeNormOfResuduals(blk,X_prev,X_tilde,lambda_x,C,Aty_tilde,AX_tilde,normb,normC,normalize,b,b_orig,AX_tilde_ORIG);
0154 
0155         balcondition = 1;
0156         if rel_norm1 > rel_norm2 && balancedinit == 4
0157             balcondition = 2;
0158         elseif  balancedinit ~= 4 && max(rel_norm1,1) > rel_norm2
0159             balcondition = 0;
0160         end
0161 
0162 
0163         while balcondition >0
0164             if balcondition == 1
0165                 theta = 2*theta;
0166             else
0167                 theta = theta/2;
0168             end
0169             lambda_x = lambda/theta;
0170             lambda_y = lambda*theta;
0171             [y,X,y_tilde,X_tilde,AX_tilde,AX_tilde_ORIG,Aty_tilde] = ...
0172                 PerformIteration(blk,y_prev,X_prev,AX_prev,lambda_y,lambda_x,lambda,theta,sigma_2,C,b,At,L, ...
0173                 normalize,adaptivelambda,1,ad_iteration,par);
0174             [rel_norm1,rel_norm2] ...
0175                 = ComputeNormOfResuduals(blk,X_prev,X_tilde,lambda_x,C,Aty_tilde,AX_tilde,normb,normC,normalize,b,b_orig,AX_tilde_ORIG);
0176             if balcondition == 1  && max(rel_norm1,1) < rel_norm2
0177                 balcondition = 1;
0178             elseif balcondition == 1 && balancedinit == 4 && rel_norm1 < rel_norm2
0179                 balcondition = 1;
0180             elseif balcondition == 2 && rel_norm1 > rel_norm2
0181                 balcondition = 2;
0182             else
0183                 balcondition = 0;
0184             end
0185         end
0186         X = X_prev;
0187         y = y_prev;
0188     end
0189 end
0190 lambda_x = lambda/theta;
0191 lambda_y = lambda*theta;
0192 
0193 %%
0194 if adaptivelambda == 1
0195     lambda_a = lambda;
0196 end
0197 
0198 objbest = [DSA_BD_trace(blk,C,X), b'*y,0,0];%[sum(sum(C.*X)), b'*y,0,0];
0199 
0200 stopCritSatisfied = 0;
0201 %%
0202 if par.print == 1
0203     fprintf(' it     secs   dual          primal      lg(rrd)    lg(rrp)   sigma    theta   lg(rrd_NORM) lg(rrp_NORM)\n');
0204 end
0205 
0206 % runhist.preprocesstime = cputime - tstart;
0207 runhist.preprocesstime = etime(clock,tstart);
0208 % tstart = cputime;
0209 tstart = clock;
0210 
0211 %%
0212 for iter = 1:par.maxit;
0213     %% Save previous variables
0214     y_prev = y;
0215     if adaptivelambda == 0 && iter > 1
0216         X_prev = X_tilde;
0217     else
0218         X_prev = X;
0219     end
0220 
0221     if iter > 1
0222         if adaptivelambda == 0
0223             AX_prev = AX_tilde;
0224         else
0225             AX_prev = (1-lambda_a/lambda)*AX_prev+(lambda_a/lambda)*AX_tilde;
0226         end
0227     else
0228         AX_prev = DSA_BD_AXfun(blk,At,X_prev,par.normalize,L);
0229     end
0230     %% Perform main iteration
0231 
0232     [y,X,y_tilde,X_tilde,AX_tilde,AX_tilde_ORIG,Aty_tilde,lambda_a] = ...
0233         PerformIteration(blk,y_prev,X_prev,AX_prev,lambda_y,lambda_x,lambda,theta,sigma_2,C,b,At,L, ...
0234         normalize,adaptivelambda,iter,ad_iteration,par);
0235     %% Obtain values of residuals
0236 
0237 
0238         [rel_norm1,rel_norm2,resNorm,eps,rel_norm1_ORIG,rel_norm2_ORIG,Z_tilde] ...
0239             = ComputeNormOfResuduals(blk,X_prev,X_tilde,lambda_x,C,Aty_tilde,AX_tilde,normb,normC,normalize,b,b_orig,AX_tilde_ORIG);
0240 
0241 
0242 
0243         %% Compute feasibility
0244             pfeas = rel_norm2;
0245             dfeas = rel_norm1;
0246             
0247             %% compute primal and dual objective.
0248             %---------------------------------------------------------------
0249             % Compute Obj Values
0250 
0251             obj = [DSA_BD_trace(blk,C,X_tilde),b'*y_tilde];%[sum(sum(C.*X_tilde)),-b'*y_tilde];
0252 
0253             if normalize == 1
0254                 obj = obj*(normC+1);
0255             end
0256 
0257             secs = etime(clock,tstart);
0258             %% check stopping criteria
0259             if max(rel_norm1, rel_norm2) < tol
0260                 fprintf( '%3.0d %8.2f %13.5e %13.5e %8.3f %8.3f  %9.6f %9.6f %8.3f %8.3f\n', iter, ...
0261                     secs, full(obj(2)), full(obj(1)), full(-log10(rel_norm1)), full(-log10(rel_norm2)), sigma, theta, full(-log10(rel_norm1_ORIG)), full(-log10(rel_norm2_ORIG)));
0262                 stopCritSatisfied = 1;
0263             end
0264 
0265 
0266             %% print current progress
0267             if par.print == 1 &&  (iter == 1 || (mod(iter,10)==0)) && stopCritSatisfied == 0;
0268 
0269                 fprintf( '%3.0d %8.2f %13.5e %13.5e %8.3f %8.3f  %9.6f %9.6f %8.3f %8.3f\n', iter, ...
0270                     secs, full(obj(2)), full(obj(1)), full(-log10(rel_norm1)), full(-log10(rel_norm2)), sigma, theta, full(-log10(rel_norm1_ORIG)), full(-log10(rel_norm2_ORIG)));
0271 
0272             end
0273 
0274             %% Calculate best so far
0275             if iter == 1 || bestNorm >= resNorm
0276                 bestNorm = resNorm;
0277                 objbest(1)=obj(1);
0278                 objbest(2)=obj(2);
0279                 objbest(3)=iter;
0280             end
0281 
0282             %% Dynamic scaling of the inner product
0283             if scalerelnorm == 1
0284                 wrel_norm2 = rel_norm2_ORIG;
0285                 wrel_norm1 = rel_norm1_ORIG;
0286             else
0287                 wrel_norm2 = rel_norm2;
0288                 wrel_norm1 = rel_norm1;
0289             end
0290             if mod(iter,dyn_scale_update_it)==0 && dynamicscaling == 1
0291                 if  wrel_norm2/wrel_norm1 > scaleratio
0292                     theta = theta/scalecorrection;
0293                 elseif  wrel_norm1/wrel_norm2 > scaleratio
0294                     theta = theta*scalecorrection;
0295 
0296                 end
0297                 lambda_x = lambda/theta;
0298                 lambda_y = lambda*theta;
0299             end
0300             %% Save history
0301             runhist.pobj(iter)  = obj(1);
0302             runhist.dobj(iter)  = obj(2);
0303             runhist.pobjbest(iter)  = objbest(1);
0304             runhist.dobjbest(iter)  = objbest(2);
0305 
0306             runhist.resnorm(iter) = resNorm;
0307 
0308             runhist.epsilon(iter) = eps;
0309             runhist.plambda(iter)  = lambda_x;
0310             runhist.dlambda(iter)  = lambda_y;
0311             runhist.cputime(iter) = secs+runhist.preprocesstime;
0312 
0313 
0314             runhist.bestresnorm(iter) = bestNorm;
0315 
0316             if adaptivelambda == 1 && mod(iter,ad_iteration) == 0
0317                 runhist.lambdaratio(iter) = (lambda_a)/(lambda);
0318             end
0319 
0320 
0321             runhist.pfeas(iter) = pfeas;%norm(AXfun3(A,X_tilde,par.normalize,L_A,St) - b);
0322             runhist.dfeas(iter) = dfeas;
0323 
0324 
0325             %% stop if criteria satisfied
0326             if stopCritSatisfied || iter>=par.maxit
0327                 runhist. maxitreached = 0;
0328                 if (iter>=par.maxit);
0329                     fprintf('max iterations reached. \n');
0330                     runhist. maxitreached = 1;
0331                 end
0332                 X = X_tilde;
0333                 if normalize ==1
0334                     y = y_tilde*(normC+1)/(1+normb+1);
0335                 elseif normalize >= 2
0336                     y = L.S*(L.R\y_tilde);
0337                 else
0338                     y = y_tilde;
0339                 end
0340 
0341                 runhist.processtime = secs;
0342                 runhist.numiterations = iter;
0343                 runhist.eachiteration = secs/iter;
0344                 runhist.totaltime = secs+runhist.preprocesstime;
0345                 runhist.finalbestresnorm = bestNorm;
0346                 runhist.bestiteration  = objbest(3);
0347 
0348                 fprintf('preprocess time: %10.3f \n',runhist.preprocesstime);
0349                 fprintf('# iterations: %d \n',runhist.numiterations);
0350                 fprintf('iterations time: %10.3f \n',runhist.processtime);
0351                 fprintf('each iteration time: %f \n',runhist.eachiteration);
0352                 fprintf('Resulting -LOG(norm): %f \n',-log10(runhist.finalbestresnorm));
0353                 fprintf('total time: %10.3f \n',runhist.totaltime);
0354                 return
0355             end
0356 
0357 end
0358 
0359 fprintf('\n\n ****End DSA-BD algorithm ****\n');
0360 %%---------------------------------------------------------------
0361 %% end of main loop
0362 %%---------------------------------------------------------------
0363 %%
0364 %% *****************************************************************************
0365 
0366 function [y,X,y_tilde,X_tilde,AX_tilde,AX_tilde_ORIG,Aty_tilde,lambda_a] = ...
0367     PerformIteration(blk,y_prev,X_prev,AX_prev,lambda_y,lambda_x,lambda,theta,sigma_2,C,b,At,L, ...
0368     normalize,adaptivelambda,iter,ad_iteration,par)
0369 %% \tilde y_k = y_{k-1} +  \lambda(A x_{k-1} - b)
0370 y_tilde = y_prev-lambda_y*(AX_prev - b);
0371 
0372 %% \tilde x_{k} = P_X(x_{k-1} - \lambda(\nabla f(x_{k-1}) - A^*\tilde y_{k}))
0373 Aty_tilde = DSA_BD_Atyfun(blk,At,y_tilde,par.normalize,L);
0374 ProjArg = DSA_BD_ops(X_prev,'calc3',C,lambda_x, Aty_tilde);
0375 X_tilde = project(blk,ProjArg);
0376 
0377 
0378 
0379 %% Get valbues of AX without normalizing if using L inverse (normaliz>=2)
0380 if normalize < 2
0381     AX_tilde = DSA_BD_AXfun(blk,At,X_tilde,par.normalize,L);
0382 else
0383     [AX_tilde,AX_tilde_ORIG] = DSA_BD_AXfun(blk,At,X_tilde,par.normalize,L);
0384 end
0385 
0386 %% Perform aggresive choice of extragradient stepsize
0387 if adaptivelambda == 0 || mod(iter,ad_iteration) ~= 0 || iter == 1
0388 
0389     X = X_tilde;
0390 
0391     y = y_tilde - lambda_y*(AX_tilde - AX_prev);
0392 
0393     lambda_a = lambda;
0394 else
0395     dif_A_x = (AX_tilde - AX_prev);
0396 
0397     v_x_norm = DSA_BD_ops(X_prev,'diffnorm',X_tilde);
0398 
0399     v_y = AX_tilde-b;
0400 
0401     gamma_k = lambda_y*((v_y)'*dif_A_x);
0402 
0403     eta_k = (v_x_norm/lambda)^2 + (norm(v_y)^2);
0404 
0405     nu_k = sigma_2*( (theta*v_x_norm)^2 + norm(y_tilde-y_prev)^2) - (lambda_y*norm(dif_A_x))^2;
0406 
0407     if nu_k < 0
0408         fprintf('discriminant = %f \n', nu_k);
0409         error('\|a\|^2-\sigma^2\|\tilde x_k - x_{k-1}\|^2 > 0');
0410     end
0411     discriminant = sqrt(gamma_k^2+ eta_k*nu_k);
0412     if gamma_k <= 0
0413         beta = (-gamma_k + discriminant)/eta_k;
0414     else
0415         beta = nu_k/(gamma_k+discriminant);
0416     end
0417     lambda_a_y = lambda_y + beta;
0418     lambda_a   = lambda_a_y /theta;
0419 
0420     X = DSA_BD_ops(X_prev,'calc7',lambda_a,lambda,X_tilde);%X = (1-lambda_a/lambda)X_prev+(lambda_a/lambda)X_tilde;
0421     y = y_prev-lambda_a_y*v_y;
0422 
0423 
0424 
0425 end
0426 
0427 function [rel_norm1,rel_norm2,resNorm,eps,rel_norm1_ORIG,rel_norm2_ORIG,p] ...
0428     = ComputeNormOfResuduals(blk,X_prev,X_tilde,lambda_x,C,Aty_tilde,AX_tilde,normb,normC,normalize,b,b_orig,AX_tilde_ORIG)
0429 %% compute p_k^x
0430 p = DSA_BD_ops(X_tilde,'calc1',X_prev,lambda_x, C,Aty_tilde);
0431 
0432 
0433 %% Compute residual vector From Theorem 5.3
0434 
0435 res1 = DSA_BD_ops(C,'calc2',Aty_tilde,p);%Aty_tilde-C+p =-C+ Aty_tilde+p %DSA_BD_ops(DSA_BD_ops(C,'+',Aty_tilde),'+',p);%(C+Aty_tilde)+p;
0436 
0437 %res1 = C+Atyfun3(At,y_tilde,n,par.normalize,L_At,S)+p;
0438 if normalize < 2
0439     res2 = b-AX_tilde;
0440 
0441     norm2 = norm(res2);
0442     norm2_ORIG = norm2;
0443     %res2 = b-AXfun3(A,X_tilde,par.normalize,L_A,St);
0444 else
0445     res2 = b_orig - AX_tilde_ORIG;
0446 
0447     norm2 = norm(res2);
0448     res2_ORIG = b-AX_tilde;
0449     norm2_ORIG = norm(res2_ORIG);
0450     %res2 = b_orig-AXfun3(A,X_tilde);
0451 end
0452 %% Compute the norm of the residual
0453 norm1 = DSA_BD_ops(res1,'norm');
0454 %         norm2 = norm(res2);
0455 
0456 
0457 %% complementarity in pointwise
0458 eps = DSA_BD_trace(blk,X_tilde,p);%sum(sum(X_tilde.*p));
0459 if normalize ~= 1
0460     rel_norm1 = norm1/(1+normC);
0461     rel_norm2 = norm2/(1+normb);
0462     rel_norm1_ORIG = rel_norm1;
0463     rel_norm2_ORIG = norm2_ORIG/(1+normb);
0464 else
0465     rel_norm1 = norm1;
0466     rel_norm2 = norm2;
0467     rel_norm1_ORIG = norm1/(1+normC);
0468     rel_norm2_ORIG = norm2_ORIG/(1+normb);
0469 end;
0470 resNorm = sqrt(rel_norm1^2 + rel_norm2^2);

Generated on Mon 19-Sep-2011 19:28:04 by m2html © 2005