Incorrect output when working with estimating normal random variables.

[SlothDragon]

Hello again, I know this is so soon after my last post, but, now, I'm facing yet another issue. It isn't with the code itself but rather the output. From my previous post, I am trying to find the parameters for mu, the mean, and sigma, the standard deviation, using the log-likelihood to find whether those parameters are the correct ones. Since the random variable I am working with is normal, therefore

$$f(x)=-\sum_{i=0}^{N-1}log(\frac{1}{\sqrt{2 \pi {\sigma}^2 }}e^{\frac{(Z(i)-\mu)^2}{2 \sigma^2}})=-\sum_{i=0}^{N-1}\frac{1}{2}log({2 \pi {\sigma}^2) }+\frac{{(Z(i)-\mu)}^2}{2\sigma^2}$$

This is because we are trying to maximize the probability of each Z(i) entry. One can compute the first and second derivatives as

$$\frac{\partial f}{\partial \sigma}=\sum_{i=0}^{N-1}\frac{{(Z(i)-\mu)}^2}{\sigma^3}-\frac{1}{\sigma}$$

$$\frac{\partial f}{\partial \mu}=\sum_{i=0}^{N-1}\frac{{(Z(i)-\mu)}}{\sigma^2}$$

$$\frac{\partial^2 f}{\partial \sigma^2}=\sum_{i=0}^{N-1}-3(\frac{{(Z(i)-\mu)}^2}{\sigma^4})+\frac{1}{\sigma^2}$$

$$\frac{\partial^2 f}{\partial \mu \partial \sigma}=\sum_{i=0}^{N-1}-2(\frac{{(Z(i)-\mu)}}{\sigma^3})$$

$$\frac{\partial^2 f}{\partial \mu^2}=\sum_{i=0}^{N-1}\frac{-1}{\sigma^2}$$

.I have implemented all of the following equations for the gradient of f and the Hessian functions; however, not only do the results remain the same after every execution, which should not occur due to the random device I am using(seen in function.cpp), but the result is nowhere close to the answer. Possible solutions I have tried to come up with rely on changing the mu_strategy and adaptive settings. However, even after I change it from adaptive to monotone, I still get an incorrect answer. Even when I initialize it with the answer as the starting points, it diverges from this answer when adaptive is chosen. On the other hand, if I were to change it to monotone, it doesn't even reach the answer and stays relatively close to the initialized points. To add one more point, it seems that in all cases, my first variable is not moved at all, regardless of what setting I choose.

When monotone is selected and the initial points are the answer of $\mu =5$ and $\sigma =2.5$, it returns $\mu =5$ and $\sigma =2.62525$. If, however, $\mu =0$ and $\sigma =1$, $\mu =0$ and $\sigma =1.03499$. If I instead choose it be adaptive instead with $\mu =5$ and $\sigma =2.5$, $\mu =5$ and $\sigma =23.9283$ is returned. Have I made some sort of error in the way I have set up my optimization?

OU_optim_cpp.txt
OU_process_cpp.txt
function_cpp.txt
function_hpp.txt
Makefile.txt
OU_optim_hpp.txt

[SlothDragon]

Even when I input Z=[5,5,5...,5,5] it does not even update mu variable, the first variable, and the standard deviation, the second variable, stays high when it should be 0.

[svigerske]

In your eval_f, eval_grad_f, and eval_h, you have Eigen::MatrixXd Z;. This overwrites the global variable Z, and the function you optimize becomes 0.

