Re: [Rd] constrOptim( ): conflict between help page and code

From: Ravi Varadhan <>
Date: Thu, 17 Jun 2010 13:22:55 -0400

One more bug in `constrOptim' that I have reported before and have also suggested a fix, but this has also not been corrected:

if (obj > obj.old) break

This is obviously wrong when the user is trying to maximize by specifying control$fnscale = -1.

The correct statement is the following:

if (obj > obj.old * sign(mu)) break


-----Original Message-----
From: Ravi Varadhan [] Sent: Thursday, June 17, 2010 10:55 AM
To: 'Ravi Varadhan'; 'Duncan Murdoch'; 'John Nolan' Cc: ''
Subject: RE: [Rd] constrOptim( ): conflict between help page and code

In `constrOptim', there is also a mistake in the reporting of function and gradient counts. These counts, as reported, correspond to that of the vary last "inner" iteration. However, they should be cumulatively summed up over each outer iteration.

I have fixed this and the convergence criterion issue that I mentioned before.

Currently, constrOptim can only use the Nelder-Mead when the analytic gradient is not specified. I have added a numerical gradient option so that it can use BFGS or other optimization functions even when the analytic gradient is not specified.

It would be desirable if Tom Lumley can commit these changes to constrOptim. In the meanwhile, I can send the function with these upgrades to anyone who is interested.


-----Original Message-----
From: [] On Behalf Of Ravi Varadhan Sent: Thursday, June 17, 2010 9:53 AM
To: 'Duncan Murdoch'; 'John Nolan'
Subject: Re: [Rd] constrOptim( ): conflict between help page and code


You are correct that there is inconsistency. The feasible region should be ui %*% theta - ci > 0, so that log(ui %*% theta - ci) is defined.

There is a more serious problem in termination criterion for iterations:

        if (is.finite(r) && is.finite(r.old) && abs(r - r.old)/(outer.eps + 
            abs(r - r.old)) < outer.eps) 

Ideally convergence is achieved when |r - r.old| is small. But according to the above code, the logical condition inside the if(.) will be TRUE only when abs(r - r.old) < (outer.eps)^2 (for small outer.eps). For example, let |r - r.old| = outer.eps. The above condition becomes: if (0.5 < outer.eps) break, which will obviously never happen for values of outer.eps < 1/2. Now, if outer.eps = 1.e-05 (default) and we let |r - r.old| <= 1.e-10, then the above condition will be satisfied.

In short, the termination criterion used is too stringent. Better termination criteria are:

        if (is.finite(r) && is.finite(r.old) && abs(r - r.old) < outer.eps) break


        if (is.finite(r) && is.finite(r.old) && abs(r - r.old)/(1 + abs(r + r.old)/2) < outer.eps) break


-----Original Message-----
From: [] On Behalf Of Duncan Murdoch Sent: Thursday, June 17, 2010 4:38 AM
To: John Nolan
Subject: Re: [Rd] constrOptim( ): conflict between help page and code

John Nolan wrote:
> There is a contradiction between what the help page says and what constrOptim actually
> does with the constraints. The issue is what happens on the boundary.

I don't know if this was a recent change, but the current help says:

"The starting value must be in the interior of the feasible region, but the minimum may be on the boundary."

The boundary is not in the interior.

Duncan Murdoch
> The help page says
> The feasible region is defined by ‘ui %*% theta - ci >= 0’,
> but the R code for constrOptim reads
> if (any(ui %*% theta - ci <= 0))
> stop("initial value not feasible")
> The following example shows that when the initial point is on the boundary of the
> feasibility region, you get the above error message and execution stops.
>> fn <- function(x) { return(sum(x)) }
>> ui <- diag(rep(1,2))
>> ci <- matrix(0,nrow=2,ncol=1)
>> constrOptim( c(0,0), fn, NULL, ui, ci )
> Error in optim(theta.old, fun, gradient, control = control, method = method, :
> function cannot be evaluated at initial parameters
>> version
> platform i386-pc-mingw32
> arch i386
> os mingw32
> system i386, mingw32
> status
> major 2
> minor 10.0
> year 2009
> month 10
> day 26
> svn rev 50208
> language R
> version.string R version 2.10.0 (2009-10-26)
> In contrast, at a different place in constrOptim - the internal function R -
> the boundary of the feasibility region is allowed: if (any(gi < 0)) return(NaN),
> and it seems to explicitly allow boundaries at another place:
> allowing gi==0 and interpreting log(gi) as -Inf.
> John
> ...........................................................................
> John P. Nolan
> Math/Stat Department
> 227 Gray Hall
> American University
> 4400 Massachusetts Avenue, NW
> Washington, DC 20016-8050
> 202.885.3140 voice
> 202.885.3155 fax
> ...........................................................................
> ______________________________________________
> mailing list
> mailing list mailing list mailing list Received on Thu 17 Jun 2010 - 17:25:26 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 17 Jun 2010 - 21:41:10 GMT.

Mailing list information is available at Please read the posting guide before posting to the list.

list of date sections of archive