I opened an issue then. By the way, indeed double forward mode is exactly what I need. The first pass will compute [dy_1/dx_1 … dy_m/dx_1] and the derivative with respect to x_1 will give the dx_1^2.
So I’ll go with this mode once I figure out how to use that with a ::Chain. Reverse-on-forward would be best If you need the full hessian.