I’m happy to report Bob and I fixed the issue.
What caused the bug: you cannot call Jacobian() inside a chain method. One run of Auto-Diff interferes with the other. In particular, once memory is freed by the Jacobian, it messes with the memory allocated for chain() or whatever gradient computation may be going on.
As a reminder, this is problematic, because we calculate the adjoint of the parameters y for the algebraic solver using Jacobians (as prescribed by the “implicit function theorem” technique, see Implicit Function Theorem Math).
The solution is to compute these Jacobians before the chain method gets called. As it is, we do this in the constructor of the vari. Jx_y (jacobian of the solution with respect to the auxiliary variables) is stored as a member of vari, and can then be used inside the chain method.
We’ll probably need to address the special case where the user passes a null vector for y (i.e. no parameters involved, all fixed data).
For completion, I’ll mention we revised the structure of the algebra_solver_vari structure. I initially followed the example of quad_form, but apparently there is some fancy memory work going on there we’ll want to revise eventually.