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

[Ravi Varadhan]

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