title: Handling box constraints in Nelder-Mead optimisation by distorting the parameter space
body: <p>I have a question on a specific implementation of a Nelder-Mead algorithm (1) that handles box contraints in an unusual way. I cannot find in anything about it in any paper (25 papers), textbook (searched 4 of them) or the internet.</p>
  
  <p>I have a typical optimisation problem: <code>min f(x)</code> with a box constraint <code>-0.25 &lt;= x_i &lt;= 250</code></p>
  
  <p>The expected approach would be using a penalty function and make sure that all instances of <code>f(x)</code> are "unattractive" when x is out of bounds. </p>
  
  <p>The algorithm works differently: the implementation in question does not touch <code>f(x)</code>. Instead it distorts the parameter space using an inverse hyperbolic tangens <code>atanh(f)</code>. Now the simplex algorithm can freely operate in a space without bounds and pick just any point. Before it gets <code>f(x)</code> in order to assess the solution at x the algorithm switches back into normal space.</p>
  
  <p>At a first glance I found the idea ingenious. This way we avoid the disadvantages of penalty functions. But now I am having doubts. The distorted space affects termination behaviour. One termination criterion is the size of the simplex. By inflating the parameter space with <code>atanh(x)</code> we also inflate the simplex size. </p>
  
  <p>Experiments with the algorithm also show that it does not work as intended. I do not yet understand how this happens, but I do get results that are out of bounds. I can say that almost half of the returned local minima are out of bounds.</p>
  
  <p>As an example, take a look at nmkb() optimising the rosenbrook function when we gradually change the width of the box constraint:</p>
  
  <pre><code>rosbkext &lt;- function(x) {
    # Extended Rosenbrock function
    n &lt;- length(x)
    sum (100*(x[1:(n-1)]^2 - x[2:n])^2 + (x[1:(n-1)] - 1)^2)
  }
  
  np &lt;- 6 #12
  for (box in c(2, 4, 12, 24, 32, 64, 128)) {
    set.seed(123)
    p0 &lt;- rnorm(np)
    p0[p0 &gt; +2] &lt;- +2 - 1E-8
    p0[p0 &lt; -2] &lt;- -2 + 1E-8
  
    ctrl &lt;- list(maxfeval = 5E4, tol = 1E-8)
    o &lt;- nmkb(fn = rosbkext, par = p0, lower = -box, upper = +box, control = ctrl)
    print(o$message)
    cat("f(", format(o$par, digits = 2), ") =", format(o$value, digits=3), "\n")
  }
  </code></pre>
  
  <p>The output shows that it claims to converge but it does not in three cases. And it does that for bounds of (-2,2) and (-12,12). I might accept that but then it also fails at (-128, 128). I also tried the same with the unconstrained dfoptim::nmk(). No trouble there. It converges perfectly.</p>
  
  <pre><code>[1] "Successful convergence"
  f( -0.99  0.98  0.97  0.95  0.90  0.81 ) = 3.97 
  [1] "Successful convergence"
  f( 1 1 1 1 1 1 ) = 4.42e-09 
  [1] "Successful convergence"
  f( -0.99  0.98  0.97  0.95  0.90  0.81 ) = 3.97 
  [1] "Successful convergence"
  f( 1 1 1 1 1 1 ) = 1.3e-08 
  [1] "Successful convergence"
  f( 1 1 1 1 1 1 ) = 4.22e-09 
  [1] "Successful convergence"
  f( 1 1 1 1 1 1 ) = 8.22e-09 
  [1] "Successful convergence"
  f( -0.99  0.98  0.97  0.95  0.90  0.81 ) = 3.97 
  </code></pre>
  
  <p>Why does the constrained algorithm have more trouble converging than the unconstrained one?</p>
  
  <hr>
  
  <p>Footnote (1):  I am referring to the Nelder-Mead implementation used in the optimx package in R. This package calls another package dfoptim with the nmkb-function.</p>