code: def _broydens_direction(s, y, x, hessinv_estimate=None, impl='first'):
      r
      assert len(s) == len(y)
  
      if hessinv_estimate is not None:
          r = hessinv_estimate(x)
      else:
          r = x.copy()
  
      for i in range(len(s)):
          if impl == 'first':
              r.lincomb(1, r, y[i].inner(r), s[i])
          elif impl == 'second':
              r.lincomb(1, r, y[i].inner(x), s[i])
          else:
              raise RuntimeError('unknown `impl`')
  
      return r
docs: r"""Compute ``Hn^-1(x)`` for Broydens method.
  
      Parameters
      ----------
      s : sequence of `LinearSpaceElement`'s'
          The ``s`` coefficients in the Broydens update, see Notes.
      y : sequence of `LinearSpaceElement`'s'
          The ``y`` coefficients in the Broydens update, see Notes.
      x : `LinearSpaceElement`
          Point in which to evaluate the product.
      hessinv_estimate : `Operator`, optional
          Initial estimate of the hessian ``H0^-1``.
      impl : {'first', 'second'}, optional
          The type of Broydens method to use.
  
      Returns
      -------
      r : ``x.space`` element
          The result of ``Hn^-1(x)``.
  
      Notes
      -----
      For ``impl = 'first'``, :math:`H_n^{-1}` is defined recursively as
  
      .. math::
          H_{n+1}^{-1} = \left(I + s_n y_n^T \right) H_{n}^{-1}
  
      and for ``impl = 'second'``:
  
      .. math::
          H_{n+1}^{-1} = H_{n}^{-1} + s_n y_n^T
  
      With :math:`H_0^{-1}` given by ``hess_estimate``.