0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
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
0035 tstart = clock;
0036
0037
0038 normb = norm(b);
0039 normC = DSA_BD_ops(C,'norm');
0040
0041
0042
0043
0044
0045
0046 normalize = par.normalize;
0047
0048
0049
0050 adaptivelambda = par.adaptivelambda;
0051
0052
0053
0054 ad_iteration = par.adaptiveiteration;
0055
0056
0057
0058 dynamicscaling = par.dynamicscaling;
0059
0060
0061
0062 dyn_scale_update_it = par.dyn_scale_updateiteration;
0063
0064
0065
0066 scaleratio = par.scaleratio;
0067
0068
0069
0070
0071 scalecorrection = par.scalecorrection;
0072
0073
0074
0075
0076
0077 scalerelnorm = par.scalerelnorm;
0078
0079
0080
0081
0082
0083
0084 balancedinit = par.balancedinit;
0085
0086
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));
0095
0096 At = DSA_BD_ops(At,'/',(normb+1));
0097
0098
0099
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);
0112 end
0113 end
0114
0115
0116
0117
0118 Lip = DSA_BD_computeLipschitzConstant(blk,At,par);
0119
0120
0121 lambda = sigma/Lip;
0122
0123
0124
0125
0126
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)));
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);
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];
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
0207 runhist.preprocesstime = etime(clock,tstart);
0208
0209 tstart = clock;
0210
0211
0212 for iter = 1:par.maxit;
0213
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
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
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
0244 pfeas = rel_norm2;
0245 dfeas = rel_norm1;
0246
0247
0248
0249
0250
0251 obj = [DSA_BD_trace(blk,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
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
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
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
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
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;
0322 runhist.dfeas(iter) = dfeas;
0323
0324
0325
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
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
0370 y_tilde = y_prev-lambda_y*(AX_prev - b);
0371
0372
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
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
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);
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
0430 p = DSA_BD_ops(X_tilde,'calc1',X_prev,lambda_x, C,Aty_tilde);
0431
0432
0433
0434
0435 res1 = DSA_BD_ops(C,'calc2',Aty_tilde,p);
0436
0437
0438 if normalize < 2
0439 res2 = b-AX_tilde;
0440
0441 norm2 = norm(res2);
0442 norm2_ORIG = norm2;
0443
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
0451 end
0452
0453 norm1 = DSA_BD_ops(res1,'norm');
0454
0455
0456
0457
0458 eps = DSA_BD_trace(blk,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);