def bfgs(x0, rho, f_df, maxiter=50, method='BFGS'):
    """
    Proximal operator for minimizing an arbitrary function using BFGS

    Uses the BFGS algorithm to find the proximal update for an arbitrary function, `f`, whose gradient is known.

    Parameters
    ----------
    x0 : array_like
        The starting or initial point used in the proximal update step

    rho : float
        Momentum parameter for the proximal step (larger value -> stays closer to x0)

    f_df : function
        The objective function and gradient

    maxiter : int, optional
        Maximum number of iterations to take (default: 50)

    method : str, optional
        Which scipy.optimize algorithm to use (default: 'BFGS')

    Returns
    -------
    theta : array_like
        The parameter vector found after running the proximal update step
    """

    # keep track of the original shape
    orig_shape = x0.shape

    # specify the objective function and gradient for the proximal operator
    def f_df_augmented(x):

        xk = x.reshape(orig_shape)

        obj, grad = f_df(xk)

        g = obj + (rho / 2.) * np.sum((xk - x0) ** 2)
        dg = (grad + rho * (xk - x0)).ravel()

        return g, dg

    # minimize via BFGS
    options = {'maxiter': maxiter, 'disp': False}
    return opt.minimize(f_df_augmented,
                        x0.ravel(),
                        method=method,
                        jac=True,
                        options=options).x.reshape(orig_shape)
Proximal operator for minimizing an arbitrary function using BFGS

    Uses the BFGS algorithm to find the proximal update for an arbitrary function, `f`, whose gradient is known.

    Parameters
    ----------
    x0 : array_like
        The starting or initial point used in the proximal update step

    rho : float
        Momentum parameter for the proximal step (larger value -> stays closer to x0)

    f_df : function
        The objective function and gradient

    maxiter : int, optional
        Maximum number of iterations to take (default: 50)

    method : str, optional
        Which scipy.optimize algorithm to use (default: 'BFGS')

    Returns
    -------
    theta : array_like
        The parameter vector found after running the proximal update step