Home > source > @penlab > eval_alddx.m

# eval_alddx

## PURPOSE

As eval_alx(), evaluate hessian of the Augmented Lagrangian and Jacobian

## SYNOPSIS

function [status] = eval_alddx(obj)

## DESCRIPTION

     As eval_alx(), evaluate hessian of the Augmented Lagrangian and Jacobian
     of the equality constraints. The results are stored in obj.ALddx, obj.eqdx
     and tickers updated.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

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

## SOURCE CODE

    0001 % As eval_alx(), evaluate hessian of the Augmented Lagrangian and Jacobian
    0002 % of the equality constraints. The results are stored in obj.ALddx, obj.eqdx
    0003 % and tickers updated.
    0004 function [status] = eval_alddx(obj)
    0005
    0006   % status TODO
    0007   status = 0;
    0008
    0009   % TODO remove!
    0010   if (true || obj.ALddxtck < obj.ticker)
    0011     starttime = cputime;
    0012
    0013     % create local copies of obj.x,obj.Y to avoid checking repetitively if they
    0014     % are 'up-to-date' with obj.xall, they get used many times in calls to
    0015     % user's functions
    0016     x=obj.x;
    0017     Y=obj.Y;
    0018     % same for obj.userdata, it might be quite expensive to get&store it in
    0019     % the object
    0020     userdata=obj.userdata;
    0021
    0022     % TODO reuse it + 'pointchanged' flag
    0023     % function inequal & equal
    0024     if (obj.NgNLN + obj.NgLIN>0)
    0025       [gx, userdata] = obj.confun(x, Y, userdata);
    0026       ineqx = obj.ineqshift + obj.ineqmlt .* gx(obj.ineqmap);
    0027       obj.eqx = obj.eqshift + gx(obj.eqmap);
    0028
    0029       [gdx, userdata] = obj.congrad(x, Y, userdata);
    0030       % because of stupid spdiags() which doesn't accept ([],0,0,0)
    0031       if (obj.Nineq>0)
    0032         ineqdx = gdx(:,obj.ineqmap) * spdiags(obj.ineqmlt,0,obj.Nineq,obj.Nineq);
    0033       else
    0034         ineqdx = [];
    0035       end
    0036       % TODO update ticker eqdx
    0037       obj.eqdx = gdx(:,obj.eqmap);
    0038     end
    0039
    0040     % prepare Lmlt for all function 2nd derivatives fddx + Lmlt*gddx
    0041     % (in user's indices, including both equal & inequal)
    0042     % AMPL wants to have the full length ... but in fact only obj.NgNLN
    0043     % should be enough.  ==> try with all and then reduce it later TODO
    0044     % TODO vectorize
    0045     % TODO this is right only if all phiind! (which will be probably true,
    0046     % as no NLN constr should go to the barrier ... why should they?)
    0047     if (obj.Nineq>0)
    0048       ind=1:obj.Nineq;
    0049       Lmlt_long=obj.ineqmlt(ind).*obj.uineq(ind).*obj.phi2_D(ineqx(ind)./obj.pineq(ind));
    0050       % gather these Lmlt belonging to the same "body" of the constraint
    0051       % (sparse sums element which belong to the same index)
    0052       Lmlt = full(sparse(obj.ineqmap(ind),1,Lmlt_long,obj.NgNLN+obj.NgLIN,1));
    0053     else
    0054       Lmlt = zeros(obj.NgNLN+obj.NgLIN,1);
    0055     end
    0056 %    Lmlt = zeros(obj.NgNLN+obj.NgLIN,1);
    0057 %    for k=1:obj.Nineq
    0058 %      kuser=obj.ineqmap(k);
    0059 %      Lmlt(kuser) = Lmlt(kuser) + obj.ineqmlt(k)*obj.uineq(k)*obj.phi2_D(ineqx(k)/obj.pineq(k));
    0060 %    end
    0061     for k=1:obj.Neq
    0062       kuser=obj.eqmap(k);
    0063       Lmlt(kuser)=obj.ueq(k);
    0064     end
    0065
    0066     % get Lagrangian directly or compose it
    0067     if (~isempty(obj.lagrhess))
    0068       [ALddx, userdata] = obj.lagrhess(x,Y,Lmlt,userdata);
    0069       if (isempty(ALddx))
    0070         ALddx=sparse(obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz);
    0071       end
    0072     else
    0073       % compose Lagrangian
    0074       [ALddx, userdata] = obj.objhess(x,Y,userdata);
    0075       if (isempty(ALddx))
    0076         ALddx=sparse(obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz);
    0077       end
    0078       for kuser=1:obj.NgNLN
    0079         [gddx, userdata] = obj.conhess(x,Y,kuser,userdata);
    0080         ALddx = ALddx + Lmlt(kuser)*gddx;
    0081       end
    0082     end
    0083
    0084     %%% Add dyadic product of function inequalitites (and box)
    0085
    0086     % box constraints
    0087     if (obj.Nxbox>0)
    0088       xboxx = obj.xboxshift + obj.xboxmlt .* obj.xall(obj.xboxmap);
    0089       % this adds only diagonal scaling as grad = (0,...,1,...0) or with -1
    0090       % so it will be either:
    0091       %   + u(i)/p(i)*phi2_D2(g(i)/p(i))   for phi/bar
    0092       % or
    0093       %   + u(i)*p(i)*phibar_D2(g(i))      for barrier
    0094       % or nothing if not bounded, don't forget on both inequalitites
    0095       diagxbox=zeros(obj.Nx+obj.NYnnz,1);
    0096       % TODO vectorize
    0097       ind=obj.xboxindbar;
    0098       if ~isempty(ind)
    0099       for k=ind
    0100         kuser=obj.xboxmap(k);
    0101         diagxbox(kuser) = diagxbox(kuser) + obj.uxbox(k)*obj.pxbox(k)*obj.phibar_D2(xboxx(k));
    0102       end
    0103       end
    0104
    0105       ind=obj.xboxindphi;
    0106       if ~isempty(ind)
    0107       for k=ind
    0108         kuser=obj.xboxmap(k);
    0109         diagxbox(kuser) = diagxbox(kuser) + obj.uxbox(k)/obj.pxbox(k)*obj.phi2_D2(xboxx(k)/obj.pxbox(k));
    0110       end
    0111       end
    0112       ALddx = ALddx + spdiags(diagxbox,0,obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz);
    0113     end
    0114
    0115     % function inequalitites
    0116     %ind=obj.ineqindbar;
    0117     %... TODO
    0118
    0119     ind=obj.ineqindphi;
    0120     if (~isempty(ind))
    0121       coef=obj.uineq(ind)./obj.pineq(ind).*obj.phi2_D2(ineqx(ind)./obj.pineq(ind));
    0122       ALddx = ALddx + ineqdx(:,ind) * spdiags(coef,0,length(ind),length(ind)) * ineqdx(:,ind)';
    0123
    0124     end
    0125 %    for k=ind
    0126       %k
    0127       %nnz(ineqdx(:,k))
    0128 %      ALddx = ALddx + ((obj.uineq(k)/obj.pineq(k)*obj.phi2_D2(ineqx(k)/obj.pineq(k)))*ineqdx(:,k))*ineqdx(:,k)';
    0129       %disp('k done')
    0130       %nnz(ALddx)
    0131 %    end
    0132
    0133     % matrix variable - log barrier (strict feasibility)
    0134     if ~isempty(obj.Yboxindbar)
    0135     for k=obj.Yboxindbar
    0136       % convert the matrix box constraint to the form:   +/-Y +/-bound >=0
    0137       pkx=obj.PYbox(k);
    0138       Ykx = Y{obj.Yboxmap(k)};
    0139       Akx=-obj.Yboxshift(k)*speye(size(Ykx)) - obj.Yboxmlt(k)*Ykx;
    0140
    0141       % can assume that ALx was computed before ALdx --> Akx must be pos. def.
    0142       chol(Akx);
    0143       invAkx = full(inv(Akx));
    0144
    0145       % compute elements of 2nd derivative of the penalty term:
    0146       %   d2/dxi dxj ... = p*trace( (A^-1) * (d/dxi A) * (A^-1) * (d/dxj A) )
    0147
    0148       mlt=-obj.Yboxmlt(k);  % +/-1
    0149       mapper=obj.vec2Ymap{obj.Yboxmap(k)};
    0150       dim=mapper.dim;
    0151       lAdep=mapper.nelem;
    0152       offset=obj.Nx + mapper.xmap(1) - 1;
    0153       irow=mapper.irow;
    0154       icol=mapper.icol;
    0155
    0156       Akdixall=cell(lAdep,1);
    0157       % store everything in a dense matrix (kernel) and copy it to the Hessian
    0158       % once when finished
    0159       Akddx=zeros(lAdep,lAdep);
    0160       for ii=1:lAdep
    0161         if (irow(ii)==icol(ii))
    0162           % diagonal element
    0163           Akdix = sparse(irow(ii),icol(ii),mlt,dim,dim);
    0164         else
    0165           % nondiag element --> add two
    0166           Akdix = sparse([irow(ii),icol(ii)],[icol(ii),irow(ii)],[mlt,mlt],dim,dim);
    0167         end
    0168
    0169         Akdixall{ii}=Akdix;
    0170
    0171         % compute all elements under the diagonal
    0172         for jj=1:ii-1
    0173           hij=0.5*pkx*mextrdsdsmat(invAkx, Akdix, invAkx, Akdixall{jj});
    0174           Akddx(ii,jj)=hij;
    0175           Akddx(jj,ii)=hij;
    0176         end
    0177
    0178         % compute diagonal element jj=ii
    0179         hij=0.5*pkx*mextrdsdsmat(invAkx, Akdix, invAkx, Akdix);
    0180         Akddx(ii,ii)=hij;
    0181
    0182       end
    0183       Adep=offset+[1:lAdep];
    0184       % copy the dense kernel back to ALddx
    0185       % ALddx=ALddx + Akddx(obj.Adep{kuser},obj.Adep{kuser}); ... need other way
    0186       Akddx_expand=sparse(obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz);
    0187       %Akddx_expand(obj.Adep{kuser},obj.Adep{kuser})=Akddx;
    0188       Akddx_expand(Adep,Adep)=Akddx;
    0189       ALddx=ALddx + Akddx_expand;
    0190
    0191     end
    0192     end
    0193
    0194     % matrix variable - pen/bar
    0195     if ~isempty(obj.Yboxindphi)
    0196     for k=obj.Yboxindphi
    0197       pkx=obj.PYbox(k);            %2*p(sdpdata.Ng+k);
    0198       Ykx = obj.Y{obj.Yboxmap(k)};
    0199       Akx=obj.Yboxshift(k)*speye(size(Ykx)) + obj.Yboxmlt(k)*Ykx;
    0200       umatk=obj.UYbox{k};
    0201
    0202       %Z=(pkx*speye(size(Akx))-Akx);
    0203       Z=(pkx*speye(size(Akx))-Akx);
    0204       %invZ=inv(Z);
    0205       invZ=full(inv(Z));
    0206       pZUZ=pkx^2*invZ*umatk*invZ;
    0207
    0208       mlt=obj.Yboxmlt(k);  % +/-1
    0209       mapper=obj.vec2Ymap{obj.Yboxmap(k)};
    0210       dim=mapper.dim;
    0211       lAdep=mapper.nelem;
    0212       offset=obj.Nx + mapper.xmap(1) - 1;
    0213       irow=mapper.irow;
    0214       icol=mapper.icol;
    0215
    0216       Akdixall=cell(lAdep,1);
    0217       % store everything in a dense matrix (kernel) and copy it to the Hessian
    0218       % once when finished
    0219       Akddx=zeros(lAdep,lAdep);
    0220       for ii=1:lAdep
    0221         if (irow(ii)==icol(ii))
    0222           % diagonal element
    0223           Akdix = sparse(irow(ii),icol(ii),mlt,dim,dim);
    0224         else
    0225           % nondiag element --> add two
    0226           Akdix = sparse([irow(ii),icol(ii)],[icol(ii),irow(ii)],[mlt,mlt],dim,dim);
    0227         end
    0228
    0229         Akdixall{ii}=Akdix;
    0230
    0231         % compute all elements under the diagonal
    0232         for jj=1:ii-1
    0233           hij=mextrdsdsmat(pZUZ, Akdix, invZ, Akdixall{jj});
    0234           Akddx(ii,jj)=hij;
    0235           Akddx(jj,ii)=hij;
    0236         end
    0237
    0238         % compute diagonal element jj=ii
    0239         hij=mextrdsdsmat(pZUZ, Akdix, invZ, Akdix);
    0240         Akddx(ii,ii)=hij;
    0241
    0242       end
    0243       Adep=offset+[1:lAdep];
    0244       % copy the dense kernel back to ALddx
    0245       % ALddx=ALddx + Akddx(obj.Adep{kuser},obj.Adep{kuser}); ... need other way
    0246       Akddx_expand=sparse(obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz);
    0247       %Akddx_expand(obj.Adep{kuser},obj.Adep{kuser})=Akddx;
    0248       Akddx_expand(Adep,Adep)=Akddx;
    0249       ALddx=ALddx + Akddx_expand;
    0250
    0251     end
    0252     end
    0253
    0254     % matrix constraints - pen/bar
    0255     if ~isempty(obj.Aindphi)
    0256     for k=obj.Aindphi
    0257       pkx=obj.PA(k);  % I used to use 2*         !!!!!!!!
    0258       % TODO need to map the matrix first! - is it correct???
    0259       kuser=obj.Amap(k);
    0260       [Akuserx, userdata] = obj.mconfun(x, Y, kuser, userdata);
    0261       Akx = obj.Ashift(k)*speye(size(Akuserx)) + obj.Amlt(k) .* Akuserx;
    0262       umatk=obj.UA{k};
    0263
    0264       Z=(pkx*speye(size(Akx))-Akx);
    0265       %invZ=inv(Z);
    0266       invZ=full(inv(Z));
    0267       %invZ=inv(full(Z));
    0268       pZUZ=pkx^2*invZ*umatk*invZ;
    0269
    0270       % MAPPING & transformation !!!
    0271       Adep=obj.Adep{kuser};
    0272       lAdep=length(Adep);
    0273       Akdixall=cell(lAdep,1);
    0274       % store everything in a dense matrix (kernel) and copy it to the Hessian
    0275       % once when finished
    0276       %Akddx=zeros(lAdep,lAdep);
    0277       Akddx_nnz=lAdep*(lAdep+1)/2;
    0278       Akddx_val=zeros(Akddx_nnz,1);
    0279       Akddx_row=zeros(Akddx_nnz,1);
    0280       Akddx_col=zeros(Akddx_nnz,1);
    0281       idx = 0;
    0282       for ii=1:lAdep
    0283         i=Adep(ii);
    0284         [Akdix, userdata] = obj.mcongrad(x,Y,kuser,i,userdata);
    0285         if (obj.Amlt(k)<0)    % TODO correct????
    0286           % must be -1, otherwise it would be 1
    0287           Akdix = -Akdix;
    0288         end
    0289         %Akdix = obj.Amlt(k)*Akdix;
    0290         Akdixall{ii}=Akdix;
    0291         Akddx_row(idx+1:idx+ii)=i;
    0292         Akddx_col(idx+1:idx+ii)=Adep(1:ii);
    0293         % compute all elements under the diagonal
    0294 %        for jj=1:ii-1
    0295           %---Akdjx=Akdixall{jj};
    0296
    0297           %AZA=Akdix*invZ*Akdjx;
    0298           %hij = trace(pZUZ*AZA);
    0299           %%hij=trdsdsmat(pZUZ, Akdix, invZ, Akdjx);
    0300           %---hij=mextrdsdsmat(pZUZ, Akdix, invZ, Akdjx);
    0301 %          hij=mextrdsdsmat(pZUZ, Akdix, invZ, Akdixall{jj});
    0302           %Akddx(ii,jj)=hij;
    0303           %Akddx(jj,ii)=hij;
    0304 %          Akddx_val(idx+jj)=hij;
    0305 %        end
    0306
    0307         % compute diagonal element jj=ii
    0308 %        hij=mextrdsdsmat(pZUZ, Akdix, invZ, Akdix);
    0309         %Akddx(ii,ii)=hij;
    0310 %        Akddx_val(idx+ii)=hij;
    0311
    0312         vec = mextrcolumn(pZUZ,Akdix,invZ,Akdixall(1:ii));
    0313         Akddx_val(idx+1:idx+ii)=vec;
    0314
    0315         idx = idx+ii;
    0316       end
    0317
    0318       % add 2nd order matrix constraints derivatives
    0319       % either by lagrangian or element by element
    0320       Akddx_lagr = [];
    0321       if (kuser<=obj.NANLN)
    0322         if (~isempty(obj.mconlagrhess))
    0323           [Akddx_lagr, userdata] = obj.mconlagrhess(x,Y,kuser,pZUZ,userdata);
    0324           if (obj.Amlt(k)~=1)
    0325             Akddx_lagr = obj.Amlt(k)*Akddx_lagr;
    0326           end
    0327         else
    0328           % lagrangian is not present, do it element by element
    0329           %Akddx_val2=zeros(Akddx_nnz,1);
    0330           idx = 0;
    0331           for ii=1:lAdep
    0332             i=Adep(ii);
    0333             for jj=1:ii-1
    0334               %[Akddijx, userdata] = obj.mconhess(x,Y,kuser,i,Adep(jj),userdata);
    0335               [Akddijx, userdata] = obj.mconhess(x,Y,kuser,Adep(jj),i,userdata);
    0336               if (~isempty(Akddijx))
    0337                 % TODO better
    0338                 % don't forget a possible transformation with Akddijx,
    0339                 % can be applied on the trace (instead of on Akddijx
    0340                 % directly)
    0341                 hij=obj.Amlt(k)*trace(Akddijx*pZUZ);
    0342                 %Akddx(ii,jj)=Akddx(ii,jj)+hij;
    0343                 %Akddx(jj,ii)=Akddx(jj,ii)+hij;
    0344                 Akddx_val(idx+jj)=Akddx_val(idx+jj)+hij;
    0345                 %Akddx_val2(idx+jj)=Akddx_val2(idx+jj)+hij;
    0346               end
    0347             end
    0348             [Akddijx, userdata] = obj.mconhess(x,Y,kuser,i,i,userdata);
    0349             if (~isempty(Akddijx))
    0350               % TODO better
    0351               hij=obj.Amlt(k)*trace(Akddijx*pZUZ);
    0352               %Akddx(ii,ii)=Akddx(ii,ii)+hij;
    0353               Akddx_val(idx+ii)=Akddx_val(idx+ii)+hij;
    0354               %Akddx_val2(idx+ii)=Akddx_val2(idx+ii)+hij;
    0355             end
    0356             idx = idx+ii;
    0357           end
    0358         end
    0359       end
    0360       % only lower triangle -> copy to full
    0361       %Akddx = tril(Akddx) + tril(Akddx,-1)';
    0362       % copy the dense kernel back to ALddx
    0363       % ALddx=ALddx + Akddx(obj.Adep{kuser},obj.Adep{kuser}); ... need other way
    0364
    0365       %%Akddx_expand=sparse(obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz);
    0366       %Akddx_expand(obj.Adep{kuser},obj.Adep{kuser})=Akddx;
    0367       %%Akddx_expand(Adep,Adep)=Akddx;
    0368       %ALddx=ALddx + Akddx_expand;
    0369       Akddx_expand = sparse(Akddx_row,Akddx_col,Akddx_val,obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz,lAdep*lAdep);
    0370       Akddx_expand = Akddx_expand + tril(Akddx_expand,-1)';
    0371       if (isempty(Akddx_lagr))
    0372         %Akddx_lagr = sparse(Akddx_row,Akddx_col,Akddx_val2,obj.Nx+obj.NYnnz,obj.Nx+obj.NYnnz,lAdep*lAdep);
    0373         %Akddx_lagr = Akddx_lagr + tril(Akddx_lagr,-1)';
    0374         ALddx=ALddx + Akddx_expand;
    0375         %ALddx=ALddx + Akddx_expand + Akddx_lagr;
    0376       else
    0377         ALddx=ALddx + Akddx_expand + Akddx_lagr;
    0378       end
    0379       % + add 2nd derivatives
    0380     end
    0381     end
    0382
    0383     % store user's data back in the object if it got changed
    0384     obj.userdata=userdata;
    0385
    0386     % update ticker
    0387     obj.ALddx = ALddx;
    0388     obj.ALddxtck = obj.ticker;
    0389
    0390     % update stats
    0391     obj.stats_ncall_alddx = obj.stats_ncall_alddx + 1;
    0392     obj.stats_time_alddx = obj.stats_time_alddx + cputime - starttime;
    0393
    0394   end
    0395