docs: def make_ggnvp(f, g=lambda x: 1./2*np.sum(x**2, axis=-1), f_argnum=0):
      """Builds a function for evaluating generalized-Gauss-Newton-vector products
      at a point. Slightly more expensive than mixed-mode."""
      @unary_to_nary
      def _make_ggnvp(f, x):
          f_vjp, f_x = _make_vjp(f, x)
          g_hvp, grad_g_x = _make_vjp(grad(g), f_x)
          f_jvp, _ = _make_vjp(f_vjp, vspace(grad_g_x).zeros())
          def ggnvp(v): return f_vjp(g_hvp(f_jvp(v)))
          return ggnvp
      return _make_ggnvp(f, f_argnum)