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 函数，但出现以下错误：

错误：梯度函数可能有误 - 检查一下！

为了检查我的梯度是否正常，我尝试使用 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、partial.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

这可能与 optimx 返回的结构与 optim 返回的结构不同：

> 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

如果您不同意允许默认梯度估计的决定，您需要将调试范围缩小到引发错误的函数：

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
        }