{'fraction_non_alphanumeric': 0.05813453413949252, 'fraction_numerical': 0.01630448531457917, 'mean_word_length': 4.04202806122449, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 36, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Linearly Convergent Douglas-Rachford Splitting Solver for Markovian Information-Theoretic Optimization Problems

In this work, we propose solving the Information bottleneck (IB) and Privacy Funnel (PF) problems with Douglas-Rachford Splitting methods (DRS). We study a general Markovian information-theoretic Lagrangian that includes IB and PF into a unified framework. We prove the linear convergence of the proposed solvers using the Kurdyka-{\L}ojasiewicz inequality. Moreover, our analysis is beyond IB and PF and applies to any convex-weakly convex pair objectives. Based on the results, we develop two types of linearly convergent IB solvers, with one improves the performance of convergence over existing solvers while the other can be independent to the relevance-compression trade-off. Moreover, our results apply to PF, yielding a new class of linearly convergent PF solvers. Empirically, the proposed IB solvers IB obtain solutions that are comparable to the Blahut-Arimoto-based benchmark and is convergent for a wider range of the penalty coefficient than existing solvers. For PF, our non-greedy solvers can characterize the privacy-utility trade-off better than the clustering-based greedy solvers.


I. INTRODUCTION
Recently, adopting information-theoretic metrics as optimization objectives drew significant attention from machine learning and data science communities. Among which, the information bottleneck (IB) methods [1], studying the complexity-relevance trade-off of representations, have been applied to representation learning, classification, clustering problems and successfully combined with deep neural networks, giving impressive performance. However, the advances of machine learning also bring new challenges. Collecting a large amount of data burdens the conventional centralized learning approach; On the other hand, the leakage of sensitive information becomes a major concern for both the end users and industries. Closely related to this, the privacy funnel (PF) [2], another information-theoretic optimization problem, has been adopted in optimization objectives as it characterizes the trade-off between the leakage of sensitive information and the utility of observations.

A. Information Bottleneck and Privacy Funnel
We start with a brief review of the IB and PF problems. Given the joint probability of observations X and targets Y , the IB methods aim to find a representation Z of the observations that is minimal in expression complexity but retains most relevance to the target. In IB, this goal is formulated as a constrained optimization problem [1]: minimize p(z|x)∈Ω I(X; Z), subject to : I(Y ; Z) > I 0 , z p(z|x) = 1, ∀x ∈ X , where the constant I 0 ≥ 0 controls the trade-off between I(Z; X) and I(Z; Y ); Ω is a compound probability simplex and the variable to optimize is the conditional probability p(Z|X). In solving the IB problem, (1) can be relaxed as the following IB Lagrangian: L IB := γI(X; Z) − I(Y ; Z), where the multiplier γ controls the trade-off between the two competing objectives. The PF problem considers the reverse scenario: denote Y as the sensitive information, X the public information and Z the observation, the goal of PF is to design an information assignment p(Z|X) such that the utility of Z to the public information is maximized but minimal sensitive information is revealed from Z. This goal satisfies the Markov chain Y → X → Z where the utility and privacy leakage are measured by the mutual information I(Z; X) and I(Z; Y ) respectively. Then to find T-H. Huang  the optimal assignment, one minimizes I(Z; Y ) but maximizes I(Z; X) over p(Z|X) with a known p(X, Y ). Through the Lagrange multipliers, this optimization problem can be written as the following PF Lagrangian: L P F := βI(Y ; Z) − I(X; Z). (3) Both IB and PF problems are known to be non-convex therefore are difficult to solve [1], [2]. A common relaxation to minimize the IB and PF Lagrangian is to fix a trade-off parameter (the multiplier) and solve the relaxed problem. Then by varying a range of the trade-off parameter, one can obtain a collection of p(Z|X), hence pairs of I(Z; X), I(Z; Y ).
To evaluate solutions obtained from both IB and PF problems, one plots each pair of I(Z; X), I(Z; Y ) as the x, y-axes respectively. The plot is called the information plane [1], [3]. On the information plane, the best I(Z; Y ) given a fixed I(Z; Y ) is revealed. For a fixed I(Z; X), a higher I(Z; Y ) is better in IB whereas a lower I(Z; Y ) is preferred in PF. Then the frontier formed from the set of best I(Z; Y ) over a range of I(Z; X) characterizes the relevance-complexity trade-off for IB and the privacy-utility trade-off for PF [1]- [4]. This frontier is known as the Pareto-frontier [5], [6] and can be used to evaluate the performance of IB and PF solvers.

B. IB and PF Solvers
The formulations of IB and PF are quite general. Hence a variety of solvers have been developed and specialized to a wide range of applications [1], [3], [7]- [14]. Inspired by the Blahut-Arimoto (BA) algorithm in rate distortion theory, a set of self-consistent equations is proposed in the seminal work for IB [1].
where M (x, γ) is a normalization function and the superscript k denotes the iteration counter. While the original formulation of the BA-based algorithm focuses on discrete random variables X, Y with known joint distribution p(X, Y ), it is later extended to a continuous special case where p(X, Y ) is jointly Gaussian [15]. The Gaussian IB has a closed-form expression and is fully characterized in terms of the structure of the latent representation Z and the associated covariance matrices.
In clustering problems, we start with a set of partitions Z formed from the realizations of X. Then, iterative greedy cluster aggregation-based algorithms merge two or more partitions, resulting in a new partition Z ′ such that the largest I(Z; Y ) − I(Z ′ ; Y ) is achieved [3], [5]. The cluster aggregation-based algorithms can be used to solve the PF problem [2], [5]. The intuition is that when minimizing (3) over p(z|x) under the restriction that the entries are either 1 or 0, the PF problem reduces to:L PF = βI(Z; Y ) − H(Z).
While the above problem is convex with respect to p(z|x), the feasible set is more restricted so only a limited number of points on the information plane are found. Following the recent innovations in non-convex optimization with splitting methods, [12] empirically shows that the IB problem can be solved with the alternating direction method of multipliers (ADMM) [16] by separating the IB Lagrangian into three terms of (conditional) Shannon entropy functions H(Z), H(Z|X) and H(Z|Y ). In our earlier work [13], we simplify the method into two terms and further prove the convergence.
The aforementioned solvers are susceptible to scalability when the cardinalities of X and Y are large. This issue has been recently addressed through variational inference techniques. For instance, recent works have proposed the deep variational IB (DVIB) [10] and a deep variational solver for the PF [14]. However, the variational inference methods are approximate in nature and will provide exact solutions for the IB or PF Lagrangian only if the variational class of distributions includes the exact posterior [17], [18]. Note that for discrete IB scenario, the variational inference-based method reduces to the BA-based algorithm (4c).

C. Contributions
In this work, we consider a general class of discrete Markovian information-theoretic optimization problems. This class includes both the IB and PF problems as special cases. We develop a framework for constructing solvers for this class of problems via the Douglas-Rachford splitting (DRS) methods [21]. We construct a unified convergence proof of the proposed set of solvers which guarantees convergence without relying on additional regularization terms; unlike earlier works for ADMMbased solvers [12], [13]. Further, we prove that the rate of convergence for the proposed solvers is locally linear. The main F :σ F -hypo-cvx. (0 < α < 2); F :σ F -str. cvx. (2 ≤ α < 4) F is Lp-smooth tool in our convergence analysis is the the Kurdyka-Łojasiewicz (KŁ) inequality [22]- [24] which is included in appendices for completeness. For the IB special case, our approach depends on decomposing the objective function into two sub-objectives via two different alternative approaches. Each approach yields a different convergent solver as described in Section IV. The first solver guarantees convergence by requiring a smaller penalty coefficient, as compared with existing ADMM based IB solvers, and exploiting the strong-convexity of sub-objective. The second solver ensures convergence by exploiting Lipschitz smoothness and weak-convexity of one of the sub-objectives, independent of the trade-off parameter γ (except for the case where X is obtained from Y through a deterministic mapping).
For the PF special case, we develop a new solver which recovers all the points on the information plane which are found by existing clustering-based greedy algorithms and adds to them more points not reachable by existing solvers. Contrary to the work reported in [14], the decoder part of our solver satisfies the Markovian condition (i.e., it is independent of the sensitive information Y ). Hence, it can explore the privacy-utility trade-off better in the discrete and Markovian settings.
Empirically, we evaluate the new IB and PF solvers on both synthetic and real-world data. For IB case, the proposed solvers are shown to be more comprehensive in exploring the information plane as compared with the benchmark BA-based solver. The new IB solvers are also convergent for a broader range of design parameters compared to existing ADMM-IB solvers. Similarly, for the PF special case, our numerical results illustrate that the new solver is able to obtain more points on the information plane than existing greedy based solvers [2], [5]. Finally, we show that our solvers are efficient in computing the variational bounds on both the IB and PF cases [10] while satisfying the Markovian condition.

II. PROBLEM FORMULATION
In this work, we study the following class of Markovian information-theoretic optimization problems with the joint probability p(x, y) assumed to be known: where Y → X → Z forms a Markov chain and ρ z , ρ z|x , ρ z|y are some constants. We focus on discrete settings and propose solving the problem with splitting methods. The IB and PF problems are special cases of (6). The IB problem corresponds to the following set of coefficients: whereas the coefficients for the PF problem are: In our earlier work, we show that the IB problem has the strongly convex-weakly convex structure [13], [25] under some mild smoothness conditions [26]- [28]. In this work, we further show that both the IB and PF problems can be expressed as a more general convex-weakly convex structure. In turns, we can apply splitting methods to solve the them, leveraging the recent advances in non-convex optimization [20], [23], [24].

A. Splitting Methods for Structured Optimization Problems
Splitting methods are based on the augmented Lagrangian methods (ALM) [29] where the constraints are added to the original Lagrangian as extra penalties. The most relevant case to ours is the following linearly constrained augmented Lagrangian (The mapping from the proposed Markovian information-theoretic optimization problem (6) to it will be explained in details in Section IV): where F, G are sub-objective functions of the original Lagrangian and p, q are primal and augmented variables. After separation, F depends on p while G depends on q only. ν is the dual variable and c > 0 is the penalty coefficient. · denotes an inner product and · is a 2-norm if not specified. The matrices A, B are subjected to certain rank regularity (e.g. full-row rank). As a remark, when the penalty vanishes, Ax * − By * ≈ 0, then a minimize (p * , q * ) to the augmented Lagrangian (9) is a minimizer to the original Lagrangian. In this sense, solving the augmented Lagrangian is easier because the equality constraints can be relaxed during optimization [29]. Instead of jointly minimizing (9) over p, q, ν, splitting methods solve (9) in an alternating fashion. More specifically, splitting methods minimize (9) by repeatedly updating p,q then ν. When both sub-objectives F, G are convex, the convergence of splitting methods are well studied [16], [30]. Furthermore, the rate of convergence of these solvers have been characterized in [31], [32]. In contrast, there are less results for the cases where F (p) + G(q) is non-convex. Recently, several works have empirically found that if one of the two sub-objectives is convex and the other weakly convex, then splitting methods can solve this problem effectively [33]. Following the discovery, convergence of splitting methods for objective with this structure has been developed in [19], [20], [25], [34], [35]. In addition to convergence, the rate of convergence can also be shown to be locally linear through the Kurdyka-Łojasiewicz inequality (see Appendix A-B for details).
Among the general class of splitting methods, a well-known solver is the alternating direction method of multiplier (ADMM). In [19], ADMM is adopted to solve non-convex quadratic programming with simplex constraint and it is shown to outperform the state-of-the-art. In [36], the ADMM algorithm is used to solve non-convex Lasso faster than known methods and with better performance. In [33] the Douglas-Rachford Splitting (DRS) [21] method, where ADMM is a special case of it, is used to find sparse solutions of a linear system [33].

B. The Proposed Solvers
We propose two solvers for the augmented Lagrangian in (9). The first solver, which will be referred to as Solver I in the sequel, is described by the set of update equation: where the superscript k denotes the iteration counter; ν 1/2 is an intermediate dual variable accounting the fixed-point relaxation [37] and Ω p , Ω q are probability simplexes for p, q respectively; α > 0 is a relaxation parameter where 1 < α < 2 is "overrelaxation" while 0 < α < 1 is "under-relaxation". We refer to [37] for details. Note that α = 1 recovers the ADMM whereas α = 2 corresponds to the Peaceman-Rachford splitting (PRS) [38]. Alternatively, by changing the updating order of (10), we obtain our second DRS solver, which will be denoted by Solver II, as described by.
Remarkably, while Solver I and Solver II seem different only in the order they optimize p, q, when implemented with gradient descent, their first-order necessary conditions are quite different. In particular, the dual ascend connects to the gradient difference of the sub-objective function F for Solver I whereas the connection between the dual variables is on the gradient of G for Solver II.
Moreover, from the perspective of parameter selection, as will be shown in the next section, the rate of convergence of the two solvers depend on the relaxation parameter. Specifically, when 1 ≤ α ≤ 2, the rates are Q-linear while for the case where 0 < α < 1, we can only show R-linear rates, which are weaker than its Q-linear counterpart [39]. However, the strong rate comes with additional assumptions on the linear constraints which limits the problems it applies to. Also, it turns out that the convergence of Solver I and II depends on the choice of the penalty coefficient c [19], [20]. Therefore, knowing the smallest penalty coefficient that assures convergence and hence the rate, is essential to the success of applying DRS solvers to non-convex optimization problems.

III. MAIN RESULTS
Our main results are the convergence and rate analysis of the two solvers under three different sets of assumptions (Table  I), which serves as a general guideline for applying DRS solvers to non-convex optimization problems beyond IB and PF. For each of the algorithm-assumption pairs listed in Table I, we characterize the rate of convergence as locally linear using the KŁ inequality. In particular, we show that the augmented Lagrangian (9) solved with the two algorithms satisfies the KŁ property with exponent θ = 1/2, which corresponds to linear rate of convergence [22], [24], [40].
Theorem 1: Suppose F is σ F -strongly convex and L p -smooth, while G is ω G -restricted weakly convex and A is positive definite. Then for 0 < α ≤ 2, the sequence {w k } k∈N , where w k := (p k , Bq k , ν k ) denotes the collective point at step k, obtained from Solver I is bounded if c > max{ω G , (L p + σ F )/(αµ 2 A )} where µ A is the smallest eigenvalue of the matrix A. Moreover, the sequence converges to a stationary point (p * , Bq * , ν * ) at a linear rate locally.
proof sketch: The full convergence analysis and the statements for enabling lemmas are deferred to Appendix A while the details for deriving this theorem are deferred to Appendix N, here we explain the key ideas.
To prove convergence, the key step is to develop a sufficient decrease lemma (Lemma 6) where the difference of the objective function L c between consecutive updates (from step k to k + 1) is lower-bounded by a combination of strictly positive squared norms.
In developing Lemma 6, we expand the consecutive steps according to the Solver I: For (12a) and (12c), due to the identity: , which links the linear constraints to the dual ascend (10a), (10c): Then in (12b), since σ F is strongly convex and L p -smooth, we can use Lemma 4 to find a lower bound that consists of ∇F (p k ) − ∇F (p k+1 ) 2 and p k − p k+1 2 . The term with the gradient difference of F connects to the dual variable as ν k − ν k+1 2 = A −T (∇F (p k )− ∇F (p k+1 )) 2 through the first-order minimizer condition (39) and A being positive definite. This balances the first negative squared norm in the r.h.s. of (13).
On the other hand, for (12d), since G is ω G -restricted weakly convex (Definition 4) w.r.t. B, we have a lower bound with an additional negative squared norm −ω G Bq k − Bq k+1 2 . This negative term is balanced by the penalty coefficient c as we end up getting (c − σ G )/2 Bq k − Bq k+1 2 . After rearranging the terms we get the sufficient decrease lemma: where L k c := L c (p k , q k , ν k ) denotes the function value at the start of step k, and From the above, we hence know the conditions in terms of the penalty coefficient c and relaxation parameter α such that δ p , δ q , δ ν are non-negative. If the coefficients are non-negative, the sufficient decrease implies convergence. Moreover, we show that the rate is locally linear.
To prove linear rate of convergence, we first show that around a local stationary point, the KŁ property [23], [40], [41] are satisfied with an exponent θ = 1/2, which gives: where K 1 > 0 is a constant. The second relation that is required is the following upper bound to the norm of gradient L k+1 where K 2 > 0 is another constant. If (14), (15) and the sufficient decrease lemma hold, then by Lemma 11, owing to [22], θ = 1/2 corresponds to the desired result, i.e. the rate of convergence is Q-linear (Appendix J). In proving (14), we separate the goal into two steps. First we find an upper bound of L k+1 c − L * c (Lemma 14) by exploiting the strong convexity of F and the restricted weak convexity of G: Then by substituting the first-order necessary conditions for minimizers (39) into the gradient of F, G, along with the relation Ap * = Bq * at a stationary point, we obtain the following relation: Second, using (39) again along with the assumption that the matrix A is positive definite, we get: Combining the two inequalities through the Cauchy-Schwarz inequality u − v 2 ≥ (1 − t) u 2 + (1 + 1/t) v 2 , t > 0, we prove the that the KŁ exponent θ = 1/2. As for (15), following the same reasons for (17), we have: where K 3 > 0. Substitute (13) into the first term in the r.h.s. of the above inequality, assuming 1 ≤ α ≤ 2, we arrive at the desired result. Note that the relaxation parameter we focus on has range 0 < α ≤ 2, and the above Q-linear rate of convergence result does not apply to 0 < α < 1. To address this, inspired by the recent work [19], we prove R-linear rate of convergence for this region. The key relation that enables this result is: where K 4 > 0. It turns out that if the sequence of the function value {L k c } k∈N converges Q-linearly, implied by the sufficient decrease lemma, and (18) holds, then we can prove R-linear rate of convergence of the sequence {w k } k∈N0 , k > N 0 ∈ N locally around the neighborhood of a stationary point w * (Lemma 13, proved in Appendix K).
In proving (18), from the sufficient decrease lemma, we can substitute c > ω G into (16) and get: Then the desired relation follows by letting K 4 := max{cµ 2 A /2, ρ * }. Note that Theorem 1 requires the function F to be strongly convex and G to be restricted weakly convex which limits the class of functions that the theorem applies to. To address this, we adopt Solver II whose convergence is also locally linear but only requires F to be convex and G to be weakly convex.
Theorem 2: Let F be convex, G be σ G -weakly convex and L q -smooth, B be positive definite and A be full row rank, then for 0 < α < 2 the sequence {w k } k∈N where w k := (Ap k , q k , ν k ), obtained from Solver II is bounded if Moreover, the sequence converges to a stationary point w * := (Ap * , q * , ν * ) at linear rate locally. proof sketch: The details are deferred to Appendix Q. The steps are similar to the proof sketch for the first theorem. For the sufficient decrease lemma (Appendix G), we again divide L k c − L k+1 c through four steps, but now according to the algorithm (11d). For (12a) and (12c), using identity (37), we have: As for (12b), by the convexity of F , first order necessary condition for minimizers (41) and the identity 2 where in the expression above, we hide q k , ν k as they are kept unchanged in p-update. On the other hand, for the q-update (12d), due to σ G -weak convexity of G, we have: To balance the negative squared norm, note that in Solver II, the q-update precedes the dual update, this connects the gradient of G to the dual variable ν as: ∇G(q k+1 ) = B T ν k+1 . This relation along with the assumption that G is L q -smooth and B is positive definite allow us to have the following relation: With (20) and combining the above results, we obtain the sufficient decrease lemma (Lemma 8): Then through elementary quadratic programming, we can determine the range in terms of the penalty coefficient c that makes the coefficient pre-multiplied by q k − q k+1 2 become positive. This proves the sufficient decrease lemma. Then we proceed to prove the linear rate of convergence. Similar to the proof sketch for Theorem 1, given the sufficient decrease lemma, the next step is again to show that the inequalities (14) and (15) hold. From the minimizer conditions and the assumption that A is full row rank, we develop Lemma 17, which gives: where K 5 , K 6 > 0. Then we can combine the above two inequalities: where K 7 > 0 and in the last inequality we consider a neighborhood around w * such that w k+1 − w * < ε, L * c < L k+1 c < L * c + δ, and using Lemma 12 that assures the existence of ∇L k+1 This proves that the Łojasiewicz exponent θ = 1/2 locally around w * and (14) is satisfied. On the other hand, for (15), since A is full row rank, from the minimizer conditions (41), we have: where K 8 > 0. Then for the last term in the above inequality, suppose 1 ≤ α < 2, we can find an upper bound for Ap k+1 − Bq k+1 2 through (19), and obtain: where K 9 > 0. This proves the relation (15). Consequently, by Lemma 11 the rate of convergence for 1 ≤ α < 2 is Q-linear locally. For the region 0 < α < 1, we again prove that the rate is R-linear instead. From (22), (20) and the sufficient decrease lemma: where K 10 > 0. This shows that the relation (18) holds, and by Lemma 13, we conclude that {w k } k>N0 , N 0 ∈ N converges R-linearly. Together, we prove that the rate of convergence is locally linear for 0 < α < 2.
Compared to Theorem 1, while Theorem 2 relaxes some required properties on the sub-objective functions F, G, it needs the matrix B to be positive definite instead of A. The change is necessary to balance the weak convexity of the function G. However, for the Markovian optimization problem (6) we considered, because the linear constraints in fact reflect the marginal or Markov relations between the (conditional) probabilities, treated as primal and augmented variables, only one of the matrices A and B is an identity matrix whereas the other one is singular with full row rank. Inspired by the recent results [35], we develop the following theorem that keeps the definiteness of the matrix A as in Theorem 1 with relaxed properties as in Theorem 2 by imposing Lipschitz continuity on G. The additional mild assumption allows us to have the reverse norm bound: q m − q n ≤ M q Bq m − Bq n without B being positive definite and therefore balancing the weak convexity of G.
Theorem 3: Suppose F is convex and L p -smooth, G is σ G -weakly convex and L q -smooth, and M q -Lipschitz continuous, and the matrix A is positive definite while B is full row rank, then for 0 < α < 2, the sequence {w k } k∈N where w k := (p k , q k , ν k ), obtained from Solver II is bounded if: where λ B denotes the largest positive singular value of the matrix B; µ BB T the smallest positive eigenvalue of the matrix BB T . Moreover, the sequence converges to a stationary point w * := (p * , q * , ν * ) at linear rate locally.
proof sketch: The details are deferred to Appendix S. The only difference of the convergence analysis from Theorem 2 lies in the linear constraints where now A is positive definite and B is full row rank. For establishing the sufficient decrease lemma, and hence the convergence (Lemma 10), we now exploit the Lipschitz continuity of G with a coefficient M q to relate the norms Bq k − Bq k+1 , q k − q k+1 . Since the q-update can be equivalently expressed as a function Ψ(u) := arg min q∈Ωq G(q) + c/2 Bq − u 2 , as in [35], we have q m − q n = Ψ(Bq m ) − Ψ(Bq n ) ≤ M q Bq m − Bq n . Replacing the terms Bq k − Bq k+1 2 , or Bq k+1 − Bq * with this relation, we can prove the local linear rate of the region 0 < α < 2 similar to the proof of Theorem 2.
As a remark, while we focus on entropy and conditional entropy functions for applications, the convergence analysis for the algorithms hold for general functions satisfying the assumptions mentioned in theorems. This is closely related to the recent strongly-weakly convex pair problems [25], [34], which are still non-convex problems. This class of functions are less explored until recently in contrast to the well-studied convex-convex counterpart. For reference and comparison purposes, we summarize the results in Table I. In the next section, we apply the results to the IB and PF problems, two information theoretic, non-convex optimization problems that are difficult to solve. Nonetheless, based on the new results, we propose new IB and PF solvers and simplify the design of existing ones [13]. Interestingly, the proposed new PF solvers are capable of exploring the information plane more than the existing greedy solvers and achieve lower privacy leakage than theirs [2], [5]. Furthermore, we note that our results apply to non-convex problems beyond the IB and PF, as long as the outlined properties and conditions are satisfied.

IV. APPLICATIONS
In the section, we apply the results stated above to practical problems. In the general Markovian Lagrangian framework (6), we focus on two specific non-convex information theoretic optimization problems, i.e., the IB and PF problems.
When applied to IB and PF, we consider vector variables whose elements are composed of the vectorized, discrete (conditional) probability mass, defined as follows: where N v := |V| denotes the cardinality of a variable V. In both IB and PF, the variable to optimize is the conditional probability p z|x . Therefore, one of the primal variables p, q in (9) must be p z|x . As for the other one, it can be assigned as p z , p z|y or formed by cascading p z and p z|y . the relation between the conditional and marginal probabilities becomes a linear penalty Ap − Bq with each row of A, B being a prior probability vector. For example, if we let p := p z and q := p z|x then A := I, B := I ⊗ p T x where ⊗ denotes the Kronecker product, then the linear penalty term penalizes the case where the marginal probability relation x p(z|x)p(x) = p(z) is violated.
Our results depend on the convexity and smoothness of the two sub-objective functions. For the convexity, the negative entropy function −H(X) is convex w.r.t. the probability p x . Similarly, for the negative conditional entropy with a known marginal, it is also convex [26]. Additionally, we use the following definition to establish smoothness conditions for the (conditional) entropy functions and find the associated Lipschitz coefficients.
Definition 2: A measure u(x) is said to be ǫ-infimal if there exists ǫ > 0, such that inf x∈X u(x) = ε. The ε-infimal assumptions are commonly adopted in density/entropy estimation problems [27], [28] for smoothness conditions that will facilitate the optimization process. A key observation used in this work is that the entropy function whose associated probability mass vector is an ε-infimal measure, is weakly convex with a coefficient proportional to 1/ε. In addition, under the Markov chain Y → X → Z, by data-processing inequality we can have a tighter bound than the weak convexity, which is defined next as the restricted weak convexity.
It is worth noting that the restricted convexity was adopted in deriving the privacy parameter in differential privacy [42]. As will be shown below, it relaxes the restrictions on the linear constraints.

A. Solvers for Information Bottleneck
As shown in Table I, our results is based on the (strongly) convex-weakly convex structure of the two sub-objective functions. Interestingly, the IB problem satisfies these conditions. Following the known result in IB, owing to [43] which states that when the trade-off parameter γ ≥ 1, the corresponding minimum loss for the IB Lagrangian (2) are trivial, e.g. I(Z; X) = I(Z; Y ) = 0. Based on this result, we can focus on the region 0 < γ < 1 accordingly, which implies that (γ − 1)H(Z) is strongly convex w.r.t. p z and −γH(Z|X) is convex w.r.t. p z|x . Therefore, for the IB problem where the associated coefficient set in (6) consists of ρ z = γ − 1, ρ z|x = −γ and ρ z|y = 1, there exist multiple combinations in terms of the variables p z , p z|x , p z|y to apply our framework. In particular, we propose two splitting methods for the IB problem that work on Solver I and Solver II respectively. By imposing p k z|y = Q x|y p k z|x as an equality constraint, we can have the following re-formulation of the IB problem to the proposed framework (9) as: where the matrix Q x|y is defined such that p(z|y) = x p(z|x)p(x|y). Specifically, if we represent the conditional probability p(x|y) as a matrix W x|y ∈ R Nx×Ny with the (i, j)-entry w ij = p(x i |y j ), then Q x|y := I Nz ⊗ W T x|y , where ⊗ denotes the Kronecker product. Similarly, we have Q x := I Nz ⊗ p T x that maintains the marginal relation between p z and p z|x . Moreover, we can further treat the Markov relation p z|y = Q x|y p z|x as an additional penalty which gives an alternative form to (23). This corresponds to the general Lagrangian (6) with the following settings: where N q := |Z| × (|Y| + 1). Interestingly, we find that the first formulation (23) satisfies all the assumptions of Theorem 1 while the second one (24) meets the assumptions of Theorem 2, we therefore have the following convergence guarantee for solving the two forms of IB with the proposed splitting methods.
Theorem 4: Suppose p z is ε z -infimal and p z|x is ε z|x -infimal, then for 0 < α ≤ 2, the IB problem formulated as in (23) satisfies: Moreover, the sequence {w k } k∈N , where w k := (p k z , Q x|y p k z|x , ν k ) converges at linear rate locally to a stationary point when solved with Solver I with a penalty coefficient: Proof: See Appendix T. Note that the definition for ǫ-infimal measures corresponds to the smoothness assumption adopted in density/entropy estimation research [27], [28]. Specifically, it assumes that the minimum of a probability mass is bounded away from zero with ε > 0. Similarly, by applying Theorem 2 to (24), we have the following result.
Theorem 5: Suppose p z , p z|y , p z|x are ε z , ε z|y , ε z|x -infimal, respectively, then for 0 < α < 2, the IB problem formulated in (24) satisfies: Moreover, the sequence {w k } k∈N , where w k := (p k z|x , q k , ν k ) converges at linear rate to a stationary point when using Solver II with a penalty coefficient satisfying: where L q := 1/ε q with ε q := min{ε z , ε z|y } and σ G := (2N z N y )/ε z|y . Proof: See Appendix U. In literature, the first form (23) is proposed in our earlier work [13], but it is limited to α = 1. While the convergence is proved therein, an additional Bregman divergence is added to the p-update to regularize it, so that the convergence is assured in theoretical analysis. By Theorem 4 we show that the convergence can be proved without additional regularization. On the other hand, the second form a new in splitting methods for IB to our knowledge. Theorem 4 and Theorem 5 allow us to compare the two DRS algorithms in terms of the smallest penalty coefficient that assures convergence.
From the formulation and convergence analysis, the advantages of each solver are clear. For Solver I, it has fewer augmented variables to optimize since it is restricted to solutions where p z|y = Q x|y p z|x holds. On the other hand, for Solver II, the smallest penalty coefficient that assures convergence c * II is independent of L p and σ F , except for the case where p(x|y) is deterministic. The independence is useful in evaluating the solutions for the IB problem on the information plane [3], [44]. A common practice to form the information plane is to vary the trade-off parameter γ over a certain range. Since varying γ does not change the weak-convexity coefficient of Solver II. This invariance therefore allows us to fix a penalty coefficient c II when collecting IB solutions from varying γ.
We further compare the convergence rate of the proposed two IB solvers to existing ones. By Theorem 4 and 5, the convergence rates are locally linear. Remarkably, the BA-based algorithm, often serves as a benchmark, is also known to be convergent with linear rate [45]. Moreover, empirically the two solvers can obtain solutions with tantamount performance on the information plane.
Remarkably, our theoretic results can extend further than solving the IB and PF Lagrangian. Inspired by the variational inference method on IB [10], we apply our solvers to a surrogate bound of the IB Lagrangian: Observe that the construction and the associated conditions are essentially the same as required in Theorem 5, but now the p z|y is viewed as a function of p z|x and the weakly convex sub-objective reduces to γH(Z) only, hence the convergence results implied from Theorem 5 apply to this related solver.

B. Solvers for Privacy Funnel
Our general framework includes the PF problem as a special case as well. We can decompose the PF problem into a convex can be shown to be weakly convex under some mild smoothness conditions. The decomposition allow us to have the following augmented Lagrangian for the PF problem: where we restrict the marginal relation p(z) = x p(z|x)p(x) to be an equality constraint. Observe that F is convex w.r.t. p z|y whereas G is weakly convex (Lemma 5) if p z|x is ε z|x -infimal. Under these conditions, we can solve the PF problem with Solver II. Moreover, we have the following convergence guarantee for the the new PF solver.
converges at linear rate locally to a stationary point when using Solver I with a penalty coefficient: Proof: See Appendix V. In literature, most PF solvers are based on the agglomerative clustering approach [3] which is restricted to deterministic mappings [2], [5]. In contrast to these works, the new proposed PF solver can recover solutions obtained by the clustering based PF solvers on the information plane. Moreover, as shown in our numerical results, we can achieve lower privacy leakage compared to them.
Similar to the application of the proposed methods to variational inference-based IB, we can solve a surrogate loss upper bound of the PF Lagrangian (3), obtained through the variational inference technique: is the variational distribution. This corresponds to the following construction, satisfying the requirements stated in Theorem 5: Note that F is a convex function w.r.t. p = p z|x because the variational decoder q(x|z) is fixed during p, q updates. In contrast to existing variational PF solvers whose decoder q(x|y, z) depends on the sensitive information Y [14], the variational decoder in (27) can be optimized without passing Y as required input to the decoder.

V. EVALUATION
In this section, we present simulation results for the proposed algorithms using synthetic and real-world datasets. We implement the algorithms with gradient descent to update the primal variables p z , p z|x , p z|y . To ensure that the updated variables remain valid probability vectors, projected is needed [39]. There are various ways to project the updated variables to a probability simplex, the one we implemented is known as the mean-subtracted gradient as we empirically find that this method is more efficient in the context of the linear penalty constraints in (6). The mean-subtracted gradient is given by: where ξ k is a sufficiently small step-size at step k. This method introduces an extra parameter, the step-size to decide where we use the standard back-tracking line-search method to decide [39].
All the proposed solvers are initialized as follows, we use Python Numpy package to randomly sample from a Unif(0, 1) source |Z| × |X | times and arrange them in to a R |Z|×|X | matrix. Then the entries are normalized. The main focus of the evaluation is the characterization of the information plane of either an IB or PF solver

A. Datasets
The synthetic conditional distribution used in our evaluation is given by: Additionally, we evaluate the performance with the following non-uniform p(X).
We set the representation dimension |Z| ≤ 4 for both the IB and PF [46]. We also evaluate the proposed methods on a real-world dataset. The dataset is named "Heart failure clinical records Data Set" [47] from the UCI Machine Learning Repository [48]. It has 299 instances with 13 attributes. Among which, we select 6 attributes including: "anaemia," "high blood pressure," "diabetes," "sex," "smoking" and "death". All selected attributes are binary. We let Y := {"sex","death"} and the rest be X , this results in |Y| = 4, |X | = 16. As for the cardinality of |Z| ≤ 17 due to [46]. The joint probability is formed by counting the 299 instances w.r.t. (Y, X) pair. We post-processed the counted results by adding 10 −3 to each entry to avoid p(x, y) = 0.

B. Privacy Funnel
The proposed PF solver is denoted as Solver II consistent with our earlier notation. The corresponding variational inference variant of this solver (28) is denoted as Solver II-V. We initialize the encoder p(z|x) as described in Section V, obtaining a feasible point in the compound simplex Ω z|x . Since p(X, Y ) is assumed to be known, computing both p(z) = x p(z|x)p(x) and p(z|y) = x p(z|x)p(x|y) is straightforward.
For the synthetic dataset, the range of the trade-off parameters is set as β ∈ [0.5, 20.0] and we generate 20 geometrically spaced grid points of the range. For a β within this range, 10 trails are performed where the termination conditions are: 1) convergent when total variation is lower than a pre-determined small threshold, i.e., set to Ap − Bq ≤ 2 × 10 −6 or 2) divergent when a maximum number of iteration is reached otherwise. In the convergent case, the resulting encoder probability p z|x can be used to compute the mutual information pair I(Z|X), I(Z|Y ). We collect the convergent cases, calculate the  Fig. 1a compares the two proposed PF solvers whereas. Fig. 1b, compares the best proposed solver with the clustering-based PF solver. In Fig. 1a, we compare the information plane of two proposed solvers on the synthetic dataset with the uniformly distributed marginal p unif (X). The first solver is denoted as Solver II and the relaxation parameter is set to α = 1.618. This solver is based on Theorem 6. The second proposed solver is the variational inference-based solver (27) where a surrogate upper bound to the PF Lagrangian (3) is solved through Solver II with α = 1.000. We denote the method as Solver II-V. As shown in Fig. 1a, under the same range of the trade-off parameter β and the same number of trials the Solver II performs better (i.e., achieves a lower information leakage). In Fig. 2a, we further consider the non-uniform marginal p non-unif (X) where it is shown that Solver II-V is better instead (refer to the I(Z; Y ) in the range I(Z; X) ∈ [0.8, 1.2] bits).
Then we evaluate the proposed solvers on the real-world dataset. In this experiment, the trade-off has a range β ∈ [1.0, 10.0] and 20 trails are performed. The results are shown in Fig. 4. In Fig. 3a, we compare the two proposed methods. For Solver II the penalty coefficient is tuned to 7000 while for the Solver II-V the penalty coefficient is tuned to 128. In this experiment setup, Solver II is found to perform better. Note that both solvers achieve the "perfect privacy" [4] region (i.e., the utility I(Z; X) > 0 while I(Z; Y ) ≈ 0).
We further examine the convergence behavior of Solver II. Our theoretical results imply that the penalty coefficient needs to be larger than a threshold to ensure convergence. In the real-world dataset, we fix a trade-off parameter β = 3.5 and sweep through a range of penalty coefficient c ∈ [1000, 8000]. For each c, we perform 2000 trials and calculate the percentage of convergent cases. The results are shown in Fig. 4a, where we compare the two cases α = 1.618 and α = 2.000. We observe that to achieve 80% of convergent cases, the smallest c for α = 2.000 is lower than that of α = 1.618, which aligns with our theoretical predictions. Lastly, we examine the rate of convergence in Fig. 4b. We compare the two cases α = 1.618 and  Fig. 3a reports the convergence behavior of the proposed solvers in a real-world dataset. α = 2.000 and fix the penalty coefficient to c = 7000. The minimum loss isL * c = −2.46. We have established that the rate of the Solver II is locally linear and Fig. 4b provides numerical evidence supporting our theoretical result.
In Fig. 1b, we compare our new solver with the state-of-the-art clustering based algorithm (referred to as Merge-Two) under the assumption of a synthetic dataset with p unif (X). We observe that in the range I(Z; X) ∈ [0, 0.9] bits, our solver obtains more points on the information plane than Merge-Two and these points have lower privacy leakage compared to Merge-Two. However, the proposed solver converges to a local minimum at I(Z; X) ≈ 0.9. This can be improved through a more optimized implementation which is the subject of our future work. In 2a we repeat the comparison with non-uniform p non-unif (X). Solver II-V is used here since it provides the best solution. Again, we observe that Solver II-V recovers the solutions of Merge-Two around I(Z; X) ≈ {0.5, 1.0, 1.3} bits and achieves lower privacy leakages other wise. Finally, Fig 3b reports the comparison with the real-world data set where the same trend is observed.

