R/lav_model_hessian.R

In lavaan: Latent Variable Analysis

# numeric approximation of the Hessian
# using an analytic gradient
lav_model_hessian <- function(lavmodel       = NULL,
                              lavsamplestats = NULL,
                              lavdata        = NULL,
                              lavoptions     = NULL,
                              lavcache       = NULL,
                              group.weight   = TRUE,
                              ceq.simple     = FALSE,
                              h              = 1e-06) {

    estimator <- lavmodel@estimator

    # catch numerical gradient
    if(lavoptions$optim.gradient == "numerical") {
        obj.f <- function(x) {
           lavmodel2 <- lav_model_set_parameters(lavmodel, x = x)
           lav_model_objective(lavmodel = lavmodel2,
               lavsamplestats = lavsamplestats, lavdata = lavdata)[1]
        }
        x <- lav_model_get_parameters(lavmodel = lavmodel)
        Hessian <- numDeriv::hessian(func = obj.f, x = x)
        return(Hessian)
    }

    # computing the Richardson extrapolation
    if(!ceq.simple && .hasSlot(lavmodel, "ceq.simple.only") &&
       lavmodel@ceq.simple.only) {
        npar <- lavmodel@nx.unco
        type.glist <- "unco"
    } else {
        npar <- lavmodel@nx.free
        type.glist <- "free"
    }
    Hessian <- matrix(0, npar, npar)
    x <- lav_model_get_parameters(lavmodel = lavmodel)
    if(!ceq.simple && .hasSlot(lavmodel, "ceq.simple.only") &&
       lavmodel@ceq.simple.only) {
        # unpack
        x <- drop(x %*% t(lavmodel@ceq.simple.K))
    }
    for(j in seq_len(npar)) {
        # FIXME: the number below should vary as a function of 'x[j]'
        h.j <- h
        x.left <- x.left2 <- x.right <- x.right2 <- x
        x.left[j]  <- x[j] - h.j; x.left2[j]  <- x[j] - 2*h.j
        x.right[j] <- x[j] + h.j; x.right2[j] <- x[j] + 2*h.j

        g.left <-
            lav_model_gradient(lavmodel       = lavmodel,
                               GLIST          = lav_model_x2GLIST(lavmodel =
                                                    lavmodel, type = type.glist,
                                                    x.left),
                               lavsamplestats = lavsamplestats,
                               lavdata        = lavdata,
                               lavcache       = lavcache,
                               type           = "free",
                               group.weight   = group.weight,
                               ceq.simple     = ceq.simple)
        g.left2 <-
            lav_model_gradient(lavmodel       = lavmodel,
                               GLIST          = lav_model_x2GLIST(lavmodel =
                                                    lavmodel, type = type.glist,
                                                    x.left2),
                               lavsamplestats = lavsamplestats,
                               lavdata        = lavdata,
                               lavcache       = lavcache,
                               type           = "free",
                               group.weight   = group.weight,
                               ceq.simple     = ceq.simple)

        g.right <-
            lav_model_gradient(lavmodel       = lavmodel,
                               GLIST          = lav_model_x2GLIST(lavmodel =
                                                    lavmodel, type = type.glist,
                                                    x.right),
                               lavsamplestats = lavsamplestats,
                               lavdata        = lavdata,
                               lavcache       = lavcache,
                               type           = "free",
                               group.weight   = group.weight,
                               ceq.simple     = ceq.simple)

        g.right2 <-
            lav_model_gradient(lavmodel       = lavmodel,
                               GLIST          = lav_model_x2GLIST(lavmodel =
                                                    lavmodel, type = type.glist,
                                                    x.right2),
                               lavsamplestats = lavsamplestats,
                               lavdata        = lavdata,
                               lavcache       = lavcache,
                               type           = "free",
                               group.weight   = group.weight,
                               ceq.simple     = ceq.simple)

        Hessian[,j] <- (g.left2 - 8*g.left + 8*g.right - g.right2)/(12*h.j)
    }

    # check if Hessian is (almost) symmetric, as it should be
    max.diff <- max(abs(Hessian - t(Hessian)))
    if(max.diff > 1e-05 * max(diag(Hessian))) {
        # hm, Hessian is not symmetric -> WARNING!
        warning("lavaan WARNING: Hessian is not fully symmetric.",
                "\n\tMax diff = ", max.diff,
                "\n\t(Max diag Hessian = ", max(diag(Hessian)), ")")
        # FIXME: use numDeriv::hessian instead?
    }
    Hessian <- ( Hessian + t(Hessian) )/2.0

    Hessian
}

# if only chol would accept a complex matrix...
lav_model_hessian_complex <- function(lavmodel       = NULL,
                                      lavsamplestats = NULL,
                                      lavdata        = NULL,
                                      lavcache       = NULL,
                                      group.weight   = TRUE) {

    gradf <- function(x) {
        GLIST <-  lav_model_x2GLIST(lavmodel = lavmodel, x = x)
        dx <- lav_model_gradient(lavmodel       = lavmodel,
                                 GLIST          = GLIST,
                                 lavsamplestats = lavsamplestats,
                                 lavdata        = lavdata,
                                 lavcache       = lavcache,
                                 type           = "free",
                                 group.weight   = group.weight)
        dx
    }

    x <- lav_model_get_parameters(lavmodel = lavmodel)
    Hessian <- lav_func_jacobian_complex(func = gradf, x = x)

    Hessian
}