Documentation for
Trust Region Subroutine Algorithm
Charles Fortin and Henry Wolkowicz
  Dept. Comb. & Opt., Univ.  of Waterloo
This file contains the documentation for the trust region subproblem algorithm.
This file is available with URL:

 http://orion.math.uwaterloo.ca/~hwolkowi/henry/software/trustreg.d/docs.html

 Variable Definitions,  Tolerances,  Initialization, 





 Call:  
  [xopt,fopt,lopt,lambdamin,muup,mulow,iter,flag] = newtrust(a,s,dgaptol) 


 solves:         min      q(x) = x'Ax - 2a'x
                such that    s^2-x'x >=0        (Lagr. mult. lopt <=0)


 which is the trust region subproblem. This can be used within a 
trust region algorithm for min f(x). It is specifically designed for
 large sparse problems. The algorithm is based on the one developed in the paper

A Semidefinite Framework for Trust Region Subproblems with Applications to
Large Scale Minimization
(Franz Rendl and Henry Wolkowicz, Math. Progr. 1997., Math. Review


                                                            
files needed: Afun.m,Dfun.m,conj_grad.m, and this file newtrust.m

 BEFORE USING NEWTRUST, THE FOLLOWING INSTRUCTION MUST BE ENTERED

                            >> global A D counter
                                                                





VARIABLES

Back to top of Documentation file

   INPUT: 

         A      is an n by n symmetric matrix. It is  assumed to be in 
                    sparse format. It is a global variable in order
                    to find the minimum eigenvalue using a Lanczos routine.
               
         a      is a column vector - it is changed to full format below.
         s>0    is a scalar. It is the radius of the trust region.
        dgaptol is the desired duality gap tolerance.
                    (The user should decrease it as one gets closer
                     to the optimum  of min f(x))
        par     is a vector that contains initial estimates if
they are available, par=[lambda lambdat t], where lambda is an eigenvalue
of A, t is the parameter in D(t) and lambdat is corresponding
eigenvalue.  Also, xparvec on input is eigenvector corresponding
to lambda while tparvec it eigenvector of D(t) corresponding to
lambdat.  ?????not done yet????


  OUTPUT: 

      xopt     approximate optimum of the quadratic (n-vector).
      fopt     approximate optimal value of the quadratic.
               ...   value is < gap(1,1)+dgaptol if assigned a value;
               ...   takes value 99999 if not assigned a value.
      lopt ... dual opt multiplier vector,
 lambdamin ... minimum eigenvalue of A
    iter   ... number of iterations (eigenvalue calls)
    flag   ... a boolean variable to know if solution is 
                   in the interior of the ball.
                   (flag ==2) and if the hard case (case 2) 
                     (or near hard case (case2)) occurs 
                   (flag ==1) and otherwise (flag == 0).

 

  OTHER VARIABLES OF INTEREST:

     iter  ... number of iterations - success occurred if < iterbnd
     gap   ... gap, optimal value is in interval [gap(1,1),gap(1,2)]
 
              on output, vector gap is:
 
              gap=[mulow muup]   
                         lower and upper bounds on mu^* (the optimum
                         of the quadratic objective )
                  ptbad(1,3) ptgood(1,3) 
                    lower and upper bounds on lambda^* (multiplier)
                  low high          lower and upper bounds on t^*
                  tarn1f t0f]
    tarn1f ... = cputime-tarn10;  
                    save initial eigenvalue cputime for statistics
       t0f ... = cputime - t0;  total computation time
      topt ... right t-value to generate opt. solution;
                      t=99999 means that lopt=0 used or 
                      0 iterations with no optimal t,
                        i.e. an unconstrained optimum was found.
                      t is parameter in eigenvalue problem.
 
 ptgood, ptbad variables: ... used for saving interpolation data
       bestxf  ... current best approximation to the x-optimum
       bestfxf  ... current best approximation to the optimum of f(x) 
                      corresponding to bestxf only (not muup)
 
       global A D counter
                                                     




 

 
              TOLERANCES                             
          Back to top of documentation file
  Tolerances - these can and should be changed by the user.
 
     The tolerance for the duality gap is for significant figures,
      i.e. it is a relative duality gap.
      The significant figures of agreement are computed for
        the duality gap as follows:
       signif = max( -log_10  |primal-dual|/(|primal|+1)  , 0  )
 
     The tolerances for ***LANCZOS** should be significantly higher
          than those used here, e.g. 4 here and 7 for Lanczos
 
  Stopping criteria:  (i)   feasibility,      ||x|| > (1+normtol)s
             (ii)  feasibility    dbestgap= (bestfxf-mulow)/(|bestfxf|+1)
                                                          > 2*dgaptol
         bestfxf is the current best opt. value corresponding to
         the current best feasible point bestxf, while muup is the current
         best upper bound of the optimal value.
              (iii)  duality gap,      |gap(2)-gap(1)|/(|gap(2)|+1)
                                                               > dgaptol
                      (iv) iteration count,  iter <= iterbnd
                      (v) interval of uncertainty for t is > zerotol
       Loop while:  ( ((i) and (ii)) or (iii))  and ((iv) and (v))
         i.e. while ( (t1&t2) | t3) & (t4 & t5),
      So stop:   (  (~t1|~t2) & ~t3) | (~t4 | ~t5)
 