C. Information Bottleneck
We adopt the same numerical setup for the synthetic dataset with both uniform and non-uniform marginal probabilities p unif (X) and p non-unif (X). The trade-off parameter γ ∈ [0.1, 1.0] and 16 geometrically-spaced grid points are evaluated. For each γ, 16 trials are performed. Each trial is initialized as described in V. The same convergence criterion for the proposed PF solvers is adopted here. For the information plane, only convergent cases are considered when characterizing the relevancecomplexity trade-off. For the proposed solvers, we denote Solver I for (23) and Solver II for (24). The proposed variational inference-based solver in (25) is denoted as Solver I-V.
We first evaluate the proposed solvers on p unif (X) in Fig. 5a with c I = 16, α I = 1.618, c II = 64, α II = 1.000 and c I-V = 64. The proposed solvers mostly obtain comparable Pareto-frontier solutions but we observe that Solver I-V converges to a local minimum at I(Z; X) ≈ 0.8bits whereas Solver I and Solver II do not. Then we evaluate the three solvers on the non-uniform p non-unif (X) in Fig. 6a. We observe that for I(Z; X) ∈ [0.2, 0.6] bits, Solver II converges to slightly sub-optimal solutions while in I(Z; X) ∈ [0.85, 1.2] bits Solver I-V converges to sub-optimal solutions, and hence, Solver I provides the best performance. We evaluate the convergence performance of the proposed solvers on the synthetic dataset with uniform p unif (X). In this simulation, the trade-off parameter is set to γ = 0.20. Each solver starts from a randomly initialized point (the same method as in the last part) for 100 trials. In Fig. 7a Solver I is evaluated with two settings α = 1.000 and α = 2.000. We observe that to reach 90% convergent percentage, α = 2.000 requires a smaller penalty coefficient compared to α = 1.000 which aligns with our theoretical convergence analysis (Theorem 4). Similar observations can be found in Fig. 7b where Solver II is configured with α = 1.000 and α = 2.000. Clearly, the case α = 2.000 requires a smaller c for to reach the same convergence percentage. This aligns with Theorem 5.
In Fig. 7c, we compare the convergence behavior of Solver I and Solver II. We fix γ = 0.20 and initialize both solvers from the same point. We run the two solvers until convergence and compare their loss decrease. We report the fastest (lowest number of iterations) configurations of each solver. In this specific case, the minimum loss among the two solvers is L * = −0.324. As shown in the figure, both solvers first explore the loss surfaces in sub-linear rate, then when the solvers operate within a neighborhood of a stationary point, then they converges to it linearly. Finally, Solver I is shown to require a smaller penalty coefficient than Solver II.
In Fig. 5b we compare the proposed solvers to BA with the synthetic dataset. The range of the trade-off parameter is γ ∈ [0.1, 1.0] and 16 geometrically-spaced grid points are generated from this range. For each γ, 16 trails are performed. In the beginning of each trail, p(z|x) is initialized as described in Section V. The convergence criteria follow Section V-B. We observe that Solver I can identify more points on the Pareto-frontier compared to in the range I(Z; X) ≈ 1 bits and I(Z; X) ∈ [0, 0.4] bits. To examine this further, we compare the two solvers in the synthetic dataset with a non-uniform  marginal p non-unif (X) in Fig. 6b where it is shown that some points obtained by BA are local minima on the information plane, whereas Solver I is not trapped at these sub-optimal solutions. Furthermore, in I(Z; X) ∈ [0, 0.4] the proposed solver again explore the relevance-complexity trade-off of this synthetic dataset better. On the other hand, in certain ranges the BA is shown to outperform the proposed solvers (slightly higher relevance). This is because the convergence of the proposed solvers will assure a feasible solution to the IB Lagrangian (2) but are not necessary the same point obtained by the BA.
In Fig. 8 compares the number of iteration to convergence versus the range of trade-off parameter γ ∈ [0.1, 1.0]. Here BA provides the best performance which is largely due to the implementation sub-optimality of the new solvers. This motivates exploring more efficient implementation, within the same splitting methods framework, which is an interesting topic for future work.
Finally, we evaluate the proposed IB solvers on the real-world dataset. The result in shown in Fig. 9. The trade-off parameter range is γ ∈ [0.01, 1.0] and 64 geometrically-spaced grid points are generated from the range. Then for a γ in the range, 16 trails are performed by each algorithm. We collect the obtained solutions and plot the convergent cases on the information plane. The proposed method is Solver I-V because we empirically found that it performs best among the proposed solvers. Compared to BA, the proposed solver can span the information plane with more points on the Pareto-frontier (observe I(Z; X) ∈ [0, 0.5] bits).

