code: def _mixed_norm_solver_bcd(
      M,
      G,
      alpha,
      lipschitz_constant,
      maxit=200,
      tol=1e-8,
      verbose=None,
      init=None,
      n_orient=1,
      dgap_freq=10,
      use_accel=True,
      K=5,
  ):
      """Solve L21 inverse problem with block coordinate descent."""
      _, n_times = M.shape
      _, n_sources = G.shape
      n_positions = n_sources // n_orient
  
      if init is None:
          X = np.zeros((n_sources, n_times))
          R = M.copy()
      else:
          X = init
          R = M - np.dot(G, X)
  
      E = []  # track primal objective function
      highest_d_obj = -np.inf
      active_set = np.zeros(n_sources, dtype=bool)  # start with full AS
  
      alpha_lc = alpha / lipschitz_constant
  
      if use_accel:
          last_K_X = np.empty((K + 1, n_sources, n_times))
          U = np.zeros((K, n_sources * n_times))
  
      # First make G fortran for faster access to blocks of columns
      G = np.asfortranarray(G)
      # Ensure these are correct for dgemm
      assert R.dtype == np.float64
      assert G.dtype == np.float64
      one_ovr_lc = 1.0 / lipschitz_constant
  
      # assert that all the multiplied matrices are fortran contiguous
      assert X.T.flags.f_contiguous
      assert R.T.flags.f_contiguous
      assert G.flags.f_contiguous
      # storing list of contiguous arrays
      list_G_j_c = []
      for j in range(n_positions):
          idx = slice(j * n_orient, (j + 1) * n_orient)
          list_G_j_c.append(np.ascontiguousarray(G[:, idx]))
  
      for i in range(maxit):
          _bcd(
              G,
              X,
              R,
              active_set,
              one_ovr_lc,
              n_orient,
              alpha_lc,
              list_G_j_c,
          )
  
          if (i + 1) % dgap_freq == 0:
              _, p_obj, d_obj, _ = dgap_l21(
                  M, G, X[active_set], active_set, alpha, n_orient
              )
              highest_d_obj = max(d_obj, highest_d_obj)
              gap = p_obj - highest_d_obj
              E.append(p_obj)
              logger.debug(
                  "Iteration %d :: p_obj %f :: dgap %f :: n_active %d"
                  % (i + 1, p_obj, gap, np.sum(active_set) / n_orient)
              )
  
              if gap < tol:
                  logger.debug(
                      "Convergence reached ! (gap: %s < %s)"
                      % (gap, tol)
                  )
                  break
  
          # using Anderson acceleration of the primal variable for faster
          # convergence
          if use_accel:
              last_K_X[i % (K + 1)] = X
  
              if i % (K + 1) == K:
                  for k in range(K):
                      U[k] = (
                          last_K_X[k + 1].ravel() - last_K_X[k].ravel()
                      )
                  C = U @ U.T
                  # at least on ARM64 we can't rely on np.linalg.solve to
                  # reliably raise LinAlgError here, so use SVD instead
                  # equivalent to:
                  # z = np.linalg.solve(C, np.ones(K))
                  u, s, _ = np.linalg.svd(C, hermitian=True)
                  if s[-1] <= 1e-6 * s[0] or not np.isfinite(s).all():
                      logger.debug(
                          "Iteration %d: LinAlg Error" % (i + 1)
                      )
                      continue
                  z = ((u * 1 / s) @ u.T).sum(0)
                  c = z / z.sum()
                  X_acc = np.sum(
                      last_K_X[:-1] * c[:, None, None], axis=0
                  )
                  _grp_norm2_acc = groups_norm2(X_acc, n_orient)
                  active_set_acc = _grp_norm2_acc != 0
                  if n_orient > 1:
                      active_set_acc = np.kron(
                          active_set_acc, np.ones(n_orient, dtype=bool)
                      )
                  p_obj = _primal_l21(
                      M, G, X[active_set], active_set, alpha, n_orient
                  )[0]
                  p_obj_acc = _primal_l21(
                      M,
                      G,
                      X_acc[active_set_acc],
                      active_set_acc,
                      alpha,
                      n_orient,
                  )[0]
                  if p_obj_acc < p_obj:
                      X = X_acc
                      active_set = active_set_acc
                      R = M - G[:, active_set] @ X[active_set]
  
      X = X[active_set]
  
      return X, active_set, E
  
imports: import functools
  from math import sqrt
  import numpy
  from ..time_frequency._stft import istft, stft, stft_norm1, stft_norm2
  from ..utils import _check_option, _get_blas_funcs, _validate_type, logger, sum_squared, verbose, warn
  from .mxne_debiasing import compute_bias
label_desc: Use image node_id 8 for calling a global function with example usage: _mixed_norm_solver_bcd(M, G, alpha, lipschitz_constant, maxit, tol, verbose, init, n_orient, dgap_freq, use_accel, K) and returns: X,  active_set,  E