Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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/// \class TQpSolverBase
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
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.
84
85 fIter = 0;
86}
87
88
89////////////////////////////////////////////////////////////////////////////////
90/// Copy constructor
91
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
154void TQpSolverBase::SteveStart(TQpProbBase * /* formulation */,
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
202void 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
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
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) {
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;
347 if (fRnorm_history) delete [] fRnorm_history;
349 if (fPhi_history) delete [] fPhi_history;
351 if (fPhi_min_history) delete [] fPhi_min_history;
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}
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
double Double_t
Definition RtypesCore.h:59
#define ClassImp(name)
Definition Rtypes.h:364
#define R__ASSERT(e)
Definition TError.h:120
@ kNOT_FINISHED
@ kINFEASIBLE
@ kMAX_ITS_EXCEEDED
@ kUNKNOWN
@ kSUCCESSFUL_TERMINATION
Mother of all ROOT objects.
Definition TObject.h:37
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition TObject.h:283
Data for the general QP formulation.
Definition TQpDataBase.h:61
virtual void Solve(TQpDataBase *prob, TQpVar *vars, TQpResidual *resids, TQpVar *step)
Solves the system for a given set of residuals.
virtual void Factor(TQpDataBase *prob, TQpVar *vars)
Sets up the matrix for the main linear system in "augmented system" form.
default general problem formulation:
Definition TQpProbBase.h:89
The Residuals class calculates and stores the quantities that appear on the right-hand side of the li...
Definition TQpResidual.h:62
Double_t GetDualityGap()
Double_t GetResidualNorm()
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...
void CalcResids(TQpDataBase *problem, TQpVar *vars)
Calculate residuals, their norms, and duality complementarity gap, given a problem and variable set.
The Solver class contains methods for monitoring and checking the convergence status of the algorithm...
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 * fPhi_min_history
TQpLinSolverBase * fSys
Double_t fDnorm
Double_t fGamma_a
virtual void DefStart(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Default starting point.
virtual void Start(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Implements a default starting-point heuristic.
Double_t * fPhi_history
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.
Double_t * fMu_history
TQpSolverBase()
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 stop_code, Int_t level)=0
Double_t fGamma_f
virtual Int_t DoStatus(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Int_t i, Double_t mu, Int_t level)
Tests for termination.
virtual Int_t DefStatus(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Int_t i, Double_t mu, Int_t level)
Default status method.
virtual Double_t FinalStepLength(TQpVar *iterate, TQpVar *step)
Implements a version of Mehrotra starting point heuristic, modified to ensure identical steps in the ...
Double_t * fRnorm_history
Double_t fMutol
virtual ~TQpSolverBase()
Deconstructor.
TQpSolverBase & operator=(const TQpSolverBase &source)
Assignment operator.
virtual void SteveStart(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Starting point algoritm according to Stephen Wright.
Double_t fArtol
Class containing the variables for the general QP formulation.
Definition TQpVar.h:60
virtual void Negate()
Perform a "negate" operation on all data vectors : x = -x.
Definition TQpVar.cxx:272
const Double_t sigma
int iterate(rng_state_t *X)
Definition mixmax.icc:34
Double_t Sqrt(Double_t x)
Definition TMath.h:691
Short_t Abs(Short_t d)
Definition TMathBase.h:120