nlm函数在解析Hessian中失效

Question

背景：nlm函数在R是一个通用的优化程序，使用牛顿的方法。为了优化一个函数，牛顿方法需要函数，以及函数的第一和第二导数(梯度向量和Hessian矩阵)。在R中，nlm函数允许您指定与梯度和Hessian计算相对应的R函数，或者可以保留这些未指定的函数，并根据数值导数(通过deriv函数)提供数值解。通过提供计算梯度和Hessian的函数可以找到更精确的解，因此它是一个有用的特性。

我的问题是:当提供解析的Hessian函数时，nlm函数速度较慢，并且常常无法在合理的时间内收敛。我猜这是底层代码中的某种错误，但我很高兴错了。有没有一种方法可以使nlm更好地使用解析的Hessian矩阵？

示例:下面的R代码使用逻辑回归示例演示了这个问题，其中

log(Pr(Y=1)/Pr(Y=0)) = b0 + Xb

其中X是维数N的多元正规的p，b是长度p的系数的向量。

library(mvtnorm)
# example demonstrating a problem with NLM
expit <- function(mu) {1/(1+exp(-mu))}
mk.logit.data <- function(N,p){
  set.seed(1232)
  U = matrix(runif(p*p), nrow=p, ncol=p)
  S = 0.5*(U+t(U)) + p*diag(rep(1,p))
  X = rmvnorm(N, mean = runif(p, -1, 1), sigma = S)  
  Design = cbind(rep(1, N), X)
  beta = sort(sample(c(rep(0,p), runif(1))))
  y = rbinom(N, 1, expit(Design%*%beta))
 list(X=X,y=as.numeric(y),N=N,p=p) 
}

# function to calculate gradient vector at given coefficient values
logistic_gr <- function(beta, y, x, min=TRUE){
  mu = beta[1] + x %*% beta[-1]
  p = length(beta)
  n = length(y)
  D = cbind(rep(1,n), x)
  gri = matrix(nrow=n, ncol=p)
  for(j in 1:p){
    gri[,j] = D[,j]*(exp(-mu)*y-1+y)/(1+exp(-mu))
  }
  gr = apply(gri, 2, sum)
  if(min) gr = -gr
  gr
}

# function to calculate Hessian matrix at given coefficient values
logistic_hess <- function(beta, y, x, min=TRUE){
  # allow to fail with NA, NaN, Inf values
  mu = beta[1] + x %*% beta[-1]
  p = length(beta)
  n = length(y)
  D = cbind(rep(1,n), x)
  h = matrix(nrow=p, ncol=p)
  for(j in 1:p){
   for(k in 1:p){
     h[j,k] = -sum(D[,j]*D[,k]*(exp(-mu))/(1+exp(-mu))^2)
   }
  }
  if(min) h = -h
  h
}

# function to calculate likelihood (up to a constant) at given coefficient values
logistic_ll <- function(beta, y,x, gr=FALSE, he=FALSE, min=TRUE){
  mu = beta[1] + x %*% beta[-1]
  lli = log(expit(mu))*y + log(1-expit(mu))*(1-y)
  ll = sum(lli)
  if(is.na(ll) | is.infinite(ll)) ll = -1e16
  if(min) ll=-ll
  # the below specification is required for using analytic gradient/Hessian in nlm function
  if(gr) attr(ll, "gradient") <- logistic_gr(beta, y=y, x=x, min=min)
  if(he) attr(ll, "hessian") <- logistic_hess(beta, y=y, x=x, min=min)
  ll
}

第一个例子，使用p=3：

dat = mk.logit.data(N=100, p=3)

glm函数估计可供参考。nlm应该给出同样的答案，考虑到由于近似而产生的小误差。

(glm.sol <- glm(dat$y~dat$X, family=binomial()))$coefficients

> (Intercept)      dat$X1      dat$X2      dat$X3 
>  0.00981465  0.01068939  0.04417671  0.01625381 

# works when correct analytic gradient is specified
(nlm.sol1 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE,  y=dat$y, x=dat$X))$estimate
> [1] 0.009814547 0.010689396 0.044176627 0.016253966

# works, but less accurate when correct analytic hessian is specified (even though the routine notes convergence is probable)
(nlm.sol2 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE, he=TRUE, y=dat$y, x=dat$X, hessian = TRUE, check.analyticals=TRUE))$estimate
> [1] 0.009827701 0.010687278 0.044178416 0.016255630

但是当p更大时，问题变得很明显，这里是10。

dat = mk.logit.data(N=100, p=10)

同样，glm解决方案可供参考。nlm应该给出同样的答案，考虑到由于近似而产生的小误差。

(glm.sol <- glm(dat$y~dat$X, family=binomial()))$coefficients
> (Intercept)      dat$X1      dat$X2      dat$X3      dat$X4      dat$X5      dat$X6      dat$X7 
> -0.07071882 -0.08670003  0.16436630  0.01130549  0.17302058  0.03821008  0.08836471 -0.16578959 
>      dat$X8      dat$X9     dat$X10 
> -0.07515477 -0.08555075  0.29119963 

