Estimates don't move far enough from poor starts in CSOLNP

I took a look at this today. It is a CSOLNP issue. NPSOL miraculously gets to the correct solution from the awful starting values provided (zeroes for everything basically). CSOLNP does not. Even given halfway decent ones with something like

Note that the correct ML estimates of the variances deviate from those calculated using var() or cov() a non-trivial amount (31/32) due to the small N of 32. This is definitely not the whole issue with CSOLNP.


> s = summary(m1)$parameters; s$Estimate = round(s$Estimate,2); s[,3:5]; diag(cov(mtcars[,manifests]))*31/32; colMeans(mtcars[,manifests])
row col Estimate
1 mpg mpg 35.19
2 disp disp 14880.81
3 gear gear 0.53
4 1 mpg 20.09
5 1 disp 230.72
6 1 gear 3.69
mpg disp gear
3.518897e+01 1.488077e+04 5.273438e-01
mpg disp gear
20.09062 230.72188 3.68750


Even given halfway decent starting values, CSOLNP fails to get var(disp) correct:

mxOption(NULL,"Default optimizer", "CSOLNP")
manifests = c("mpg", "disp", "gear")

# =============================
# = simple independence model =
# =============================
m1 <- mxModel("ind", type = "RAM",
manifestVars = manifests,
mxPath(from = manifests, arrows = 2, values=diag(cov(mtcars[,manifests]))*1.05),
mxPath(from = "one", to = manifests, values=colMeans(mtcars[,manifests])*1.1),
mxData(mtcars[,manifests], type="raw")
)

m1 = mxRun(m1)

# ===========================
# = Our parameter estiamtes =
# ===========================
> s = summary(m1)$parameters; s$Estimate = round(s$Estimate,2); s[,3:5];
row col Estimate
1 mpg mpg 35.19
2 disp disp 16128.81
3 gear gear 0.53
4 1 mpg 20.09
5 1 disp 230.72
6 1 gear 3.69


Although the estimates are off, the -2lnL is only off by .2 or so, this is likely due to flattish gradients on account of the lack of data. But it should do better.

Useful test case, should probably go in models/failing in a CSOLNPmtcars.R


mxOption(NULL,"Default optimizer", "CSOLNP")
manifests = c("mpg", "disp", "gear")

# =============================
# = simple independence model =
# =============================
m1 <- mxModel("ind", type = "RAM",
manifestVars = manifests,
mxPath(from = manifests, arrows = 2, values=diag(cov(mtcars[,manifests]))*1.05),
mxPath(from = "one", to = manifests, values=colMeans(mtcars[,manifests])*1.1),
mxData(mtcars[,manifests], type="raw")
)

m1 = mxRun(m1)

# ===========================
# = Our parameter estiamtes =
# ===========================
s = summary(m1)$parameters; s$Estimate = round(s$Estimate,2); s[,3:5]; summary(m1); diag(cov(mtcars[,manifests]))*31/32; colMeans(mtcars[,manifests])


This is not a raw data issue. The problem is the same with covariance data.


# =============================
m1 <- mxModel("ind", type = "RAM",
manifestVars = manifests,
mxPath(from = manifests, arrows = 2),
mxPath(from = "one", to = manifests),
mxData(cov(mtcars[,manifests]), means=colMeans(mtcars[,manifests]), numObs=nrow(mtcars), type="cov")
)

m1 = mxRun(m1)
#Running ind
s = summary(m1)$parameters; s$Estimate = round(s$Estimate,2); s[,3:5]
# row col Estimate
#1 mpg mpg 695012.58
#2 disp disp 108907335.59
#3 gear gear 0.54
#4 1 mpg 20.09
#5 1 disp 230.59
#6 1 gear 3.69


Scaling is almost certainly the major problem. When free parameters span five orders of magnitude (1e-1 to 1e+4), only the most robust optimizers will be able to handle that on their own without blinking. It really screws with the magnitude of the gradient because you have to make tiny steps in some directions (for the 1e-1 parameters) and giant steps in other directions (for the 1e+4 parameters). Physicists suggest that you re-parameterize the problem to avoid parameter estimation at very different orders of magnitude.

Another part of the problem is definitely the starting values. Zero is just about the worst starting value for a variance. With that being said, CSOLNP does wiggle off from zero starting values. It just doesn't wiggle far enough.

Re-scaling the data (ala physics) and letting CSOLNP wiggle the starting values off zero resolves the behavioral issue.


library(OpenMx)
manifests = c("mpg", "disp", "gear")

mtcars$mpg <- mtcars$mpg/10
mtcars$disp <- mtcars$disp/300
mtcars$gear <- mtcars$gear

# =============================
# = simple independence model =
# =============================
m1 <- mxModel("ind", type = "RAM",
manifestVars = manifests,
mxPath(from = manifests, arrows = 2, values=0),
mxPath(from = "one", to = manifests),
mxData(mtcars[,manifests], type="raw")
)

m1 = mxRun(m1)

# ===========================
# = Our parameter estiamtes =
# ===========================
s = summary(m1)$parameters; s$Estimate = round(s$Estimate,2); s[,3:5]
# row col Estimate
#1 mpg mpg 35.19
#2 disp disp 16.53
#3 gear gear 52.73
#4 1 mpg 20.09
#5 1 disp 7.69
#6 1 gear 36.87

round(diag(cov(mtcars[,manifests])),2)
# mpg disp gear
#36.32 17.07 54.44


Of course, we would probably want CSOLNP to more robustly find the correct solution. Does anyone know if NPSOL rescales the gradient similar to CSOLNP?

Oh, we use our own omx routine for that, which reduces the number of required libraries.
E.g.,

omxCheckCloseEnough(growthCurveFit$output$estimate[["meani"]], 9.930, 0.01)

so a translation would be:

omxCheckCloseEnough(obtainedVars[1], mlexpected[1], epsilon = 0.01)
omxCheckCloseEnough(obtainedVars[2], mlexpected[2], epsilon = 0.05)
omxCheckCloseEnough(obtainedVars[3], mlexpected[3], epsilon = 0.01)

This passes with NPSOL but not with CSOLNP.

However, omxCheckCloseEnough() only works on absolute differences, whereas testthat() can do relative as well. I like the flexibility of testthat() though -- not sure if it's worth changing over to using it more broadly, a question for the developers generally.

OpenMx has omxCheckCloseEnough() for absolute differences and omxCheckWithinPercentError() for relative (percent) differences. These functions are defined in trunk/R/MxUnitTesting.R.

OpenMx has


omxCheckIdentical # uses R's identical() function
omxCheckEquals # uses R's == function
omxCheckSetEquals # checks for equality of sets (order independent)
omxCheckTrue # similar to boolean test
omxCheckCloseEnough # approximate equality on absolute scale
omxCheckWithinPercentError # approximate equality on relative scale
omxCheckWarning # verify existence of warning and the right message is printed
omxCheckError # verify existence of error and the right message is printed


As far as I can tell we've had all of these functions for 3+ years.