This work presents efficient formulas for the computation of Jacobian matrices arising from entropy stable summation-by-parts schemes. Competing methods for computing Jacobian matrices include finite difference, automatic differentiation, graph coloring, and Jacobian-free Newton-Krylov methods. In contrast to these methods the formulas proposed provide a sparsity-informed method for computing Jacobian matrices that are free of truncation error. Computational timings confirm that the proposed formulas scale very robustly with respect to the size of the system and easily outperform existing methods on denser Jacobian matrices. Two applications of Jacobian matrices, two-derivative and implicit time stepping, are also explored in numerical experiments with Burgers' and the shallow water equations.