title: Error in optim/constrOptim "initial value in 'vmmin' is not finite"
body: <p>My problem (and what I think could help to solve it) is explained until the line "FOR REPRODUCTION". After that I just posted my code, just in case reproducing might help to solve it.</p>
  
  <p>I use optim and constrOptim.nl to solve an optimization problem with constraints within the function g (see below) I wrote. </p>
  
  <p>I know that the initial values used below are not ideal, but I chose them 
  because they cause the problem I face in a shorter program. I use this program to calibrate model parameters to data and there this problem also occurs for better inital values, higher tolerances etc. </p>
  
  <p><strong>The Error</strong></p>
  
  <p>I call the function get_par I wrote with:</p>
  
  <pre><code>v&lt;-c(0.12504710,0.09329359,0.06778733, 0.04883216, 0.04187344,0.02886261,0.02332951,0.02178576,0.02282214,0.02956336,0.03478598)
  Ti=1/12
  x&lt;-log(cbind(0.8,0.85,0.9,0.95,0.975,1,1.025,1.05,1.1,1.15,1.2))
  
  g(par2=c(-5,5),v=v,Ti=Ti,x=x)
  </code></pre>
  
  <p>Then I get</p>
  
  <pre><code>Error in optim: inital value 'vmmin' is not finite.
  </code></pre>
  
  <p><strong>What I have observed so far</strong></p>
  
  <p>So I started to debug my code to find out where exactly this error occurs. The error occurs in the function g (see below) in the line ( with the values sigma=5,m=-5,y=(x-m)/sigma,vtilde=v/12)</p>
  
  <pre><code>#print(paste("vW: sigma: ",sigma,"mv:",mv))
  
  argmin&lt;-constrOptim.nl(par=c(3*sigma,sigma,mv/2),fn=f,hin.jac=hinv.jac,
  hin=hinv,heq.jac=heqv.jac,heq=heqv,control.outerlist(trace=T),
  control.optim=list(abstol=10^(-10)),y=y,vtilde=vtilde,sigma=sigma)
  </code></pre>
  
  <p>The Trace of the funciton constrOptim.nl displays</p>
  
  <pre><code>Outer iteration:  18 
  Min(hin):  1.026858e-19 Max(abs(heq)):  0 
  par:  10 9.99998 1.02686e-19 
  fval =   6399 
  </code></pre>
  
  <p>for the last iteration. I guess that there is some sort of a numerical problem with 1.02686e-19 appearing in the last iteration.</p>
  
  <p>I had a look into the function constrOptim.nl and albama (with debug() ), and the error exactly occurs in the line </p>
  
  <pre><code>theta.old &lt;- theta
  atemp &lt;- optim(par = theta, fn = fun, gr = grad, control = control.optim, 
      method = "BFGS", hessian = TRUE, ...)
  </code></pre>
  
  <p>where theta=theta.old has the value</p>
  
  <pre><code>Browse[2]&gt; theta.old
  [1]  1.000002e+01  9.999985e+00 -3.349452e-20
  </code></pre>
  
  <p>Hence it has an entry that is just below zero (its absolute value is even smaller than machine accuracy, isn't it?).</p>
  
  <p>When you look at the function fun you realize that it calls the function </p>
  
  <pre><code>R:
  function (theta, theta.old, ...) 
  {
  gi &lt;- hin(theta, ...)
  if (any(gi &lt; 0)) 
      return(NaN)
  gi.old &lt;- hin(theta.old, ...)
  hjac &lt;- hin.jac(theta.old, ...)
  bar &lt;- sum(gi.old * log(gi) - hjac %*% theta)
  if (!is.finite(bar)) 
      bar &lt;- -Inf
  fn(theta, ...) - mu * bar
  }
  </code></pre>
  
  <p>hin(theta,...)=hinv(theta,...) returns a vector with a negative entry, thus this function returns NaN. I suppose that this should cause the error message: "Error in optim: inital value 'vmmin' is not finite". My question is now:</p>
  
  <p>How can I fix that? I thought of forcing the program to terminate somehow when such small values occur, but I have not managed to do that yet. What do yo suggest?</p>
  
  <p>Many thanks in advance,</p>
  
  <p>FOR REPRODUCTION:</p>
  
  <p><strong>Here is my program:</strong> </p>
  
  <p>The functions hinv, hinv.jac, heq and heq.jac are just for the constraints. The function where I optimize is g. </p>
  
  <pre><code>library(alabama)
  library(dfoptim)
  
  #function f, par = (c,d,atilde)
  f&lt;-function(par3,y,vtilde,sigma){
  sum((par3[3]+par3[2]*y+par3[1]*sqrt(y^2+1)-vtilde)^2)
  }
  
  #Equality/Inequality constraints
  
  heqv&lt;-function(par3,y,vtilde,sigma){
  
  J1&lt;-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2) 
  J2&lt;-matrix(0,nrow=3,ncol=3)
  J2[1:2,1:2]&lt;-J1
  J2[3,3]&lt;-1
  
  j&lt;-J2%*%par3
  j[2]-2*sqrt(2)*sigma
  }
  
  #Jacobian-matrix
  hinv.jac&lt;-function(par3,y,vtilde,sigma){
  
  #J1, J2: Drehungen für die constraints
  J1&lt;-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2) 
  J2&lt;-matrix(0,nrow=3,ncol=3)
  J2[1:2,1:2]&lt;-J1
  J2[3,3]&lt;-1
  
  hjac&lt;-matrix(cbind(1,-1,0,0,0,0,0,0,0,0,1,-1),nrow=4)%*%J2
  hjac
  }
  
  hinv&lt;-function(par3,y,vtilde,sigma){
  
  #J1, J2: Drehungen für die constraints
  J1&lt;-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2) 
  J2&lt;-matrix(0,nrow=3,ncol=3)
  J2[1:2,1:2]&lt;-J1
  J2[3,3]&lt;-1
  
  j&lt;-J2%*%par3
  h&lt;-rep(NA,4)
  h[1]&lt;- j[1]
  h[2]&lt;- sqrt(2)*2*sigma-j[1]
  h[3]&lt;-j[3]
  h[4]&lt;-max(vtilde)-j[3]
  h
  }
  
  #Jacobian-matrix
  heqv.jac&lt;-function(par3,y,vtilde,sigma){
  #J1, J2: Drehungen für die constraints
  J1&lt;-matrix(1/2*cbind(sqrt(2),sqrt(2),-sqrt(2),sqrt(2)),nrow=2,ncol=2) 
  J2&lt;-matrix(0,nrow=3,ncol=3)
  J2[1:2,1:2]&lt;-J1
  J2[3,3]&lt;-1
  
  cbind(J2[2,1],J2[2,2],0)
  }
  
  #function g input: par2= (m,sigma): optimization of function f
  
  g&lt;-function(par2,v,Ti,x){
  
  
  #definition of parameters being used
  m&lt;-par2[1]
  sigma&lt;-par2[2]
  y&lt;-(x-m)/sigma #Transformation von x zu y gemäß paper
  vtilde&lt;-Ti*v
  mv&lt;-max(vtilde)
  
  #print(paste("vW: sigma: ",sigma,"mv:",mv))
  argmin&lt;-constrOptim.nl(par=c(3*sigma,sigma,mv/2),fn=f,hin.jac=hinv.jac,hin=hinv,heq.jac=heqv.jac,heq=heqv,control.outer=list(trace=F),control.optim=list(abstol=10^(-10)),y=y,vtilde=vtilde,sigma=sigma)
  argmin$par
  }
  </code></pre>