docs: def pgm(X, prox_f, step_f, accelerated=False, relax=None, e_rel=1e-6, max_iter=1000, traceback=None):
      """Proximal Gradient Method
  
      Adapted from Combettes 2009, Algorithm 3.4.
      The accelerated version is Algorithm 3.6 with modifications
      from Xu & Yin (2015).
  
      Args:
          X: initial X, will be updated
          prox_f: proxed function f (the forward-backward step)
          step_f: step size, < 1/L with L being the Lipschitz constant of grad f
          accelerated: If Nesterov acceleration should be used
          relax: (over)relaxation parameter, 0 < relax < 1.5
          e_rel: relative error of X
          max_iter: maximum iteration, irrespective of residual error
          traceback: utils.Traceback to hold variable histories
  
      Returns:
          converged: whether the optimizer has converged within e_rel
          error: X^it - X^it-1
      """
  
      # init
      stepper = utils.NesterovStepper(accelerated=accelerated)
  
      if relax is not None:
          assert relax > 0 and relax < 1.5
  
      if traceback is not None:
          traceback.update_history(0, X=X, step_f=step_f)
          if accelerated:
              traceback.update_history(0, omega=0)
          if relax is not None:
              traceback.update_history(0, relax=relax)
  
      for it in range(max_iter):
  
          # use Nesterov acceleration (if omega > 0), automatically incremented
          omega = stepper.omega
          if omega > 0:
              _X = X + omega*(X - X_)
          else:
              _X = X
          # make copy for convergence test and acceleration
          X_ = X.copy()
  
          # PGM step
          X[:] = prox_f(_X, step_f)
  
          if relax is not None:
              X += (relax-1)*(X - X_)
  
          if traceback is not None:
              traceback.update_history(it+1, X=X, step_f=step_f)
              if accelerated:
                  traceback.update_history(it+1, omega=omega)
              if relax is not None:
                  traceback.update_history(it+1, relax=relax)
  
          # test for fixed point convergence
          converged = utils.l2sq(X - X_) <= e_rel**2*utils.l2sq(X)
          if converged:
              break
  
      logger.info("Completed {0} iterations".format(it+1))
      if not converged:
          logger.warning("Solution did not converge")
  
      return converged, X-X_