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
57
58////////////////////////////////////////////////////////////////////////////////
59/// Default constructor
60
62{
63 fSys = nullptr;
64
65 fDnorm = 0.;
66
67 // define parameters associated with the step length heuristic
68 fMutol = 1.0e-8;
69 fArtol = 1.0e-8;
70 fGamma_f = 0.99;
71 fGamma_a = 1.0/(1.0-fGamma_f);
72
73 fPhi = 0.0;
74
75 fMaxit = 100;
76
77 // allocate space to track the sequence of complementarity gaps,
78 // residual norms, and merit functions.
83
84 fIter = 0;
85}
86
87
88////////////////////////////////////////////////////////////////////////////////
89/// Copy constructor
90
95
96
97////////////////////////////////////////////////////////////////////////////////
98/// Deconstructor
99
101{
102 if (fSys) { delete fSys; fSys = nullptr; }
103
104 if (fMu_history) { delete [] fMu_history; fMu_history = nullptr; }
105 if (fRnorm_history) { delete [] fRnorm_history; fRnorm_history = nullptr; }
106 if (fPhi_history) { delete [] fPhi_history; fPhi_history = nullptr; }
107 if (fPhi_min_history) { delete [] fPhi_min_history; fPhi_min_history = nullptr; }
108}
109
110
111////////////////////////////////////////////////////////////////////////////////
112/// Implements a default starting-point heuristic. While interior-point theory
113/// places fairly loose restrictions on the choice of starting point, the choice
114/// of heuristic can significantly affect the robustness and efficiency of the
115/// algorithm.
116
122
123
124////////////////////////////////////////////////////////////////////////////////
125/// Default starting point
126
129{
133
134 iterate->InteriorPoint(a,b);
135 resid->CalcResids(prob,iterate);
136 resid->Set_r3_xz_alpha(iterate,0.0);
137
139 fSys->Solve(prob,iterate,resid,step);
140 step->Negate();
141
142 // Take the full affine scaling step
143 iterate->Saxpy(step,1.0);
144 // resid.CalcResids(prob,iterate); // Calc the resids if debugging.
145 Double_t shift = 1.e3+2*iterate->Violation();
146 iterate->ShiftBoundVariables(shift,shift);
147}
148
149
150////////////////////////////////////////////////////////////////////////////////
151/// Starting point algorithm according to Stephen Wright
152
153void TQpSolverBase::SteveStart(TQpProbBase * /* formulation */,
155 TQpResidual *resid,TQpVar *step)
156{
158 Double_t a = 0.0;
159 Double_t b = 0.0;
160
161 iterate->InteriorPoint(a,b);
162
163 // set the r3 component of the rhs to -(norm of data), and calculate
164 // the residuals that are obtained when all values are zero.
165
166 resid->Set_r3_xz_alpha(iterate,-sdatanorm);
167 resid->CalcResids(prob,iterate);
168
169 // next, assign 1 to all the complementary variables, so that there
170 // are identities in the coefficient matrix when we do the solve.
171
172 a = 1.0; b = 1.0;
173 iterate->InteriorPoint(a,b);
175 fSys->Solve (prob,iterate,resid,step);
176 step->Negate();
177
178 // copy the "step" into the current vector
179
180 iterate = step;
181
182 // calculate the maximum violation of the complementarity
183 // conditions, and shift these variables to restore positivity.
184 Double_t shift = 1.5*iterate->Violation();
185 iterate->ShiftBoundVariables(shift,shift);
186
187 // do Mehrotra-type adjustment
188
189 const Double_t mutemp = iterate->GetMu();
190 const Double_t xsnorm = iterate->Norm1();
191 const Double_t delta = 0.5*iterate->fNComplementaryVariables*mutemp/xsnorm;
192 iterate->ShiftBoundVariables(delta,delta);
193}
194
195
196////////////////////////////////////////////////////////////////////////////////
197/// Alternative starting point heuristic: sets the "complementary" variables to a large
198/// positive value (based on the norm of the problem data) and the remaining variables
199/// to zero .
200
201void TQpSolverBase::DumbStart(TQpProbBase * /* formulation */,
202 TQpVar *iterate,TQpDataBase * /* prob */,
203 TQpResidual * /* resid */,TQpVar * /* step */)
204{
205 const Double_t sdatanorm = fDnorm;
206 const Double_t a = 1.e3;
207 const Double_t b = 1.e5;
208 const Double_t bigstart = a*sdatanorm+b;
209 iterate->InteriorPoint(bigstart,bigstart);
210}
211
212
213////////////////////////////////////////////////////////////////////////////////
214/// Implements a version of Mehrotra starting point heuristic,
215/// modified to ensure identical steps in the primal and dual variables.
216
218{
221 const Double_t maxAlpha = iterate->FindBlocking(step,primalValue,primalStep,
223 Double_t mufull = iterate->MuStep(step,maxAlpha);
224
225 mufull /= fGamma_a;
226
227 Double_t alpha = 1.0;
228 switch (firstOrSecond) {
229 case 0:
230 alpha = 1; // No constraints were blocking
231 break;
232 case 1:
234 break;
235 case 2:
237 break;
238 default:
239 R__ASSERT(0 && "Can't get here");
240 break;
241 }
242 // make it at least fGamma_f * maxStep
243 if (alpha < fGamma_f*maxAlpha) alpha = fGamma_f*maxAlpha;
244 // back off just a touch
245 alpha *= .99999999;
246
247 return alpha;
248}
249
250
251////////////////////////////////////////////////////////////////////////////////
252/// Monitor progress / convergence aat each interior-point iteration
253
256 Int_t stop_code,Int_t level)
257{
258 this->DefMonitor(data,vars,resids,alpha,sigma,i,mu,stop_code,level);
259}
260
261
262////////////////////////////////////////////////////////////////////////////////
263/// Tests for termination. Unless the user supplies a specific termination
264/// routine, this method calls another method defaultStatus, which returns
265/// a code indicating the current convergence status.
266
268 Int_t i,Double_t mu,Int_t level)
269{
270 return this->DefStatus(data,vars,resids,i,mu,level);
271}
272
273
274////////////////////////////////////////////////////////////////////////////////
275/// Default status method
276
279{
281
282 const Double_t gap = TMath::Abs(resids->GetDualityGap());
283 const Double_t rnorm = resids->GetResidualNorm();
284
285 Int_t idx = iterate-1;
286 if (idx < 0 ) idx = 0;
287 if (idx >= fMaxit) idx = fMaxit-1;
288
289 // store the historical record
290 fMu_history[idx] = mu;
291 fRnorm_history[idx] = rnorm;
292 fPhi = (rnorm+gap)/fDnorm;
293 fPhi_history[idx] = fPhi;
294
295 if (idx > 0) {
297 if (fPhi < fPhi_min_history[idx]) fPhi_min_history[idx] = fPhi;
298 } else
299 fPhi_min_history[idx] = fPhi;
300
301 if (iterate >= fMaxit) {
303 }
304 else if (mu <= fMutol && rnorm <= fArtol*fDnorm) {
306 }
307 if (stop_code != kNOT_FINISHED) return stop_code;
308
309 // check infeasibility condition
310 if (idx >= 10 && fPhi >= 1.e-8 && fPhi >= 1.e4*fPhi_min_history[idx])
312 if (stop_code != kNOT_FINISHED) return stop_code;
313
314 // check for unknown status: slow convergence first
315 if (idx >= 30 && fPhi_min_history[idx] >= .5*fPhi_min_history[idx-30])
317
318 if (rnorm/fDnorm > fArtol &&
319 (fRnorm_history[idx]/fMu_history[idx])/(fRnorm_history[0]/fMu_history[0]) >= 1.e8)
321
322 return stop_code;
323}
324
325
326////////////////////////////////////////////////////////////////////////////////
327/// Assignment operator
328
330{
331 if (this != &source) {
333
334 fSys = source.fSys;
335 fDnorm = source.fDnorm;
336 fMutol = source.fMutol;
337 fArtol = source.fArtol;
338 fGamma_f = source.fGamma_f;
339 fGamma_a = source.fGamma_a;
340 fPhi = source.fPhi;
341 fIter = source.fIter;
342
343 if (fMaxit != source.fMaxit) {
344 if (fMu_history) delete [] fMu_history;
346 if (fRnorm_history) delete [] fRnorm_history;
348 if (fPhi_history) delete [] fPhi_history;
350 if (fPhi_min_history) delete [] fPhi_min_history;
352 }
353
354 fMaxit = source.fMaxit;
355 memcpy(fMu_history,source.fMu_history,fMaxit*sizeof(Double_t));
356 memcpy(fRnorm_history,source.fRnorm_history,fMaxit*sizeof(Double_t));
357 memcpy(fPhi_history,source.fPhi_history,fMaxit*sizeof(Double_t));
358 memcpy(fPhi_min_history,source.fPhi_min_history,fMaxit*sizeof(Double_t));
359 }
360 return *this;
361}
#define b(i)
Definition RSha256.hxx:100
#define a(i)
Definition RSha256.hxx:99
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
@ kNOT_FINISHED
@ kINFEASIBLE
@ kMAX_ITS_EXCEEDED
@ kUNKNOWN
@ kSUCCESSFUL_TERMINATION
Mother of all ROOT objects.
Definition TObject.h:41
TObject & operator=(const TObject &rhs) noexcept
TObject assignment operator.
Definition TObject.h:299
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
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
~TQpSolverBase() override
Deconstructor.
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
TQpSolverBase & operator=(const TQpSolverBase &source)
Assignment operator.
virtual void SteveStart(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Starting point algorithm 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:271
const Double_t sigma
int iterate(rng_state_t *X)
Definition mixmax.icc:34
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124