code: def _bfgs_direction(s, y, x, hessinv_estimate=None):
      r
      assert len(s) == len(y)
  
      r = x.copy()
      alphas = np.zeros(len(s))
      rhos = np.zeros(len(s))
  
      for i in reversed(range(len(s))):
          rhos[i] = 1.0 / y[i].inner(s[i])
          alphas[i] = rhos[i] * (s[i].inner(r))
          r.lincomb(1, r, -alphas[i], y[i])
  
      if hessinv_estimate is not None:
          r = hessinv_estimate(r)
  
      for i in range(len(s)):
          beta = rhos[i] * (y[i].inner(r))
          r.lincomb(1, r, alphas[i] - beta, s[i])
  
      return r
docs: r"""Compute ``Hn^-1(x)`` for the L-BFGS method.
  
      Parameters
      ----------
      s : sequence of `LinearSpaceElement`
          The ``s`` coefficients in the BFGS update, see Notes.
      y : sequence of `LinearSpaceElement`
          The ``y`` coefficients in the BFGS update, see Notes.
      x : `LinearSpaceElement`
          Point in which to evaluate the product.
      hessinv_estimate : `Operator`, optional
          Initial estimate of the hessian ``H0^-1``.
  
      Returns
      -------
      r : ``x.space`` element
          The result of ``Hn^-1(x)``.
  
      Notes
      -----
      :math:`H_n^{-1}` is defined recursively as
  
      .. math::
          H_{n+1}^{-1} =
          \left(I - \frac{ s_n y_n^T}{y_n^T s_n} \right)
          H_{n}^{-1}
          \left(I - \frac{ y_n s_n^T}{y_n^T s_n} \right) +
          \frac{s_n s_n^T}{y_n^T \, s_n}
  
      With :math:`H_0^{-1}` given by ``hess_estimate``.