R help - Apr 2022 - hessian in solnp

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 




	[[alternative HTML version deleted]]

? 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

Thanks. Yes, there is a Lagrange multiplier (though really small for the
example). However, AFAIK the hessian is the square matrix of second-order
partial derivatives and, given the simple constraint in the example, the second
(partial) derivatives of the lagrangian function should simplify to the second
(partial) derivatives of the loss function. Or do I miss something? Moreover, I
would expect smaller standard errors of the estimated parameters when additional
constraints are included, and it is not the case in the given example.


De: "Ivan Krylov" <krylov.r00t at gmail.com> 
?: "CAMARDA Carlo Giovanni via R-help" <r-help at r-project.org>
Cc: "Carlo Giovanni Camarda" <carlo-giovanni.camarda at ined.fr>
Envoy?: Jeudi 28 Avril 2022 06:30:00 
Objet: Re: [R] hessian in solnp 

? 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 


	[[alternative HTML version deleted]]