title: Computing hessian with pydrake autodiff
body: <p>One of Drake's selling points is the easy availability of gradients via AutoDiff, but I'm struggling to see how to easily compute second-order derivatives in pydrake.</p>
  <p>Given a function f(x), I know of two ways to compute the Jacobian. The first way uses the <code>forwarddiff.jacobian</code> helper function, e.g.:</p>
  <pre><code>def f(x):
      return x.T@x
  
  x = np.array([1,2,3])
  fx = jacobian(f,x)     # = [2,4,6]
  </code></pre>
  <p>The second way uses the <code>autodiffutils</code> bindings more directly:</p>
  <pre><code>x = InitializeAutoDiff([1,2,3])
  y = f(x)
  fx = ExtractGradient(y)   # = [2,4,6]
  </code></pre>
  <p>Are there similar ways to get the Hessian? Nested calls to the <code>jacobian</code> helper function don't work, since the second argument can't be an <code>AutoDiffXd</code> type. But perhaps there is some way more analogous to the second method above?</p>