1 #include <../src/tao/complementarity/impls/ssls/ssls.h>
2 /*
3 Context for ASXLS
4 -- active-set - reduced matrices formed
5 - inherit properties of original system
6 -- semismooth (S) - function not differentiable
7 - merit function continuously differentiable
8 - Fischer-Burmeister reformulation of complementarity
9 - Billups composition for two finite bounds
10 -- infeasible (I) - iterates not guaranteed to remain within bounds
11 -- feasible (F) - iterates guaranteed to remain within bounds
12 -- linesearch (LS) - Armijo rule on direction
13
14 Many other reformulations are possible and combinations of
15 feasible/infeasible and linesearch/trust region are possible.
16
17 Basic theory
18 Fischer-Burmeister reformulation is semismooth with a continuously
19 differentiable merit function and strongly semismooth if the F has
20 lipschitz continuous derivatives.
21
22 Every accumulation point generated by the algorithm is a stationary
23 point for the merit function. Stationary points of the merit function
24 are solutions of the complementarity problem if
25 a. the stationary point has a BD-regular subdifferential, or
26 b. the Schur complement F'/F'_ff is a P_0-matrix where ff is the
27 index set corresponding to the free variables.
28
29 If one of the accumulation points has a BD-regular subdifferential then
30 a. the entire sequence converges to this accumulation point at
31 a local q-superlinear rate
32 b. if in addition the reformulation is strongly semismooth near
33 this accumulation point, then the algorithm converges at a
34 local q-quadratic rate.
35
36 The theory for the feasible version follows from the feasible descent
37 algorithm framework. See {cite}`billups:algorithms`, {cite}`deluca.facchinei.ea:semismooth`,
38 {cite}`ferris.kanzow.ea:feasible`, {cite}`fischer:special`, and {cite}`munson.facchinei.ea:semismooth`.
39 */
40
TaoSetUp_ASFLS(Tao tao)41 static PetscErrorCode TaoSetUp_ASFLS(Tao tao)
42 {
43 TAO_SSLS *asls = (TAO_SSLS *)tao->data;
44
45 PetscFunctionBegin;
46 PetscCall(VecDuplicate(tao->solution, &tao->gradient));
47 PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
48 PetscCall(VecDuplicate(tao->solution, &asls->ff));
49 PetscCall(VecDuplicate(tao->solution, &asls->dpsi));
50 PetscCall(VecDuplicate(tao->solution, &asls->da));
51 PetscCall(VecDuplicate(tao->solution, &asls->db));
52 PetscCall(VecDuplicate(tao->solution, &asls->t1));
53 PetscCall(VecDuplicate(tao->solution, &asls->t2));
54 PetscCall(VecDuplicate(tao->solution, &asls->w));
55 asls->fixed = NULL;
56 asls->free = NULL;
57 asls->J_sub = NULL;
58 asls->Jpre_sub = NULL;
59 asls->r1 = NULL;
60 asls->r2 = NULL;
61 asls->r3 = NULL;
62 asls->dxfree = NULL;
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
Tao_ASLS_FunctionGradient(TaoLineSearch ls,Vec X,PetscReal * fcn,Vec G,void * ptr)66 static PetscErrorCode Tao_ASLS_FunctionGradient(TaoLineSearch ls, Vec X, PetscReal *fcn, Vec G, void *ptr)
67 {
68 Tao tao = (Tao)ptr;
69 TAO_SSLS *asls = (TAO_SSLS *)tao->data;
70
71 PetscFunctionBegin;
72 PetscCall(TaoComputeConstraints(tao, X, tao->constraints));
73 PetscCall(VecFischer(X, tao->constraints, tao->XL, tao->XU, asls->ff));
74 PetscCall(VecNorm(asls->ff, NORM_2, &asls->merit));
75 *fcn = 0.5 * asls->merit * asls->merit;
76 PetscCall(TaoComputeJacobian(tao, tao->solution, tao->jacobian, tao->jacobian_pre));
77
78 PetscCall(MatDFischer(tao->jacobian, tao->solution, tao->constraints, tao->XL, tao->XU, asls->t1, asls->t2, asls->da, asls->db));
79 PetscCall(VecPointwiseMult(asls->t1, asls->ff, asls->db));
80 PetscCall(MatMultTranspose(tao->jacobian, asls->t1, G));
81 PetscCall(VecPointwiseMult(asls->t1, asls->ff, asls->da));
82 PetscCall(VecAXPY(G, 1.0, asls->t1));
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85
TaoDestroy_ASFLS(Tao tao)86 static PetscErrorCode TaoDestroy_ASFLS(Tao tao)
87 {
88 TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
89
90 PetscFunctionBegin;
91 PetscCall(VecDestroy(&ssls->ff));
92 PetscCall(VecDestroy(&ssls->dpsi));
93 PetscCall(VecDestroy(&ssls->da));
94 PetscCall(VecDestroy(&ssls->db));
95 PetscCall(VecDestroy(&ssls->w));
96 PetscCall(VecDestroy(&ssls->t1));
97 PetscCall(VecDestroy(&ssls->t2));
98 PetscCall(VecDestroy(&ssls->r1));
99 PetscCall(VecDestroy(&ssls->r2));
100 PetscCall(VecDestroy(&ssls->r3));
101 PetscCall(VecDestroy(&ssls->dxfree));
102 PetscCall(MatDestroy(&ssls->J_sub));
103 PetscCall(MatDestroy(&ssls->Jpre_sub));
104 PetscCall(ISDestroy(&ssls->fixed));
105 PetscCall(ISDestroy(&ssls->free));
106 PetscCall(KSPDestroy(&tao->ksp));
107 PetscCall(PetscFree(tao->data));
108 PetscFunctionReturn(PETSC_SUCCESS);
109 }
110
TaoSolve_ASFLS(Tao tao)111 static PetscErrorCode TaoSolve_ASFLS(Tao tao)
112 {
113 TAO_SSLS *asls = (TAO_SSLS *)tao->data;
114 PetscReal psi, ndpsi, normd, innerd, t = 0;
115 PetscInt nf;
116 TaoLineSearchConvergedReason ls_reason;
117
118 PetscFunctionBegin;
119 /* Assume that Setup has been called!
120 Set the structure for the Jacobian and create a linear solver. */
121
122 PetscCall(TaoComputeVariableBounds(tao));
123 PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, Tao_ASLS_FunctionGradient, tao));
124 PetscCall(TaoLineSearchSetObjectiveRoutine(tao->linesearch, Tao_SSLS_Function, tao));
125 PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
126
127 PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
128
129 /* Calculate the function value and fischer function value at the
130 current iterate */
131 PetscCall(TaoLineSearchComputeObjectiveAndGradient(tao->linesearch, tao->solution, &psi, asls->dpsi));
132 PetscCall(VecNorm(asls->dpsi, NORM_2, &ndpsi));
133
134 tao->reason = TAO_CONTINUE_ITERATING;
135 while (1) {
136 /* Check the converged criteria */
137 PetscCall(PetscInfo(tao, "iter %" PetscInt_FMT ", merit: %g, ||dpsi||: %g\n", tao->niter, (double)asls->merit, (double)ndpsi));
138 PetscCall(TaoLogConvergenceHistory(tao, asls->merit, ndpsi, 0.0, tao->ksp_its));
139 PetscCall(TaoMonitor(tao, tao->niter, asls->merit, ndpsi, 0.0, t));
140 PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
141 if (TAO_CONTINUE_ITERATING != tao->reason) break;
142
143 /* Call general purpose update function */
144 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
145 tao->niter++;
146
147 /* We are going to solve a linear system of equations. We need to
148 set the tolerances for the solve so that we maintain an asymptotic
149 rate of convergence that is superlinear.
150 Note: these tolerances are for the reduced system. We really need
151 to make sure that the full system satisfies the full-space conditions.
152
153 This rule gives superlinear asymptotic convergence
154 asls->atol = min(0.5, asls->merit*sqrt(asls->merit));
155 asls->rtol = 0.0;
156
157 This rule gives quadratic asymptotic convergence
158 asls->atol = min(0.5, asls->merit*asls->merit);
159 asls->rtol = 0.0;
160
161 Calculate a free and fixed set of variables. The fixed set of
162 variables are those for the d_b is approximately equal to zero.
163 The definition of approximately changes as we approach the solution
164 to the problem.
165
166 No one rule is guaranteed to work in all cases. The following
167 definition is based on the norm of the Jacobian matrix. If the
168 norm is large, the tolerance becomes smaller. */
169 PetscCall(MatNorm(tao->jacobian, NORM_1, &asls->identifier));
170 asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier);
171
172 PetscCall(VecSet(asls->t1, -asls->identifier));
173 PetscCall(VecSet(asls->t2, asls->identifier));
174
175 PetscCall(ISDestroy(&asls->fixed));
176 PetscCall(ISDestroy(&asls->free));
177 PetscCall(VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed));
178 PetscCall(ISComplementVec(asls->fixed, asls->t1, &asls->free));
179
180 PetscCall(ISGetSize(asls->fixed, &nf));
181 PetscCall(PetscInfo(tao, "Number of fixed variables: %" PetscInt_FMT "\n", nf));
182
183 /* We now have our partition. Now calculate the direction in the
184 fixed variable space. */
185 PetscCall(TaoVecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1));
186 PetscCall(TaoVecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2));
187 PetscCall(VecPointwiseDivide(asls->r1, asls->r1, asls->r2));
188 PetscCall(VecSet(tao->stepdirection, 0.0));
189 PetscCall(VecISAXPY(tao->stepdirection, asls->fixed, 1.0, asls->r1));
190
191 /* Our direction in the Fixed Variable Set is fixed. Calculate the
192 information needed for the step in the Free Variable Set. To
193 do this, we need to know the diagonal perturbation and the
194 right-hand side. */
195
196 PetscCall(TaoVecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1));
197 PetscCall(TaoVecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2));
198 PetscCall(TaoVecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3));
199 PetscCall(VecPointwiseDivide(asls->r1, asls->r1, asls->r3));
200 PetscCall(VecPointwiseDivide(asls->r2, asls->r2, asls->r3));
201
202 /* r1 is the diagonal perturbation
203 r2 is the right-hand side
204 r3 is no longer needed
205
206 Now need to modify r2 for our direction choice in the fixed
207 variable set: calculate t1 = J*d, take the reduced vector
208 of t1 and modify r2. */
209
210 PetscCall(MatMult(tao->jacobian, tao->stepdirection, asls->t1));
211 PetscCall(TaoVecGetSubVec(asls->t1, asls->free, tao->subset_type, 0.0, &asls->r3));
212 PetscCall(VecAXPY(asls->r2, -1.0, asls->r3));
213
214 /* Calculate the reduced problem matrix and the direction */
215 PetscCall(TaoMatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type, &asls->J_sub));
216 if (tao->jacobian != tao->jacobian_pre) {
217 PetscCall(TaoMatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub));
218 } else {
219 PetscCall(MatDestroy(&asls->Jpre_sub));
220 asls->Jpre_sub = asls->J_sub;
221 PetscCall(PetscObjectReference((PetscObject)asls->Jpre_sub));
222 }
223 PetscCall(MatDiagonalSet(asls->J_sub, asls->r1, ADD_VALUES));
224 PetscCall(TaoVecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree));
225 PetscCall(VecSet(asls->dxfree, 0.0));
226
227 /* Calculate the reduced direction. (Really negative of Newton
228 direction. Therefore, rest of the code uses -d.) */
229 PetscCall(KSPReset(tao->ksp));
230 PetscCall(KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub));
231 PetscCall(KSPSolve(tao->ksp, asls->r2, asls->dxfree));
232 PetscCall(KSPGetIterationNumber(tao->ksp, &tao->ksp_its));
233 tao->ksp_tot_its += tao->ksp_its;
234
235 /* Add the direction in the free variables back into the real direction. */
236 PetscCall(VecISAXPY(tao->stepdirection, asls->free, 1.0, asls->dxfree));
237
238 /* Check the projected real direction for descent and if not, use the negative
239 gradient direction. */
240 PetscCall(VecCopy(tao->stepdirection, asls->w));
241 PetscCall(VecScale(asls->w, -1.0));
242 PetscCall(VecBoundGradientProjection(asls->w, tao->solution, tao->XL, tao->XU, asls->w));
243 PetscCall(VecNorm(asls->w, NORM_2, &normd));
244 PetscCall(VecDot(asls->w, asls->dpsi, &innerd));
245
246 if (innerd >= -asls->delta * PetscPowReal(normd, asls->rho)) {
247 PetscCall(PetscInfo(tao, "Gradient direction: %5.4e.\n", (double)innerd));
248 PetscCall(PetscInfo(tao, "Iteration %" PetscInt_FMT ": newton direction not descent\n", tao->niter));
249 PetscCall(VecCopy(asls->dpsi, tao->stepdirection));
250 PetscCall(VecDot(asls->dpsi, tao->stepdirection, &innerd));
251 }
252
253 PetscCall(VecScale(tao->stepdirection, -1.0));
254 innerd = -innerd;
255
256 /* We now have a correct descent direction. Apply a linesearch to
257 find the new iterate. */
258 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
259 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &psi, asls->dpsi, tao->stepdirection, &t, &ls_reason));
260 PetscCall(VecNorm(asls->dpsi, NORM_2, &ndpsi));
261 }
262 PetscFunctionReturn(PETSC_SUCCESS);
263 }
264
265 /*MC
266 TAOASFLS - Active-set feasible linesearch algorithm for solving complementarity constraints
267
268 Options Database Keys:
269 + -tao_ssls_delta - descent test fraction
270 - -tao_ssls_rho - descent test power
271
272 Level: beginner
273
274 Note:
275 See {cite}`billups:algorithms`, {cite}`deluca.facchinei.ea:semismooth`,
276 {cite}`ferris.kanzow.ea:feasible`, {cite}`fischer:special`, and {cite}`munson.facchinei.ea:semismooth`.
277
278 .seealso: `Tao`, `TaoType`, `TAOASILS`
279 M*/
TaoCreate_ASFLS(Tao tao)280 PETSC_EXTERN PetscErrorCode TaoCreate_ASFLS(Tao tao)
281 {
282 TAO_SSLS *asls;
283 const char *armijo_type = TAOLINESEARCHARMIJO;
284
285 PetscFunctionBegin;
286 PetscCall(PetscNew(&asls));
287 tao->data = (void *)asls;
288 tao->ops->solve = TaoSolve_ASFLS;
289 tao->ops->setup = TaoSetUp_ASFLS;
290 tao->ops->view = TaoView_SSLS;
291 tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
292 tao->ops->destroy = TaoDestroy_ASFLS;
293 tao->subset_type = TAO_SUBSET_SUBVEC;
294 asls->delta = 1e-10;
295 asls->rho = 2.1;
296 asls->fixed = NULL;
297 asls->free = NULL;
298 asls->J_sub = NULL;
299 asls->Jpre_sub = NULL;
300 asls->w = NULL;
301 asls->r1 = NULL;
302 asls->r2 = NULL;
303 asls->r3 = NULL;
304 asls->t1 = NULL;
305 asls->t2 = NULL;
306 asls->dxfree = NULL;
307 asls->identifier = 1e-5;
308
309 PetscCall(TaoParametersInitialize(tao));
310
311 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
312 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
313 PetscCall(TaoLineSearchSetType(tao->linesearch, armijo_type));
314 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
315 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
316
317 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
318 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
319 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
320 PetscCall(KSPSetFromOptions(tao->ksp));
321
322 /* Override default settings (unless already changed) */
323 PetscObjectParameterSetDefault(tao, max_it, 2000);
324 PetscObjectParameterSetDefault(tao, max_funcs, 4000);
325 PetscObjectParameterSetDefault(tao, gttol, 0);
326 PetscObjectParameterSetDefault(tao, grtol, 0);
327 PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1.0e-6 : 1.0e-16);
328 PetscObjectParameterSetDefault(tao, fmin, PetscDefined(USE_REAL_SINGLE) ? 1.0e-4 : 1.0e-8);
329 PetscFunctionReturn(PETSC_SUCCESS);
330 }
331