docs: def bsdmm(X, proxs_f, steps_f_cb, proxs_g=None, steps_g=None, Ls=None, update_order=None, steps_g_update='steps_f', max_iter=1000, e_rel=1e-6, e_abs=0, traceback=None):
      """Block-Simultaneous Method of Multipliers.
  
      This method is an extension of the linearized SDMM, i.e. ADMM for multiple
      constraints, for functions f with several arguments.
      The function f needs to be proper convex in every argument.
      It performs a block optimization for each argument while propagating the
      changes to other arguments.
  
      Args:
          X: list of initial Xs, will be updated
          proxs_f: proxed function f
              Signature prox(X,step, j=None, Xs=None) -> X'
          steps_f_cb: callback function to compute step size for proxs_f[j]
              Signature: steps_f_cb(j, Xs) -> Reals
          proxs_g: list of proxed functions
              [[prox_X0_0, prox_X0_1...],[prox_X1_0, prox_X1_1,...],...]
          steps_g: specific value of step size for proxs_g (experts only!)
              If set, needs to have same format as proxs_g.
              By default, set to the maximum value of step_f_j * ||L_i||_s^2.
          Ls: linear operators of the argument of g_i.
              If set, needs to have same format as proxs_g.
              Matrices can be numpy.array, scipy.sparse, or None (for identity).
          update_order: list of components to update in desired order.
          steps_g_update: relation between steps_g and steps_f (experts only!)
              'steps_f':  update steps_g as required by most conservative limit
              'fixed':    never update initial value of steps_g
              'relative': update initial values of steps_g propertional to changes
                          of steps_f
          e_rel: relative error threshold for primal and dual residuals
          e_abs: absolute error threshold for primal and dual residuals
          max_iter: maximum iteration number, irrespective of current residuals
          traceback: utils.Traceback to hold variable histories
  
      Returns:
          converged: whether the optimizer has converged within e_rel
          error: X^it - X^it-1
  
      Warning:
          Because of the potentially large list of optimization variables,
          setting traceback may exhaust memory. It should thus be run
          with a sufficiently small max_iter.
  
      See also:
          algorithms.sdmm
  
      Reference:
          Moolekamp & Melchior, Algorithm 3 (arXiv:1708.09066)
      """
  
      # Set up
      N = len(X)
      if proxs_g is None:
          proxs_g = [None] * N
      assert len(proxs_g) == N
      assert steps_g_update.lower() in ['steps_f', 'fixed', 'relative']
  
      if np.isscalar(e_rel):
          e_rel = [e_rel] * N
      if np.isscalar(e_abs):
          e_abs = [e_abs] * N
      steps_f = [None] * N
  
      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
  
      if steps_g_update.lower() == 'steps_f':
          if steps_g is not None:
              logger.debug("Setting steps_g = None for update strategy 'steps_f'.")
              steps_g = None
      if steps_g_update.lower() in ['fixed', 'relative']:
          if steps_g is None:
              logger.debug("Ignoring steps_g update strategy '%s' because steps_g is None." % steps_g_update)
              steps_g_update = 'steps_f'
  
      # if steps_g / Ls are None or single: create N duplicates
      if not hasattr(steps_g, "__iter__"):
          steps_g = [steps_g] * N
      if not hasattr(Ls, "__iter__"):
          Ls = [Ls] * N
      # check for cases in which a list was given
      assert len(steps_g) == N
      assert len(Ls) == N
  
      M = [0] * N
      for j in range(N):
          if proxs_g[j] is not None:
              if not hasattr(proxs_g[j], "__iter__"):
                  proxs_g[j] = [proxs_g[j]]
              M[j] = len(proxs_g[j])
              if not hasattr(steps_g[j], "__iter__"):
                  steps_g[j] = [steps_g[j]] * M[j]
              if not hasattr(Ls[j], "__iter__"):
                  Ls[j] = [Ls[j]] * M[j]
              assert len(steps_g[j]) == M[j]
              assert len(Ls[j]) == M[j]
  
      # need container for current-iteration steps_g and matrix adapters
      steps_g_ = []
      _L = []
      for j in range(N):
          if proxs_g[j] is None:
              steps_g_.append(None)
              _L.append(utils.MatrixAdapter(None))
          else:
              steps_g_.append([[None] for i in range(M[j])])
              _L.append([ utils.MatrixAdapter(Ls[j][m]) for m in range(M[j])])
  
      # Initialization
      Z, U = [],[]
      LX, R, S = [None] * N, [None] * N, [None] * N
      for j in range(N):
          Zj, Uj = utils.initZU(X[j], _L[j])
          Z.append(Zj)
          U.append(Uj)
  
      # containers
      converged, errors = [None] * N, [None] * N
      slack = [1.] * N
      it = 0
  
      if traceback is not None:
          for j in update_order:
              if M[j]>0:
                  _S = [np.zeros(X[j].shape, dtype=X[j].dtype) for n in range(M[j])]
              else:
                  _S = np.zeros(X[j].shape, dtype=X[j].dtype)
              traceback.update_history(it, j=j, X=X[j], steps_f=steps_f[j])
              traceback.update_history(it, j=j, M=M[j], steps_g=steps_g_[j], Z=Z[j], U=U[j],
                                R=U[j],
                                S=[np.zeros(X[j].shape, dtype=X[j].dtype) for n in range(M[j])])
  
      while it < max_iter:
  
          # iterate over blocks X_j
          for j in update_order:
              proxs_f_j = partial(proxs_f, j=j, Xs=X)
              steps_f_j = steps_f_cb(j, X) * slack[j]
  
              # update steps_g relative to change of steps_f ...
              if steps_g_update.lower() == 'relative':
                  for i in range(M[j]):
                      steps_g[j][i] *= steps_f_j / steps_f[j]
              steps_f[j] = steps_f_j
              # ... or update them as required by the most conservative limit
              if steps_g_update.lower() == 'steps_f':
                  for i in range(M[j]):
                      steps_g_[j][i] = utils.get_step_g(steps_f[j], _L[j][i].spectral_norm, N=N, M=M[j])
  
              # update the variables
              LX[j], R[j], S[j] = utils.update_variables(X[j], Z[j], U[j], proxs_f_j, steps_f[j], proxs_g[j], steps_g_[j], _L[j])
  
              # convergence criteria, adapted from Boyd 2011, Sec 3.3.1
              converged[j], errors[j] = utils.check_constraint_convergence(X[j], _L[j], LX[j], Z[j], U[j],
                  R[j], S[j], steps_f[j],steps_g_[j],e_rel[j], e_abs[j])
              # Optionally update the new state
              if traceback is not None:
                  traceback.update_history(it+1, j=j, X=X[j], steps_f=steps_f[j])
                  traceback.update_history(it+1, j=j, M=M[j], steps_g=steps_g_[j], Z=Z[j], U=U[j], R=R[j], S=S[j])
  
          if all(converged):
              break
          it += 1
  
      logger.info("Completed {0} iterations".format(it+1))
      if not all(converged):
          logger.warning("Solution did not converge")
  
      return converged, errors