Hi,
Yes we follow a slightly different semantic.
The usual definition for multivariate functions is given a 1D input and 1D output, you end up with as 2D Jacobian.
What sympy is doing here is to linearize the inputs from a (2, 2) Tensor into a (4,) Tensor. And so the Jacobian is (2, 4).
Because we usually allow multiple Tensors, people very often multiple dimensions and linearizing the Tensors can be confusing if you don’t know how they are layed out in memory, we made a different choice here (similar to what tensorflow or jax are doing).
Here we don’t linearize Tensors and so you can index your Jacobian the same way you would index your original Tensor.