# Testing Derivatives

#### 2023-04-05

It is fairly easy to use OpenMx to compare numerical and analytic function derivatives. This vignette shows how to do it. The main tool that we are going to use are two custom compute plans called aPlan and nPlan.

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(  #analytic
mxComputeReportDeriv()))

nPlan <- mxComputeSequence(list(  #numerical
mxComputeNumericDeriv(analytic = FALSE, hessian=FALSE, checkGradient = FALSE),
mxComputeReportDeriv()))

Now that we have the plans ready, we can use them to debug a fitfunction. Here’s a fitfunction from the test suite that is somewhat contrived, 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")
mv1 <- mxModel("mv1", mat1, obj, grad,
mxFitFunctionAlgebra("obj", gradient="grad"))

Since we are not very good at calculus, the gradient function contains 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.0826322 2.165264

The optimizer is not run so we get the results immediately, even for large 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 analytic even when there are no errors. We can try another point in the parameter space.

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 -1.066184 -0.1323689

The results do not correspond very closely so we look for math errors. Indeed, there are errors in the gradient function. We replace it with 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 -1.066184 -1.066184

Wow, looks much better! Still, it is prudent to check at a few more points.

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.5334338 -0.5334338