Logo ROOT   6.13/01
Reference Guide
TGondzioSolver.cxx
Go to the documentation of this file.
1 // @(#)root/quadp:$Id$
2 // Author: Eddy Offermann May 2004
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /*************************************************************************
13  * Parts of this file are copied from the OOQP distribution and *
14  * are subject to the following license: *
15  * *
16  * COPYRIGHT 2001 UNIVERSITY OF CHICAGO *
17  * *
18  * The copyright holder hereby grants you royalty-free rights to use, *
19  * reproduce, prepare derivative works, and to redistribute this software*
20  * to others, provided that any changes are clearly documented. This *
21  * software was authored by: *
22  * *
23  * E. MICHAEL GERTZ gertz@mcs.anl.gov *
24  * Mathematics and Computer Science Division *
25  * Argonne National Laboratory *
26  * 9700 S. Cass Avenue *
27  * Argonne, IL 60439-4844 *
28  * *
29  * STEPHEN J. WRIGHT swright@cs.wisc.edu *
30  * Computer Sciences Department *
31  * University of Wisconsin *
32  * 1210 West Dayton Street *
33  * Madison, WI 53706 FAX: (608)262-9777 *
34  * *
35  * Any questions or comments may be directed to one of the authors. *
36  * *
37  * ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES OF *
38  * ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT, AND *
39  * OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A CONTRACT *
40  * WITH THE DEPARTMENT OF ENERGY. *
41  *************************************************************************/
42 
43 //////////////////////////////////////////////////////////////////////////
44 // //
45 // TGondzioSolver //
46 // //
47 // Derived class of TQpSolverBase implementing Gondzio-correction //
48 // version of Mehrotra's original predictor-corrector algorithm. //
49 // //
50 //////////////////////////////////////////////////////////////////////////
51 
52 
53 #include "Riostream.h"
54 #include "TMath.h"
55 #include "TGondzioSolver.h"
56 #include "TQpLinSolverDens.h"
57 
58 ClassImp(TGondzioSolver);
59 
60 ////////////////////////////////////////////////////////////////////////////////
61 /// Default constructor
62 
64 {
65  fPrintlevel = 0;
66  fTsig = 0.0;
69 
70  fStepFactor0 = 0.0;
71  fStepFactor1 = 0.0;
72  fAcceptTol = 0.0;
73  fBeta_min = 0.0;
74  fBeta_max = 0.0;
75 
76  fCorrector_step = 0;
77  fStep = 0;
78  fCorrector_resid = 0;
79  fFactory = 0;
80 }
81 
82 
83 ////////////////////////////////////////////////////////////////////////////////
84 /// Constructor
85 
87 {
88  fFactory = of;
89  fStep = fFactory->MakeVariables(prob);
92 
93  fPrintlevel = verbose;
94  fTsig = 3.0; // the usual value for the centering exponent (tau)
95 
96  fMaximum_correctors = 3; // maximum number of Gondzio correctors
97 
99 
100  // the two StepFactor constants set targets for increase in step
101  // length for each corrector
102  fStepFactor0 = 0.08;
103  fStepFactor1 = 1.08;
104 
105  // accept the enhanced step if it produces a small improvement in
106  // the step length
107  fAcceptTol = 0.005;
108 
109  //define the Gondzio correction box
110  fBeta_min = 0.1;
111  fBeta_max = 10.0;
112 }
113 
114 
115 ////////////////////////////////////////////////////////////////////////////////
116 /// Copy constructor
117 
119 {
120  *this = another;
121 }
122 
123 
124 ////////////////////////////////////////////////////////////////////////////////
125 /// Solve the quadratic programming problem as formulated through prob, store
126 /// the final solution in iterate->fX . Monitor the residuals during the iterations
127 /// through resid . The status is returned as defined in TQpSolverBase::ETerminationCode .
128 
130 {
131  Int_t status_code;
132  Double_t alpha = 1;
133  Double_t sigma = 1;
134 
135  fDnorm = prob->DataNorm();
136 
137  // initialization of (x,y,z) and factorization routine.
138  fSys = fFactory->MakeLinSys(prob);
139  this->Start(fFactory,iterate,prob,resid,fStep);
140 
141  fIter = 0;
143  Double_t mu = iterate->GetMu();
144 
145  Int_t done = 0;
146  do {
147  fIter++;
148  // evaluate residuals and update algorithm status:
149  resid->CalcResids(prob,iterate);
150 
151  // termination test:
152  status_code = this->DoStatus(prob,iterate,resid,fIter,mu,0);
153  if (status_code != kNOT_FINISHED ) break;
154  if (fPrintlevel >= 10)
155  this->DoMonitor(prob,iterate,resid,alpha,sigma,fIter,mu,status_code,0);
156 
157  // *** Predictor step ***
158 
159  resid->Set_r3_xz_alpha(iterate,0.0);
160 
161  fSys->Factor(prob,iterate);
162  fSys->Solve(prob,iterate,resid,fStep);
163  fStep->Negate();
164 
165  alpha = iterate->StepBound(fStep);
166 
167  // calculate centering parameter
168  Double_t muaff = iterate->MuStep(fStep,alpha);
169  sigma = TMath::Power(muaff/mu,fTsig);
170 
171  if (fPrintlevel >= 10)
172  this->DoMonitor(prob,iterate,resid,alpha,sigma,fIter,mu,status_code,2);
173 
174  // *** Corrector step ***
175 
176  // form right hand side of linear system:
177  resid->Add_r3_xz_alpha(fStep,-sigma*mu);
178 
179  fSys->Solve(prob,iterate,resid,fStep);
180  fStep->Negate();
181 
182  // calculate distance to boundary along the Mehrotra
183  // predictor-corrector step:
184  alpha = iterate->StepBound(fStep);
185 
186  // prepare for Gondzio corrector loop: zero out the
187  // corrector_resid structure:
189 
190  // calculate the target box:
191  Double_t rmin = sigma*mu*fBeta_min;
192  Double_t rmax = sigma*mu*fBeta_max;
193 
194  Int_t stopCorrections = 0;
196 
197  // enter the Gondzio correction loop:
198  if (fPrintlevel >= 10)
199  std::cout << "**** Entering the correction loop ****" << std::endl;
200 
202  alpha < 1.0 && !stopCorrections) {
203 
204  // copy current variables into fcorrector_step
206 
207  // calculate target steplength
208  Double_t alpha_target = fStepFactor1*alpha+fStepFactor0;
209  if (alpha_target > 1.0) alpha_target = 1.0;
210 
211  // add a step of this length to corrector_step
212  fCorrector_step->Saxpy(fStep,alpha_target);
213 
214  // place XZ into the r3 component of corrector_resids
216 
217  // do the projection operation
218  fCorrector_resid->Project_r3(rmin,rmax);
219 
220  // solve for corrector direction
222 
223  // add the current step to corrector_step, and calculate the
224  // step to boundary along the resulting direction
226  Double_t alpha_enhanced = iterate->StepBound(fCorrector_step);
227 
228  // if the enhanced step length is actually 1, make it official
229  // and stop correcting
230  if (alpha_enhanced == 1.0) {
232  alpha = alpha_enhanced;
234  stopCorrections = 1;
235  }
236  else if(alpha_enhanced >= (1.0+fAcceptTol)*alpha) {
237  // if enhanced step length is significantly better than the
238  // current alpha, make the enhanced step official, but maybe
239  // keep correcting
241  alpha = alpha_enhanced;
243  stopCorrections = 0;
244  }
245  else {
246  // otherwise quit the correction loop
247  stopCorrections = 1;
248  }
249  }
250 
251  // We've finally decided on a step direction, now calculate the
252  // length using Mehrotra's heuristic.x
253  alpha = this->FinalStepLength(iterate,fStep);
254 
255  // alternatively, just use a crude step scaling factor.
256  // alpha = 0.995 * iterate->StepBound(fStep);
257 
258  // actually take the step (at last!) and calculate the new mu
259 
260  iterate->Saxpy(fStep,alpha);
261  mu = iterate->GetMu();
262  } while (!done);
263 
264  resid->CalcResids(prob,iterate);
265  if (fPrintlevel >= 10)
266  this->DoMonitor(prob,iterate,resid,alpha,sigma,fIter,mu,status_code,1);
267 
268  return status_code;
269 }
270 
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Print information about the optimization process and monitor the convergence
274 /// status of thye algorithm
275 
276 void TGondzioSolver::DefMonitor(TQpDataBase* /* data */,TQpVar* /* vars */,
277  TQpResidual *resid,
278  Double_t alpha,Double_t sigma,Int_t i,Double_t mu,
279  Int_t status_code,Int_t level)
280 {
281  switch (level) {
282  case 0 : case 1:
283  {
284  std::cout << std::endl << "Duality Gap: " << resid->GetDualityGap() << std::endl;
285  if (i > 1) {
286  std::cout << " Number of Corrections = " << fNumberGondzioCorrections
287  << " alpha = " << alpha << std::endl;
288  }
289  std::cout << " *** Iteration " << i << " *** " << std::endl;
290  std::cout << " mu = " << mu << " relative residual norm = "
291  << resid->GetResidualNorm()/fDnorm << std::endl;
292 
293  if (level == 1) {
294  // Termination has been detected by the status check; print
295  // appropriate message
296  if (status_code == kSUCCESSFUL_TERMINATION) {
297  std::cout << std::endl
298  << " *** SUCCESSFUL TERMINATION ***"
299  << std::endl;
300  }
301  else if (status_code == kMAX_ITS_EXCEEDED) {
302  std::cout << std::endl
303  << " *** MAXIMUM ITERATIONS REACHED *** " << std::endl;
304  }
305  else if (status_code == kINFEASIBLE) {
306  std::cout << std::endl
307  << " *** TERMINATION: PROBABLY INFEASIBLE *** "
308  << std::endl;
309  }
310  else if (status_code == kUNKNOWN) {
311  std::cout << std::endl
312  << " *** TERMINATION: STATUS UNKNOWN *** " << std::endl;
313  }
314  }
315  } break;
316  case 2:
317  std::cout << " *** sigma = " << sigma << std::endl;
318  break;
319  }
320 }
321 
322 
323 ////////////////////////////////////////////////////////////////////////////////
324 /// Deconstructor
325 
327 {
328  if (fCorrector_step) { delete fCorrector_step; fCorrector_step = 0; }
329  if (fStep) { delete fStep; fStep = 0; }
331 }
332 
333 
334 ////////////////////////////////////////////////////////////////////////////////
335 /// Assignment operator
336 
338 {
339  if (this != &source) {
340  TQpSolverBase::operator=(source);
341 
342  fPrintlevel = source.fPrintlevel;
343  fTsig = source.fTsig ;
346 
347  fStepFactor0 = source.fStepFactor0;
348  fStepFactor1 = source.fStepFactor1;
349  fAcceptTol = source.fAcceptTol;
350  fBeta_min = source.fBeta_min;
351  fBeta_max = source.fBeta_max;
352 
353  if (fCorrector_step) delete fCorrector_step;
354  if (fStep) delete fStep;
356 
357  fCorrector_step = new TQpVar(*source.fCorrector_step);
358  fStep = new TQpVar(*source.fStep);
360  fFactory = source.fFactory;
361  }
362  return *this;
363 }
virtual Double_t DataNorm()=0
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
virtual TQpVar * MakeVariables(const TQpDataBase *data)=0
virtual void DoMonitor(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Double_t alpha, Double_t sigma, Int_t i, Double_t mu, Int_t stop_code, Int_t level)
Monitor progress / convergence aat each interior-point iteration.
void CalcResids(TQpDataBase *problem, TQpVar *vars)
Calculate residuals, their norms, and duality complementarity gap, given a problem and variable set...
TGondzioSolver()
Default constructor.
virtual void DefMonitor(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Double_t alpha, Double_t sigma, Int_t i, Double_t mu, Int_t status_code, Int_t level)
Print information about the optimization process and monitor the convergence status of thye algorithm...
virtual Double_t MuStep(TQpVar *step, Double_t alpha)
Compute the complementarity gap resulting from a step of length "alpha" along direction "step"...
Definition: TQpVar.cxx:210
void Project_r3(Double_t rmin, Double_t rmax)
Perform the projection operation required by Gondzio algorithm: replace each component r3_i of the co...
Double_t fAcceptTol
Int_t fNumberGondzioCorrections
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:734
virtual void Solve(TQpDataBase *prob, TQpVar *vars, TQpResidual *resids, TQpVar *step)
Solves the system for a given set of residuals.
virtual Double_t FinalStepLength(TQpVar *iterate, TQpVar *step)
Implements a version of Mehrotra starting point heuristic, modified to ensure identical steps in the ...
virtual Int_t DoStatus(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Int_t i, Double_t mu, Int_t level)
Tests for termination.
TQpSolverBase & operator=(const TQpSolverBase &source)
Assignment operator.
virtual TQpResidual * MakeResiduals(const TQpDataBase *data)=0
virtual Double_t GetMu()
compute complementarity gap, obtained by taking the inner product of the complementary vectors and di...
Definition: TQpVar.cxx:191
Double_t fBeta_max
static Double_t StepBound(TVectorD &v, TVectorD &dir, Double_t maxStep)
Find the maximum stepsize of v in direction dir before violating the nonnegativity constraints...
Definition: TQpVar.cxx:347
Double_t fStepFactor0
virtual void Start(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Implements a default starting-point heuristic.
Double_t GetResidualNorm()
Definition: TQpResidual.h:108
void Add_r3_xz_alpha(TQpVar *vars, Double_t alpha)
Modify the "complementarity" component of the residuals, by adding the pairwise products of the compl...
Double_t fDnorm
Definition: TQpSolverBase.h:78
virtual Int_t Solve(TQpDataBase *prob, TQpVar *iterate, TQpResidual *resid)
Solve the quadratic programming problem as formulated through prob, store the final solution in itera...
TQpResidual * fCorrector_resid
Double_t GetDualityGap()
Definition: TQpResidual.h:109
TQpLinSolverBase * fSys
Definition: TQpSolverBase.h:76
virtual TQpLinSolverBase * MakeLinSys(const TQpDataBase *data)=0
virtual void Negate()
Perform a "negate" operation on all data vectors : x = -x.
Definition: TQpVar.cxx:271
Int_t fMaximum_correctors
void Set_r3_xz_alpha(TQpVar *vars, Double_t alpha)
Set the "complementarity" component of the residuals to the pairwise products of the complementary va...
virtual void Saxpy(TQpVar *b, Double_t alpha)
Perform a "saxpy" operation on all data vectors : x += alpha*y.
Definition: TQpVar.cxx:231
virtual void Factor(TQpDataBase *prob, TQpVar *vars)
Sets up the matrix for the main linear system in "augmented system" form.
TGondzioSolver & operator=(const TGondzioSolver &source)
Assignment operator.
void Clear_r1r2()
set the noncomplementarity components of the residual (the terms arising from the linear equalities i...
Double_t fBeta_min
Definition: TQpVar.h:59
Double_t fStepFactor1
virtual ~TGondzioSolver()
Deconstructor.
TQpProbBase * fFactory
TQpVar * fCorrector_step