title: Efficient way to compute Jacobian x Jacobian.T
body: <p>Assume <code>J</code> is the Jacobian of some function <code>f</code> with respect to some parameters. Are there efficient ways (in PyTorch or perhaps Jax) to have a function that takes two inputs (<code>x1</code> and <code>x2</code>) and computes <code>J(x1)*J(x2).transpose()</code> <em>without</em> instantiating the entire <code>J</code> matrices in memory?</p>
  <p>I have come across something like <code>jvp(f, input, v=vjp(f, input))</code> but don't quite understand it and not sure is what I want.</p>