If you remove these, it becomes more interesting:

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.0127506e+04 0.00e+00 1.01e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  3.5392535e+04 0.00e+00 1.37e+00  -1.0 5.00e+01    -  1.94e-02 2.50e-01f  3
   2  3.5270879e+04 0.00e+00 1.39e+00  -1.0 1.69e-01    -  1.00e+00 1.00e+00f  1
   3  3.5146324e+04 0.00e+00 1.40e+00  -1.0 1.71e-01    -  1.00e+00 1.00e+00f  1
   4  3.5018735e+04 0.00e+00 1.42e+00  -1.0 1.73e-01    -  1.00e+00 1.00e+00f  1
   5  3.3909657e+04 0.00e+00 1.58e+00  -1.0 1.42e+00    -  1.00e+00 1.00e+00f  1
   6  3.2523171e+04 0.00e+00 1.80e+00  -1.0 1.58e+00    -  1.00e+00 1.00e+00f  1
   7  3.0693730e+04 0.00e+00 2.12e+00  -1.0 1.80e+00    -  1.00e+00 1.00e+00f  1
   8  2.8083247e+04 0.00e+00 2.62e+00  -1.0 2.12e+00    -  1.00e+00 1.00e+00f  1
   9  2.4200875e+04 0.00e+00 2.68e+00  -1.0 2.61e+00    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.3543825e+04 0.00e+00 3.15e+00  -1.0 2.66e+00    -  1.00e+00 5.00e-01f  2
  11  2.3527839e+04 0.00e+00 1.80e+00  -1.0 3.09e+00    -  6.90e-01 2.50e-01f  3
  12  2.3333226e+04 0.00e+00 7.78e-01  -1.0 2.78e-01    -  1.00e+00 1.00e+00f  1
  13  2.3309218e+04 0.00e+00 5.09e-01  -1.7 2.25e-01    -  1.00e+00 1.00e+00f  1
  14  2.3300113e+04 0.00e+00 7.10e-02  -1.7 8.71e-02    -  1.00e+00 1.00e+00f  1
  15  2.3299851e+04 0.00e+00 5.78e-03  -2.5 1.17e-02    -  1.00e+00 1.00e+00f  1
  16  2.3299849e+04 0.00e+00 8.94e-05  -3.8 1.12e-03    -  1.00e+00 1.00e+00f  1
  17  2.3299849e+04 0.00e+00 3.11e-05  -5.7 9.92e-06    -  1.00e+00 1.00e+00f  1
  18  2.3299849e+04 0.00e+00 8.20e-07  -5.7 9.86e-06    -  1.00e+00 1.00e+00f  1
  19  2.3299849e+04 0.00e+00 7.37e-10  -7.8 2.50e-07    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 19

                                   (scaled)                 (unscaled)
Objective...............:   4.4909906661436580e+01    2.3299848989762690e+04
Dual infeasibility......:   7.3710322821838571e-10    3.8241882880826270e-07
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.7522695211499044e-08    9.0910042498962874e-06
Overall NLP error.......:   1.7522695211499044e-08    9.0910042498962874e-06


Number of objective function evaluations             = 37
Number of objective gradient evaluations             = 20
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 0
Total seconds in IPOPT                               = 0.419

EXIT: Optimal Solution Found.


Solution of the primal variables, x
x[0] = 5.01436
x[1] = 2.48755

This is with Hessian approximation, because there seems to be a mistake in you Hessian implementation:

Starting derivative checker for second derivatives with obj_factor or lambda[i] set to 1.5.

*             obj_hess[    0,    0] =  4.0009331360077982e+02 v  ~  6.0014002200027755e+02  [ 3.333e-01]
*             obj_hess[    0,    1] =  1.4701659452422311e+03 v  ~  2.2052486474905518e+03  [ 3.333e-01]
*             obj_hess[    1,    0] =  1.4701659452422311e+03 v  ~  2.2052496316132783e+03  [ 3.333e-01]
*             obj_hess[    1,    1] = -7.6467101639868133e+03 v  ~  5.9225746215838053e+03  [ 2.291e+00]

[SlothDragon]

Thank you once again. Overwriting the variable Z was the issue, and changing it solved this. How do I begin fixing my Hessian implementation? Is it a problem with the code I wrote, since it uses a for loop, or is it due to other settings I would need to change? After some testing, I found that when I reduce the size of Z, no issues appear during my optimization; however, when it exceeds 1000 in length, the optimization fails. I am guessing this is where my improper Hessian implementation comes into play.

[svigerske]

Turn on the Hessian check in Ipopt (derivative_test second-order) and check your formulas. The errors reported by Ipopt should be small or vanish.
One thing is that you missed multiplying by obj_factor, but also after fixing that, there seems to be still an issue with $$\frac{\partial^2 f}{\partial \mu^2}$$:

*             obj_hess[    1,    1] = -1.1397872195406137e+04 v  ~  5.8792543659224957e+03  [ 2.939e+00]

[SlothDragon]

Thank you. I have solved this issue. The issue I had with $\frac{\partial^2 f}{\partial \mu^2}$ was that I forgot a negative sign when using the chain rule in the implementation.