Error with calculation of hessian, NaN detected

[Gabriele Fadini]

Hello,

First of all thanks to the developers for your work in the creation of this very powerful optimization tool, it's proving a very interesting library in the problems I'm trying to tackle. I really appreciate its versatility.

I'm using casadi for a problem in  which an external library computes the symbolic update in the state given by the x_next = Phi(x, u) [SX, SX] -> [SX], where Phi is of casadi.Function type.
The notation is the following:
- x is the complete state including positions q and velocities v
- u is the command torque at the joints (nu = 2)
Additionally I'm using the casadi.Opti class to express/solve the OCP

I've already tested some OCP on simple cases but I came up with a problem when trying to solve a problem for a floating base system.
The state includes a quaternion that ideally can be properly integrated by such external library so that the transition between one state and the following is well determined (through the transition law Phi), but it makes the position and the velocity vectors have different dimension (nq = 9, nv = 8) (so internally Phi is already taking care of the proper integration of the quaternion and the unity constraint).

This said, when trying to launch the solver I have the following error which I'm not sure on how to interpret but I suspect is related to the introduction of the quaternion, since the entry [3,3] is indeed connected to the position of the first element of the quaternion.

The error is the following:
[...]
CasADi - 2021-02-11 15:44:47 WARNING("solver:nlp_hess_l failed: NaN detected for output hess_gamma_x_x, at nonzero index 9 (row 3, col 3).") [.../casadi/core/oracle_function.cpp:265]
Number of Iterations....: 0
[ ...]
RuntimeError: Error in Opti::solve [OptiNode] at .../casadi/core/optistack.cpp:167:
.../casadi/core/optistack_internal.cpp:999: Assertion "return_success(accept_limit)" failed:
Solver failed. You may use opti.debug.value to investigate the latest values of variables. return_status is 'Invalid_Number_Detected'

I'd very very grateful if someone could explain me what the function hess_gamma_x_x is about and if you have some suggestion on how to debug/solve this problem

[gabf...@gmail.com]

EDIT:

trying with the IPOPT option 'hessian_approximation':'limited-memory' makes the problem disappear and correctly launches the solver.
So I guess that the problem is connected to the exact computation of the hessian, the differentiation there seems to fail for some reason.

I've inspected the derivatives of the jacobian and it seems that if I compute it with the casadi.jacobian method it works as expected, while with the iteration of casadi.gradient I can observe some nans in the quaterion block.
I fear that it possible that the same thing repeats also for the hessian, so I'm checking that as well

Can someone of the dev team help me in debugging this problem?
Where is the check of the hess_gamma_x_x launched? What part of the problem is it checking for consistency: has it to do with the dynamics constraints or cost function?

[Joel Andersson]

Hi,

"hess_gamma_x_x, at nonzero index 9 (row 3, col 3)"

Means that the Hessian of the Lagrangian function evaluates to NaN at that particular entry (row 3, column 3). You want to check how your third decision variable enters into the NLP objective and constraints.

Joel

[gabf...@gmail.com]

Hi Joel,

Thanks for your reply! So if I understand correctly it should be possible to obtain such matrix with:
the evaluation of casadi.hessian(opti.f, opti.x) for the init value of x, is this correct?
Maybe it is also accessible directly with opti.debug?

I have also an additional question: performing some tests I found that in the case of the dynamics that I consider, computing the jacobian of a SX vector or iterating the gradient for each component don't yield the same result and when using the gradient some NaNs appear.
What would be the difference in these two approaches?

Thank you again,
Gabriele

[Joel Andersson]

It's not the Hessian of the objective, it's the Hessian of the Lagrangian function. So you have to include a linear combination of the constraints, weighted by the corresponding multipliers. If you have the particular IPOPT instance, you should be able to extract the Hessian function using the command.

hfcn = solver.get_function('nlp_hess_l')

If you're using the Opti abstraction, you probably need to replace "solver" with "opti.debug().casadi_solver().

Not sure what you mean with your second question. But you should use a different post for it, to keep with one topic per thread.

Joel

[gabf...@gmail.com]

Oh, I see my error, my bad for that I misread your message.
Yes you're right, the way to recover you suggested is correct, thanks for pointing that out Joel.
I will open a new issue for the second question then, like you suggested.

Gabriele

[sep]

Hi, 

I have the same issue, ("WARNING("solver:nlp_hess_l failed: NaN detected for output hess_gamma_x_x, at nonzero index 728 (row 96, col 96).") [.../casadi/core/oracle_function.cpp:265]" 

But I am encountering something super odd whilst debugging: If I create the hessian of the Lagrangian function manually (i.e. via "cs.hessian(opti.f+cs.dot(opti.lam_g,opti.g),opti.x)[0]"), then evaluate the value of the hessian (using "opti.debug.value(...)"), there are no NaNs. 

However, if I check the value of the hessian calculated by the 'nlp_hess_l' function, i.e.

"hess_gamma_x_x_func = opti.debug.casadi_solver.get_function('nlp_hess_l')

hess_gamma_x_x_SX = hess_gamma_x_x_func(x=opti.x, p=opti.p, lam_f=1, lam_g=opti.lam_g)['hess_gamma_x_x']

debug_hess_gamma_x_x = opti.debug.value(hess_gamma_x_x_SX).toarray()", 

then there appear to be NaNs. 

I am stumped as to why the hessian I calculate manually is fine, but the one obtained from the function in the solver object is not. In particular, I cannot understand what "gamma" refers to.  

In my Opti stack, I do expand MX to SX, could this be related to my issue? I would greatly appreciate any guidance.

Regards, 

Sepehr 

[Joel Andersson]

Gamma here refers to a linear combination of the objective and the constraints. That is gamma := lam_f*f + dot(lam_g, g)

It might indeed be an SX vs. MX issue. It's possible that the SX version simplifies the code in a way that eliminates the NaNs.

For more question, I suggest using a separate post, since this is deviating from the original topic.

Joel

[sep]

thanks Joel! No more questions, but for others who stumble onto this thread, I was able to solve the issue by overriding the automatically generated nlp_hess_l: 

""""

lag = opti.f+cs.dot(opti.lam_g,opti.g)
temp_lam_f = cs.MX.sym("temp_lam_f")

H_lag = cs.triu(cs.hessian(lag, opti.x)[0], True)
hessLag = cs.Function('nlp_hess_l',{'x':opti.x, 'p': opti.p, 'lam_f':temp_lam_f, 'lam_g': opti.lam_g, 'hess_gamma_x_x':temp_lam_f*H_lag},
             ['x','p','lam_f','lam_g'], ['hess_gamma_x_x']).expand() 

...

opts["hess_lag"] = hessLag

opti.solver('ipopt', opts)

"""