code: def bfgs(x0, rho, f_df, maxiter=50, method='BFGS'):
      
  
      # 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)
docs: 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