Home > source > @penlab > eval_alx.m

# eval_alx

## PURPOSE

Evaluate Augmented Lagrangian at the current 'xall' point with the current

## SYNOPSIS

function [status] = eval_alx(obj)

## DESCRIPTION

     Evaluate Augmented Lagrangian at the current 'xall' point with the current
     penalty parameters and Lagrangian multipliers. The result will be stored
     in 'obj.ALx' and its ticker updated. All equalities 'obj.eqx' get updated
     as well. Return nonzero if evaluation failed (matrix not pos def, etc.) and
     in this case, the value ALx will be set to Inf.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- penlab PenLab (Pennon Laboratory for Matlab, previously PennonM)

## SOURCE CODE

    0001 % Evaluate Augmented Lagrangian at the current 'xall' point with the current
    0002 % penalty parameters and Lagrangian multipliers. The result will be stored
    0003 % in 'obj.ALx' and its ticker updated. All equalities 'obj.eqx' get updated
    0004 % as well. Return nonzero if evaluation failed (matrix not pos def, etc.) and
    0005 % in this case, the value ALx will be set to Inf.
    0006 function [status] = eval_alx(obj)
    0007
    0008   % TODO ... set status when cannot evaluate or something like that
    0009   status=0;
    0010
    0011   % work only if the value is out of date
    0012   if (obj.ALxtck < obj.ticker)
    0013     starttime = cputime;
    0014
    0015     % create local copies of obj.x,obj.Y to avoid checking repetitively if they
    0016     % are 'up-to-date' with obj.xall, they get used many times in calls to
    0017     % user's functions
    0018     x=obj.x;
    0019     Y=obj.Y;
    0020     % same for obj.userdata, it might be quite expensive to get&store it in
    0021     % the object
    0022     userdata=obj.userdata;
    0023
    0024     [fx, userdata] = obj.objfun(x, Y, userdata);
    0025     obj.objx=fx;
    0026
    0027     % function inequal & equal
    0028     %ineqx = zeros(obj.Nineq,1);   % probably not necessary, rather = []?
    0029     ineqx=[];
    0030     %obj.eqx = zeros(obj.Neq,1);
    0031     if (obj.NgNLN + obj.NgLIN>0)
    0032       [gx, userdata] = obj.confun(x, Y, userdata);
    0033       ineqx = obj.ineqshift + obj.ineqmlt .* gx(obj.ineqmap);
    0034       obj.ineqx = ineqx;
    0035       obj.eqx = obj.eqshift + gx(obj.eqmap);
    0036     end
    0037
    0038     % box constraints
    0039     xboxx = [];
    0040     if (obj.Nxbox>0)
    0041       xboxx = obj.xboxshift + obj.xboxmlt .* obj.xall(obj.xboxmap);
    0042       obj.xboxx=xboxx;
    0043     end
    0044
    0045     ALx=fx;
    0046
    0047     %  inequalitites: phi or bar?
    0048     %  L = L + u'*(p.*phi2(g./p));
    0049     %  L = L + u(i)*p(i)*phi2(g(i)/p(i));
    0050
    0051     % box constraints
    0052     ind=obj.xboxindbar;
    0053     if (~isempty(ind))
    0054       ALx=ALx + obj.uxbox(ind)'*(obj.pxbox(ind).*obj.phibar(xboxx(ind)));
    0055     end
    0056
    0057     ind=obj.xboxindphi;
    0058     if (~isempty(ind))
    0059       ALx=ALx + obj.uxbox(ind)'*(obj.pxbox(ind).*obj.phi2(xboxx(ind)./obj.pxbox(ind)));
    0060     end
    0061
    0062     % function inequalitites
    0063     %ind=obj.ineqindbar;
    0064     %if (~isempty(ind))
    0065     %  ALx=ALx + obj.uineq(ind)'*(obj.pineq(ind).*phibar(ineqx(ind)));
    0066     %end
    0067
    0068     ind=obj.ineqindphi;
    0069     if (~isempty(ind))
    0070       ALx=ALx + obj.uineq(ind)'*(obj.pineq(ind).*obj.phi2(ineqx(ind)./obj.pineq(ind)));
    0071     end
    0072
    0073     % (function) equalities
    0074     if (obj.Neq>0)
    0075       ALx = ALx + obj.ueq'*obj.eqx;
    0076     end
    0077
    0078     % matrix variable - log barrier (strict feasibility)
    0079     if ~isempty(obj.Yboxindbar)
    0080     for k=obj.Yboxindbar
    0081       % convert the matrix box constraint to the form:   +/-Y +/-bound >=0
    0082       pkx=obj.PYbox(k);
    0083       Ykx = Y{obj.Yboxmap(k)};
    0084       Akx=-obj.Yboxshift(k)*speye(size(Ykx)) - obj.Yboxmlt(k)*Ykx;
    0085       % check that it is pos. def.
    0086       [R,iii]=chol(Akx);
    0087       if (iii~=0)
    0088         ALx=Inf;
    0089         % jump out and even avoid anything other cycle
    0090         break;
    0091       end
    0092
    0093       % compute penalty term: -p*log(det(A)) = -p*log( det(R)^2 ) =
    0094       %    = -2*p*log( prod(diag(R)) ) = -2*p*sum(log( diag(R) ))
    0095       ALx = ALx -2*pkx*sum(log(diag(R)));
    0096
    0097     end
    0098     end
    0099
    0100
    0101     % matrix variable - pen/bar
    0102     if ~isempty(obj.Yboxindphi)
    0103     for k=obj.Yboxindphi
    0104       pkx=obj.PYbox(k);            %2*p(sdpdata.Ng+k);
    0105       Ykx = Y{obj.Yboxmap(k)};
    0106       Akx=obj.Yboxshift(k)*speye(size(Ykx)) + obj.Yboxmlt(k)*Ykx;
    0107       umatk=obj.UYbox{k};
    0108
    0109       %Z=(pkx*speye(size(Akx))-Akx);
    0110       Z=(pkx*speye(size(Akx))-Akx);
    0111       [R,iii]=chol(Z);   % to double check that it is pos. def!
    0112       if (iii~=0)
    0113         ALx=Inf;
    0114         % jump out and even avoid anything other cycle
    0115         break;
    0116       end
    0117
    0118       invZ=inv(Z);
    0119       % ... TODO ...
    0120       ALx = ALx + trace(pkx^2*umatk*invZ-pkx*umatk);
    0121     end
    0122     end
    0123
    0124     % matrix constraints - pen/bar
    0125     if ~isempty(obj.Aindphi)
    0126     for k=obj.Aindphi
    0127       pkx=obj.PA(k);  % I used to use 2*         !!!!!!!!
    0128       % TODO need to map the matrix first! - is it correct???
    0129       kuser=obj.Amap(k);
    0130       [Akuserx, userdata] = obj.mconfun(x, Y, kuser, userdata);
    0131       Akx = obj.Ashift(k)*speye(size(Akuserx)) + obj.Amlt(k) .* Akuserx;
    0132       umatk=obj.UA{k};
    0133
    0134       %Z=(pkx*speye(size(Akx))-Akx);
    0135       Z=(pkx*speye(size(Akx))-Akx);
    0136       [R,iii]=chol(Z);   % to double check that it is pos. def!
    0137       if (iii~=0)
    0138         ALx=Inf;
    0139         break;
    0140       end
    0141       invZ=inv(Z);
    0142       % ... TODO ...
    0143       %ALx = ALx + trace(pkx^2*umatk*invZ-pkx*umatk);
    0144       ALx = ALx + pkx^2*(umatk(:)'*invZ(:))-pkx*trace(umatk);
    0145     end
    0146     end
    0147
    0148     % store user's data back in the object if it got changed
    0149     obj.userdata=userdata;
    0150
    0151     % update ticker
    0152     obj.ALx = ALx;
    0153     obj.ALxtck = obj.ticker;
    0154
    0155     % update stats
    0156     obj.stats_ncall_alx = obj.stats_ncall_alx + 1;
    0157     obj.stats_time_alx = obj.stats_time_alx + cputime - starttime;
    0158   end
    0159