Home > source > feasm.m

# feasm

## PURPOSE

Compute a perturbation needed for a given symmetric matrix to be pos. def.

## SYNOPSIS

function [pert] = feasm(M, lb, ub, pstart, pstop)

## DESCRIPTION

     Compute a perturbation needed for a given symmetric matrix to be pos. def.
     I.e., compute a lower estimate of the smallest (negative) eigenvalue (if any).
     Input:
       M ... symmetric matrix to find 'pert>=0' so that M+pert*I>0 (pos.def.)
       lb ... (>=0) lower bound on pert, i.e., find pert>=lb>=0,
              e.g., current known feasibility of other matrix constraints
       ub ... (>0, Inf if not known), upper bound on pert which is known that
              M+ub*I>0, typically the penalty parameter
       pstart ... if lb==0, pstart is the first perturbation to try if M not posdef
       pstop ... tolerance of the requested perturbation, the exact perturbation
              will be in interval [pert,pert+ptol]
     Output:
       pert ... estimated perturbation within the given tolerance unless
              M+lb*I>=0 is already pos def

     TODO
       how to make pstart&ptol dynamic to save some time...

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

## SOURCE CODE

    0001 % Compute a perturbation needed for a given symmetric matrix to be pos. def.
    0002 % I.e., compute a lower estimate of the smallest (negative) eigenvalue (if any).
    0003 % Input:
    0004 %   M ... symmetric matrix to find 'pert>=0' so that M+pert*I>0 (pos.def.)
    0005 %   lb ... (>=0) lower bound on pert, i.e., find pert>=lb>=0,
    0006 %          e.g., current known feasibility of other matrix constraints
    0007 %   ub ... (>0, Inf if not known), upper bound on pert which is known that
    0008 %          M+ub*I>0, typically the penalty parameter
    0009 %   pstart ... if lb==0, pstart is the first perturbation to try if M not posdef
    0010 %   pstop ... tolerance of the requested perturbation, the exact perturbation
    0011 %          will be in interval [pert,pert+ptol]
    0012 % Output:
    0013 %   pert ... estimated perturbation within the given tolerance unless
    0014 %          M+lb*I>=0 is already pos def
    0015 %
    0016 % TODO
    0017 %   how to make pstart&ptol dynamic to save some time...
    0018 function [pert] = feasm(M, lb, ub, pstart, pstop)
    0019
    0020   nfactor=0;
    0021
    0022   % Get ordering if M really sparse & shuffle M to avoid recomputing
    0023   [n m] = size(M);
    0024   Missparse = n>10 && issparse(M) && nnz(M)<0.15*n*n;
    0025      % nnz() used directly and behind issparse() because otherwise it
    0026      % counts nonzero elements even in the dense matrix
    0027
    0028   if (Missparse)
    0029     perm=amd(M);
    0030     M=M(perm,perm);
    0031     I=speye(n,n);
    0032   else
    0033     % it is usually faster to compute with dense matrices in dense format
    0034     M=full(M);
    0035     I=eye(n,n);
    0036   end
    0037
    0038   % first quick try
    0039   pert=max(0.0, lb);
    0040   [R,p] = chol(M+pert*I);
    0041   %pert
    0042   nfactor=nfactor+1;
    0043   if (p==0)
    0044     % first match -> go
    0045     return;
    0046   end
    0047
    0048   % guess a reasonable step how to increase the perturbation
    0049   if (ub<Inf)
    0050     % perturbation with ub should be pos. def; with 'pert' wasnt --> pert<ub
    0051     %pstep=min(20,max((ub-pert)/4, 2));
    0052     pstep=min(50,max(1.01*realpow(ub/pert,1/4),2));
    0053       % in 4 steps jump above the known perturbation
    0054   else
    0055     pstep=10;
    0056   end
    0057   %pstep
    0058   pertlowstart = pert;
    0059   pert=max(pert,pstart/pstep);
    0060
    0061   while (p~=0)
    0062     pertlow = pert;
    0063     pert=pert*pstep;
    0064     [R,p] = chol(M+pert*I);
    0065     %disp(sprintf('up   pert=%e (%i)',pert,p));
    0066     nfactor=nfactor+1;
    0067   end
    0068   pertup=pert;
    0069
    0070   if (nfactor==2)
    0071     % this is probably not necessary but it the upper loop finishes after one
    0072     % step, perlow might be actually higher than the one which failed in the
    0073     % "first quick try"
    0074     pertlow=pertlowstart;
    0075   end
    0076   %nfactor
    0077
    0078   pstopn = max(pstop,norm(pertlow)*pstop);
    0079
    0080   while (pertup-pertlow>pstopn)
    0081     %pert=pertlow + (pertup-pertlow)/3;
    0082     pert=(pertlow+pertup)/2;
    0083     %pert=realsqrt(pertlow*pertup);
    0084     [R,p] = chol(M+pert*I);
    0085     %disp(sprintf('down pert=%e (%i)',pert,p));
    0086     nfactor=nfactor+1;
    0087     if (p==0)
    0088       pertup=pert;
    0089     else
    0090       pertlow=pert;
    0091     end
    0092   end
    0093   pert = pertup;   % return upper bound on the perturbation (=overestimate)
    0094   % (=lower bound on the (negative) eigenvalue)
    0095
    0096   %nfactor
    0097 end
    0098