from pytensor.gradient import jacobian import pytensor.tensor as pt x = pt.vector("x", shape=(3,)) y = pt.outer(x[1:], x[2:]) assert y.type.shape == (2, 1) try: jacobian(y, x) except Exception as exc: print(exc) # jacobian expects a 1 dimensional variable as `expression`. If not use flatten to make it a vector jac_y = jacobian(y.ravel(), x).reshape((*y.shape, *x.shape)) assert jac_y.type.shape == (2, 1, 3)