Home > source > @penlab > mpen_update.m

# mpen_update

## PURPOSE

Penalty for matrix constraints update

## SYNOPSIS

function [] = mpen_update(obj)

## DESCRIPTION

     Penalty for matrix constraints update

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SUBFUNCTIONS

- function [p] = p_update(p, pfactor, pborder, pmin)
- function [p, nfactor] = p_check(M, p, pold)

## SOURCE CODE

    0001 % Penalty for matrix constraints update
    0002 function  [] = mpen_update(obj)
    0003
    0004   % add barrier penalty update as well...
    0005   idx=obj.Yboxindbar;
    0006   pnew=obj.PYbox;   % to get the right shape
    0007   pnew(idx)=p_update(obj.PYbox(idx), obj.allopts.mpenalty_update, obj.allopts.mpenalty_border, obj.allopts.mpenalty_min);
    0008
    0009   if ~isempty(obj.Yboxindbar)
    0010   for k=obj.Yboxindbar
    0011     Ykx = obj.Y{obj.Yboxmap(k)};
    0012     Akx=obj.Yboxshift(k)*speye(size(Ykx)) + obj.Yboxmlt(k)*Ykx;
    0013
    0014     %pkx=pnew(k); %obj.PA(k);
    0015     [pnew(k), nfactor] = p_check(-Akx, pnew(k), obj.PYbox(k));
    0016   end
    0017   end
    0018   obj.PYbox(idx) = max(pnew(idx));
    0019
    0020   idx=obj.Yboxindphi;
    0021   pnew=obj.PYbox;   % to get the right shape
    0022   pnew(idx)=p_update(obj.PYbox(idx), obj.allopts.mpenalty_update, obj.allopts.mpenalty_border, obj.allopts.mpenalty_min);
    0023
    0024   if ~isempty(obj.Yboxindphi)
    0025   for k=obj.Yboxindphi
    0026     Ykx = obj.Y{obj.Yboxmap(k)};
    0027     Akx=obj.Yboxshift(k)*speye(size(Ykx)) + obj.Yboxmlt(k)*Ykx;
    0028
    0029     %pkx=pnew(k); %obj.PA(k);
    0030     [pnew(k), nfactor] = p_check(-Akx, pnew(k), obj.PYbox(k));
    0031   end
    0032   end
    0033   obj.PYbox(idx) = max(pnew(idx));
    0034
    0035
    0036   idx=obj.Aindphi;
    0037   %pold=obj.PA;
    0038   pnew=obj.PA;   % to get the right shape
    0039   pnew(idx)=p_update(obj.PA(idx), obj.allopts.mpenalty_update, obj.allopts.mpenalty_border, obj.allopts.mpenalty_min);
    0040      % p, pfactor, pborder, pmin
    0041
    0042   % check that each matrix constraints is still feasible with the new penalty parameter
    0043   if ~isempty(obj.Aindphi)
    0044   for k=obj.Aindphi
    0045     kuser=obj.Amap(k);
    0046     [Akuserx, obj.userdata] = obj.mconfun(obj.x, obj.Y, kuser, obj.userdata);
    0047     Akx = obj.Ashift(k)*speye(size(Akuserx)) + obj.Amlt(k) .* Akuserx;
    0048
    0049     %pkx=pnew(k); %obj.PA(k);
    0050     [pnew(k), nfactor] = p_check(-Akx, pnew(k), obj.PA(k));
    0051   end
    0052   end
    0053
    0054   % use the minima for all
    0055   obj.PA(idx) = max(pnew(idx));
    0056
    0057 end
    0058
    0059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0060 % penalty vector update according to the rules
    0061 %   p - vector of original penalty parameters
    0062 %   pfactor - update factor for "normal" operation
    0063 %   pborder - where to switch to slower update
    0064 %   pmin - minimal penalty which will be used
    0065 % OR do it up there??
    0066 function [p] = p_update(p, pfactor, pborder, pmin)
    0067   idxnormal = p>pborder;
    0068
    0069   p(idxnormal) = p(idxnormal)*pfactor;
    0070   p(~idxnormal) = p(~idxnormal)*0.9;
    0071
    0072   % don't go below pmin
    0073   p = max(p,pmin);
    0074
    0075 end
    0076
    0077 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0078 % check if the new penalty parameter is usable for the matrix
    0079 % constraints, i.e., check that M+p*I >=0 (positive semidefinite)
    0080 %   M - matrix of the constraint (M>=0) at the current point
    0081 %   p - new penalty parameter
    0082 %   pold - previous penalty parameter (to backtrack)
    0083 % returns
    0084 %   p - valid penalty argument
    0085 %   nfactor - how many times it was necessary to factorize
    0086 %     (nfactor>1 ==> p had to be shortened)
    0087 function [p, nfactor] = p_check(M, p, pold)
    0088
    0089   rFactor=0.75;   % if need to refactorize, prefer new penalty parameter
    0090   nfactor=0;
    0091
    0092   [n m] = size(M);
    0093   Missparse = n>10 && issparse(M) && nnz(M)<0.15*n*n;
    0094
    0095   if (Missparse)
    0096     perm=amd(M);
    0097     M=M(perm,perm);
    0098     I=speye(n,n);
    0099   else
    0100     % it is usually faster to compute with dense matrices in dense format
    0101     M=full(M);
    0102     I=eye(n,n);
    0103   end
    0104
    0105   [R,k] = chol(M+p*I);
    0106   nfactor=nfactor+1;
    0107   if (k==0)
    0108     % first match -> go
    0109     return;
    0110   end
    0111
    0112   while (k~=0)
    0113     p=rFactor*p + (1-rFactor)*pold;
    0114     [R,k] = chol(M+p*I);
    0115     %disp(sprintf('up   pert=%e (%i)',p,k));
    0116     nfactor=nfactor+1;
    0117     % add save bounds??
    0118   end
    0119
    0120 end
    0121
    0122