It is fairly easy to use OpenMx to compare numerical and analyticfunction derivatives. This vignette shows how to do it. The main toolthat we are going to use are two custom compute plans calledaPlan andnPlan.
library(OpenMx)## To take full advantage of multiple cores, use:## mxOption(key='Number of Threads', value=parallel::detectCores()) #now## Sys.setenv(OMP_NUM_THREADS=parallel::detectCores()) #before library(OpenMx)aPlan<-mxComputeSequence(list(#analyticmxComputeOnce('fitfunction',c('gradient')),mxComputeReportDeriv()))nPlan<-mxComputeSequence(list(#numericalmxComputeNumericDeriv(analytic =FALSE,hessian=FALSE,checkGradient =FALSE),mxComputeReportDeriv()))Now that we have the plans ready, we can use them to debug afitfunction. Here’s a fitfunction from the test suite that is somewhatcontrived, but can serve our pedagogical needs.
mat1<-mxMatrix("Full",rnorm(1),free=TRUE,nrow=1,ncol=1,labels="m1",name="mat1")obj<-mxAlgebra(-.5* (log(2*pi)+log(2)+ (mat1[1,1])^2/2),name ="obj")grad<-mxAlgebra(-(mat1[1,1])+2,name ="grad",dimnames=list("m1",NULL))mv1<-mxModel("mv1", mat1, obj, grad,mxFitFunctionAlgebra("obj",gradient="grad"))Since we are not very good at calculus, the gradient functioncontains some errors.
nu<-mxRun(mxModel(mv1, nPlan),silent =TRUE)an<-mxRun(mxModel(mv1, aPlan),silent =TRUE)cbind(numerical=nu$output$gradient,analytic=an$output$gradient)## numerical analytic## m1 -0.2401617 1.519677The optimizer is not run so we get the results immediately, even forlarge complex models. The function also does not need to be(approximately) convex. Any function will do.
The numerical approximation can be pretty different from the analyticeven when there are no errors. We can try another point in the parameterspace.
p1<-runif(2,-10,10)mv1<-omxSetParameters(mv1,labels ='m1',values=p1)nu<-mxRun(mxModel(mv1, nPlan),silent =TRUE)an<-mxRun(mxModel(mv1, aPlan),silent =TRUE)cbind(numerical=nu$output$gradient,analytic=an$output$gradient)## numerical analytic## m1 -4.86523 -7.730459The results do not correspond very closely so we look for matherrors. Indeed, there are errors in the gradient function. We replace itwith the correct gradient,
grad<-mxAlgebra(-(mat1[1,1])/2,name ="grad",dimnames=list("m1",NULL))mv2<-mxModel(mv1, grad)Let’s check the correspondence now.
nu<-mxRun(mxModel(mv2, nPlan),silent =TRUE)an<-mxRun(mxModel(mv2, aPlan),silent =TRUE)cbind(numerical=nu$output$gradient,analytic=an$output$gradient)## numerical analytic## m1 -4.86523 -4.86523Wow, looks much better! Still, it is prudent to check at a few morepoints.
mv2<-omxSetParameters(mv2,labels ='m1',values=rnorm(1))nu<-mxRun(mxModel(mv2, nPlan),silent =TRUE)an<-mxRun(mxModel(mv2, aPlan),silent =TRUE)cbind(numerical=nu$output$gradient,analytic=an$output$gradient)## numerical analytic## m1 0.2396829 0.2396829