{'fraction_non_alphanumeric': 0.05607452112306481, 'fraction_numerical': 0.019653634216741014, 'mean_word_length': 3.9591411841249187, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical optimality and stability of tangent transform algorithms in logit models

A systematic approach to finding variational approximation in an otherwise intractable non-conjugate model is to exploit the general principle of convex duality by minorizing the marginal likelihood that renders the problem tractable. While such approaches are popular in the context of variational inference in non-conjugate Bayesian models, theoretical guarantees on statistical optimality and algorithmic convergence are lacking. Focusing on logistic regression models, we provide mild conditions on the data generating process to derive non-asymptotic upper bounds to the risk incurred by the variational optima. We demonstrate that these assumptions can be completely relaxed if one considers a slight variation of the algorithm by raising the likelihood to a fractional power. Next, we utilize the theory of dynamical systems to provide convergence guarantees for such algorithms in logistic and multinomial logit regression. In particular, we establish local asymptotic stability of the algorithm without any assumptions on the data-generating process. We explore a special case involving a semi-orthogonal design under which a global convergence is obtained. The theory is further illustrated using several numerical studies.


Introduction
Variational Inference (VI) has gained substantial momentum in recent years as an efficient way of performing approximate Bayesian inference. VI seeks to minimize a divergence measure between a tractable family of probability distributions and the posterior distribution, utilizing optimization based techniques to arrive at a minima. In many high dimensional examples where sampling based techniques such as the Markov chain Monte Carlo require expert vigilance and care for scalability, VI provides a viable answer with relatively lower computational cost. Some notable application areas include graphical models [Jordan et al., 1999, Wainwright et al., 2008, hidden markov models [MacKay, 1997], latent class models [Blei et al., 2003], neural networks [Graves, 2011] to name a few. Refer to Chapter 10 in Bishop [2006] and Blei et al. [2017] for excellent reviews on the topic.
The empirical success of VI has prompted researchers to investigate their theoretical properties.
Two distinct directions of research seem to have emerged over the last few years. One line of research concerns the statistical aspects of variational estimators [Alquier and Ridgway, 2020, Alquier et al., 2016, Chérief-Abdellatif and Alquier, 2018, Pati et al., 2018, Wang and Blei, 2019a,b, Yang et al., 2020, Zhang and Gao, 2020 in a general setting, delineating sufficient conditions on the data generation mechanism and the variational family under which the variational estimators have optimal first or second-order statistical properties. Motivated by the robustness properties of a fractional likelihood Ridgway, 2020, Bhattacharya et al., 2019], Yang et al. [2020] proposed a simple modification to the variational objective function, deemed as the α-Variational Bayes (α-VB), that only requires the variational family to be sufficiently flexible and the prior density to be appropriately concentrated around the true parameter to obtain optimal risk bounds.
The other line of research studies convergence of the algorithms employed to arrive at the variational optimizer. In this aspect, the coordinate ascent variational inference (CAVI) algorithm for mean-field VI (refer to Chapter 10 of Bishop [2006]) has arguably received the most attention due to its simplicity and generality. An early result on algorithmic convergence (and lack thereof) of CAVI in Gaussian mixture models appears in Wang and Titterington [2006]. Mukherjee et al. [2018], Zhang and Zhou [2020] analyzed CAVI for stochastic block models, a popular model for networks belonging to cEXP. Yin et al. [2020] obtained convergence of cluster labels in a stochastic block model by considering a structured variational family which was not possible using mean field VI. Ghorbani et al. [2018] noted instability of naive mean-field VI in latent Dirichlet allocation and provided a remedy by optimizing a different type of free energy (TAP) instead of the standard variational objective. Campbell and Li [2019], Locatello et al. [2018] analyzed convergence of a more flexible class of boosting algorithms which aim to approximate the target class by a mixture of Gaussians rather than a single Gaussian or a product distribution.
Our goal in this article is to explore a popular class of variational approximations outside cEXP, called the tangent-transform approach [Jaakkola, 1997, Jaakkola andJordan, 2000]. The tangent transform approach is an example of a structured variational approximation, lying on the spectrum between the two extremes given by the restrictive mean-field inference and the highly flexible variational boosting. In this specific instance, the structure exploited is convex duality [Jordan et al., 1999, Wainwright et al., 2005, Wainwright and Jordan, 2003] to minorize the log-likelihood function and provide sharp bounds for the log-partition functions in exponential family models.
Assume p(x, θ) is an exponential family on a discrete space X indexed by parameter θ ∈ Θ, p(x; θ) = exp{ θ, t(x) − B(θ)}, B(θ) = log x∈X exp{ θ, t(x) } . Ideas related to the tangent-transform have found widespread applications ranging from approximate inference in graphical models [Jordan et al., 1999], low-rank approximations [Srebro and Jaakkola, 2003], inference in large scale generalized linear models [Nickisch and Seeger, 2009], non-conjugate latent Gaussian models [Emtiyaz Khan et al., 2013] to more recently in sparse kernel machines [Shi and Yu, 2019], hierarchical relevance determination [Hirose et al., 2020], online prediction [Konagayoshi and Watanabe, 2019] among others. Jaakkola and Jordan [2000] exploits convex duality to minorize the marginal likelihood by introducing a variational parameter that allows the minorant to be arbitrarily close to the marginal likelihood. Logistic and multinomial logit regression models are notable examples where a clever use of this idea results in a straightforward Expectation-Maximization (EM) algorithm to compute the variational parameters.
In this article, we investigate both the statistical and algorithmic aspects of the tangent transform algorithm in logit models. Despite its widespread usage, statistical properties of the point estimate of the regression coefficients resulting from a tangent transform algorithm has not been previously studied. One possible reason is that unlike mean-field VI, where the global objective is to minimize the Kullback-Leibler (or another) divergence between a product distribution and the posterior distribution, the tangent transform algorithm is defined locally, without a clear global objective function that is being minimized. A key observation underlying our statistical analysis expresses any stationary point of the EM algorithm as a minimzer of a suitably chosen global variational objective function. This observation allows us to extend some previously developed variational risk bounds for mean-field VI [Pati et al., 2018, Yang et al., 2020 to the present setting with some non-trivial modifications. We show that with minimal assumptions on the data generating process and the prior density on the regression coefficients, the variational risk bound is minimax optimal (up to logarithmic terms). Moreover, the assumption on the data generating process can be completely relaxed by raising the standard logistic likelihood by a fractional exponent [Bhattacharya et al., 2019].
Next, we investigate the convergence of the EM algorithm to the fixed point of the EM iterations.
There has been some previous efforts to shed more light into the EM sequence of tangent-transform algorithms. Hunter and Lange [2004] studied connections between minimization-majorization (MM) in case of logistic likelihood to argue convergence of the updates. Durante and Rigon [2019] drew a connection with the Pólya-Gamma data augmentation technique [Polson et al., 2013] to provide a probabilistic interpretation of the EM updates and showed that the optimal evidence lower bound of the tangent transformation approach coincides with the same obtained in a bonafide variational inference with a suitably defined conditionally conjugate exponential family. However, statistical and algorithmic convergence guarantees of the tangent transform itself continue to be an open problem. It may appear on the surface that the EM algorithm underlying tangent-transforms can be analyzed using the general sufficient conditions for convergence of the EM (we refer to the recent article [Balakrishnan et al., 2017] and the references therein for more on this topic), a careful inspection however reveals that these general-purpose conditions pose significant difficulty to verify for the present EM iterates and demand stringent conditions on the design matrix and other data generating parameters. Our approach, on the other hand, is to directly analyze the EM sequence without resorting to any high-level results.
By viewing the EM updates as iterations in a discrete time autonomous dynamical system, we show that the EM updates converge to the desired fixed point under suitable initialization, a phenomenon known as local asymptotic stability. While stability is typically a weaker statement compared to global convergence as it only ensures convergence if the system is initialized in a neighborhood around the fixed point, our stability result is essentially assumption-free -it does not require any assumption on the design matrix, on the sparsity of the coefficients, and on the dimension p and the sample size n. Although the notion of such convergence is local, to the best of our knowledge, this is the first assumption-free result on the stability of a variational algorithm.
The main technical contribution is to show that the spectral radius of the Jacobian matrix of the linearized operator of the EM sequence is strictly smaller than one at the fixed point. In the special case when the design matrix is semi-orthogonal, we show that the EM sequence is globally convergent with an exponential rate of convergence (logarithmic run time) independent of the initialization. We also provide a straightforward extension of this result to the case of multinomial logit model.

Tangent Transformation Approach
Denote the data by X and the likelihood conditioned on parameter θ ∈ Θ by p(X | θ), where Θ is the parameter space. For a prior density p(·) on Θ, the goal of VI is to approximate the posterior p(θ | X) ∝ p(θ) p(X | θ) by a member of a tractable family Γ of densities on Θ with respect to the Kullback-Leibler (KL) divergence. 1 Notationally, VI seeks to find q = arg min q∈Γ D(q p(· | X)), (2.1) which is equivalent to maximizing the evidence lower bound (ELBO), L(q) = q(θ) log{p(X, θ)/q(θ)} dθ with respect to q ∈ Γ. Using a component-wise product structure on Γ popularly known as the mean field family [Parisi, 1988], closed-form updates of a coordinate ascent algorithm (CAVI) can be generally derived in conditionally conjugate exponential families [Blei et al., 2017]. However, many non-conjugate models such as logistic regression, multinomial logit regression, graphical and, topic models, do not lead to closed-form CAVI updates, necessitating various specialized techniques [Blei et al., 2017, Jordan et al., 1999. One such approach is to introduce variational parameters to minorize the log-marginal likelihood by a tractable family, which when combined with an appropriate prior enjoys conjugate inference. For Bayesian logistic regression models, Jaakkola and Jordan [2000] introduced a tangent-transform of the logistic function using convex duality. By a standard result in convex analysis [Rockafellar, 1970], a convex function f (·) on R d can be represented via a conjugate or dual function f * as, Geometrically, the evaluation of a convex function at any point x can be viewed as the maxima of The usage of duality is not restricted to linear approximations, i.e., hyperplanes. In fact, Jaakkola and Jordan [2000] used a quadratic bound for the logistic function that induces conjugacy with Gaussian priors. In the following subsection, we discuss the salient features of the tangent transform approach.

Convex minorant construction for logistic likelihood
We discuss a slightly general version of the tangent transform approach where we raise the usual logistic likelihood to a power α ∈ (0, 1] before combining with the prior. Variational Bayes procedures with fractional likelihoods have been recently considered in Alquier and Ridgway [2020], Alquier et al. [2016], Yang et al. [2020]. The case α = 1 recovers the usual tangent transform. Assuming we observe binary responses y i corresponding to fixed covariates x i ∈ R p (i = 1, . . . , n), consider the usual logistic regression model, Denote by X the n × p covariate matrix with ith row x T i (i = 1, 2, . . . , n). Consider a Gaussian prior β ∼ N p (µ β , Σ β ), denoted by π(β).
Since p α l (y, β | X, ξ) serves as a lower bound to p α (y, β | X) for any ξ ∈ R n , similar to Jaakkola and Jordan [2000] we use an empirical Bayes approach to estimate the variational parameters ξ by maximizing p α l (y | X, ξ) = p α l (y, β | X, ξ)dβ with respect to ξ. The true posterior distribution of β in (2.4) is not available in closed form. However, assuming (2.5) to be a working (pseudo)likelihood of y, β given X, ξ, it is straightforward to see that the corresponding conditional posterior Treating β as latent variables and augmenting with y to get the complete data, one obtains the E-step, where tr(A) denotes the trace of a matrix A. Upon differentiating the above expression with respect to ξ t+1 and using the fact that C T (x) = −x 2 A T (x), we get the M-step, (2.8) The square operation in the above display is to be interpreted elementwise. We assume convergence when the increment in p α l (y | X, ξ t ) is negligible which implies convergence of ξ t by virtue of EM algorithm. The EM sequence in (2.8) is recognized to be a fixed point iteration corresponding to the fixed point equation given by, (2.9) Assuming (2.8) converges to a fixed point ξ * , µ α (ξ * ) gives the variational estimate of β.

Statistical optimality of the variational estimate
In this section we develop a rigorous framework to obtain frequentist risk bounds of the variational approximation obtained in (2.6) at any fixed point ξ * of (2.8). Throughout the section, we assume that the data is generated from a logistic regression model It is not immediately clear whether the empirical likelihood based inference of ξ as discussed in Section 2 falls into the framework of variational inference in the sense of (2.1). In the following, we propose an objective function whose minimizer satisfies the fixed point iteration (2.8). Let our working model be It is important to note here that p α l (y | β, X, ξ) is not a probability density, even when α = 1. Let F be the set of densities on R p . Define a mapping from L : F × R n to R as where p α l (y, β | X, ξ) is defined in (2.5). Observe that L(q, ξ) is the negative of the evidence lower bound obtained in a variational inference with (3.2) as the working likelihood, N p (µ β , Σ β ) the prior on β, and variational family F × {δ ξ : ξ ∈ R n } where δ ξ is the Dirac delta measure on ξ ∈ R n . In Lemma 3.1, we show that the tangent transform algorithm maximizes −L(q, ξ).
We first set the gradient of L(q, ξ) with respect to ξ to zero holding q fixed. As the second term in (3.5) is independent of ξ, this is equivalent to setting the gradient of E q log p α l (y, β | ξ, X) with respect to ξ to be zero, By an application of Fubini, (3.6) is equivalent to For fixed ξ, to maximize (3.3), we simply apply Lemma A.1 in the Appendix. This leads to the optimal choice of q(β) being the conditional distribution p α l (β | y, ξ, X) which is N p (µ α (ξ), Σ α (ξ)/α).
This when combined with (3.7) yields To show that the solution of (3.8) satisfies (2.9), recall that the first-order stationarity condition Although (3.3) is reminiscent of the α-variational objective function of Yang et al. [2020], we note a couple of key differences : (a) p α l (y | β, ξ, X) is not a valid probability density, but it is a lower bound to p α (y | β, X), (b) The latent variables ξ lack a probabilistic interpretation as in Yang et al.
[2020], where one recovers the original likelihood after marginalization over the latent variables.
Here, the latent variables instead correspond to tuning parameters appearing from convex duality.
The usage of fractional likelihood for α ∈ (0, 1) results in only minor changes from a methodological and implementation perspective. However, from a theoretical perspective, like Yang et al.

Variational Risk Bounds
In the following, we develop risk bounds for the variational estimator separately for the case α ∈ (0, 1) and α = 1. In the former case, to quantify the discrepancy between the variational estimate and the true parameter, we use an α-Rényi divergence Refer to Bhattacharya et al. [2019] for more on posterior risk bounds under the α-Rényi divergence.
The factor (1/n) is used to measure average discrepancy per observation. We can further simplify The next theorem derives an upper bound to the risk obtained by integrating the α-Rényi divergence with respect to the optimal variational solution. Denote by The proof of Theorem 3.1 can be found in §A.1 in the Appendix.
Remark 3.1. Setting ε = (p log n/n) 1/2 , the risk bound for discrepancy D α is p/n upto logarithmic terms which is minimax optimal. The explicit bound is non-asymptotic and depends on prior parameters, the covariate matrix X and the true data generating density.
Next, we separately deal with the case α = 1. In doing so, we work with a limiting metric of α-Rényi divergence as α tends to 1. Let a(t) = log(1 + e t ) and a (1) and a (2) denote the first and second derivatives. a satisfies a(t + h) ≥ a(t) + h a (1) (t) + r(|h|) a (2) (t)/2 for all t, h, where The proof of Theorem 3.2 can be found in §A.2 in the Appendix.
Remark 3.2. To obtain a risk bound, we keep γ to be a fixed number in (0, 1) and set ε = (p log n/n) 1/2 . Then KL divergence risk is p/n upto logarithmic terms which is again minimax optimal. As opposed to Theorem 3.1, Theorem 3.2 requires the eigen values of X T W X/n to be bounded from below and the eigen values of X T X/n to be bounded from above.
Now we conduct a numerical study to empirically support the conclusions of the theorems above.
For fixed (n, p) we construct a n × p design matrix X where x T i (i = 1, . . . , n) are independently drawn from N p 0, (0.5 I p + 0.5 1 p 1 T p ) . We then normalize each row of X by √ p. We fix β * to be We place a zero-mean Gaussian prior β ∼ N p (0 p , Σ β ) and set Σ β = 5 2 I p . Given a dataset (y, X) and fixed α ∈ {0.50, 0.65, 0.80, 0.95, 1.00} we calculate the fixed point solution ξ * using (2.8). We use p(y | µ(ξ * ), X) as the final estimated density and calculate the discrepency along with D(β * , µ(ξ * )) that corresponds to α = 1. In panel (b) we plot the 2 norm between µ(ξ * ) and β * . We repeat this process for 500 independent samples with (n = 100, p = 4) and (n = 200, p = 4). Clearly, increasing the sample size leads to improved estimation as seen from either panel of Figure 1. Also, D α (µ(ξ * ), β * ) slightly increases as α increases to 1. Since α-Rényi divergence counterbalances the effect of the misspecified likelihood and reinforces concentration around the truth, this behavior is expected. 4 Stability and convergence of tangent transform algorithm

Preliminaries
We provide a brief review of stability of dynamical systems here; a more detailed review and relevant references can be found in §B of the Appendix. Consider the following discrete-time autonomous system, where f : R n → R n (or, f : D → R n , D ⊆ R n ) is a twice continuously differentiable function. Any ψ * ∈ R n satisfying ψ * = f (ψ * ) is called a fixed point for this system. A fixed point ψ * of (4.1) is called locally asymptotically stable if given any > 0, there exists δ = δ( ) such that whenever The following well-known result is instrumental to show that a fixed point is locally asymptotically stable. Denote by ρ(J) the spectral radius of a square matrix J, the largest eigenvalue of J in absolute value.
Lemma 4.1. Let ψ * be a fixed point solution to the discrete-time autonomous system given by Then, ψ * is locally asymptotically stable if ρ(J) is less than 1.

Asymptotic stability of tangent transform EM
In this subsection, we study the EM sequence of ξ from equation (2.8) viewed as a discrete time dynamical system in ξ 2 . As noted above, the convergence and stability aspects of the system depends crucially on the properties of the Jacobian of the map. Since the function A(ξ) = −tanh(ξ/2)/4ξ is symmetric around 0, and Σ α (ξ) and µ α (ξ) are dependent on ξ through A(·), only the magnitude of ξ is relevant and hence we will discuss the nature of EM iterates on R + . The properties of the function A(·) play a crucial role in such an analysis. In Proposition C.1 in the Appendix, we in particular show that the function A : R + → R − is monotonically increasing and twice continuously differentiable with A(0) = −1/8 and A(ξ) + ξA T (ξ) < 0 for all ξ ∈ R + .
In Theorem 4.1 below, we show that the EM sequence in equation (2.8) is locally asymptotically stable.
Theorem 4.1. Suppose the design matrix X does not have any row equal to the zero vector. For any α ∈ (0, 1] and positive definite Σ β , any fixed point solution ξ * of the EM sequence in (2.8) is locally asymptotically stable.
In light of Lemma 4.1, one needs to check the spectral radius of the Jacobian of the system at the fixed point to prove Theorem 4.1. We present an outline of the proof here; refer to §C.1 in the Appendix for a complete proof. Given positive semi-definite matrices A, B of the same dimension, one can conclude ξ * i > 0 from (2.9). Next, we show that the Jacobian matrix evaluated at the fixed point ξ * can be analytically expressed as where • denotes the Hadamard (or, elementwise) product, µ α (ξ), Σ α (ξ) are defined in (2.6), and D = diag {A T (ξ * )/ξ * }, with the / operation interpreted elementwise. By similarity, J α and have the same set of eigenvalues. Clearly,J α is real symmetric and positive semi-definite by the Schur product theorem. Therefore,J α , and hence J α , have non-negative eigenvalues. Hence, the spectral radius ρ(J α ) is simply the largest eigenvalue of J α , which we proceed to bound next.
Using the fact that, . Then the matrix on the right hand side of the (4.3) can be written as, A result from Horn and Johnson [1994], . Lemma C.3 shows that the largest eigenvalue of M −1/2 N M −1/2 is strictly less than 1 where N ≺ M and M, N are positive definite and positive semi-definite matrices respectively. This delivers the proof that ρ(J α ) < 1.
In the special case when p = 1, we can make substantial simplifications and show that (see §C.2 in the Appendix for details), where the first inequality follows from the fact that the second term in (4.4) is positive as A T (x)/x > 0 for all x ∈ R. The second inequality follows from the fact that A(x) + xA T (x) < 0 for all x ∈ R.
It is important to note that Theorem 4.1 places minimal restriction on the design matrix X and its conclusion remains true for any p and n. We conduct a replicated numerical study to empirically demonstrate some of these features. We use the same simulation design corresponding to Figure 1 except now for a fixed (n, p), we provide a sufficiently flat prior Σ β = 10 2 I p while fixing the first p/2 (resp. [p/2] + 1) entries of β * to be −4, and the remaining p/2 (resp. [p/2]) to be 4 when p is even (resp. odd). To remain faithful to the assumptions of Theorem 4.1, we do not normalize X with √ p. We compute the spectral radius ρ ≡ ρ(J α ) of the Jacobian matrix J α for α ∈ {0.50, 1.00} at the fixed point ξ * for different values of (n, p) over 500 independent replicates, with summary boxplots shown in Figure 3. In panel (a), we fix n = 150 and vary p ∈ {2, 5, 10, 20}. In panel (b), we fix p at 15 and vary n ∈ {5, 10, 50, 100}. It is evident that ρ remains less than 1 for all combinations of (n, p). Observe also that the first two cases in panel (b) correspond to p > n, and as predicted by the theory, the spectral radius continues to be smaller than 1. It can be seen from either panel that on an average ρ at α = 0.5 is higher than the corresponding value at α = 1. Both the plots are produced with 500 replications using the same data (y, X). It can be clearly seen that, for both the α ∈ {0.50, 1.00} the spectral radius is strictly less than 1 irrespective of n and p. . . , 500}, we observe the number of iterations (logscale) required by the algorithm for convergence and calculate the ρ at the fixed point solution for α = 1. We plot (ρ i ,log(iter i )) for all i ∈ {1, 2, . . . , 500} and fitted with LOWESS line to explore the relationship between these two dependent variables. In (a), we generated the data (y, X) with (n = 150, p = 2). In (b), data (y, X) is generated with (n = 150, p = 20). It can be seen that number of iteration grows almost exponentially with the increasing ρ. Also, for fixed n, bigger p leads to higher ρ and as a consequence more iteration are required for convergence.
It is worth noting that local asymptotic stability does not provide any information other than existence of a δ -neighborhood around ξ * such that, if the system is initialized in that region the iterates converge to ξ * as t → ∞. Also, the definition does not say anything about the rate of convergence. In the following, we provide a heuristic argument to connect the notion of rate of convergence with the spectral radius.
For simplicity, consider the one-dimensional system x t+1 = g(x t ) for some function g : R → R which is twice continuously differentiable. If x * is a fixed point of this system, using Taylor's theorem we have for some Recall that the linear rate of convergence [Romero et al., 2019] is given by, γ = lim t→∞ x t+1 − x * / x t − x * , provided the limit exists. In the above scenario, the iterates converge when g(x * ) < 1 and the rate of convergence is g(x * ). For a general d-dimensional linear system α (t+1) = Aα (t) with fixed point and · 2 is the Euclidean norm. Hence ρ(A) acts as a rate of convergence for this case. Figure   4 is an illustration of the number of iterations needed for the system given by (2.8) to converge to the fixed point as a function of ρ(J α=1 ). It is evident that the number of iterations increases exponentially as ρ(J α=1 ) tends to 1.

A special case of semi-orthogonal design
In this section, we shall consider a simple hierarchical logistic regression model given by, p(y ij = 1 | β) = 1/{1 + exp (−β j )} (i = 1, 2, . . . , n; j = 1, 2, . . . , p), (4.5) We assume a prior β ∼ N p (0, σ 2 β I p ). In this case, the results of Section 4.2 can be strengthened to obtain a global convergence rate of the EM sequence (2.8). One key advantage here is the ability to decouple the EM sequence into independent coordinate-wise updates. This is illustrated in Lemma 4.2.
Lemma 4.2. The EM updates for the model (4.5) can be simplified to, where,ȳ j = n i=1 y ij /n, for all j = 1, 2, . . . , p and ζ t is the update at the t th iteration.
It is important to distinguish between the EM update ξ in (2.8) and ζ in (4.6). In the general setting (2.3), the variational parameter ξ are introduced for each individual i ∈ {1, 2, . . . , n}, whereas ζ is introduced here for different groups j ∈ {1, 2, . . . , p}. Though we used similar techniques to get the updates, they have different interpretation. Figure 3 and Figure 5 are not comparable in that sense.
The parallelization of the updates of ζ makes the posterior of β independent. Also, since the updates are independent and identical for all j = 1, 2, . . . , p given the initial point, it suffices to study the stability of a single coordinate. The following theorem assures the global asymptotic stability of the EM sequence in (4.6).
It can be seen from (4.6) that the updates of ζ j depend on the (y, X) throughȳ j , which is a sufficient statistic for β j . Therefore if for some j = j T we haveȳ j =ȳ j T , the sequences {ζ t j } and {ζ t j T } converge to the same limit. Figure 5a shows the global convergence of the EM sequence for differentȳ j when σ β = 1 and n = 50 with arbitrary initializations. Numerically we assumed convergence when |ζ t+1 − ζ t | < 10 −10 . Interestingly, it is observed that convergence is slower when the data becomes more imbalanced, i.e. |ȳ j − 0.5| → ±0.5. A similar behavior for the mixing time of the Pólya-Gamma data augmentation Gibbs sampling in Bayesian logistic regression is observed in Johndrow et al. [2019], which is all the more interesting given the connection between Pólya-Gamma augmentation and tangent transforms established by Durante and Rigon [2019].  Figure 5c shows that for fixed z, h T 1,n (z) is an increasing function of n and less than 1. Lemma C.5 proves this fact and in addition shows that for fixed z, lim n→∞ h T 1,n (z) < 1. It is important to note is that h T σ β ,n (z) is dependent on σ β and for large σ β and fixed z, h T σ β ,n (z) may not be an increasing function of n. Finally, Figure 5d shows h T u,σ β ,10 ∞ increases as σ β increases. It can be easily verified that for fixed u and n, h T u,σ β ,n (z) is an increasing function of σ β and also for fixed n and σ β , an increasing function of |u|. Numerically it can be seen that h T 0.5,σ β ,10 ∞ ≥ 1 when σ β ≥ 12.894. Overall, as the data gets more imbalanced, a flatter prior β increasingly hurts the convergence.
Lemma 5.1. The EM updates to the above multinomial logit regression are given by, Proof. We begin with the log-fractional likelihood,  [Lawrence et al., 2004], Jensen's inequality [Blei and Lafferty, 2007], quadratic approximation [Bouchard, 2008, Jebara andChoromanska, 2012]. In the following, we use a very recent method by Titsias [2016] which has been numerically shown to outperform the others which is based on the following inequality, In an ideal scenario, if y i = l both the terms in the inequality above should be ≈ log(1 + exp(x T i β l )).
This inequality is not too loose in that sense. Using (5.3) in (5.2), we get a lower bound to p α (y | X, β) given by, Now we use the quadratic bound proposed by Jaakkola and Jordan [2000] on the right hand side of the above inequality. This leads to a lower bound to log p α (y | X, β) similar to (2.5) where Call ξ * = χ * j and ξ t = χ t j . Then, (5.1) and (5.5) reduces to (2.8) and (2.9) respectively. Now, we directly apply Theorem 4.1 to conclude that the updates in (5.1) are locally asymptotically stable.

A Proof of Statistical optimality results in Section 3
In the following, we first provide the proofs of Theorems 3.1 and 3.2 in §A.1 and A.2 respectively and then provide the proofs of some of the auxiliary results used in subsequent §A.3.

A.1 Proof of Theorem 3.1
The proof consists of two major steps.
Risk majorization. In this first step, we obtain an upper bound to the integrated risk in terms of easily controllable quantities. We denote E β * as taking expectation under (3.1). From the definition of the α-Renyi divergence and the fact that p l lower bounds p(y | β, X) Thus, for any ε ∈ (0, 1), we have Integrating both side of this inequality with respect to the prior π β and interchanging the integrals using Fubini's theorem, we obtain Now, recall the variational inequality for a probability measure µ and for h such that e h is integrable, If we choose ρ = q * β ≡ φ p β; µ α (ξ * ), Σ α (ξ * ) as the variational approximation and set ξ = ξ * E β * exp α log p l (y | β, ξ * , X) p(y | β * , X) By applying Markov's inequality, we further obtain that with P β * probability at least (1 − ε), Optimizing the majorized risk. Our second step consists of optimizing the term obtained in (A.2) by choosing suitable candidates for q and ξ. We refer to them asq andξ. The idea is to choosẽ q andξ so thatq places almost all its mass into a small neighborhood around truth Xβ * , so that the first term in the right hand side of (A.2) becomes small; on the other hand, the neighborhood is large enough so that the second regularization term n −1 D(q || π β ) is not too large. We chooseq first andξ later. Letqq be the restriction of the prior density π β into the KL neighborhood B n (β * , ε) around β * with radius ε defined as where for two non-negative functions f, g, D(f g) = f | log(f /g)| and V(f, g) : Note thatD(f | g) is an extension of the usual KL distance for probability measures to positive functions which may not integrate to one. With this substitution, the second term in (A.2) becomes the negative log prior mass [n (1 − α)] −1 log π β [B n (β * , ε)] −1 and it remains to provide a high-probability bound for the first term and an upper bound for the log-prior concentration term log π β [B n (β * , ε)] .

A.3 Auxiliary results for proofs in Section 3
Lemma A.1. Let u and v denote two continuous random vectors with joint density function p(u, v).
The maximum value of q(u) log p(u, v) q(u) du over all density functions q is attained by q * (u) = p(u | v).

B Review of Dynamical Systems & Notion of Stability
Dynamical systems theory is a classical technique that deals with stability and convergence of complex iterative methods. We call a dynamical system to be discrete-time if the system is observed on discrete time points {t 0 , t 1 , t 2 , . . .}. Usually, we consider the time-points to be evenly placed, i.e. t j+1 = t j + h for some h > 0. Moreover, a system is considered autonomous if the function is independent of time and non-autonomous otherwise. In this section, we will discuss the notion of stability for discrete time autonomous systems. Let us consider the following discrete-time autonomous system given by, where f : R n → R n (or, f : D → R n , D ⊆ R n ) is a diffeomorphism, i.e. a smooth function with smooth inverse and ψ t ∈ R n . ψ * ∈ R n is called a fixed point to this system if ψ * = f (ψ * ). We recall the following definition from Bof et al. [2018].
(c) locally unstable if it is not locally stable.
The locality in the definition is used to denote the fact that we are initializing the system in a δ-ball around the fixed point. We say the stability is global if the system converges to the fixed point independent of the initialization, i.e. we can initialize at any point in the function domain.
Lemma B.1. Consider a system x t+1 = g(x t ) where g : D → R (D ⊆ R) with a fixed point x * such that, |g T (x)| ≤ δ for all x ∈ D − {x * }, for some δ < 1. Then, x * is globally asymptotically stable.
Proof. Given a fixed point x * = g(x * ), use the Mean Value Theorem to get, Let β k be the k-th coordinate of a vector β ∈ R n . Consider the system in (B.1) with a fixed point ψ * . Using generalized Taylor's theorem we get, where ∇f k (ψ) is the gradient vector with i th entry given by ∂f k (ψ)/∂ψ i and h : R n → R such that, lim ψ→ψ * h(ψ) = 0. If ψ t is close to ψ * , the convergence of the system depends on ∇f k (ψ * ) by the following approximation, where J is the n × n Jacobian matrix evaluated at ψ * with i th row given by ∇f T i (ψ * ). Thus the behavior of the dynamical system (B.1) around a small neighbourhood of ψ * is exactly same as that of the linearization in (B.2). This is formalized in the Hartman-Grobman theorem.
Definition B.2. A fixed point ψ * , for a map ψ → f (ψ), ψ ∈ R n is called hyperbolic if none of the eigenvalues of J has magnitude 1.
(b) ψ * is locally unstable if at least one eigenvalue of J is greater than one in absolute value.
Lemma B.2 along with Theorem B.1 provides sufficient conditions for the local convergence of a system. Consider the linear system given by α (t+1) = Aα (t) , with a fixed point α * = 0. Let us consider, Aν i = λ i ν i (i = 1, 2, . . . , n) and |λ 1 | ≥ |λ 2 | ≥ . . . ≥ |λ n |. Suppose A has a complete set of eigenvectors, i.e. the set of eigenvectors {ν 1 , ν 2 , . . . , ν n } form a basis of R n . Then it can be easily seen that a solution to the system is α (t) = c 1 λ t 1 ν 1 + c 2 λ t 2 ν 2 + . . . + c n λ t n ν n for some arbitrary constants c 1 , c 2 , . . . , c n . Also, c 1 λ t 1 ν 1 + c 2 λ t 2 ν 2 + . . . + c n λ t n ν n → 0 iff |λ 1 | < 1. This illustrates the Lemma above, in the most simplistic scenario, can be extended to the case where A does not have a complete set of eigenvectors using Jordan Canonical form of A (refer to Wood and O'Neill [2003]). In the following, we first provide the proof of Theorems 4.1 in §C.1, calculation of spectral radius for p = 1 in §C.2 and then provide the proofs of some of the auxiliary results used in subsequent §C.3.

C.1 Proof of Theorem 4.1
For some fixed α ∈ (0, 1], the update equation in (2.8) can be rewritten as, We calculate the partial derivatives in order to get the Jacobian matrix, Then Jacobian Matrix(J α ) at ξ = ξ t is given by, where • denotes the Hadamard Product. Let us denote the maximum eigenvalue of a matrix A by λ 1 (A). Our objective is to show that λ 1 (J α )| ξ=ξ * < 1. We call D = diag (A T (ξ)/ξ). By Lemma C.1 has the same set of eigenvalues as with J α | ξ=ξ * . XΣ α (ξ * )X T and X {Σ α (ξ * )/α + 2µ α (ξ * )µ T α (ξ * )} X T are positive semi-definite matrices as well as symmetric. Since the eigenvalues of a real symmetric positive semi-definite matrix are real and non-negative, the eigenvalues of J α | ξ=ξ * are real and non-negative. We denote ξ * by ξ in the following discussion for notational simplicity.
By rearranging the terms it follows that λ 1 < 1.

C.3 Auxiliary results for the proofs in Section 4
Lemma C.1. For a symmetric matrix M n×n and a invertible diagonal matrix N n×n , the set of eigenvalues of MN and N 1/2 MN 1/2 are the same.