# works when correct analytic gradient is specified
(nlm.sol1 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE,  y=dat$y, x=dat$X))$estimate
> [1] -0.07071879 -0.08670005  0.16436632  0.01130550  0.17302057  0.03821009  0.08836472
> [8] -0.16578958 -0.07515478 -0.08555076  0.29119967

# fails to converge in 5000 iterations when correct analytic hessian is specified
(nlm.sol2 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE, he=TRUE,   y=dat$y, x=dat$X, hessian = TRUE, iterlim=5000, check.analyticals=TRUE))$estimate

> [1]  0.31602065 -0.06185190  0.10775381 -0.16748897  0.05032156  0.34176104  0.02118631
> [8] -0.01833671 -0.20364929  0.63713991  0.18390489

编辑:我还应该补充说，我已经通过多种不同的方法确认了正确的Hessian矩阵。

Answer 1

我尝试了这些代码，但一开始它似乎使用的是与我在CRAN上找到的不同的rmvnorm。我在dae软件包中找到了一个rmvnorm，然后在mvtnorm包中找到了一个。后者是可以使用的。

nlm()是关于上述发布时间的补丁。我目前正在尝试验证这些补丁，现在它似乎正常工作了。请注意，我是许多R的优化代码的作者，包括optim()中的3/5。

uottawa.ca的nashjc

密码在下面。

Answer 2

经修订的守则：

# example demonstrating a problem with NLM
expit <- function(mu) {1/(1+exp(-mu))}
mk.logit.data <- function(N,p){
  set.seed(1232)
  U = matrix(runif(p*p), nrow=p, ncol=p)
  S = 0.5*(U+t(U)) + p*diag(rep(1,p))
  X = rmvnorm(N, mean = runif(p, -1, 1), sigma = S)  
  Design = cbind(rep(1, N), X)
  beta = sort(sample(c(rep(0,p), runif(1))))
  y = rbinom(N, 1, expit(Design%*%beta))
  list(X=X,y=as.numeric(y),N=N,p=p) 
}

# function to calculate gradient vector at given coefficient values
logistic_gr <- function(beta, y, x, min=TRUE){
  mu = beta[1] + x %*% beta[-1]
  p = length(beta)
  n = length(y)
  D = cbind(rep(1,n), x)
  gri = matrix(nrow=n, ncol=p)
  for(j in 1:p){
    gri[,j] = D[,j]*(exp(-mu)*y-1+y)/(1+exp(-mu))
  }
  gr = apply(gri, 2, sum)
  if(min) gr = -gr
  gr
}

# function to calculate Hessian matrix at given coefficient values
logistic_hess <- function(beta, y, x, min=TRUE){
  # allow to fail with NA, NaN, Inf values
  mu = beta[1] + x %*% beta[-1]
  p = length(beta)
  n = length(y)
  D = cbind(rep(1,n), x)
  h = matrix(nrow=p, ncol=p)
  for(j in 1:p){
    for(k in 1:p){
      h[j,k] = -sum(D[,j]*D[,k]*(exp(-mu))/(1+exp(-mu))^2)
    }
  }
  if(min) h = -h
  h
}

# function to calculate likelihood (up to a constant) at given coefficient values
logistic_ll <- function(beta, y,x, gr=FALSE, he=FALSE, min=TRUE){
  mu = beta[1] + x %*% beta[-1]
  lli = log(expit(mu))*y + log(1-expit(mu))*(1-y)
  ll = sum(lli)
  if(is.na(ll) | is.infinite(ll)) ll = -1e16
  if(min) ll=-ll
  # the below specification is required for using analytic gradient/Hessian in nlm function
  if(gr) attr(ll, "gradient") <- logistic_gr(beta, y=y, x=x, min=min)
  if(he) attr(ll, "hessian") <- logistic_hess(beta, y=y, x=x, min=min)
  ll
}

##!!!! NOTE: Must have this library loaded
library(mvtnorm)
dat = mk.logit.data(N=100, p=3)
(glm.sol <- glm(dat$y~dat$X, family=binomial()))$coefficients
# works when correct analytic gradient is specified
(nlm.sol1 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE,  y=dat$y, x=dat$X))$estimate
# works, but less accurate when correct analytic hessian is specified (even though the routine notes convergence is probable)
(nlm.sol2 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE, he=TRUE, y=dat$y, x=dat$X, hessian = TRUE, check.analyticals=TRUE))$estimate

dat = mk.logit.data(N=100, p=10)

# Again, glm solution for reference. nlm should give the same answer, allowing for small errors due to approximation.

(glm.sol <- glm(dat$y~dat$X, family=binomial()))$coefficients

# works when correct analytic gradient is specified
(nlm.sol1 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE,  y=dat$y, x=dat$X))$estimate

# fails to converge in 5000 iterations when correct analytic hessian is specified
(nlm.sol2 <- nlm(p=runif(dat$p+1), f=logistic_ll, gr=TRUE, he=TRUE,   y=dat$y, x=dat$X, hessian = TRUE, iterlim=5000, check.analyticals=TRUE))$estimate