R - 有解析梯度的数值误差？

Question

我有下面的代码：

theta=0.05 
n=1000 
m=200 
r=rnorm(2000) 

#ER check function 
nu=Vectorize(function(a,tau){return(abs(tau-(a<0))*a^2)}) 

#Selecting 10 lowest sum values (lowest10 function returns indices) 
lowest10=function(x){ 
    values=sort(x)[1:min(10,length(x))] 
    indices=match(values,x) 
    return(indices) 
} 
sym.expectile=function(beta,e,abs.r){return(beta[1]+beta[2]*e+beta[3]*abs.r)} 

ERsum=function(beta,tau,start,end){ 
    y=r[(start+1):end] 
    X1=rep(1,n-1) 
    X3=abs(r[start:(end-1)]) 
    X2=c() 
    X2[1]=e.sym.optimal[start-m] 
    for (i in 2:(n-1)){ 
    X2[i]=sym.expectile(beta,X2[i-1],X3[i-1]) 
    } 
    X=matrix(c(X1,X2,X3),ncol=3) 
    res=y-X%*%beta 
    sum.nu=mean(nu(res,tau)) 
    return(sum.nu) 
} 

ERsum.gr=function(beta,tau,start,end){ 
    y=r[(start+1):end] 
    X1=rep(1,n-1) 
    X3=abs(r[start:(end-1)]) 
    X2=c() 
    X2[1]=e.sym.optimal[start-m] 
    for (i in 2:(n-1)){ 
    X2[i]=sym.expectile(beta,X2[i-1],X3[i-1]) 
    } 
    X=matrix(c(X1,X2,X3),ncol=3) 
    partial.beta0=c() 
    for (i in 1:(n-1)){partial.beta0[i]=-(1-beta[2]^(i))/(1-beta[2])} 
    gr.beta0=2/T*sum(abs(tau-(y<X%*%beta))*(y-X%*%beta)*partial.beta0)/1000 
    partial.beta1=c() 
    partial.beta1[1]=-X2[1] 
    for (i in 2:(n-1)){partial.beta1[i]=partial.beta1[i-1]*beta[2]-X2[i]} 
    gr.beta1=2/T*sum(abs(tau-(y<X%*%beta))*(y-X%*%beta)*partial.beta1)/1000 
    partial.beta2=c() 
    partial.beta2[1]=-X3[1] 
    for (i in 2:(n-1)){partial.beta2[i]=partial.beta2[i-1]*beta[2]-X3[i]} 
    gr.beta2=2/T*sum(abs(tau-(y<X%*%beta))*(y-X%*%beta)*partial.beta2)/1000 
    c(gr.beta0,gr.beta1,gr.beta2) 
} 

beta=matrix(nrow=1e4,ncol=3) 
beta[,1]=runif(1e4,-1,0)#beta0 
beta[,2]=runif(1e4,0,1)#beta1 
beta[,3]=runif(1e4,-1,0)#beta2 

e.sym.optimal=c() 
tau.found.sym.optim=0.02234724 
library('expectreg') 
e.sym.optimal[1]=expectile(r[1:m],tau.found.sym.optim) 

ERsums.sym=c() 
for (i in 1:nrow(beta)){ 
    ERsums.sym[i]=ERsum(beta[i,],tau.found.sym.optim,m+1,m+n) 
} 

initialbeta.esym=beta[lowest10(ERsums.sym),] 

intermedietebeta.esym=matrix(ncol=3,nrow=10) 
for (i in 1:10){ 
    intermedietebeta.esym[i,]=optim(initialbeta.esym[i,],ERsum, 
            gr=ERsum.gr,tau=tau.found.sym.optim, 
            start=m+1,end=m+n, 
            method="BFGS")$par 
} 

我试图取代optimx在Optim功能，而且得到了以下错误：

Error: Gradient function might be wrong - check it!

要检查我的渐变是确定我试着使用numDeriv中的函数grad评估梯度函数的值，并直接调用我的ERsum.gr函数。对于样品矢量

beta 
[1] -0.8256490 0.7146256 -0.4945032 

我得到以下结果：

>grad(function(beta) ERsum(c(beta[1],beta[2],beta[3]),tau.found.sym.optim,m+1,m+n),beta) 
[1] -0.6703170 2.8812666 -0.5573101 
> ERsum.gr2(beta,tau.found.sym.optim,m+1,m+n) 
[1] -0.6696467 2.8783853 -0.5567527 

因此，这里是我的问题：是有可能，这些差异只是造成舍去了一些partial.beta0数值误差，部分.beta1，partial.beta2哪些只是代表渐变的组成部分？我认为是这样，因为如果我的梯度解析公式错过某些东西，差异可能会大得多，但我怎么能确定呢？如果这是一种情况，还有其他方法可以获得更准确的渐变值吗？

Answer 1

即使你解决了这个问题是否真的是一个合适的梯度，我认为这个问题太复杂了。如果你拿出gr参数，并尝试以仅optimx代替optim运行，您可以：

Error in intermedietebeta.esym[i, ] <- optimx(initialbeta.esym[i, ], ERsum, : 
    number of items to replace is not a multiple of replacement length 

这可能涉及到的事实，如由optim返回optimx不会返回相同的结构：

> optimx(initialbeta.esym[i,],ERsum, 
+         tau=tau.found.sym.optim, 
+         start=m+1,end=m+n, 
+         method="BFGS")$par 
NULL 
> optimx(initialbeta.esym[i,],ERsum, 
+         tau=tau.found.sym.optim, 
+         start=m+1,end=m+n, 
+         method="BFGS") # leave out `$par` 
      p1  p2   p3  value fevals gevals niter convcode kkt1 kkt2 xtimes 
BFGS -1.0325 0.2978319 0.04921863 0.09326904 102 100 NA  1 TRUE FALSE 3.366 

如果您不同意的决定，允许默认坡度估计，HTEN您需要将您的调试缩小到引发错误的功能：

Error: Gradient function might be wrong - check it! 
> traceback() 
3: stop("Gradient function might be wrong - check it! \n", call. = FALSE) 
2: optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, 
     upper, hessian, optcfg$ctrl, have.bounds = optcfg$have.bounds, 
     usenumDeriv = optcfg$usenumDeriv, ...) 
1: optimx(initialbeta.esym[i, ], ERsum, gr = ERsum.gr, tau = tau.found.sym.optim, 
     start = m + 1, end = m + n, method = "BFGS") 

然后查看文档（没有帮助页面）和代码optimx:::optimx.check。这是检查的代码部分：

if (!is.null(ugr) && !usenumDeriv) { 
     gname <- deparse(substitute(ugr)) 
     if (ctrl$trace > 0) 
      cat("Analytic gradient from function ", gname, 
       "\n\n") 
     fval <- ufn(par, ...) 
     gn <- grad(func = ufn, x = par, ...) 
     ga <- ugr(par, ...) 
     teps <- (.Machine$double.eps)^(1/3) 
     if (max(abs(gn - ga))/(1 + abs(fval)) >= teps) { 
      stop("Gradient function might be wrong - check it! \n", 
       call. = FALSE) 
      optchk$grbad <- TRUE 
     }