Home > source > @penlab > init.m

# init

## PURPOSE

Initialize lagrangian multipliers, penalty parameters etc.

## SYNOPSIS

function [] = init(obj, forceupdate)

## DESCRIPTION

     Initialize lagrangian multipliers, penalty parameters etc.
     If called without argument (or with forceupdate=true), restores default
     settings for all (not only for mis-shaped or empty items).

     TODO: if not forced and wrong --> warning?
     TODO: force transposing if column vector?
     problem with vectors 0x1 --> changed it, if empty then 0x0 size

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

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

## SUBFUNCTIONS

- function [p, nfactor] = p_check2(M, p)

## SOURCE CODE

    0001 % Initialize lagrangian multipliers, penalty parameters etc.
    0002 % If called without argument (or with forceupdate=true), restores default
    0003 % settings for all (not only for mis-shaped or empty items).
    0004 %
    0005 % TODO: if not forced and wrong --> warning?
    0006 % TODO: force transposing if column vector?
    0007 % problem with vectors 0x1 --> changed it, if empty then 0x0 size
    0008 function [] = init(obj, forceupdate)
    0009
    0010   if (nargin<2)
    0011     forceupdate = true;
    0012   end
    0013
    0014   % function constraints: Lagr. multipliers & penalty
    0015   if (forceupdate || ~isvector(obj.uxbox) || length(obj.uxbox)~=obj.Nxbox)
    0016     if (obj.Nxbox>0)
    0017       obj.uxbox = obj.allopts.uinit_box*ones(obj.Nxbox,1);
    0018     else
    0019       obj.uxbox = [];
    0020     end
    0021   end
    0022   obj.uxboxnew = ones(obj.Nxbox,1);
    0023
    0024   if (forceupdate || ~isvector(obj.uineq) || length(obj.uineq)~=obj.Nineq)
    0025     if (obj.Nineq>0)
    0026       obj.uineq = obj.allopts.uinit*ones(obj.Nineq,1);
    0027     else
    0028       obj.uineq = [];
    0029     end
    0030   end
    0031   obj.uineqnew = ones(obj.Nineq,1);
    0032
    0033   if (forceupdate || ~isvector(obj.ueq) || length(obj.ueq)~=obj.Neq)
    0034     if (obj.Neq>0)
    0035       obj.ueq = obj.allopts.uinit_eq*ones(obj.Neq,1);
    0036     else
    0037       obj.ueq = [];
    0038     end
    0039   end
    0040
    0041   if (forceupdate || ~isvector(obj.pxbox) || length(obj.pxbox)~=obj.Nxbox)
    0042     if (obj.Nxbox>0)
    0043       obj.pxbox = obj.allopts.pinit*ones(obj.Nxbox,1);
    0044     else
    0045       obj.pxbox = [];
    0046     end
    0047   end
    0048
    0049   if (forceupdate || ~isvector(obj.pineq) || length(obj.pineq)~=obj.Nineq)
    0050     if (obj.Nineq>0)
    0051       obj.pineq = obj.allopts.pinit*ones(obj.Nineq,1);
    0052     else
    0053       obj.pineq = [];
    0054     end
    0055   end
    0056
    0057   % matrix constraints: NYbox, NA
    0058   % MATRIX variables (Y) TODO   don't forget obj.PY...!!
    0059   if (forceupdate || ~isvector(obj.PYbox) || length(obj.PYbox)~=obj.NYbox)
    0060     if (obj.NYbox>0)
    0061       % same as below, should be length(Ybox used in Phi, not barrier)
    0062       pnew=1;  % TODO get me from the option settings!!
    0063       for k=1:obj.NYbox   %obj.Yboxindphi
    0064         kuser=obj.Yboxmap(k);
    0065         Ykuserx=obj.Y{kuser};
    0066         Ykx = obj.Yboxshift(k)*speye(size(Ykuserx)) + obj.Yboxmlt(k) .* Ykuserx;
    0067
    0068         %pkx=pnew(k); %obj.PA(k);
    0069         [pnew, nfactor] = p_check2(-Ykx, pnew);
    0070       end
    0071
    0072       % use the minima for all  TODO shoudln't it be maxima????
    0073       obj.PYbox=pnew*ones(obj.NYbox,1);
    0074
    0075     else
    0076       obj.PYbox = [];
    0077     end
    0078   end
    0079
    0080   % perhpas in addition to these, it should be tested if PA is sufficiently large and increase it as necessary...
    0081   if (forceupdate || ~isvector(obj.PA) || length(obj.PA)~=obj.NA)
    0082     if (obj.NA>0)
    0083       % TODO this should be infact length(obj.Aindphi)>0 rather than obj.NA>0
    0084       % needs to be set up big enough that the point is feasible
    0085       % start with one from option settins and assign the biggest one to all
    0086       pnew=1;  % TODO get me from the option settings!!
    0087       for k=obj.Aindphi
    0088         kuser=obj.Amap(k);
    0089         [Akuserx, obj.userdata] = obj.mconfun(obj.x, obj.Y, kuser, obj.userdata);
    0090         Akx = obj.Ashift(k)*speye(size(Akuserx)) + obj.Amlt(k) .* Akuserx;
    0091
    0092         %pkx=pnew(k); %obj.PA(k);
    0093         [pnew, nfactor] = p_check2(-Akx, pnew);
    0094       end
    0095
    0096       % use the minima for all  TODO shoudln't it be maxima????
    0097       obj.PA=pnew*ones(obj.NA,1);
    0098     else
    0099       obj.PA=[];
    0100     end
    0101   end
    0102
    0103   % force update always... (just for now)
    0104   if (forceupdate || isempty(obj.UA))
    0105     obj.UA = [];
    0106     obj.UAnew = [];
    0107     for k=1:obj.NA
    0108       kuser=obj.Amap(k);
    0109       [Akuserx, obj.userdata] = obj.mconfun(obj.x, obj.Y, kuser, obj.userdata);
    0110       obj.UA{k} = eye(size(Akuserx));
    0111       % does this get used???
    0112       obj.UAnew{k} = eye(size(Akuserx));
    0113     end
    0114   end
    0115
    0116   if (forceupdate || isempty(obj.UYbox))
    0117     obj.UYbox = [];
    0118     obj.UYboxnew = [];
    0119     for k=1:obj.NYbox
    0120       kuser=obj.Yboxmap(k);
    0121       obj.UYbox{k} = eye(size(obj.Y{kuser}));
    0122       obj.UYboxnew{k} = eye(size(obj.Y{kuser}));
    0123     end
    0124   end
    0125
    0126 end
    0127
    0128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0129 % TODO this is the same as from mpen_update ... put extra if needed!
    0130 % the difference is that there is no parameter to backtrack ...
    0131 %  ... probably just double it every time??
    0132 %  ... original Pennon adjusts rUpdate to control how fast it goes down
    0133 %
    0134 % check if the new penalty parameter is usable for the matrix
    0135 % constraints, i.e., check that M+p*I >=0 (positive semidefinite)
    0136 %   M - matrix of the constraint (M>=0) at the current point
    0137 %   p - new penalty parameter
    0138 %   xxxxpold - previous penalty parameter (to backtrack)
    0139 % returns
    0140 %   p - valid penalty argument
    0141 %   nfactor - how many times it was necessary to factorize
    0142 %     (nfactor>1 ==> p had to be shortened)
    0143 function [p, nfactor] = p_check2(M, p)
    0144
    0145   rSafetyFactor = 1.2;  % add 20% to the initial penalty to be on the safe side
    0146   %rSafetyFactor = 1;
    0147   rFactor=2;   % if need to refactorize, try rFactor*p
    0148   nfactor=0;
    0149
    0150   [n m] = size(M);
    0151   Missparse = n>10 && issparse(M) && nnz(M)<0.15*n*n;
    0152
    0153   if (Missparse)
    0154     perm=amd(M);
    0155     M=M(perm,perm);
    0156     I=speye(n,n);
    0157   else
    0158     % it is usually faster to compute with dense matrices in dense format
    0159     M=full(M);
    0160     I=eye(n,n);
    0161   end
    0162
    0163   p = p/rSafetyFactor;
    0164
    0165   [R,k] = chol(M+p*I);
    0166   nfactor=nfactor+1;
    0167   if (k==0)
    0168     % first match -> go
    0169     return;
    0170   end
    0171
    0172   while (k~=0)
    0173     p=rFactor*p;
    0174     [R,k] = chol(M+p*I);
    0175     %disp(sprintf('up   pert=%e (%i)',p,k));
    0176     nfactor=nfactor+1;
    0177     % add save bounds??
    0178   end
    0179
    0180   p = p*rSafetyFactor;
    0181
    0182 end
    0183
    0184
    0185
    0186