def broydens_method(f, x, line_search=1.0, impl='first', maxiter=1000,
                    tol=1e-15, hessinv_estimate=None,
                    callback=None):
    r"""Broyden's first method, a quasi-Newton scheme.

    Notes
    -----
    This is a general and optimized implementation of Broyden's method,
    a quasi-Newton method for solving a general unconstrained optimization
    problem

    .. math::
        \min f(x)

    for a differentiable function
    :math:`f: \mathcal{X}\to \mathbb{R}` on a Hilbert space
    :math:`\mathcal{X}`. It does so by finding a zero of the gradient

    .. math::
        \nabla f: \mathcal{X} \to \mathcal{X}

    using a Newton-type update scheme with approximate Hessian.

    The algorithm is described in [Bro1965] and [Kva1991], and in a
    `Wikipedia article
    <https://en.wikipedia.org/wiki/Broyden's_method>`_.

    Parameters
    ----------
    f : `Functional`
        Functional with ``f.gradient``
    x : ``f.domain`` element
        Starting point of the iteration
    line_search : float or `LineSearch`, optional
        Strategy to choose the step length. If a float is given, uses it as a
        fixed step length.
    impl : {'first', 'second'}, optional
        What version of Broydens method to use. First is also known as Broydens
        'good' method, while the second is known as Broydens 'bad' method.
    maxiter : int, optional
        Maximum number of iterations.
        ``tol``.
    tol : float, optional
        Tolerance that should be used for terminating the iteration.
    hessinv_estimate : `Operator`, optional
        Initial estimate of the inverse of the Hessian operator. Needs to be an
        operator from ``f.domain`` to ``f.domain``.
        Default: Identity on ``f.domain``
    callback : callable, optional
        Object executing code per iteration, e.g. plotting each iterate.

    References
    ----------
    [Bro1965] Broyden, C G. *A class of methods for solving nonlinear
    simultaneous equations*. Mathematics of computation, 33 (1965),
    pp 577--593.

    [Kva1991] Kvaalen, E. *A faster Broyden method*. BIT Numerical
    Mathematics 31 (1991), pp 369--372.
    """
    grad = f.gradient
    if x not in grad.domain:
        raise TypeError('`x` {!r} is not in the domain of `grad` {!r}'
                        ''.format(x, grad.domain))

    if not callable(line_search):
        line_search = ConstantLineSearch(line_search)

    impl, impl_in = str(impl).lower(), impl
    if impl not in ('first', 'second'):
        raise ValueError('`impl` {!r} not understood'
                         ''.format(impl_in))

    ss = []
    ys = []

    grad_x = grad(x)
    for i in range(maxiter):
        # find step size
        search_dir = -_broydens_direction(ss, ys, grad_x,
                                          hessinv_estimate, impl)
        dir_deriv = search_dir.inner(grad_x)
        if np.abs(dir_deriv) == 0:
            return  # we found an optimum

        step = line_search(x, search_dir, dir_deriv)

        # update x
        x_update = step * search_dir
        x += x_update

        # compute new gradient
        grad_x, grad_x_old = grad(x), grad_x
        delta_grad = grad_x - grad_x_old

        # update hessian.
        # TODO: reuse from above
        v = _broydens_direction(ss, ys, delta_grad, hessinv_estimate,
                                impl)
        if impl == 'first':
            divisor = x_update.inner(v)

            # Test for convergence
            if np.abs(divisor) < tol:
                if grad_x.norm() < tol:
                    return
                else:
                    # Reset if needed
                    ys = []
                    ss = []
                    continue
            u = (x_update - v) / divisor
            ss.append(u)
            ys.append(x_update)
        elif impl == 'second':
            divisor = delta_grad.inner(delta_grad)

            # Test for convergence
            if np.abs(divisor) < tol:
                if grad_x.norm() < tol:
                    return
                else:
                    # Reset if needed
                    ys = []
                    ss = []
                    continue
            u = (x_update - v) / divisor
            ss.append(u)
            ys.append(delta_grad)

        if callback is not None:
            callback(x)
r"""Broyden's first method, a quasi-Newton scheme.

    Notes
    -----
    This is a general and optimized implementation of Broyden's method,
    a quasi-Newton method for solving a general unconstrained optimization
    problem

    .. math::
        \min f(x)

    for a differentiable function
    :math:`f: \mathcal{X}\to \mathbb{R}` on a Hilbert space
    :math:`\mathcal{X}`. It does so by finding a zero of the gradient

    .. math::
        \nabla f: \mathcal{X} \to \mathcal{X}

    using a Newton-type update scheme with approximate Hessian.

    The algorithm is described in [Bro1965] and [Kva1991], and in a
    `Wikipedia article
    <https://en.wikipedia.org/wiki/Broyden's_method>`_.

    Parameters
    ----------
    f : `Functional`
        Functional with ``f.gradient``
    x : ``f.domain`` element
        Starting point of the iteration
    line_search : float or `LineSearch`, optional
        Strategy to choose the step length. If a float is given, uses it as a
        fixed step length.
    impl : {'first', 'second'}, optional
        What version of Broydens method to use. First is also known as Broydens
        'good' method, while the second is known as Broydens 'bad' method.
    maxiter : int, optional
        Maximum number of iterations.
        ``tol``.
    tol : float, optional
        Tolerance that should be used for terminating the iteration.
    hessinv_estimate : `Operator`, optional
        Initial estimate of the inverse of the Hessian operator. Needs to be an
        operator from ``f.domain`` to ``f.domain``.
        Default: Identity on ``f.domain``
    callback : callable, optional
        Object executing code per iteration, e.g. plotting each iterate.

    References
    ----------
    [Bro1965] Broyden, C G. *A class of methods for solving nonlinear
    simultaneous equations*. Mathematics of computation, 33 (1965),
    pp 577--593.

    [Kva1991] Kvaalen, E. *A faster Broyden method*. BIT Numerical
    Mathematics 31 (1991), pp 369--372.