share/contrib/augmented_lagrangian.mac

   1 /* augmented Lagrangian method for constrained optimization
   2  *
   3  * copyright Robert Dodier, October 2005
   4  * Released under the terms of the GNU General Public License
   5  *
   6  * See http://www-fp.mcs.anl.gov/otc/Guide/OptWeb/continuous/constrained/nonlinearcon/auglag.html
   7  * and http://www.cs.ubc.ca/spider/ascher/542/chap10.pdf
   8  *
   9  * This implementation minimizes the augmented Lagrangian by
  10  * applying LBFGS (a quasi-Newton algorithm).
  11  *
  12  * Function: augmented_lagrangian_method (FOM, xx, constraints, yy, [optional_args])
  13  *
  14  * FOM           = figure of merit expression
  15  * xx            = list of variables over which to minimize
  16  * constraints   = list of expressions to be held equal to zero
  17  * yy            = list of initial guesses for xx
  18  *
  19  * optional_args = more arguments of the form symbol = value
  20  *                 optional arguments recognized are:
  21  *
  22  *                 niter           = number of iterations of the augmented Lagrangian algorithm
  23  *                 lbfgs_tolerance = tolerance supplied to LBFGS
  24  *                 iprint          = IPRINT parameter supplied to LBFGS
  25  *                 %lambda         = initial value of %lambda
  26  *
  27  * LBFGS (to solve grad L = 0) has to be loaded before calling augmented_lagrangian_method.
  28  *
  29  * Example 1:
  30  *
  31    load (lbfgs);
  32    load (augmented_lagrangian);
  33    FOM: x^2 + 2*y^2;
  34    xx: [x, y];
  35    C: [x + y - 1];
  36    yy: [1, 1];
  37    augmented_lagrangian_method (FOM, xx, C, yy, niter = 20);
  38  *
  39  *  => [[x = 0.6661347183787, y = 0.33306735918935], %lambda = [- 1.331471514325457]]
  40  *
  41  * Example 2:
  42  * same as Example 1, but with gradient computed by Lisp functions
  43  * Construct Lisp functions suitable for this example via translate
  44  *
  45    load (lbfgs);
  46    load (augmented_lagrangian);
  47    FOM: x^2 + 2*y^2;
  48    xx: [x, y];
  49    C: [x + y - 1];
  50    yy: [1, 1];
  51    F(x, y) := ''FOM;
  52    F1(x, y) := ''(diff(F(x, y), x));
  53    F2(x, y) := ''(diff(F(x, y), y));
  54    translate (F, F1, F2);
  55    :lisp (mremprop '|$f| 'mexpr)
  56    :lisp (mremprop '|$f1| 'mexpr)
  57    :lisp (mremprop '|$f2| 'mexpr)
  58    augmented_lagrangian_method ('[F(x, y), [F1(x, y), F2(x, y)]], xx, C, yy, niter = 20);
  59  *
  60  *  => [[x = 0.6661347183787, y = 0.33306735918935], %lambda = [- 1.331471514325457]]
  61  */
  62
  63 /* fboundp detects various kinds of Maxima functions */
  64 fboundp (a) :=
  65     ?fboundp (a) # false
  66  or ?mget (a, ?mexpr) # false
  67  or ?get (a, ?mfexpr\*) # false
  68  or ?mget (a, ?mmacro) # false;
  69
  70 if not fboundp (lbfgs) then load (lbfgs);
  71
  72 with_parameters ([L]) ::= buildq ([a : subst (":", "=",  ev (L [1])), e : rest (L)], block (a, splice (e)));
  73
  74 augmented_lagrangian_method (FOM, xx, constraints, yy, [optional_args]) :=
  75
  76   block ([n, augmented_lagrangian, %lambda, %nu, FOM_augment,
  77           niter : 10, lbfgs_tolerance : 1e-4, iprint : [1, 1]],
  78
  79   nconstraints: length (constraints),
  80
  81   FOM_augment :
  82       apply ("+", makelist (%lambda[i] * constraints[i], i, 1, nconstraints))
  83     + (1/2) * apply ("+", makelist (%nu[i] * constraints[i]^2, i, 1, nconstraints)),
  84
  85   if listp (FOM)
  86     then block ([FOM_proper, FOM_grad],
  87         FOM_proper : first (FOM) + FOM_augment,
  88         FOM_grad : second (FOM) + map (lambda ([x], diff (FOM_augment, x)), xx),
  89         augmented_lagrangian : [FOM_proper, FOM_grad])
  90     else augmented_lagrangian: FOM + FOM_augment,
  91
  92   %lambda: makelist (1, i, 1, nconstraints),
  93
  94   %nu: makelist (1, i, 1, nconstraints),
  95
  96   with_parameters (optional_args,
  97
  98     for i:1 thru niter do
  99
 100       (soln: lbfgs (augmented_lagrangian, xx, yy, lbfgs_tolerance, iprint),
 101
 102       yy: map (rhs, soln),
 103
 104       %lambda: %lambda + apply ("+", %nu * map (lambda ([c], subst (soln, c)), constraints)),
 105
 106       %nu : 2 * %nu),
 107
 108     [map ("=", xx, yy), '%lambda = %lambda]));