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;
149: TaoInt nn, nf;
150: TaoTerminateReason reason;
151: TaoTruth flag;
153: TaoFunctionBegin;
155: // Assume that Setup has been called!
157: // Start by obtaining the solution vector, Jacobian,
158: // and problem size
159: info = TaoGetSolution(tao, &x); CHKERRQ(info);
160: info = TaoGetJacobian(tao, &J); CHKERRQ(info);
161: info = x->GetDimension(&nn); CHKERRQ(info);
163: // Evaluate the bounds
164: info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
165: info = TaoEvaluateVariableBounds(tao, l, u); CHKERRQ(info);
167: // Create vector for derivative wrt smoothing parameter
168: info = x->Clone(&(isls->dm)); CHKERRQ(info);
170: // Create temporary vectors and index sets (for active set)
171: info = x->CreateIndexSet(&FreeVariableSet); CHKERRQ(info);
172: info = x->CreateIndexSet(&FixedVariableSet); CHKERRQ(info);
173: info = x->Clone(&r1); CHKERRQ(info);
174: info = x->Clone(&r2); CHKERRQ(info);
175: info = x->Clone(&r3); CHKERRQ(info);
177: // Create submatrix
178: info = J->CreateReducedMatrix(FreeVariableSet, FreeVariableSet, &Jsub); CHKERRQ(info);
180: // Initialize the merit function
181: // Used in the linesearch
182: // Termination based on different function
183: info = TaoSetMeritFunction(tao, Tao_ISLS_SFunction, Tao_ISLS_SFunctionGradient,
184: TAO_NULL, TAO_NULL, isls); CHKERRQ(info);
186: // Initialize the value for the smoothing parameter
187: isls->mu = isls->mu_init;
189: #if 0
190: // Choice of mu:
191: // Simple search
193: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
195: psival = fabs(isls->merit_eqn - nn*isls->merit_mu);
196: printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
197: isls->mu, isls->merit_eqn, nn*isls->merit_mu);
199: double fact;
200: if (isls->merit_eqn < nn*isls->merit_mu) {
201: fact = 0.5;
202: } else
203: {
204: fact = 2.0;
205: }
207: isls->mu *= fact;
208: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
209: printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
210: isls->mu, isls->merit_eqn, nn*isls->merit_mu);
211:
212: while (fabs(isls->merit_eqn - nn*isls->merit_mu) < psival) {
213: psival = fabs(isls->merit_eqn - nn*isls->merit_mu);
215: isls->mu *= fact;
216: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
217: printf("mu: %5.4e merit_eqn: %5.4e merit_mu: %5.4e\n",
218: isls->mu, isls->merit_eqn, nn*isls->merit_mu);
219: }
221: isls->mu /= fact;
222: printf("mu_choice : %5.4e\n", isls->mu);
223: #endif
225: // DONE:
226: // 1. Eliminated smoothing parameter vector. Do not know how to deal
227: // with mixture of differentiable and nondifferentiable functions.
228: // Maybe let the smoothing parameter be a vector later.
229: // 2. Rewrote smoothed functions to provide the following information:
230: // - derivative of the function with respect to smoothing parameter
231: // - gives nonsmooth function when smoothing parameter == 0
232: // 3. Rewrote the function and function/gradient routines to add
233: // the smoothing constraint (e^{mu}-1 = 0) and the appropriate
234: // derivatives.
235: // 4. Modified the right-hand side of the direction calculation to
236: // account for the smoothing constraint. I.e. the direction in
237: // the smoothing parameter is substituted out of the problem.
238: // 5. Modified descent direction test to deal with the augmented
239: // problem with the smoothing parameter.
240: // 6. Hacked in a monotone linesearch to deal with the augmented system.
241: // 7. When smoothing parameter == 0, make sure we get the same
242: // behavior as asils.
243: // TODO:
244: // 8. Add method to optimally choose initial value for mu.
246:
247: // Calculate the fischer function at the initial iterate
248: Tao_ASLS_FunctionGradient(tao, x, &psi, isls->dpsi, isls);
249: info = isls->dpsi->Norm2(&ndpsi); CHKERRQ(info);
251: while (1) {
253: // Assert:
254: // Fischer function for terminatation
255: // Evaluated at current iterate
256: // Norm gradient of Fischer function calculated
258: printf("nd_merit : %5.4e\n", isls->merit);
259: printf("nd_ndpsi : %5.4e\n", ndpsi);
261: // Check the termination criteria
263: info = TaoMonitor(tao, iter++, isls->merit, ndpsi, 0.0, t, &reason); CHKERRQ(info);
264: if (TAO_CONTINUE_ITERATING != reason) break;
266: // Compute the smoothed Fischer function.
267: // This forms the basis for the interior-point method.
269: info = TaoComputeMeritFunctionGradient(tao, x, &psi, isls->dpsi); CHKERRQ(info);
271: printf("merit : %5.4e\n", isls->merit);
272: printf("merit_eqn: %5.4e\n", isls->merit_eqn);
273: printf("merit_mu : %5.4e\n", isls->merit_mu);
274: printf("mu : %5.4e\n", isls->mu);
275: printf("mucon : %5.4e\n", isls->mucon);
276: printf("d_mucon : %5.4e\n", isls->d_mucon);
277: printf("g_mucon : %5.4e\n", isls->g_mucon);
279: #if 0
280: // Compute tolerances for linear system solver
282: // We want to set tolerances so that we maintain a superlinear
283: // asymptotic rate of convergence. Note: the tolerances set are
284: // for the reduced system. We really need to make sure that the
285: // full system satisfies the termination conditions.
287: // Superlinear asymptotic convergence:
288: // isls->atol = TaoMin(0.5, isls->merit*sqrt(isls->merit));
289: // isls->rtol = 0.0;
291: // Quadratic asymptotic convergence:
292: isls->atol = TaoMin(0.5, isls->merit*isls->merit);
293: isls->rtol = 0.0;
295: // Set the convergence tolerances.
296: info = TaoGetLinearSolver(tao, &lsolver); CHKERRQ(info);
297: // info = lsolver->SetTolerances(isls->rtol, isls->atol, 1e30, 10*nn); CHKERRQ(info);
298: // info = lsolver->SetTolerances(0.0, 0.0, 1e30, 10*nn); CHKERRQ(info);
299: // info = lsolver->SetTolerances(0.0, 0.0, 1e30, 10*nn); CHKERRQ(info);
300: #endif
302: // Determine whether we are taking a predictor or corrector step
303: // Calculate the direction in the smoothing parameter
304: // Substitute out of the problem
306: if ((isls->merit_eqn <= 1e-10) ||
307: (isls->merit_eqn <= nn*isls->merit_mu)) {
308: // Only take a predictor step when ||eqn|| <= const*||mu||
309: printf("Predictor step.\n");
310: step_type = 1;
311: isls->dmu = isls->mucon / isls->d_mucon;
312: printf("dmu : %5.4e\n", isls->dmu);
313: }
314: else {
315: // Otherwise take a correct step
316: printf("Corrector step.\n");
317: step_type = 0;
318: isls->dmu = 0;
319: }
320: info = isls->w->Waxpby(1.0, isls->ff, -isls->dmu, isls->dm);
322: // Assert: w vector now contains right-hand side
323: // We solve the system:
324: // [DPhi DMu ][dx ] = Phi
325: // [ 0 DRho][dmu] = Rho
326: // where Rho is the forcing function.
328: // Now calculate a free and fixed set of variables. The fixed set of
329: // variables are those for which d_b is approximately equal to zero.
330: // The definition of approximately changes as we approact the
331: // solution.
333: // No one rule is guaranteed to work in all cases. The following
334: // definition is based on the norm of the Jacobian matrix. If the
335: // norm is large, the tolerance becomes smaller.
337: info = J->Norm1(&(isls->identifier)); CHKERRQ(info);
338: isls->identifier = TaoMin(isls->merit, 1e-2) / (1 + isls->identifier);
340: info = isls->t1->SetToConstant(-isls->identifier); CHKERRQ(info);
341: info = isls->t2->SetToConstant( isls->identifier); CHKERRQ(info);
343: info = isls->t1->SetToConstant(0.0); CHKERRQ(info);
344: info = isls->t2->SetToConstant(0.0); CHKERRQ(info);
346: info = FixedVariableSet->WhichBetweenOrEqual(isls->t1, isls->db, isls->t2); CHKERRQ(info);
347: info = FreeVariableSet->ComplementOf(FixedVariableSet); CHKERRQ(info);
349: info = FixedVariableSet->GetSize(&nf);
350: if (nf) PetscPrintf(((PetscObject)tao)->comm,"Fixed size: %d\n", nf);
352: // Assert:
353: // Partition created.
354: // Now calculate the direction in the fixed variable space.
355: info = r1->SetReducedVec(isls->w, FixedVariableSet); CHKERRQ(info);
356: info = r2->SetReducedVec(isls->da, FixedVariableSet); CHKERRQ(info);
357: info = r1->PointwiseDivide(r1, r2); CHKERRQ(info);
359: // Throw the fixed variable direction into a global direction
360: info = isls->d->SetToZero(); CHKERRQ(info);
361: info = isls->d->ReducedXPY(r1, FixedVariableSet); CHKERRQ(info);
363: // Our direction in the FixedVariableSet is fixed. Calculate the
364: // information needed for the step in the FreeVariableSet. To
365: // do this, we need to know the diagonal perturbation and the
366: // right hand side.
367: info = r1->SetReducedVec(isls->da, FreeVariableSet); CHKERRQ(info);
368: info = r2->SetReducedVec(isls->w, FreeVariableSet); CHKERRQ(info);
369: info = r3->SetReducedVec(isls->db, FreeVariableSet); CHKERRQ(info);
371: info = r1->PointwiseDivide(r1, r3); CHKERRQ(info);
372: info = r2->PointwiseDivide(r2, r3); CHKERRQ(info);
374: // Assert:
375: // r1 is the diagonal perturbation
376: // r2 is the right hand side
377: // r3 is no longer needed
378: // Now need to modify r2 for our direction choice in the fixed
379: // variable set: calculate t1 = J*d, take the reduced vector
380: // of t1 and modify r2.
381: info = J->Multiply(isls->d, isls->t1); CHKERRQ(info);
382: info = r3->SetReducedVec(isls->t1, FreeVariableSet); CHKERRQ(info);
383: info = r2->Axpy(-1.0, r3);
385: // Calculate the reduced problem matrix
386: info = Jsub->SetReducedMatrix(J, FreeVariableSet, FreeVariableSet); CHKERRQ(info);
387: info = Jsub->AddDiagonal(r1); CHKERRQ(info);
389: // Assert:
390: // r1 is no longer needed
391: // Calculate the reduced direction. (Really negative of Newton
392: // direction. Therefore, rest of the code uses -d.)
393: info = TaoPreLinearSolve(tao, Jsub); CHKERRQ(info);
394: info = TaoLinearSolve(tao, Jsub, r2, r1, &flag); CHKERRQ(info);
396: // Add the direction in the free variables back into the real direction.
397: info = isls->d->ReducedXPY(r1, FreeVariableSet); CHKERRQ(info);
399: // Calculate the norm of the relative residual in the real space
400: info = isls->t1->PointwiseMultiply(isls->da, isls->d); CHKERRQ(info);
401: info = J->Multiply(isls->d, isls->t2); CHKERRQ(info);
402: info = isls->t2->PointwiseMultiply(isls->db, isls->t2); CHKERRQ(info);
403: info = isls->t1->Axpy( 1.0, isls->t2); CHKERRQ(info);
404: info = isls->t1->Axpy(-1.0, isls->w); CHKERRQ(info);
405: info = isls->t1->Norm2(&two_norm); CHKERRQ(info);
407: printf("eqn_resid: %5.4e\n", two_norm);
409: // Check direction for descent.
410: info = isls->d->Norm2squared(&normd); CHKERRQ(info);
411: normd += isls->dmu*isls->dmu;
412: normd = sqrt(normd);
414: info = isls->d->Dot(isls->dpsi, &innerd); CHKERRQ(info);
415: innerd += isls->g_mucon*isls->dmu;
417: printf("normd : %5.4e\n", normd);
418: printf("innerd : %5.4e\n", innerd);
420: if (innerd <= isls->delta*pow(normd, isls->rho)) {
421: // Descent direction test failed, use gradient direction
422: info = PetscInfo1(tao, "TaoSolve_ISILS: %d: newton direction not descent\n", iter); CHKERRQ(info);
423: printf("Gradient direction\n");
425: info = isls->d->CopyFrom(isls->dpsi); CHKERRQ(info);
426: if (1 == step_type) {
427: // Predictor step gradient
428: isls->dmu = isls->g_mucon;
429: }
430: else {
431: // Corrector step gradient
432: isls->dmu = 0;
433: }
435: info = isls->d->Dot(isls->dpsi, &innerd); CHKERRQ(info);
436: innerd += isls->g_mucon*isls->dmu;
438: printf("g_innerd : %5.4e\n", innerd);
439: }
441: // We have a descent direction. We need to take the negative to
442: // correct the sign of the direction.
443: info = isls->d->Negate(); CHKERRQ(info);
444: isls->dmu = -isls->dmu;
445: innerd = -innerd;
447: // Apply an armijo linesearch to find the new iterate.
448: innerd *= sigma;
449: muval = isls->mu;
450: t = 1;
452: while (t >= minstep) {
453: // Calculate iterate
454: info = isls->w->Waxpby(1.0, x, t, isls->d); CHKERRQ(info);
455: isls->mu = muval + t*isls->dmu;
456:
457: // Calculate function at new iterate
458: info = TaoComputeMeritFunction(tao, isls->w, &psival); CHKERRQ(info);
460: printf("psival: %5.4e psi+t*innerd: %5.4e\n", psival, psi + t*innerd);
462: // Check descent condition
463: if (psival <= psi + t*innerd) {
464: break;
465: }
466: t *= beta;
467: }
469: x->CopyFrom(isls->w);
471: printf("t : %5.4e\n", t);
472: printf("psi : %5.4e\n", psi);
473: printf("psival : %5.4e\n", psival);
474: printf("diff : %5.4e\n", psival - psi);
476: // Have an improvement for the merit function. Calculate the function
477: // values for the global merit function.
479: printf("ip_merit: %5.4e\n", isls->merit);
481: Tao_ASLS_FunctionGradient(tao, x, &psi, isls->dpsi, isls);
482: info = isls->dpsi->Norm2(&ndpsi); CHKERRQ(info);
483: }
485: info = TaoVecDestroy(isls->dm); CHKERRQ(info);
486: isls->dm = TAO_NULL;
488: info = TaoMatDestroy(Jsub); CHKERRQ(info);
489: info = TaoVecDestroy(r1); CHKERRQ(info);
490: info = TaoVecDestroy(r2); CHKERRQ(info);
491: info = TaoVecDestroy(r3); CHKERRQ(info);
493: info = TaoIndexSetDestroy(FreeVariableSet); CHKERRQ(info);
494: info = TaoIndexSetDestroy(FixedVariableSet); CHKERRQ(info);
495: TaoFunctionReturn(0);
496: }
498: /* ---------------------------------------------------------- */
502: int TaoCreate_ISILS(TAO_SOLVER tao)
503: {
504: TAO_SSLS *isls;
505: int info;
507: TaoFunctionBegin;
509: info = TaoNew(TAO_SSLS,&isls); CHKERRQ(info);
510: info = PetscLogObjectMemory(tao, sizeof(TAO_SSLS)); CHKERRQ(info);
512: isls->delta = 1e-10;
513: isls->rho = 2.1;
514: isls->mu_init = 0.1;
516: isls->identifier = 1e-5;
518: info=TaoSetTaoSolveRoutine(tao,TaoSolve_ISILS,(void*)isls); CHKERRQ(info);
519: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_SSLS,TaoSetDown_SSLS); CHKERRQ(info);
520: info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_SSLS); CHKERRQ(info);
521: info=TaoSetTaoViewRoutine(tao,TaoView_SSLS); CHKERRQ(info);
523: info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
524: info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);
526: info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
527: info = TaoSetGradientTolerances(tao,1.0e-16,0.0,0.0); CHKERRQ(info);
528: info = TaoSetFunctionLowerBound(tao,1.0e-8); CHKERRQ(info);
530: TaoFunctionReturn(0);
531: }