qpqc <- function (yold, a, b, c, itmax = 100, eps = 1e-10, analytic = TRUE, verbose = FALSE) { hfunc <- ifelse (analytic, hfgh, hfghnum) hval <- nnnewton ( yold, hfunc, a = a, b = b, c = c, eps = eps, verbose = verbose, itmax = itmax ) fval <- primal (hval\$solution, a, b, c) return ( list ( h = hval\$value, f = fval\$value, xmin = fval\$solution, multipliers = hval\$solution, constraints = hval\$gradient ) ) } qfunc <- function (x, a, b, c) { return (c + sum (b * x) + sum (x * (a %*% x)) / 2) } hfgh <- function (x, a, b, c) { m <- length (x) n <- nrow (a[, , 1]) asum <- a[, , 1] bsum <- b [, 1] csum <- c[1] for (j in 1:m) { asum <- asum + x[j] * a[, , j + 1] bsum <- bsum + x[j] * b[, j + 1] csum <- csum + x[j] * c[j + 1] } vinv <- solve (asum) xmin <- -drop (vinv %*% bsum) h <- csum + sum (bsum * xmin) / 2 dh <- rep (0, m) dg <- matrix (0, n, m) for (j in 1:m) { dh[j] <- qfunc (xmin, a[, , j + 1], b[, j + 1], c[j + 1]) dg[, j] <- drop (b[, j + 1] + a[, , j + 1] %*% xmin) } d2h <- -crossprod (dg, vinv %*% dg) return (list ( h = h, dh = dh, d2h = d2h )) } hfghnum <- function (x, a, b, c) { hfunc <- function (x) { return (hfgh (x, a, b, c)\$h) } return (list (h = hfunc (x), dh = grad (hfunc, x), d2h = hessian (hfunc, x))) } primal <- function (x, a, b, c) { m <- length (x) asum <- a[, , 1] bsum <- b [, 1] csum <- c[1] for (j in 1:m) { asum <- asum + x[j] * a[, , j + 1] bsum <- bsum + x[j] * b[, j + 1] csum <- csum + x[j] * c[j + 1] } xmin <- -drop (solve (asum, bsum)) fmin <- qfunc (xmin, a[, , 1], b[, 1], c[1]) return (list (solution = xmin, value = fmin)) }