Some recent MCMC samplers such as "Randomize-then-optimize" exploit least squares structure in the log-likelihood; i.e., the log-density has the form $f(x) = \sum_j f_j(x)^2$. Could we add this kind of derivative information into the API? As a proposal, I would suggest something like
residual(ℓ, x)
jacobian(ℓ, x)
which is similar to the NLSolversBase.jl API. Non-allocating versions would be helpful as well, e.g., writing into the last argument
residual!(ℓ, x, res)
jacobian!(ℓ, x, jac)