Logo ROOT   6.13/01
Reference Guide
TQpSolverBase.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 // TSolverBase //
46 // //
47 // The Solver class contains methods for monitoring and checking the //
48 // convergence status of the algorithm, methods to determine the step //
49 // length along a given direction, methods to define the starting point,//
50 // and the solve method that implements the interior-point algorithm //
51 // //
52 //////////////////////////////////////////////////////////////////////////
53 
54 #include "TMath.h"
55 #include "TQpSolverBase.h"
56 
57 ClassImp(TQpSolverBase);
58 
59 ////////////////////////////////////////////////////////////////////////////////
60 /// Default constructor
61 
63 {
64  fSys = 0;
65 
66  fDnorm = 0.;
67 
68  // define parameters associated with the step length heuristic
69  fMutol = 1.0e-8;
70  fArtol = 1.0e-8;
71  fGamma_f = 0.99;
72  fGamma_a = 1.0/(1.0-fGamma_f);
73 
74  fPhi = 0.0;
75 
76  fMaxit = 100;
77 
78  // allocate space to track the sequence of complementarity gaps,
79  // residual norms, and merit functions.
80  fMu_history = new Double_t[fMaxit];
81  fRnorm_history = new Double_t[fMaxit];
82  fPhi_history = new Double_t[fMaxit];
83  fPhi_min_history = new Double_t[fMaxit];
84 
85  fIter = 0;
86 }
87 
88 
89 ////////////////////////////////////////////////////////////////////////////////
90 /// Copy constructor
91 
92 TQpSolverBase::TQpSolverBase(const TQpSolverBase &another) : TObject(another)
93 {
94  *this = another;
95 }
96 
97 
98 ////////////////////////////////////////////////////////////////////////////////
99 /// Deconstructor
100 
102 {
103  if (fSys) { delete fSys; fSys = 0; }
104 
105  if (fMu_history) { delete [] fMu_history; fMu_history = 0; }
106  if (fRnorm_history) { delete [] fRnorm_history; fRnorm_history = 0; }
107  if (fPhi_history) { delete [] fPhi_history; fPhi_history = 0; }
108  if (fPhi_min_history) { delete [] fPhi_min_history; fPhi_min_history = 0; }
109 }
110 
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 /// Implements a default starting-point heuristic. While interior-point theory
114 /// places fairly loose restrictions on the choice of starting point, the choice
115 /// of heuristic can significantly affect the robustness and efficiency of the
116 /// algorithm.
117 
119  TQpResidual *resid,TQpVar *step)
120 {
121  this->DefStart(formulation,iterate,prob,resid,step);
122 }
123 
124 
125 ////////////////////////////////////////////////////////////////////////////////
126 /// Default starting point
127 
129  TQpDataBase *prob,TQpResidual *resid,TQpVar *step)
130 {
131  Double_t sdatanorm = TMath::Sqrt(fDnorm);
132  Double_t a = sdatanorm;
133  Double_t b = sdatanorm;
134 
135  iterate->InteriorPoint(a,b);
136  resid->CalcResids(prob,iterate);
137  resid->Set_r3_xz_alpha(iterate,0.0);
138 
139  fSys->Factor(prob,iterate);
140  fSys->Solve(prob,iterate,resid,step);
141  step->Negate();
142 
143  // Take the full affine scaling step
144  iterate->Saxpy(step,1.0);
145  // resid.CalcResids(prob,iterate); // Calc the resids if debugging.
146  Double_t shift = 1.e3+2*iterate->Violation();
147  iterate->ShiftBoundVariables(shift,shift);
148 }
149 
150 
151 ////////////////////////////////////////////////////////////////////////////////
152 /// Starting point algoritm according to Stephen Wright
153 
154 void TQpSolverBase::SteveStart(TQpProbBase * /* formulation */,
155  TQpVar *iterate,TQpDataBase *prob,
156  TQpResidual *resid,TQpVar *step)
157 {
158  Double_t sdatanorm = TMath::Sqrt(fDnorm);
159  Double_t a = 0.0;
160  Double_t b = 0.0;
161 
162  iterate->InteriorPoint(a,b);
163 
164  // set the r3 component of the rhs to -(norm of data), and calculate
165  // the residuals that are obtained when all values are zero.
166 
167  resid->Set_r3_xz_alpha(iterate,-sdatanorm);
168  resid->CalcResids(prob,iterate);
169 
170  // next, assign 1 to all the complementary variables, so that there
171  // are identities in the coefficient matrix when we do the solve.
172 
173  a = 1.0; b = 1.0;
174  iterate->InteriorPoint(a,b);
175  fSys->Factor(prob,iterate);
176  fSys->Solve (prob,iterate,resid,step);
177  step->Negate();
178 
179  // copy the "step" into the current vector
180 
181  iterate = step;
182 
183  // calculate the maximum violation of the complementarity
184  // conditions, and shift these variables to restore positivity.
185  Double_t shift = 1.5*iterate->Violation();
186  iterate->ShiftBoundVariables(shift,shift);
187 
188  // do Mehrotra-type adjustment
189 
190  const Double_t mutemp = iterate->GetMu();
191  const Double_t xsnorm = iterate->Norm1();
192  const Double_t delta = 0.5*iterate->fNComplementaryVariables*mutemp/xsnorm;
193  iterate->ShiftBoundVariables(delta,delta);
194 }
195 
196 
197 ////////////////////////////////////////////////////////////////////////////////
198 /// Alternative starting point heuristic: sets the "complementary" variables to a large
199 /// positive value (based on the norm of the problem data) and the remaining variables
200 /// to zero .
201 
202 void TQpSolverBase::DumbStart(TQpProbBase * /* formulation */,
203  TQpVar *iterate,TQpDataBase * /* prob */,
204  TQpResidual * /* resid */,TQpVar * /* step */)
205 {
206  const Double_t sdatanorm = fDnorm;
207  const Double_t a = 1.e3;
208  const Double_t b = 1.e5;
209  const Double_t bigstart = a*sdatanorm+b;
210  iterate->InteriorPoint(bigstart,bigstart);
211 }
212 
213 
214 ////////////////////////////////////////////////////////////////////////////////
215 /// Implements a version of Mehrotra starting point heuristic,
216 /// modified to ensure identical steps in the primal and dual variables.
217 
219 {
220  Int_t firstOrSecond;
221  Double_t primalValue; Double_t primalStep; Double_t dualValue; Double_t dualStep;
222  const Double_t maxAlpha = iterate->FindBlocking(step,primalValue,primalStep,
223  dualValue,dualStep,firstOrSecond);
224  Double_t mufull = iterate->MuStep(step,maxAlpha);
225 
226  mufull /= fGamma_a;
227 
228  Double_t alpha = 1.0;
229  switch (firstOrSecond) {
230  case 0:
231  alpha = 1; // No constraints were blocking
232  break;
233  case 1:
234  alpha = (-primalValue+mufull/(dualValue+maxAlpha*dualStep))/primalStep;
235  break;
236  case 2:
237  alpha = (-dualValue+mufull/(primalValue+maxAlpha*primalStep))/dualStep;
238  break;
239  default:
240  R__ASSERT(0 && "Can't get here");
241  break;
242  }
243  // make it at least fGamma_f * maxStep
244  if (alpha < fGamma_f*maxAlpha) alpha = fGamma_f*maxAlpha;
245  // back off just a touch
246  alpha *= .99999999;
247 
248  return alpha;
249 }
250 
251 
252 ////////////////////////////////////////////////////////////////////////////////
253 /// Monitor progress / convergence aat each interior-point iteration
254 
256  Double_t alpha,Double_t sigma,Int_t i,Double_t mu,
257  Int_t stop_code,Int_t level)
258 {
259  this->DefMonitor(data,vars,resids,alpha,sigma,i,mu,stop_code,level);
260 }
261 
262 
263 ////////////////////////////////////////////////////////////////////////////////
264 /// Tests for termination. Unless the user supplies a specific termination
265 /// routine, this method calls another method defaultStatus, which returns
266 /// a code indicating the current convergence status.
267 
269  Int_t i,Double_t mu,Int_t level)
270 {
271  return this->DefStatus(data,vars,resids,i,mu,level);
272 }
273 
274 
275 ////////////////////////////////////////////////////////////////////////////////
276 /// Default status method
277 
278 Int_t TQpSolverBase::DefStatus(TQpDataBase * /* data */,TQpVar * /* vars */,
279  TQpResidual *resids,Int_t iterate,Double_t mu,Int_t /* level */)
280 {
281  Int_t stop_code = kNOT_FINISHED;
282 
283  const Double_t gap = TMath::Abs(resids->GetDualityGap());
284  const Double_t rnorm = resids->GetResidualNorm();
285 
286  Int_t idx = iterate-1;
287  if (idx < 0 ) idx = 0;
288  if (idx >= fMaxit) idx = fMaxit-1;
289 
290  // store the historical record
291  fMu_history[idx] = mu;
292  fRnorm_history[idx] = rnorm;
293  fPhi = (rnorm+gap)/fDnorm;
294  fPhi_history[idx] = fPhi;
295 
296  if (idx > 0) {
297  fPhi_min_history[idx] = fPhi_min_history[idx-1];
298  if (fPhi < fPhi_min_history[idx]) fPhi_min_history[idx] = fPhi;
299  } else
300  fPhi_min_history[idx] = fPhi;
301 
302  if (iterate >= fMaxit) {
303  stop_code = kMAX_ITS_EXCEEDED;
304  }
305  else if (mu <= fMutol && rnorm <= fArtol*fDnorm) {
306  stop_code = kSUCCESSFUL_TERMINATION;
307  }
308  if (stop_code != kNOT_FINISHED) return stop_code;
309 
310  // check infeasibility condition
311  if (idx >= 10 && fPhi >= 1.e-8 && fPhi >= 1.e4*fPhi_min_history[idx])
312  stop_code = kINFEASIBLE;
313  if (stop_code != kNOT_FINISHED) return stop_code;
314 
315  // check for unknown status: slow convergence first
316  if (idx >= 30 && fPhi_min_history[idx] >= .5*fPhi_min_history[idx-30])
317  stop_code = kUNKNOWN;
318 
319  if (rnorm/fDnorm > fArtol &&
320  (fRnorm_history[idx]/fMu_history[idx])/(fRnorm_history[0]/fMu_history[0]) >= 1.e8)
321  stop_code = kUNKNOWN;
322 
323  return stop_code;
324 }
325 
326 
327 ////////////////////////////////////////////////////////////////////////////////
328 /// Assignment operator
329 
331 {
332  if (this != &source) {
333  TObject::operator=(source);
334 
335  fSys = source.fSys;
336  fDnorm = source.fDnorm;
337  fMutol = source.fMutol;
338  fArtol = source.fArtol;
339  fGamma_f = source.fGamma_f;
340  fGamma_a = source.fGamma_a;
341  fPhi = source.fPhi;
342  fIter = source.fIter;
343 
344  if (fMaxit != source.fMaxit) {
345  if (fMu_history) delete [] fMu_history;
346  fMu_history = new Double_t[fMaxit];
347  if (fRnorm_history) delete [] fRnorm_history;
348  fRnorm_history = new Double_t[fMaxit];
349  if (fPhi_history) delete [] fPhi_history;
350  fPhi_history = new Double_t[fMaxit];
351  if (fPhi_min_history) delete [] fPhi_min_history;
352  fPhi_min_history = new Double_t[fMaxit];
353  }
354 
355  fMaxit = source.fMaxit;
356  memcpy(fMu_history,source.fMu_history,fMaxit*sizeof(Double_t));
357  memcpy(fRnorm_history,source.fRnorm_history,fMaxit*sizeof(Double_t));
358  memcpy(fPhi_history,source.fPhi_history,fMaxit*sizeof(Double_t));
359  memcpy(fPhi_min_history,source.fPhi_min_history,fMaxit*sizeof(Double_t));
360  }
361  return *this;
362 }
Double_t * fRnorm_history
Definition: TQpSolverBase.h:90
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
TQpSolverBase()
Default constructor.
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...
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
virtual void SteveStart(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Starting point algoritm according to Stephen Wright.
virtual void DefStart(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Default starting point.
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 Double_t Violation()
The amount by which the current variables violate the non-negativity constraints. ...
Definition: TQpVar.cxx:573
virtual Double_t GetMu()
compute complementarity gap, obtained by taking the inner product of the complementary vectors and di...
Definition: TQpVar.cxx:191
virtual ~TQpSolverBase()
Deconstructor.
Double_t fGamma_a
Definition: TQpSolverBase.h:84
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
virtual void DumbStart(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Alternative starting point heuristic: sets the "complementary" variables to a large positive value (b...
Double_t fGamma_f
Definition: TQpSolverBase.h:83
Double_t fDnorm
Definition: TQpSolverBase.h:78
virtual void ShiftBoundVariables(Double_t alpha, Double_t beta)
Add alpha to components of (u,t,v,w) and beta to components of (lambda,pi,phi,gamma) ...
Definition: TQpVar.cxx:614
virtual Int_t DefStatus(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Int_t i, Double_t mu, Int_t level)
Default status method.
Double_t * fMu_history
Definition: TQpSolverBase.h:89
Double_t GetDualityGap()
Definition: TQpResidual.h:109
TQpLinSolverBase * fSys
Definition: TQpSolverBase.h:76
virtual void Negate()
Perform a "negate" operation on all data vectors : x = -x.
Definition: TQpVar.cxx:271
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 DefMonitor(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)=0
static Double_t FindBlocking(TVectorD &w, TVectorD &wstep, TVectorD &u, TVectorD &ustep, Double_t maxStep, Double_t &w_elt, Double_t &wstep_elt, Double_t &u_elt, Double_t &ustep_elt, int &first_or_second)
See other FindBlocking function.
Definition: TQpVar.cxx:445
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.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t fArtol
Definition: TQpSolverBase.h:81
virtual void InteriorPoint(Double_t alpha, Double_t beta)
Sets components of (u,t,v,w) to alpha and of (lambda,pi,phi,gamma) to beta.
Definition: TQpVar.cxx:533
virtual Double_t Norm1()
Return the sum of the vector-norm1&#39;s.
Definition: TQpVar.cxx:675
Definition: TQpVar.h:59
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Int_t fNComplementaryVariables
Definition: TQpVar.h:88
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
Double_t fPhi
Definition: TQpSolverBase.h:85
Double_t fMutol
Definition: TQpSolverBase.h:80
Double_t * fPhi_history
Definition: TQpSolverBase.h:91
Double_t * fPhi_min_history
Definition: TQpSolverBase.h:93