In the commented snippet below, a behavior of the solnp() I cannot fully explain: when a constraint is added, hessian matrix is obviously changing, but in a way I don't understand. I know that the model below could be estimated differently and in a much efficient way. However I wanted to present the simplest example I could think.
In particular, I encountered this issue because I aim to analytically compute variance-covariance matrix in a (non-linear) model with equality constraints. Any practical hint to find a solution to this more general issue would be also welcome.
Thanks, GC
## simulating a simple linear model
n <- 100
x <- sort(runif(n))
X <- cbind(1,x)
k <- ncol(X)
y <- X%*%c(3,4) + rnorm(n, sd=0.5)
## estimating with lm()
fitlm <- lm(y~x)
## using solnp
library(Rsolnp)
## objecting function
fn1 <- function(par){
y.hat <- X%*%par
res <- y-y.hat
RSS <- t(res)%*%res/2
return(c(RSS))
}
## using solnp without any constraint
solUR <- solnp(c(3.2, 3.5), fun = fn1)
## getting the same fitted objects as in lm()
## coefficients
cbind(solUR$pars, fitlm$coef)
## and standard errors by variance-covariance matrix computed
## with hessian internally estimated
(H.UR <- solUR$hessian)
y.UR <- X%*%solUR$pars
RSS.UR <- t(y-y.UR)%*%(y-y.UR)
sigma2.UR <- RSS.UR/(n-k)
V.UR <- c(sigma2.UR)*solve(H.UR)
se.UR <- sqrt(diag(V.UR))
cbind(se.UR,summary(fitlm)$coef[,2])
## adding equality constraint
eqn1 <- function(par) sum(par)
c <- sum(solUR$pars)
## I use the sum of the previously estimated coefficients
solR <- solnp(c(3.2, 3.5), fun = fn1, eqfun = eqn1, eqB = c)
## as expected we get the same coefficients
cbind(solR$pars,fitlm$coef)
## but the hessian is rather different
(H.R <- solR$hessian)
## which leads to completely different (much larger here)
## standard errors
y.R <- X%*%solR$pars
RSS.R <- t(y-y.R)%*%(y-y.R)
sigma2.R <- RSS.R/(n-k)
V.R <- c(sigma2.R)*solve(H.R)
se.R <- sqrt(diag(V.R))
cbind(se.R,summary(fitlm)$coef[,2])
## I noticed that differences in the hessian
## is ~constant over both dimensions
## and it depends (also) on sample size
H.R-H.UR
hessian in solnp
2 messages · CAMARDA Carlo Giovanni, Ivan Krylov
? Thu, 28 Apr 2022 00:16:05 +0200 (CEST) CAMARDA Carlo Giovanni via R-help <r-help at r-project.org> ?????:
when a constraint is added, hessian matrix is obviously changing, but in a way I don't understand.
Isn't it the point of augmented Lagrange multiplier, to solve the constrained optimisation problem by modifying the loss function and optimising the result in an unconstrained manner? Apologies if I misunderstood your question. Starting on <https://github.com/cran/Rsolnp/blob/4b56bb5cd7c5d1096d1ba2f3946df7afa9af4201/R/subnp.R#L282>, we can see how the functions are called: both the loss and the constraint function are concatenated into the `obm` vector (if there's no constraint, the function returns NULL, which is eaten by concatenation), which form the vectors `g` (seems to be the gradient) and `p`, which, in turn, form the matrix `hessv`. My reading of the code could be wrong (and so could be my understanding of augmented Lagrangian methods). Contacting maintainer('Rsolnp') could be an option; maybe there's some documentation for the original MATLAB version of the code at <https://web.stanford.edu/~yyye/Col.html>?
Best regards, Ivan