VI. CONCLUSIONS
In this work, we considered a general discrete rate distortion Lagrangian following a three letter Markov chain. The general framework includes the IB and PF problems as special cases. We proposed solving the general problem with splitting methods that are capable of handling large-scale problems which include important applications in multi-view learning [49], [50] and multi-source privacy problems [51]- [54].
Our convergence analysis is general for any objective function that can be decomposed to a convex-weakly convex pair. We further proved that our proposed algorithms are linearly convergent. Based on these theoretical insights, we developed optimized new solvers for both the IB and PF problems. For the two classes of the developed IB solvers, the first class has fewer variables to optimize by restricting the Markov relation to hold strictly while the second class is convergent independent of the selection of the trade-off parameter controlling the relevance-compression trade-off (except for one special case). In the PF case, our new solvers are shown to outperform the state of the art clustering-based solvers. Our empirical evaluations include synthetic and real-world data sets and explored both uniform and non-uniform priors.
For future work, we plan to extend the proposed framework to the continuous setting which is still an open challenge where only special cases are known [10], [55]. Another direction is multi-view learning via deep neural networks [56]- [60] where splitting methods can shed light on solver architectures with better parallelism and efficiency [61], [62].

APPENDIX A CONVERGENCE ANALYSIS
In this section, we prove the convergence and the corresponding rates for Solver I and Solver II. We start with the preliminaries including definitions and properties that will be used in the following proofs.

A. Preliminaries
where L > 0 is the Lipschitz coefficient.
Note that if f ∈ C 1 and ∇f (x) is L-Lipschitz continuous, then the function f is said to be a L-smooth function.
Definition 2: A measure u(x) is said to be ǫ-infimal if there exists ǫ > 0, such that inf x∈X u(x) = ε. In other words, if a distribution is ǫ-infimal then the smallest mass is strictly bounded away from zero by a positive constant ǫ. The infimal measure is commonly assumed in non-parametric entropy/density estimation for smoothness of the estimators [27], [28].
Lemma 1: let f (u) = |µ| i=1 µ i log µ i be the negative entropy function where two distinct measures µ, ν are ǫ-infimal. Then f is | log ǫ|-Lipschitz continuous and 1/ǫ-smooth Proof: The Lipschitz continuity follows as: As for smoothness: where the inequality is due to the following identity and the fact that log x < x − 1 for x > 0: We can establish similar smoothness condition for the conditional entropy. Corollary 1: Let p x be given, p z|x be ǫ-infimal, then the conditional entropy H(Z|X) = − x p(x) z p(z|x) log p(z|x) is | log ǫ|-Lipschitz continuous and 1/ǫ-smooth.
Proof: Following lemma 1, for two measures u, v ∈ Ω z|x , where Ω z|x denotes a compound simplex for the conditional probability p(z|x), the Lipschitz continuity follows as: On the other hand, to prove the smoothness, similar to the r.h.s. of the inequality (32), we have:

Definition 3:
A differentiable function f : R n → [0, ∞) is said to be σ-hypoconvex, σ ∈ R if the following holds: If σ = 0, (33) reduces to the definition of convex function; σ > 0 corresponds to strong convexity whereas when σ < 0, it is known as the weak convexity [20], [25], [34]. A well-known example is the negative entropy function, which is 1-strongly convex in 1-norm [26], and consequently in 2-norm. Another example is the conditional entropy, which is weakly convex if the corresponding conditional probability mass is ε-infimal as shown in the follow lemma.
Lemma 2: Let G(p z|y ) = H(Z|Y ). If p z|y is an ǫ z|y -infimal measure, then the function G is (2N z N y /ǫ z|y )-weakly convex, where N z = |Z|, N y = |Y| denote the cardinalities of the random variables Z, Y , respectively.
Proof: See Appendix C. A closely related concept to hypoconvexity that we called restricted weakly convexity is defined as follows: The restricted-weak convexity property is adopted in our earlier work [13] to prove the convergence of an ADMM solver for IB. We further extend the application of restricted weak convexity to prove the locally linear rate of convergence for nonconvex splitting methods. This is based on the observation that if a function f is composed of a combination of the negative conditional entropy and positive marginal/conditional entropy functions where the conditional measure corresponding to the negative conditional entropy function is the primal variable and the marginal/Markovian conditional measure induced by it are augmented ones, then the deviation of f from a convex function can be lower bounded by the total variation of the augmented measures. Lemma 3: Assume p z|x is ǫ z|x -infimal. Let G(p z|x ) := −γH(Z|X) + H(Z|Y ) and Y → X → Z forms a Markov chain. If 0 < γ < 1, then for two p m z|x , p n z|x ∈ Ω z|x , where Ω z|x := {p(z|x)| z p(z|x) = 1, ∀z ∈ Z, x ∈ X }, G(p z|x ) is ω G -restricted weakly convex.
Proof: see Appendix D. Beyond smooth functions, if in addition, convexity applies, then we have the following descent lemma, commonly used in first-order optimization methods [24], [63], [64].
Lemma 4 (Theorem 2.1.12 [64]): If f : R n → [0, +∞) is σ-strongly convex and L-smooth, then for any x, y, the following holds: A recent result generalized the above to σ-hypoconvex f which can be found in the reference therein [20]. Under ǫ-infimality, the following results show that the positive entropy function is weakly convex. Lemma 5: Given p x , let G(q) be defined as in (26).
Proof: See Appendix B. The following elementary identities are useful for the convergence proof. We list them for completeness.
Lastly, by "linear" rate of convergence, we refer to the definition in [39].
Definition 5: Let {w k } be a sequence in R n that converges to a stationary point w * when k > K 0 ∈ N. If it converges Q-linearly, then ∃Q ∈ (0, 1) such that On the other hand, the convergence of the sequence is R-linear if there is Q-linearly convergent sequence {µ k }, ∀k ∈ N, µ k ≥ 0 such that:

B. Kurdyka-Łojasiewicz Inequality
In the main contents, we assume that p z is ǫ z -infimal while p z|x is ǫ z|x -infimal. Following this, we can adopt the standard alternating direction method of multiplier (ADMM) [16] to prove the convergence of the proposed solvers by showing that the corresponding augmented Lagrangian satisfies the Kurdyka Łojasiewicz (KŁ) property. [22]- [24]. Moreover, the rate of convergence can be determined in terms of the Łojasiewicz exponent of the augmented Lagrangian. This section is provided as a brief review of this tool.
The convergence for convex splitting methods are well-studied and has been applied to a variety of algorithms [37], [64]. Interestingly, recent works found that splitting methods are convergent in solving a rather broad class of non-convex functions and give remarkable performance [35], [62]. However, compared to their convex counterpart, the fundamental understanding for non-convex splitting methods is less addressed until recently [20], [23], [34], [35], [65]- [67].
The KŁ inequality is a generalization of the well-known Łojasiewicz inequality to potentially non-smooth functions. But even in the smooth objective function class, the convergence of splitting methods is characterized through the Łojasiewicz inequality [22].
Definition 6: A function f (x) : R |X | → R is said to satisfy the Łojasweicz inequality if there exists an exponent θ ∈ [0, 1), δ > 0 and a critical point x * ∈ Ω * with a constant C > 0, and a neighborhood x − x * ≤ ε such that: Leveraging the main result in [22], the rate of convergence of splitting methods can be determined immediately by the associated Łojasiewicz exponents θ if each of the sub-objective function satisfies the KŁ property.
where ∂f denotes the sub-gradient of f (·) for non-smooth functions and reduces to gradient ∇f for smooth functions. Clearly, if ϕ(s) = Cs 1−θ , then (38) reduces to the Łojasiewicz inequality. The attribution to Kurdyka is due to the discovery of a variety class of practical functions satisfying the KŁ property, where the class of functions is said to have the o-minimal structure [41], i.e., sub-analytic and semi-algebraic functions. Once satisfying the KŁ property and knowing the exponent θ, the objective function, if solved with splitting methods, has the corresponding rate of convergence characterized depending on the value of θ. In particular, the most relevant case in the sequel, if the exponent θ = 1/2 then the rate of convergence is locally linear around a neighborhood of a stationary point [22].
While the Łojasiewicz exponents of the o-minimal function class are often easy to calculate [23], [68], for more general functions, the exponents are difficulty to determine [40], [69].
Recently, since the application of the KŁ inequality for convergence analysis of splitting methods in non-convex settings [23], [24], a wealth of optimization mathematics research has devoted to characterizing the convergence conditions under a assumed structure of the non-convex objective function. Among which, the one that is most relevant to ours is the (strongly) convexweakly convex structure [19], [20], [25], [34]. The main discovery is that under this structure the convergence is assured if the penalty coefficient c is sufficiently large, characterized by the Lipschitz smoothness and properties of the operators in the linear constraints. We refer to [35] for a summary of convergence conditions of non-convex ADMM and [20] for Douglas-Rachford splitting (DRS) [21] and the references therein for recent advances.

C. Proof of Convergence
In proving the convergence of the two algorithms, we consider three different sets of assumptions. We start with the most restricted one paired with Solver I: Assumption A: • There exists a stationary point w * := (p * , Bq * , ν * ) that belongs to a set Ω * := {w|w ∈ Ω, ∇L c = 0}.
is L q -smooth and ω G -restricted weakly convex.
• A is positive definite.
• The penalty coefficient c > c min , where c min is defined as: A ]/α}. We consider first-order optimization methods for (10), which gives the following minimizer conditions: Note that at a stationary point (p * , q * , ν * ), the above reduces to: With the minimizer conditions, we present a sufficient decrease lemma for Solver I. Lemma 6 (Sufficient Decrease I): Let L c be defined as in (9) and Assumption A is satisfied, then with Solver I, we have: where the coefficients δ p , δ q , δ ν are defined as: where µ A denotes the largest eigenvalue of the positive definite matrix A.
Proof: See Appendix E. By Lemma 6, the conditions that assure sufficient decrease are equivalent to the range of the penalty coefficient c and the relaxation parameter α such that δ p , δ q , δ ν are non-negative. When the conditions are satisfied, the sufficient decrease lemma implies the convergence of Solver I. Lemma 7 (Convergence I): Suppose Assumption A is satisfied and 0 < α ≤ 2. Define the collective point at step k as w k := (p k , Bq k , ν k ), then the sequence {w k } k∈N obtained from Solver I is convergent to a stationary point w * ∈ Ω * .
Proof: See Appendix F. As a remark, convergence is not point-wise. This can be observed as q in the collective point is pre-multiplied by the matrix B. In practice, take IB for example, this corresponds to the symmetry of solutions [70], [71]. Nonetheless, point-wise convergence is not necessary as the mutual information, the metric typically involved in the information-theoretic optimization problems, is symmetric.
Observe that in Assumption A, the function F is required to be strongly convex while G is restricted weakly convex, which limits the class of objective functions that our results can apply to. To relax this assumption, we consider Solver II instead and develop a sufficient decrease lemma with relaxed assumptions: Assumption B: • There exist stationary points w * := (Ap * , q * , ν * ) that belong to a set Ω * := {w|w ∈ Ω, ∇L c = 0}, • The function F (p) is L p -smooth and convex while G(q) is L q -smooth and σ G -weakly convex.
• B is positive definite and A is full row rank.
• The penalty coefficient c satisfies: The corresponding first-order minimizer conditions are: Lemma 8 (Sufficient Decrease II): Let L c be defined as in (9) and Assumption B is satisfied, then using Solver II, we have: Proof: See Appendix G. In parallel to Lemma 7, we have the following convergence result for Solver II.
Lemma 9 (Convergence II): Suppose Assumption B is satisfied and 0 < α < 2. Define w k := (Ap k , q k , ν k ) the collective point at step k. Then the sequence {w k } k∈N obtained from Solver II is convergent to a stationary point w * ∈ Ω * .
Proof: See Appendix H. Note that the convergence of Solver II requires no strong convexity for the sub-objective function F (p). Moreover, the assumption for G(q) is more relaxed than that of Solver I, and hence the results apply to wider class of functions.
Another major difference between Assumption A and B lies in the linear constraints. In Assumption A, A is positive definite while B is positive definite in Assumption B. In the Markovian information theoretic optimization problem we considered (6), the linear constraints Ap − Bq are in essence the marginal/Markov relations of (conditional) probabilities. Therefore, only one of the two matrices A, B is identity, while the other will be singular. Then for problems such as PF, whose convex sub-objective function is not strongly convex, with A being positive definite instead of B, neither the assumptions mentioned above hold. Inspired by [35], when F and G are further assumed to be Lipschitz continuous, we can relax Assumption A but keep A to be positive definite as in Assumption B.
Assumption C: • There exists a stationary point w * := (p * , q * , ν * ) that belongs to a set Ω * := {w|w ∈ Ω, ∇L c = 0}, • The function F (p) is L p -smooth and convex while G(q) is L q -smooth and σ G -weakly convex, • In addition, G(q) is M q -Lipschitz continuous, • A is positive definite and B is full row rank, • The penalty coefficient c satisfies: When the above assumptions are imposed on Solver II, which reuses the minimization conditions (41), we have the following sufficient decrease lemma.
Lemma 10 (Convergence III): Suppose Assumption C is satisfied and 0 < α < 2. Define w k := (p k , q k , ν k ) the collective point at step k, then the sequence {w k } k∈N obtained from Solver II is convergent to a stationary point w * ∈ Ω * .
Proof: See Appendix I.

D. Rate of Convergence Analysis
In this part, we show that the rates of convergence of the algorithms, under the three sets of assumptions discussed in the previous part, are all locally linear. Specifically, the linear convergence is independent of initialization and the sequence obtained from the two corresponding algorithms converges to local minimizers when the current update of the variables lies around their neighborhood [35], [65]. The results are based on the KŁ inequality that recently applied to characterize the rate of convergence for splitting methods in non-convex problems [22]- [24], [63]. The analysis consists of two steps. First we show that (9), solved either with Solver I or Solver II, satisfies the KŁ property with a Łojasiewicz exponent θ = 1/2. Then due to the following result, owing to [22], [24], [40], we prove the linear convergence rate.
Lemma 11 (Theorem 2 [22]): Assume that a function L c (p, q, ν) satisfies the KŁ property, define w k the collective point at step k, and let {w k } k∈N be a sequence generated by either Solver I or Solver II. Suppose {w k } k∈N is bounded and the following relation holds: ∇L k c ≤ C * w k − w k−1 , where L k c := L c (p k , q k , ν k ) and C * > 0 is some constant. Denote the Łojasiewicz exponent of L c with {w ∞ } as θ. Then the following holds: (i) If θ = 0, the sequence {w k } k∈N converges in a finite number of steps, (ii) If θ ∈ (0, 1/2] then there exist τ > 0 and Q ∈ [0, 1) such that If θ ∈ (1/2, 1) then there exists τ > 0 such that Proof: See Appendix J. We only prove the case corresponding to θ = 1/2 as it is relevant to the following discussion. For the proof for other scenarios, we refer the reader to [22]. The above result characterizes the rate of convergence in terms of the KŁ exponent, but except for certain types of functions, the calculation of the KŁ exponent is difficult. The following key result, due to [40], is useful in calculating the KŁ exponent of (9) and is included for completeness.
In literature, the KŁ inequality has been successfully adopted to find the rate of convergence for alternating algorithms such as ADMM and recently PRS or DRS with α = (1+ √ 5)/2. For more general DRS methods in terms of the relaxation parameter α, we find that proving locally linear rate through the KŁ inequality only holds for 1 ≤ α ≤ 2. As for 0 < α < 1, inspired by the recent results that show locally R-linear rate of convergence for the primal ADMM [19], we adopt and extend the approach to Solver I and Solver II under the three sets of assumptions. Combining the two methods, we therefore theoretically prove that the rates are locally linear for 0 < α ≤ 2.
Lemma 13: Let L c be defined as in (9) and let the sequence {w k } k∈N obtained through either Solver I or Solver II be bounded. Denote L k c := L c (p k , q k , ν k ). Suppose the following holds for some K * > 0: and there exists a neighborhood around a stationary point w * , such that w − w * < ǫ, L * c < L c < L * c + δ with δ, ǫ > 0. Then {L k c } k∈N is Q-linearly convergent and {w k } k∈N converges R-linearly to w * around the neighborhood. Proof: See Appendix K. Remarkably, the rate of convergence with KŁ inequality is Q-linear, or in other words, monotonic convergence in terms of the error between variables w k − w k−1 in consecutive steps is guaranteed, while the R-linear rate is non-monotonic, hence a weaker rate. However, the weaker R-linear rate comes with milder assumptions imposed on the linear constraints, in particular, the full row rank assumptions are lifted.
In the rest of this part, we aim to prove that the sequence {w k } k∈N obtained from any of the proposed two algorithms satisfies the KŁ property. The results are based on the following lemmas. We start with a lemma developed for Solver I. Lemma 14: Let L c be defined as in (9). For the sequence {w k } k∈N obtained from Solver I where w k := (p k , Bq k , ν k ), if it is bounded and converges to a stationary point w * satisfying (40), then we have: and ∇L c (w k+1 ) ≥ (c 2 µ 2 A + 1) Ap k+1 − Bq k+1 , where L k c := L c (p k , q k , ν k ); λ A , µ A denote the largest and smallest eigenvalue of a positive definite matrix A, respectively. Proof: See Appendix L. Lemma 15: Suppose Assumption A is satisfied, if the augmented Lagrangian (9) is solved with Solver I, then it satisfies the KŁ property with an exponent θ = 1/2.
Proof: See Appendix M. Given the exponent θ = 1/2, by mapping the exponent according to Lemma 11, we show the linear rate of convergence, which extends the convergence (Lemma 7) to the following result.
Theorem 1: Suppose Assumption A is satisfied. For 0 < α ≤ 2, define w k := (p k , Bq k , ν k ) the collective point at step k. Then the sequence {w k } k∈N obtained from Solver I is bounded. Moreover, the sequence converges to a stationary point w * at linear rate locally.
Proof: See Appendix N. Similarly, for Solver II, we show that the Łojasiewicz exponent of the corresponding augmented Lagrangian is θ = 1/2. However, this requires an additional assumption that A be full row rank. We later show that this additional assumption is not necessary to prove locally linear convergence rate in an alternative approach.
Lemma 16: Let L c be defined as in (9). For the sequence {w k } k∈N obtained from Solver II, where w k := (Ap k , q k , ν k ), if the sequence is bounded and converges to a stationary point w * satisfying (40), 0 < α < 2, then we have: where L k c := L c (p k , q k , ν k ). Moreover, if AA T ≻ 0, then: where µ W denotes the smallest positive eigenvalue of a matrix W . Proof: See Appendix O. Lemma 17: Suppose Assumption B is satisfied and the matrix A is full row rank. For 0 < α < 2, if the augmented Lagrangian (9) is solved with Solver II, then it satisfies the KŁ inequality with an exponent θ = 1/2.
Proof: See Appendix P. Observe that in Lemma 17, an additional full row rank assumption is imposed on the matrix A. This is necessary to prove Q-linear rate of convergence with KŁ inequality. It turns out that we can relax this condition by showing local linear rate of convergence without assuming A to be full row rank, which is due Lemma 13.
The above lemma shows that the sequence {L k c } k∈N is locally Q-linear convergent, which in turns allows us to show locally R-linear rate of convergence of the sequence {w k } k∈N obtained from Solver II.
Theorem 2: Suppose Assumption B is satisfied and the sequence {w k } k∈N with w k := (Ap k , q k , ν k ) obtained from Solver II is bounded, then {w k } is R-linearly convergent to a stationary point w * locally around a neighborhood w − w * 2 < ǫ and L * c < L c < L * c + η for some ǫ, η > 0. Proof: See Appendix Q. Lastly, when imposing Assumption C on Solver II, we prove that the Łojasiewicz exponent in solving the augmented Lagrangian (9) is θ = 1/2 and apply the KŁ inequality to prove its linear rate of convergence. The key difference of the Assumption C is that G is required to be Lipschitz continuous, which allows us to have inequalities such as q m − q n ≤ M q Bq m − Bq n with B not necessarily being positive definite [35]. We first adopt this result to show θ = 1/2. Lemma 18: Suppose Assumption C is satisfied and the matrix B is full row rank. For 0 < α < 2, if the augmented Lagrangian (9) is solved with Solver II, then it satisfies the KŁ property with an exponent θ = 1/2.
Proof: See Appendix R. Theorem 3: Suppose Assumption C is satisfied. For 0 < α < 2, define w k := (p k , q k , ν k ) the collective point at step k. Then the sequence {w k } k∈N obtained from Solver II is bounded. Moreover, the sequence converges to a stationary point w * at linear rate locally.
Proof: See Appendix S.
Then for p-update, due to the σ F -strong convexity and using Lemma 4, we have: where the first inequality is due to σ F -strong convexity; the second is due to A being positive definite. Then, for the dual update, we have: Combining (46) and (48) using the identity (37), we get: Lastly, for the q-update, since G is ω G -restricted weakly convex w.r.t. the matrix B: Summing (47), (49) and (50), and using (39), we have: Then by positive definiteness of A, we have: Ap k − Ap k+1 ≥ µ 2 A p k − p k+1 , where µ A denotes the smallest eigenvalue of A. Substitute this into (51), then we complete the proof.

APPENDIX F PROOF OF LEMMA 7
By Assumption A, the coefficients δ p , δ q , δ ν defined in Lemma 7 are non-negative, so the next step is to show {L k c } k∈N is finite. Denote L k c := L c (p k , q k , ν k ) for simplicity. From the above, assume a penalty coefficient c * satisfying Assumption A, we have: where C * = min{δ p , δ q , δ ν } > 0. Define the collective point at step k as w k := (p k , Bq k , ν k ), then since there exist stationary points w * , the l.h.s. of (52) is lower semi-continuous. Let N → ∞ and denote the limit point w ∞ , since L 1 c − L ∞ c is finite, the r.h.s. of (52) is finite. This implies w k − w k+1 2 → 0 as k → ∞, since ∞ w k − w k+1 2 is a Cauchy sequence. From this, we know that w ∞ ∈ Ω * , or equivalently, for k > N 0 ∈ N sufficiently large, w k → w * as k → ∞, which proves that {w k } k∈N is convergent to w * . APPENDIX G PROOF OF LEMMA 8 First, by assumption, F is convex, hence: where the last equality is due to the minimizer conditions (41). Then for the relaxation step (11b): On the other hand, by assumption, G is σ G -weakly convex, so we have the following lower bound for q-update (11c): Lastly, for the dual ascend (11d): Combining (54) and (56) using the identity (37), we get: Summing (53), (55) and (57), we get: which completes the proof.
APPENDIX H PROOF OF LEMMA 9 By assumption, B is positive definite, denote its smallest eigenvalue µ B and L k c := L c (p k , q k , ν k ) for simplicity, we have: Note that B −1 = µ B . On the other hand, for the dual variable, we have: combining the two results above, we have the following lower bound to Lemma 8: Then, we would like to make the coefficient of the squared norm q k − q k+1 2 be positive. Observe that this is simply an elementary quadratic programming, and we therefore have the desired range of the penalty coefficient c as: which is satisfied as listed in Assumption B. Then, for any c satisfying (59), define c * := c/2, we have: Then, consider the following: Define the collective point at step k as w k := (Ap k , q k , ν k ), since there exists stationary points w * := (Ap * , q * , ν * ) by assumption, the l.h.s. of (60) is lower semi-continuous. Note that the the r.h.s of (60) does not depend on the dual variable ν, we can further define a condensed point at step k as z k := (Ap k , q k ). Then, by letting N → ∞ and denoting the limit point z ∞ , since L 1 c − L ∞ c is finite, the r.h.s. of (60) is finite. This implies z k − z k+1 2 → 0 as k → ∞, since ∞ z k − z k+1 2 is a Cauchy sequence. From this we know that z ∞ = z * . Moreover, due to (20), L 2 q µ 2 B q k − q k+1 2 ≥ ν k − ν k+1 2 and hence ν k − ν k+1 2 → 0 as k → ∞. Therefore ν ∞ = ν * . So, together we have w k → w * as k → ∞ which proves that {w k } k∈N is convergent to w * .

APPENDIX I PROOF OF LEMMA 10
Following the steps (53), (55) and (57) in Appendix G, we start from (58). Define L k c := L c (p k , q k , ν k ) the function value evaluated with variables at step k for simplicity: where in the last inequality, the first term is by A being positive definite, and for the second term, we follow [35] and use Lipschitz continuity of G to have q k − q k+1 ≤ M q Bq k − Bq k+1 ; we denote λ B := B as the largest positive singular value of a matrix B and µ B for the smallest positive eigenvalue of B; For ν k − ν k+1 , since B is full row rank and G is L q -smooth, we have: (62) From elementary quadratic programming, the range in terms of the penalty coefficient c that assures that the second term of the last inequality in (61) is positive: Then by assumption, c satisfies the above. Rewrite the coefficients as τ p , τ q > 0 for simplicity, then there exists a τ * := min{τ p , τ q } such that: Then denote w k := (p k , q k , ν k ) the collective point at step k; L k c := L c (w k ) the function value evaluated with w k . Summing both sides of the inequality (62), we have: By assumption the l.h.s. of the above inequality is lower semi-continuous and therefore is finite. So as N → ∞, L 1 c −L ∞ c < +∞. This implies the r.h.s. is finite and therefore p k − p k+1 2 → 0 and q k − q k+1 2 → 0 as k → ∞. Due to (62), we know that ν k − ν k+1 2 → 0 as well. Given the results, denote the limit points as w ∞ := (p ∞ , q ∞ , ν ∞ ), since w k − w k+1 2 → 0 as k → ∞, w ∞ = w * which proves that {w k } k∈N is convergent to w * .
APPENDIX J PROOF OF LEMMA 11 L c (p, q, ν) satisfies the KŁ property with an exponent θ = 1/2. Denote L k c := L c (p k , q k , ν k ), without loss of generality let L * c = 0, and define a concave function Φ(s) := C 0 s 1−θ with C 0 , s > 0. For k > N 0 ∈ N sufficiently large, by the concavity of Φ (Note that the gradient is evaluated at w k ): where the second inequality is due to Lemma 7 and the last inequality is due to Lemma 15. Then, by assumption, for some constant C * > 0, we have: Substitute the above into (63), define C 1 := C ′ /C * (1 − θ), we get: Substitute the above into (64), we have: for some K p , K q , K ν > 0 and 0 < Q < 1. Combine the above together, we have: ρ w n − w n+1 2 ≤KQ n , whereK = K p +K q +K ν andρ = min {ρ p , ρ q , ρ ν }. Now, for the sequence {w n } n∈N\[N0] around w * , by taking m > n ≥ N 0 , we have: Since the above is a Cauchy sequence, by taking limit with m → ∞, which gives w m → w * as m → ∞, we get: and therefore prove that {w n } n>N0 is R-linearly convergent.
APPENDIX L PROOF OF LEMMA 14 By the definition in (9) along with the properties of F and G, using Solver I with the first order minimizer conditions (39), and denote L k c := L c (p k , q k , ν k ) for simplicity, we have: where the first inequality is due to Lemma 4 and the restricted-weak convexity of G. Then, by letting w = w * , which gives Ap * − Bq * = 0, and using identity (36) for the second inner product in the last line of (67), we get: where λ A denotes the largest eigenvalue of the matrix A. Then, for the second part, consider: where the last equality follows from (39). By showing that ∇L k+1 c 2 ≥ K Ap k+1 − Bq k+1 2 with K := c 2 µ 2 A + 1 where µ A is the smallest eigenvalue of the matrix A, we complete the proof.

APPENDIX M PROOF OF LEMMA 15
By assumption, c > ω G and by Cauchy-Schwarz inequality u − v 2 ≥ (1 − 1/t) u 2 + (1 − t) v 2 , t > 1, consider: Denote L c (p k , q k , ν k ) := L k c , then for some t > 1, the following holds: Then define c G := max{0, c(1 − 1/t) + ω G /t − 2σ F L p /[λ 2 A (L p + σ G )]}/2, we have c G ≥ 0. Substitute the above into Lemma 14, we have: where K 1 := c 2 µ 2 A + 1 > 0; the first inequality is due to Lemma 14 and w k+1 − w * < ǫ around the neighborhood of w * ; the last inequality follows from Lemma 12. By taking square root of both sides we complete the proof. Substitute (71) into (70), using identities (36) and (37), and let w := w * , which gives Ap * = Bq * , we have: By identity (37) and the minimizer conditions (41): Substitute the above into (72), we have: and we complete the proof for the first part. For the second part, consider: Denote the smallest positive eigenvalue of a matrix W as µ W , by assumption, since AA T ≻ 0, we have: where in the last inequality, we use the minimizer condition (41) and L q -smoothness of G.