Skip to content

hessian in solnp

2 messages · CAMARDA Carlo Giovanni, Ivan Krylov

#
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
#
? Thu, 28 Apr 2022 00:16:05 +0200 (CEST)
CAMARDA Carlo Giovanni via R-help <r-help at r-project.org> ?????:
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>?