# The problem in the hessian matrix in optimizing function

**URL:** <https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407>\
**Category:** General\
**Created:** [November 2, 2017, 1:21pm UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407 "2017-11-02T13:21:50Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![bistoon](https://avatars.discourse-cdn.com/v4/letter/b/6f9a4e/32.png) [@bistoon](https://discourse.mc-stan.org/u/bistoon)\
**Post date:** [November 2, 2017, 1:21pm UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/1 "2017-11-02T13:21:50Z")

</div>

[sd\_in\_optimizing\_3.r](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/69eebd21218df4c5370a923a0a5cce5446ed513f.r) (1.5 KB)

Dear all,  
Hi

I think that the optimizing function in Stan is fast and better than optim function in R, but the hessian matrix results (and, of course, the standard deviation (sd) of the MLEs) are different.  
If X~N(mu, sigma), we can obtain the exact value of MLEs and their sds. By running the attached file, we can see that optim and sampling (in Stan) functions have good sd, but the optimizing function results is unacceptable (of course, only for the sd of the sigma parameter).  
Thank you for checking the file and specifying my mistake  
.  
Best,  
-Bistoon

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [November 2, 2017, 2:30pm UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/2 "2017-11-02T14:30:24Z")

</div>

```
set.seed(100)

n = 2000

x = rnorm(n, 0, 2)
loglike = function(x, par) -sum(log(dnorm(x,par[1],par[2])))

fit1 = optim(par=c(mean(x), sd(x)), loglike, x=x , hessian = TRUE)
fit1
$par
[1] 0.02095186 2.01174852

$value
[1] 4235.386

$counts
function gradient 
      35 NA 

$convergence
[1] 0

$message
NULL

$hessian
         [,1] [,2]
[1,] 494.1771 0.0000
[2,] 0.0000 987.6151

fit2 = optimizing(mystanmodel, hessian = TRUE, data=list(n=n, x=x), init = list(mu = mean(x), sigma = sd(x)))
Initial log joint probability = -2397.51
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance
fit2
$par
        mu sigma 
0.02095186 2.01124552 

$value
[1] -2397.508

$return_code
[1] 0

$hessian
                 mu sigma
mu -4.944243e+02 2.247603e-11
sigma 2.247603e-11 -4.000003e+03

```

I’m not a Stan developer. I don’t see very much difference in parameter estimate in mean and sd of fit1 and fit2, the hessian estimate differs in sigma/sigma. Is it that optimization is calculated at the **unconstrained** parameter space?

See: " hessian: ‘TRUE’ or ‘FALSE’ (the default): flag indicating whether to  
calculate the Hessian (via numeric differentiation of the  
gradient function in the unconstrained parameter space)."

---

<div class="post-metadata">

**Author:** ![bistoon](https://avatars.discourse-cdn.com/v4/letter/b/6f9a4e/32.png) [@bistoon](https://discourse.mc-stan.org/u/bistoon)\
**Post date:** [November 2, 2017, 4:43pm UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/3 "2017-11-02T16:43:24Z")

</div>

[sd\_in\_optimizing\_3.r](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/69eebd21218df4c5370a923a0a5cce5446ed513f.r) (1.5 KB)

in optim: fit1$hessian[2,2] = 987.6151  
in optimizing: fit2$hessian[2,2] = -4000

I changed fit1 as follows and attached new file here:

### optim function in R (stats)

loglike = function(x, par) -sum(log(dnorm(x,par[1],par[2])))  
lower = c(-Inf, 0)  
upper = c(Inf, Inf)  
fit1 = optim(par=c(mean(x), sd(x)), loglike, x=x , hessian = TRUE, method = “L-BFGS-B”, lower=lower, upper=upper)

Here, I think fit1 and fit2 have similar parameter space and optimizing method. The optim function seems to be a good result (in sd) and I want to get the same results with the optimizing function, but I have not been able to do this yet.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [November 6, 2017, 12:13am UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/4 "2017-11-06T00:13:14Z")

</div>

> [@Andre\_Pfeuffer](#):
>
> Is it that optimization is calculated at the unconstrained parameter space?

The Hessian’s being computed on the unconstrained scale so that’d be `log(sigma)`.

The accuracy with finite differences of gradients isn’t great, but it should be in the ballpark to several decimal places.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [November 6, 2017, 12:15am UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/5 "2017-11-06T00:15:08Z")

</div>

Here’s the doc for `optimizing` option `hessian`:

```
hessian: ‘TRUE’ or ‘FALSE’ (the default): flag indicating whether to
          calculate the Hessian (via numeric differentiation of the
          gradient function in the unconstrained parameter space).

```

---

<div class="post-metadata">

**Author:** ![Andre\_Pfeuffer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andre_pfeuffer/32/4939_2.png) [@Andre\_Pfeuffer](https://discourse.mc-stan.org/u/Andre_Pfeuffer)\
**Post date:** [November 6, 2017, 3:56am UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/6 "2017-11-06T03:56:18Z")

</div>

Bistoon has a point, the hessian looks not ok.

> in optim: fit1$hessian[2,2] = 987.6151  
> in optimizing: fit2$hessian[2,2] = -4000

Exactly in this entry [2,2] of the matrix and the “end”-element too.  
Bob, could it be some for(0 … (n-1)) or for(1 … n) coping thing?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [November 6, 2017, 5:26am UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/7 "2017-11-06T05:26:36Z")

</div>

No, it is an unconstrained parameterization thing.

---

<div class="post-metadata">

**Author:** ![bistoon](https://avatars.discourse-cdn.com/v4/letter/b/6f9a4e/32.png) [@bistoon](https://discourse.mc-stan.org/u/bistoon)\
**Post date:** [November 6, 2017, 6:29am UTC](https://discourse.mc-stan.org/t/the-problem-in-the-hessian-matrix-in-optimizing-function/2407/8 "2017-11-06T06:29:50Z")

</div>

I’ve taken my answer in this topic (Optimizing function: incorrect estimation of the standard deviation of MLE). This code is useful (thanks Dr. Goodrich):

fit2 = optimizing(Stan\_Normal,hessian=TRUE,data=list(n=n,x=x), init=list(mu=mean(x),sigma=sd(x)),draws=20000)  
optimizing\_mle = fit2$par;  
optimizing\_sd\_mle = sqrt(diag(cov(fit2$theta\_tilde)))  
optimizing\_sd\_mle

-Bistoon