Logo ROOT  
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/// \class 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 <iostream>
54#include "TMath.h"
55#include "TGondzioSolver.h"
56#include "TQpLinSolverDens.h"
57
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
77 fStep = 0;
79 fFactory = 0;
80}
81
82
83////////////////////////////////////////////////////////////////////////////////
84/// Constructor
85
87{
88 fFactory = of;
92
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
276void TGondzioSolver::DefMonitor(TQpDataBase* /* data */,TQpVar* /* vars */,
277 TQpResidual *resid,
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{
329 if (fStep) { delete fStep; fStep = 0; }
331}
332
333
334////////////////////////////////////////////////////////////////////////////////
335/// Assignment operator
336
338{
339 if (this != &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
354 if (fStep) delete fStep;
356
358 fStep = new TQpVar(*source.fStep);
360 fFactory = source.fFactory;
361 }
362 return *this;
363}
ClassImp(TGondzioSolver)
@ kNOT_FINISHED
Definition: TQpSolverBase.h:66
@ kINFEASIBLE
Definition: TQpSolverBase.h:68
@ kMAX_ITS_EXCEEDED
Definition: TQpSolverBase.h:67
@ kUNKNOWN
Definition: TQpSolverBase.h:69
@ kSUCCESSFUL_TERMINATION
Definition: TQpSolverBase.h:65
Derived class of TQpSolverBase implementing Gondzio-correction version of Mehrotra's original predict...
Double_t fAcceptTol
Double_t fBeta_min
Int_t fMaximum_correctors
TQpVar * fCorrector_step
Int_t fNumberGondzioCorrections
Double_t fStepFactor0
TQpProbBase * fFactory
TGondzioSolver()
Default constructor.
Int_t Solve(TQpDataBase *prob, TQpVar *iterate, TQpResidual *resid) override
Solve the quadratic programming problem as formulated through prob, store the final solution in itera...
Double_t fBeta_max
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) override
Print information about the optimization process and monitor the convergence status of thye algorithm...
~TGondzioSolver() override
Deconstructor.
TGondzioSolver & operator=(const TGondzioSolver &source)
Assignment operator.
Double_t fStepFactor1
TQpResidual * fCorrector_resid
Data for the general QP formulation.
Definition: TQpDataBase.h:61
virtual Double_t DataNorm()=0
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
virtual TQpLinSolverBase * MakeLinSys(const TQpDataBase *data)=0
virtual TQpResidual * MakeResiduals(const TQpDataBase *data)=0
virtual TQpVar * MakeVariables(const TQpDataBase *data)=0
The Residuals class calculates and stores the quantities that appear on the right-hand side of the li...
Definition: TQpResidual.h:62
void Clear_r1r2()
set the noncomplementarity components of the residual (the terms arising from the linear equalities i...
Double_t GetDualityGap()
Definition: TQpResidual.h:109
Double_t GetResidualNorm()
Definition: TQpResidual.h:108
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.
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...
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...
The Solver class contains methods for monitoring and checking the convergence status of the algorithm...
Definition: TQpSolverBase.h:73
TQpLinSolverBase * fSys
Definition: TQpSolverBase.h:76
Double_t fDnorm
Definition: TQpSolverBase.h:78
virtual void Start(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Implements a default starting-point heuristic.
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.
virtual Int_t DoStatus(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Int_t i, Double_t mu, Int_t level)
Tests for termination.
virtual Double_t FinalStepLength(TQpVar *iterate, TQpVar *step)
Implements a version of Mehrotra starting point heuristic, modified to ensure identical steps in the ...
TQpSolverBase & operator=(const TQpSolverBase &source)
Assignment operator.
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
virtual void Saxpy(TQpVar *b, Double_t alpha)
Perform a "saxpy" operation on all data vectors : x += alpha*y.
Definition: TQpVar.cxx:232
const Double_t sigma
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719