Inspecting NLP constraint Jacobian in C++

[Dave Greer]

Hi,

I'm hitting a problem with detection of a NaN or Inf in my constraint Jacobian when trying to solve an NLP with IPOPT.  I'm using the C++ interface rather than python.  The constraints are manually defined, but modified from code which previously worked with IPOPT and ADOL-C.  I'd like to inspect the Jacobian so I can figure out which constraints are causing the problem.  I figured out the following code

FX J = solver.jacG();

J.init();

J.setInput(x_init);

J.evaluate();

cout << J.output() << endl;

to output the Jacobian but it's printing a dense representation.  My problem is fairly large (20000 constraints) and redirecting the output gives me a file >1GB in size.  I tried grepping for NaN and Inf but didn't find anything that way and it's hard to do any more detailed inspection with such a large output.  Is there a way to either generate a sparse output, or iterate through the non-zero elements in code?

Best wishes,

Dave

[Joel Andersson]

Hello again Dave!

Sure. "J.output()" is a sparse matrix and if you call "J.output().data()" you get a reference to a std::vector<double> containing the nonzeros of the matrix. We have a long-standing issue of implementing iterators for the this matrix class (CasADi::Matrix<>), see issue https://github.com/casadi/casadi/issues/57. At the moment what you can do is to access the nonzero entries and the sparsity pattern (which is currently compressed row storage, but will change to compressed column storage in the next release).

Also note that the Jacobian of the constraints that you (correctly) queried above with the command "FX J = solver.jacG();" is something that you can easily generate yourself from the callback function you passed to the IpoptSolver constructor : "FX J = nlp.jacobian(NL_X,NL_G);". If your "nlp" data structure is an instance of the MXFunction of SXFunction classes, you can also generate symbolic expressions for the for Jacobian (by calling "nlp.jac(NL_X,NL_G)", which might help to analyze what goes wrong.

Good luck!

Joel

[Dave Greer]

Hi Joel,

thanks a lot for responding so quickly!  

The problem turned out to be a dopy mistake on my part (some parameters of my model are loaded from disk and when I moved the code to Cygwin I didn't test they were working ok with the different file system, which led to some divisions by zero).  I got the code working with CasADi and got a speedup of about 100x in the function evaluations over my previous ADOL-C implementation :)

One question - the number of nonzero elements being reported seems to be about 10% lower with CasADi than ADOL-C.  Is this anything to worry about or does it just mean that CasADi is cleverer / less conservative at identifying zero elements?

Thanks again for all your advice, you have been really helpful.

Dave

[Joel Andersson]

Hello Dave!

Great to hear that it is working well!

I do not know exactly how ADOL-C/ColPack identifies nonzeros. What I can imagine is that if you have an expression such as :"f(x) = x-x", CasADi would simplify the expression to "0", whereas ADOL-C/ColPack might not (or might not in a more advanced cases) and conclude that f(x) depends on x. But that is just my speculation.

If you are worried that the output is not correct, you can of course check that the corresponding values of the Jacobian from ADOL-C are numerically zero.

Joel

P.S: If speed is an issue, generating C-code usually gives a speedup of around 5 times over the virtual machine. See the examples/cplusplus/nlp_codegen.cpp for an example.

[Dave Greer]

Hi Joel,

I had guessed that it might be something along those lines.  The output looks sensible so I think there must just have been some non-trivial zeros which ADOL-C was less rigorous at finding.

I was already aware of the code generation option, but the cost of the function evaluations which was previously dominating the run time is now negligible compared to the linear solver in IPOPT, so I think there are other things I can try first which will have more impact, such as reducing iteration count in the solver by improving my problem scaling, initial guess etc.

Thanks,

Dave