docs: def sdmm(X, prox_f, step_f, proxs_g=None, steps_g=None, Ls=None, e_rel=1e-6, e_abs=0, max_iter=1000, traceback=None):
      """Simultaneous-Direction Method of Multipliers
  
      This method is an extension of the linearized ADMM for multiple constraints.
  
      Args:
          X: initial X, will be updated
          prox_f: proxed function f
          step_f: step size for prox_f
          proxs_g: list of proxed functions
          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 * ||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).
          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
  
      See also:
          algorithms.admm
  
      Reference:
          Moolekamp & Melchior, Algorithm 2 (arXiv:1708.09066)
      """
  
      # fall-back to simple ADMM
      if proxs_g is None or not hasattr(proxs_g, '__iter__'):
          return admm(X, prox_f, step_f, prox_g=proxs_g, step_g=steps_g, L=Ls, e_rel=e_rel, max_iter=max_iter, traceback=traceback)
  
      # from here on we know that proxs_g is a list
      M = len(proxs_g)
  
      # if steps_g / Ls are None or single: create M duplicates
      if not hasattr(steps_g, "__iter__"):
          steps_g = [steps_g] * M
      if not hasattr(Ls, "__iter__"):
          Ls = [Ls] * M
      # check for cases in which a list was given
      assert len(steps_g) == M
      assert len(Ls) == M
  
      # get/check compatible step sizes for g
      # use matrix adapter for convenient & fast notation
      _L = []
      for i in range(M):
          _L.append(utils.MatrixAdapter(Ls[i]))
          # get/check compatible step size for g
          if steps_g[i] is None:
              steps_g[i] = utils.get_step_g(step_f, _L[i].spectral_norm, M=M)
  
      # Initialization
      Z,U = utils.initZU(X, _L)
      it, omega = 0, 0
  
      if traceback is not None:
          traceback.update_history(it, X=X, step_f=step_f, omega=omega)
          traceback.update_history(it, M=M, Z=Z, U=U, R=U, S=[np.zeros(X.shape, dtype=X.dtype) for n in range(M)], steps_g=steps_g)
  
      while it < max_iter:
  
          # update the variables
          LX, R, S = utils.update_variables(X, Z, U, prox_f, step_f, proxs_g, steps_g, _L)
  
          if traceback is not None:
              traceback.update_history(it+1, X=X, step_f=step_f, omega=omega)
              traceback.update_history(it+1, M=M, Z=Z, U=U, R=R, S=S, steps_g=steps_g)
  
          # convergence criteria, adapted from Boyd 2011, Sec 3.3.1
          converged, errors = utils.check_constraint_convergence(X, _L, LX, Z, U, R, S, step_f, steps_g,
                                                                   e_rel, e_abs)
  
          if converged:
              break
  
          it += 1
  
          # if X and primal residual does not change: decrease step_f and step_g, and restart
          if it > 1:
              if (X == X_).all() and all([(R[i] == R_[i]).all() for i in range(M)]):
                  step_f /= 2
                  for i in range(M):
                      steps_g[i] /= 2
  
                  # re-init
                  it = 0
  
                  Z,U  = utils.initZU(X, _L)
                  if traceback is not None:
                      traceback.reset()
                      traceback.update_history(it, X=X, step_f=step_f)
                      traceback.update_history(it, M=M, Z=Z, U=U, R=U,
                                        S=[np.zeros(X.shape, dtype=X.dtype) for n in range(M)], steps_g=steps_g)
                  logger.info("Restarting with step_f = %.3f" % step_f)
  
          R_ = R
          X_ = X.copy()
  
      logger.info("Completed {0} iterations".format(it+1))
      if not converged:
          logger.warning("Solution did not converge")
  
      return converged, errors