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

From: Ravi Varadhan <rvaradhan_at_jhmi.edu>
Date: Thu, 17 Jun 2010 09:52:31 -0400

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 Received on Thu 17 Jun 2010 - 13:54:48 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 - 16:01:04 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive