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
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 ) : 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 : 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 not hasattr ( steps_g , "__iter__" ) : steps_g = [ steps_g ] * N if not hasattr ( Ls , "__iter__" ) : Ls = [ Ls ] * N 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 ] 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 ] ) ] ) 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 ) 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 : 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 ] 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 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 ] ) 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 ] ) 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 ] ) 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
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)