Actual source code: isils.c
1: #include "src/complementarity/impls/ssls/ssls.h"
3: // Context for ISXLS
4: // -- interior-point - smoothed fisched function
5: // -- active-set - reduced matrices formed
6: // - inherit properties of original system
7: // -- semismooth (S) - function not differentiable
8: // - merit function continuously differentiable
9: // - Fischer-Burmeister reformulation of complementarity
10: // - Billups composition for two finite bounds
11: // -- infeasible (I) - iterates not guaranteed to remain within bounds
12: // -- feasible (F) - iterates guaranteed to remain within bounds
13: // -- linesearch (LS) - Armijo rule on direction
14: //
15: // Many other reformulations are possible and combinations of
16: // feasible/infeasible and linesearch/trust region are possible.
17: //
18: // Basic theory
19: // Fischer-Burmeister reformulation is semismooth with a continuously
20: // differentiable merit function and strongly semismooth if the F has
21: // lipschitz continuous derivatives.
22: //
23: // Every accumulation point generated by the algorithm is a stationary
24: // point for the merit function. Stationary points of the merit function
25: // are solutions of the complementarity problem if
26: // a. the stationary point has a BD-regular subdifferential, or
27: // b. the Schur complement F'/F'_ff is a P_0-matrix where ff is the
28: // index set corresponding to the free variables.
29: //
30: // If one of the accumulation points has a BD-regular subdifferential then
31: // a. the entire sequence converges to this accumulation point at
32: // a local q-superlinear rate
33: // b. if in addition the reformulation is strongly semismooth near
34: // this accumulation point, then the algorithm converges at a
35: // local q-quadratic rate.
36: //
37: // The theory for the feasible version follows from the feasible descent
38: // algorithm framework.
39: //
40: // References:
41: // Billups, "Algorithms for Complementarity Problems and Generalized
42: // Equations," Ph.D thesis, University of Wisconsin - Madison, 1995.
43: // De Luca, Facchinei, Kanzow, "A Semismooth Equation Approach to the
44: // Solution of Nonlinear Complementarity Problems," Mathematical
45: // Programming, 75, pages 407-439, 1996.
46: // Ferris, Kanzow, Munson, "Feasible Descent Algorithms for Mixed
47: // Complementarity Problems," Mathematical Programming, 86,
48: // pages 475-497, 1999.
49: // Fischer, "A Special Newton-type Optimization Method," Optimization,
50: // 24, pages 269-284, 1992
51: // Munson, Facchinei, Ferris, Fischer, Kanzow, "The Semismooth Algorithm
52: // for Large Scale Complementarity Problems," Technical Report 99-06,
53: // University of Wisconsin - Madison, 1999.
55: /*------------------------------------------------------------*/
58: int Tao_ISLS_SFunction(TAO_SOLVER tao, TaoVec *X, double *fcn, void *solver)
59: {
60: TAO_SSLS *isls = (TAO_SSLS *)solver;
61: TaoVec *l, *u;
62: int info;
64: info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
66: info = TaoComputeConstraints(tao, X, isls->f); CHKERRQ(info);
67: info = isls->ff->SFischer(X, isls->f, l, u, isls->mu); CHKERRQ(info);
68: info = isls->ff->Norm2squared(&isls->merit_eqn); CHKERRQ(info);
70: // Add in the constraint for mu
71: isls->mucon = exp(isls->mu) - 1.0;
72: isls->merit_mu = isls->mucon*isls->mucon;
74: isls->merit = isls->merit_eqn + isls->merit_mu;
76: // Calculate the objective and merit function
77: *fcn = 0.5*isls->merit;
78: isls->merit = sqrt(isls->merit);
79: return 0;
80: }
82: /*------------------------------------------------------------*/
85: int Tao_ISLS_SFunctionGradient(TAO_SOLVER tao, TaoVec *X, double *fcn,
86: TaoVec *G, void *solver)
87: {
88: TAO_SSLS *isls = (TAO_SSLS *)solver;
89: TaoVec *l, *u;
90: TaoMat *J;
91: int info;
93: info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
94: info = TaoGetJacobian(tao, &J); CHKERRQ(info);
96: info = TaoComputeConstraints(tao, X, isls->f); CHKERRQ(info);
97: info = isls->ff->SFischer(X, isls->f, l, u, isls->mu); CHKERRQ(info);
98: info = isls->ff->Norm2squared(&isls->merit_eqn); CHKERRQ(info);
100: // Add in the constraint for mu
101: isls->mucon = exp(isls->mu) - 1.0;
102: isls->merit_mu = isls->mucon*isls->mucon;
104: isls->merit = isls->merit_eqn + isls->merit_mu;
106: // Calculate the objective and merit function
107: *fcn = 0.5*isls->merit;
108: isls->merit = sqrt(isls->merit);
110: info = TaoComputeJacobian(tao, X, J); CHKERRQ(info);
111: info = J->D_SFischer(X, isls->f, l, u, isls->mu,
112: isls->t1, isls->t2,
113: isls->da, isls->db, isls->dm); CHKERRQ(info);
114: isls->d_mucon = exp(isls->mu);
116: info = isls->t1->PointwiseMultiply(isls->ff, isls->db);
117: info = J->MultiplyTranspose(isls->t1, G); CHKERRQ(info);
118: info = isls->t1->PointwiseMultiply(isls->ff, isls->da); CHKERRQ(info);
119: info = G->Axpy(1.0, isls->t1); CHKERRQ(info);
120: info = isls->dm->Dot(isls->ff, &isls->g_mucon); CHKERRQ(info);
121: isls->g_mucon += exp(isls->mu)*isls->mucon;
122: return 0;
123: }
125: /*------------------------------------------------------------*/
128: static int TaoSolve_ISILS(TAO_SOLVER tao, void *solver)
129: {
130: TAO_SSLS *isls = (TAO_SSLS *)solver;
132: TaoVec *x; // Solution vector
133: TaoVec *l, *u; // Problem bounds
134: TaoMat *J; // Jacobian matrix
136: TaoIndexSet *FreeVariableSet; // Free variables
137: TaoIndexSet *FixedVariableSet;// Fixed variables
139: TaoMat *Jsub; // Reduced Jacobian
140: TaoVec *r1, *r2, *r3; // Temporary vectors
142: // TaoLinearSolver *lsolver;
144: double psi, ndpsi, normd, innerd, t=0;
145: double two_norm;
146: double sigma = 1e-4, beta = 0.5, minstep = TAO_EPSILON;
147: double muval, psival;
148: int iter=0, step_type, info, nn, nf;
149: TaoTerminateReason reason;
150: TaoTruth flag;
152: TaoFunctionBegin;
154: // Assume that Setup has been called!
156: // Start by obtaining the solution vector, Jacobian,
157: // and problem size
158: info = TaoGetSolution(tao, &x); CHKERRQ(info);
159: info = TaoGetJacobian(tao, &J); CHKERRQ(info);
160: info = x->GetDimension(&nn); CHKERRQ(info);
162: // Evaluate the bounds
163: info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
164: info = TaoEvaluateVariableBounds(tao, l, u); CHKERRQ(info);
166: // Create vector for derivative wrt smoothing parameter
167: info = x->Clone(&(isls->dm)); CHKERRQ(info);
169: // Create temporary vectors and index sets (for active set)
170: info = x->CreateIndexSet(&FreeVariableSet); CHKERRQ(info);
171: info = x->CreateIndexSet(&FixedVariableSet); CHKERRQ(info);
172: info = x->Clone(&r1); CHKERRQ(info);
173: info = x->Clone(&r2); CHKERRQ(info);
174: info = x->Clone(&r3); CHKERRQ(info);
176: // Create submatrix
177: info = J->CreateReducedMatrix(FreeVariableSet, FreeVariableSet, &Jsub); CHKERRQ(info);
179: // Initialize the merit function
180: // Used in the linesearch
181: // Termination based on different function
182: info = TaoSetMeritFunction(tao, Tao_ISLS_SFunction, Tao_ISLS_SFunctionGradient,
183: TAO_NULL, TAO_NULL, isls); CHKERRQ(info);
185: // Initialize the value for the smoothing parameter
186: isls->mu = isls->mu_init;
188: #if 0
189: // Choice of mu:
190: // Simple search
192: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
194: psival = fabs(isls->merit_eqn - nn*isls->merit_mu);
195: printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
196: isls->mu, isls->merit_eqn, nn*isls->merit_mu);
198: double fact;
199: if (isls->merit_eqn < nn*isls->merit_mu) {
200: fact = 0.5;
201: } else
202: {
203: fact = 2.0;
204: }
206: isls->mu *= fact;
207: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
208: printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
209: isls->mu, isls->merit_eqn, nn*isls->merit_mu);
210:
211: while (fabs(isls->merit_eqn - nn*isls->merit_mu) < psival) {
212: psival = fabs(isls->merit_eqn - nn*isls->merit_mu);
214: isls->mu *= fact;
215: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
216: printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
217: isls->mu, isls->merit_eqn, nn*isls->merit_mu);
218: }
220: isls->mu /= fact;
221: printf("mu_choice : %5.4e\n", isls->mu);
222: #endif
224: // DONE:
225: // 1. Eliminated smoothing parameter vector. Do not know how to deal
226: // with mixture of differentiable and nondifferentiable functions.
227: // Maybe let the smoothing parameter be a vector later.
228: // 2. Rewrote smoothed functions to provide the following information:
229: // - derivative of the function with respect to smoothing parameter
230: // - gives nonsmooth function when smoothing parameter == 0
231: // 3. Rewrote the function and function/gradient routines to add
232: // the smoothing constraint (e^{mu}-1 = 0) and the appropriate
233: // derivatives.
234: // 4. Modified the right-hand side of the direction calculation to
235: // account for the smoothing constraint. I.e. the direction in
236: // the smoothing parameter is substituted out of the problem.
237: // 5. Modified descent direction test to deal with the augmented
238: // problem with the smoothing parameter.
239: // 6. Hacked in a monotone linesearch to deal with the augmented system.
240: // 7. When smoothing parameter == 0, make sure we get the same
241: // behavior as asils.
242: // TODO:
243: // 8. Add method to optimally choose initial value for mu.
245:
246: // Calculate the fischer function at the initial iterate
247: Tao_ASLS_FunctionGradient(tao, x, &psi, isls->dpsi, isls);
248: info = isls->dpsi->Norm2(&ndpsi); CHKERRQ(info);
250: while (1) {
252: // Assert:
253: // Fischer function for terminatation
254: // Evaluated at current iterate
255: // Norm gradient of Fischer function calculated
257: printf("nd_merit : %5.4e\n", isls->merit);
258: printf("nd_ndpsi : %5.4e\n", ndpsi);
260: // Check the termination criteria
262: info = TaoMonitor(tao, iter++, isls->merit, ndpsi, 0.0, t, &reason); CHKERRQ(info);
263: if (TAO_CONTINUE_ITERATING != reason) break;
265: // Compute the smoothed Fischer function.
266: // This forms the basis for the interior-point method.
268: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
270: printf("merit : %5.4e\n", isls->merit);
271: printf("merit_eqn: %5.4e\n", isls->merit_eqn);
272: printf("merit_mu : %5.4e\n", isls->merit_mu);
273: printf("mu : %5.4e\n", isls->mu);
274: printf("mucon : %5.4e\n", isls->mucon);
275: printf("d_mucon : %5.4e\n", isls->d_mucon);
276: printf("g_mucon : %5.4e\n", isls->g_mucon);
278: #if 0
279: // Compute tolerances for linear system solver
281: // We want to set tolerances so that we maintain a superlinear
282: // asymptotic rate of convergence. Note: the tolerances set are
283: // for the reduced system. We really need to make sure that the
284: // full system satisfies the termination conditions.
286: // Superlinear asymptotic convergence:
287: // isls->atol = TaoMin(0.5, isls->merit*sqrt(isls->merit));
288: // isls->rtol = 0.0;
290: // Quadratic asymptotic convergence:
291: isls->atol = TaoMin(0.5, isls->merit*isls->merit);
292: isls->rtol = 0.0;
294: // Set the convergence tolerances.
295: info = TaoGetLinearSolver(tao, &lsolver); CHKERRQ(info);
296: // info = lsolver->SetTolerances(isls->rtol, isls->atol, 1e30, 10*nn); CHKERRQ(info);
297: // info = lsolver->SetTolerances(0.0, 0.0, 1e30, 10*nn); CHKERRQ(info);
298: // info = lsolver->SetTolerances(0.0, 0.0, 1e30, 10*nn); CHKERRQ(info);
299: #endif
301: // Determine whether we are taking a predictor or corrector step
302: // Calculate the direction in the smoothing parameter
303: // Substitute out of the problem
305: if ((isls->merit_eqn <= 1e-10) ||
306: (isls->merit_eqn <= nn*isls->merit_mu)) {
307: // Only take a predictor step when ||eqn|| <= const*||mu||
308: printf("Predictor step.\n");
309: step_type = 1;
310: isls->dmu = isls->mucon / isls->d_mucon;
311: printf("dmu : %5.4e\n", isls->dmu);
312: }
313: else {
314: // Otherwise take a correct step
315: printf("Corrector step.\n");
316: step_type = 0;
317: isls->dmu = 0;
318: }
319: info = isls->w->Waxpby(1.0, isls->ff, -isls->dmu, isls->dm);
321: // Assert: w vector now contains right-hand side
322: // We solve the system:
323: // [DPhi DMu ][dx ] = Phi
324: // [ 0 DRho][dmu] = Rho
325: // where Rho is the forcing function.
327: // Now calculate a free and fixed set of variables. The fixed set of
328: // variables are those for which d_b is approximately equal to zero.
329: // The definition of approximately changes as we approact the
330: // solution.
332: // No one rule is guaranteed to work in all cases. The following
333: // definition is based on the norm of the Jacobian matrix. If the
334: // norm is large, the tolerance becomes smaller.
336: info = J->Norm1(&(isls->identifier)); CHKERRQ(info);
337: isls->identifier = TaoMin(isls->merit, 1e-2) / (1 + isls->identifier);
339: info = isls->t1->SetToConstant(-isls->identifier); CHKERRQ(info);
340: info = isls->t2->SetToConstant( isls->identifier); CHKERRQ(info);
342: info = isls->t1->SetToConstant(0.0); CHKERRQ(info);
343: info = isls->t2->SetToConstant(0.0); CHKERRQ(info);
345: info = FixedVariableSet->WhichBetweenOrEqual(isls->t1, isls->db, isls->t2); CHKERRQ(info);
346: info = FreeVariableSet->ComplementOf(FixedVariableSet); CHKERRQ(info);
348: info = FixedVariableSet->GetSize(&nf);
349: if (nf) printf("Fixed size: %d\n", nf);
351: // Assert:
352: // Partition created.
353: // Now calculate the direction in the fixed variable space.
354: info = r1->SetReducedVec(isls->w, FixedVariableSet); CHKERRQ(info);
355: info = r2->SetReducedVec(isls->da, FixedVariableSet); CHKERRQ(info);
356: info = r1->PointwiseDivide(r1, r2); CHKERRQ(info);
358: // Throw the fixed variable direction into a global direction
359: info = isls->d->SetToZero(); CHKERRQ(info);
360: info = isls->d->ReducedXPY(r1, FixedVariableSet); CHKERRQ(info);
362: // Our direction in the FixedVariableSet is fixed. Calculate the
363: // information needed for the step in the FreeVariableSet. To
364: // do this, we need to know the diagonal perturbation and the
365: // right hand side.
366: info = r1->SetReducedVec(isls->da, FreeVariableSet); CHKERRQ(info);
367: info = r2->SetReducedVec(isls->w, FreeVariableSet); CHKERRQ(info);
368: info = r3->SetReducedVec(isls->db, FreeVariableSet); CHKERRQ(info);
370: info = r1->PointwiseDivide(r1, r3); CHKERRQ(info);
371: info = r2->PointwiseDivide(r2, r3); CHKERRQ(info);
373: // Assert:
374: // r1 is the diagonal perturbation
375: // r2 is the right hand side
376: // r3 is no longer needed
377: // Now need to modify r2 for our direction choice in the fixed
378: // variable set: calculate t1 = J*d, take the reduced vector
379: // of t1 and modify r2.
380: info = J->Multiply(isls->d, isls->t1); CHKERRQ(info);
381: info = r3->SetReducedVec(isls->t1, FreeVariableSet); CHKERRQ(info);
382: info = r2->Axpy(-1.0, r3);
384: // Calculate the reduced problem matrix
385: info = Jsub->SetReducedMatrix(J, FreeVariableSet, FreeVariableSet); CHKERRQ(info);
386: info = Jsub->AddDiagonal(r1); CHKERRQ(info);
388: // Assert:
389: // r1 is no longer needed
390: // Calculate the reduced direction. (Really negative of Newton
391: // direction. Therefore, rest of the code uses -d.)
392: info = TaoPreLinearSolve(tao, Jsub); CHKERRQ(info);
393: info = TaoLinearSolve(tao, Jsub, r2, r1, &flag); CHKERRQ(info);
395: // Add the direction in the free variables back into the real direction.
396: info = isls->d->ReducedXPY(r1, FreeVariableSet); CHKERRQ(info);
398: // Calculate the norm of the relative residual in the real space
399: info = isls->t1->PointwiseMultiply(isls->da, isls->d); CHKERRQ(info);
400: info = J->Multiply(isls->d, isls->t2); CHKERRQ(info);
401: info = isls->t2->PointwiseMultiply(isls->db, isls->t2); CHKERRQ(info);
402: info = isls->t1->Axpy( 1.0, isls->t2); CHKERRQ(info);
403: info = isls->t1->Axpy(-1.0, isls->w); CHKERRQ(info);
404: info = isls->t1->Norm2(&two_norm); CHKERRQ(info);
406: printf("eqn_resid: %5.4e\n", two_norm);
408: // Check direction for descent.
409: info = isls->d->Norm2squared(&normd); CHKERRQ(info);
410: normd += isls->dmu*isls->dmu;
411: normd = sqrt(normd);
413: info = isls->d->Dot(isls->dpsi, &innerd); CHKERRQ(info);
414: innerd += isls->g_mucon*isls->dmu;
416: printf("normd : %5.4e\n", normd);
417: printf("innerd : %5.4e\n", innerd);
419: if (innerd <= isls->delta*pow(normd, isls->rho)) {
420: // Descent direction test failed, use gradient direction
421: info = PetscLogInfo((tao, "TaoSolve_ISILS: %d: newton direction not descent\n", iter)); CHKERRQ(info);
422: printf("Gradient direction\n");
424: info = isls->d->CopyFrom(isls->dpsi); CHKERRQ(info);
425: if (1 == step_type) {
426: // Predictor step gradient
427: isls->dmu = isls->g_mucon;
428: }
429: else {
430: // Corrector step gradient
431: isls->dmu = 0;
432: }
434: info = isls->d->Dot(isls->dpsi, &innerd); CHKERRQ(info);
435: innerd += isls->g_mucon*isls->dmu;
437: printf("g_innerd : %5.4e\n", innerd);
438: }
440: // We have a descent direction. We need to take the negative to
441: // correct the sign of the direction.
442: info = isls->d->Negate(); CHKERRQ(info);
443: isls->dmu = -isls->dmu;
444: innerd = -innerd;
446: // Apply an armijo linesearch to find the new iterate.
447: innerd *= sigma;
448: muval = isls->mu;
449: t = 1;
451: while (t >= minstep) {
452: // Calculate iterate
453: info = isls->w->Waxpby(1.0, x, t, isls->d); CHKERRQ(info);
454: isls->mu = muval + t*isls->dmu;
455:
456: // Calculate function at new iterate
457: info = TaoComputeMeritFunction(tao, isls->w, &psival); CHKERRQ(info);
459: printf("psival: %5.4e psi+t*innerd: %5.4e\n", psival, psi + t*innerd);
461: // Check descent condition
462: if (psival <= psi + t*innerd) {
463: break;
464: }
465: t *= beta;
466: }
468: x->CopyFrom(isls->w);
470: printf("t : %5.4e\n", t);
471: printf("psi : %5.4e\n", psi);
472: printf("psival : %5.4e\n", psival);
473: printf("diff : %5.4e\n", psival - psi);
475: // Have an improvement for the merit function. Calculate the function
476: // values for the global merit function.
478: printf("ip_merit: %5.4e\n", isls->merit);
480: Tao_ASLS_FunctionGradient(tao, x, &psi, isls->dpsi, isls);
481: info = isls->dpsi->Norm2(&ndpsi); CHKERRQ(info);
482: }
484: info = TaoVecDestroy(isls->dm); CHKERRQ(info);
485: isls->dm = TAO_NULL;
487: info = TaoMatDestroy(Jsub); CHKERRQ(info);
488: info = TaoVecDestroy(r1); CHKERRQ(info);
489: info = TaoVecDestroy(r2); CHKERRQ(info);
490: info = TaoVecDestroy(r3); CHKERRQ(info);
492: info = TaoIndexSetDestroy(FreeVariableSet); CHKERRQ(info);
493: info = TaoIndexSetDestroy(FixedVariableSet); CHKERRQ(info);
494: TaoFunctionReturn(0);
495: }
497: /* ---------------------------------------------------------- */
501: int TaoCreate_ISILS(TAO_SOLVER tao)
502: {
503: TAO_SSLS *isls;
504: int info;
506: TaoFunctionBegin;
508: info = TaoNew(TAO_SSLS,&isls); CHKERRQ(info);
509: info = PetscLogObjectMemory(tao, sizeof(TAO_SSLS)); CHKERRQ(info);
511: isls->delta = 1e-10;
512: isls->rho = 2.1;
513: isls->mu_init = 0.1;
515: isls->identifier = 1e-5;
517: info=TaoSetTaoSolveRoutine(tao,TaoSolve_ISILS,(void*)isls); CHKERRQ(info);
518: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_SSLS,TaoSetDown_SSLS); CHKERRQ(info);
519: info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_SSLS); CHKERRQ(info);
520: info=TaoSetTaoViewRoutine(tao,TaoView_SSLS); CHKERRQ(info);
522: info = TaoCreateArmijoLineSearch(tao); CHKERRQ(info);
524: info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
525: info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);
527: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
528: info = TaoSetGradientTolerances(tao,1.0e-16,0.0,0.0); CHKERRQ(info);
529: info = TaoSetFunctionLowerBound(tao,1.0e-8); CHKERRQ(info);
531: TaoFunctionReturn(0);
532: }