Friday, March 10, 2017

ifthen adventures in GAMS (and R)

GAMS has a function ifthen which can be used in equations. For example:

e.. y =e= ifthen(x<1,(1+x)/2,sqrt(x));

This function implements:

\[y=\begin{cases}\displaystyle\frac{1+x}{2}&\text{if $x<1$}\\
                         \sqrt{x}&\text{if $x\ge 1$}\end{cases}\]

This is a perfectly smooth function: it is continuous and differentiable.

image_thumb[1]image_thumb

However when we use it in GAMS we are in for a surprise:

variable x,y;
equation
e;

e.. y =e= ifthen(x<1,(1+x)/2,sqrt(x));

x.fx = -1;

model m /all/
;
solve
m minimizing x using dnlp;

Here we fix \(x=-1\) so we just want to evaluate at that point. We see:

**** Exec Error at line 4: sqrt: FUNC DOMAIN: x < 0

The underlying reason is that GAMS always evaluates both branches. Obviously we would much prefer it would not do that.

If we replace fixing the variable by:

x.lo = -1;
x.l = 2;

i.e. start at \(x=2\) and move to the optimal value \(x=-1\), we see:

   C O N O P T 3   version 3.17A
   Copyright (C)   ARKI Consulting and Development A/S
                   Bagsvaerdvej 246 A
                   DK-2880 Bagsvaerd, Denmark


  Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK
     0   0        1.4142135624E+00 (Input point)

                  Pre-triangular equations:   0
                  Post-triangular equations:  1

     1   0        0.0000000000E+00 (After pre-processing)
     2   0        0.0000000000E+00 (After scaling)

** Feasible solution. Value of objective =    2.00000000000

  Iter Phase Ninf     Objective     RGmax    NSB   Step InItr MX OK
    11   3        6.5613264487E-03 1.0E+00     1 3.9E-03      F  F
    21   3        1.6608134823E-07 1.0E+00     1 6.0E-08      F  F
    31   3        3.2205059005E-13 1.0E+00     1 1.5E-11      F  T
    41   3        8.5912908646E-17 1.0E+00     1 2.2E-16      F  T
    51   3        5.0871607672E-23 1.0E+00     1 3.4E-21      F  T
    52   3        5.0871607672E-23 1.0E+00     1

** Feasible solution. Convergence too slow. The change in objective
   has been less than 3.0000E-12 for 20 consecutive iterations

The solver gets stuck at \(x=0\). This is obviously incorrect. The only hint of the source of the problem is listed here:

RESOURCE USAGE, LIMIT          0.062      1000.000
ITERATION COUNT, LIMIT        52    2000000000
EVALUATION ERRORS            803             0

This behavior is really a bug: we should see good feedback where the problem with domain errors originates from. Note that setting the limit of allowed evaluation errors (using option domlim) does not help in this case.

R

Interestingly R has something similar. Besides the standard if statement it has a vectorized ifelse function. The standard if is not vectorized which we can see from:

> curve(if (x<1) (1+x)/2 else sqrt(x),-5,5)

Warning message:

In if (x < 1) (1 + x)/2 else sqrt(x) :

  the condition has length > 1 and only the first element will be used

This only plots \(f(x)=\displaystyle\frac{1+x}{2}\):

image

Basically this if statement is interpreted as:

if (x[1]<1) (1+x)/2 else sqrt(x)

With ifelse we get better results:

> curve(ifelse(x<1,(1+x)/2,sqrt(x)),-5,5)
Warning message:
In sqrt(x) : NaNs produced

image

The plot is correct, but from the warning message we see that ifelse is also evaluating the \(\sqrt{x}\) for all \(x\). The ifelse function evaluates both branches.

To be precise this is only for vector arguments with some but not all \(x_i<1\). For scalar arguments (vectors of length one) or more general if all elements of the vector are either less than one or all are \(\ge 1\) things are ok:

image

Basically this is due to lazy but vector “all-or-nothing” evaluation.

Workaround

Of course the workaround is easy in this particular case. We can use something like:

e.. y =e= ifthen(x<1,(1+x)/2,sqrt(abs(x)));

Now the solver has no problems with this model, and can solve it very quickly (again with starting point \(x=2\) and optimal solution \(x=-1\)):

   C O N O P T 3   version 3.17A
   Copyright (C)   ARKI Consulting and Development A/S
                   Bagsvaerdvej 246 A
                   DK-2880 Bagsvaerd, Denmark


  Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK
     0   0        1.4142135624E+00 (Input point)

                  Pre-triangular equations:   0
                  Post-triangular equations:  1

     1   0        0.0000000000E+00 (After pre-processing)
     2   0        0.0000000000E+00 (After scaling)

** Feasible solution. Value of objective =    2.00000000000

  Iter Phase Ninf     Objective     RGmax    NSB   Step InItr MX OK
     4   3       -1.0000000000E+00 0.0E+00     0

** Optimal solution. There are no superbasic variables.

The same trick can be used in the R code:

curve(ifelse(x<1,(1+x)/2,sqrt(abs(x))),-5,5)

does not give any errors or warnings.