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

[Ravi Varadhan]

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.

Best,
Ravi.

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

Nolan,

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) 
            break

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        

or 

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

Best,
Ravi.

-----Original Message-----
From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org] On Behalf Of Duncan Murdoch
Sent: Thursday, June 17, 2010 4:38 AM
To: John Nolan
Cc: r-devel at r-project.org
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
>
>  jpnolan at american.edu
>  202.885.3140 voice
>  202.885.3155 fax
>  http://academic2.american.edu/~jpnolan
>  ...........................................................................
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
R-devel at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel