docs: def bpgm(X, proxs_f, steps_f_cb, update_order=None, accelerated=False, relax=None, max_iter=1000, e_rel=1e-6, traceback=None):
      """Block Proximal Gradient Method.
  
      Also know as Alternating Proximal Gradient Method, it performs Proximal
      gradient (forward-backward) updates on each block in alternating fashion.
  
      Args:
          X: list of initial Xs, will be updated
          proxs_f: proxed function f
              Signature prox(X,step, j=None, Xs=None) -> X'
          step_f: step size, < 1/L with L being the Lipschitz constant of grad f
          steps_f_cb: callback function to compute step size for proxs_f[j]
              Signature: steps_f_cb(j, Xs) -> Reals
          update_order: list of components to update in desired order.
          accelerated: If Nesterov acceleration should be used for all variables
          relax: (over)relaxation parameter for each variable, < 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
      """
      # Set up
      N = len(X)
      if np.isscalar(e_rel):
          e_rel = [e_rel] * N
  
      if relax is not None:
          assert relax > 0 and relax < 1.5
  
      if update_order is None:
          update_order = range(N)
      else:
          # we could check that every component is in the list
          # but one can think of cases when a component is *not* to be updated.
          #assert len(update_order) == N
          pass
  
      # init
      X_ = [None] * N
      stepper = utils.NesterovStepper(accelerated=accelerated)
  
      if traceback is not None:
          for j in update_order:
              traceback.update_history(0, j=j, X=X[j], steps_f=None)
              if accelerated:
                  traceback.update_history(0, j=j, omega=0)
              if relax is not None:
                  traceback.update_history(0, j=j, relax=relax)
  
      for it in range(max_iter):
          # use Nesterov acceleration (if omega > 0), automatically incremented
          omega = stepper.omega
  
          # iterate over blocks X_j
          for j in update_order:
  
              # tell prox the state of other variables
              proxs_f_j = partial(proxs_f, j=j, Xs=X)
              steps_f_j = steps_f_cb(j, X)
  
              # acceleration?
              # check for resizing: if resize ocurred, temporily skip acceleration
              if omega > 0 and X[j].shape == X_[j].shape:
                  _X = X[j] + omega*(X[j] - X_[j])
              else:
                  _X = X[j]
  
              # keep copy for convergence test (and acceleration)
              X_[j] = X[j].copy()
  
              # PGM step, force inline update
              X[j][:] = proxs_f_j(_X, steps_f_j)
  
              if relax is not None:
                  X[j] += (relax-1)*(X[j] - X_[j])
  
              if traceback is not None:
                  traceback.update_history(it+1, j=j, X=X[j], steps_f=steps_f_j)
                  if accelerated:
                      traceback.update_history(it+1, j=j, omega=omega)
                  if relax is not None:
                      traceback.update_history(it+1, j=j, relax=relax)
  
          # test for fixed point convergence
          # allowing for transparent resizing of X: need to check shape of X_
          errors = [X[j] - X_[j] if X[j].shape == X_[j].shape else X[j] for j in range(N)]
          converged = [utils.l2sq(errors[j]) <= e_rel[j]**2*utils.l2sq(X[j]) for j in range(N)]
          if all(converged):
              break
  
      logger.info("Completed {0} iterations".format(it+1))
      if not all(converged):
          logger.warning("Solution did not converge")
  
      return converged, errors