Inverse Interpolation Starts

 if side>0,     last iterate was on the GOOD SIDE   


  if  sum(ptgood(:,4)) >=2,  two t-iterates with data good side  

         find interpolant using two points - then check for only
        one point type interpolant using point from bad side
      if (lambda<0)& (min(ptgood(1:2,2)) > (10^(side-1)*zerotol) ),  
       try to avoid too many on good side??????
        correct sign of multiplier
        lambda(t) is concave and nondecreasing in t
      ???what about >0?? interpol using lambda data to
       get to <=0 quickly?????
       INVERSE INTERPOLATION - interpolate on good side
       Let t(psi) be the inverse function of 
      psi(t), then we linearly interpolate
       the value t(0), which is t^*, using the two points 
      from the good side.
     tint=ptgood(1:2,1);
     yint=ptgood(1:2,2);
     f1=-1/yint(1)   + 1/ystar;
     f2=-1/yint(2)   + 1/ystar;
        Slope on the good side should always guarantee f1-f2 not =0.
     slope=(tint(1)-tint(2))/(f1-f2);
     tt=tint(1) - f1* slope;      linear inverse interpolation
     if newlow <=tt&tt<=newhigh,       roundoff error problem check  

          We now decide whether to chose the new point
        t=tt;
        tind=[tind 3];       indicator for t value - for debugging
     end
      elseif side>4,      We need a point on the bad side if too many ???
     if side>4&(min(ptgood(1:2,2)) < (10^(side-1)*zerotol) ),  

         too many iterates on good side and it is near bad case
        t=.95*newlow+.05*newhigh;
         ????use weights from k' if indexb>0?????
        tind=[tind 4];       indicator for t value - for debugging 
     end                 
  elseif sum(ptgood(:,4)) >=1 & indexb>0 >= 1,
     if ptgood(1,2)>zerotol,     not near hard case
         interpolate with one point from each side if not hard case
       tintg=ptgood(1,1);
       tintb=ptbad(1,1);
       psig=-1/ptgood(1,2)+sqrt(ss1);
       psib=-1/ptbad(1,2)+sqrt(ss1);
       tt=(tintg*psib-tintb*psig)/(psib-psig); 
       if newlow <= tt & tt <= newhigh,
          t=tt;
          tind=[tind 31];
       end
     else
        use spline on k(t) - to do still
          tind=[tind 32];
     end
  end     end of if sum(ptgood) >=2 else >=1




 Top of Documentation page.