Inference for Generalized Linear Models via Alternating Directions and Bethe Free Energy Minimization
Abstract
Generalized Linear Models (GLMs), where a random vector is observed through a noisy, possibly nonlinear, function of a linear transform , arise in a range of applications in nonlinear filtering and regression. Approximate Message Passing (AMP) methods, based on loopy belief propagation, are a promising class of approaches for approximate inference in these models. AMP methods are computationally simple, general, and admit precise analyses with testable conditions for optimality for large i.i.d. transforms . However, the algorithms can diverge for general . This paper presents a convergent approach to the generalized AMP (GAMP) algorithm based on direct minimization of a largesystem limit approximation of the Bethe Free Energy (LSLBFE). The proposed method uses a doubleloop procedure, where the outer loop successively linearizes the LSLBFE and the inner loop minimizes the linearized LSLBFE using the Alternating Direction Method of Multipliers (ADMM). The proposed method, called ADMMGAMP, is similar in structure to the original GAMP method, but with an additional leastsquares minimization. It is shown that for strictly convex, smooth penalties, ADMMGAMP is guaranteed to converge to a local minimum of the LSLBFE, thus providing a convergent alternative to GAMP that is stable under arbitrary transforms. Simulations are also presented that demonstrate the robustness of the method for nonconvex penalties as well.
I Introduction
Consider the problem of estimating a random vector from observations as shown in Fig. 1. The unknown vector is assumed to have a prior density of the form and the observations are described by a likelihood function of the form for some known transform . In statistics, this model is a special case of a generalized linear model (GLM) [1, 2] and arises in a range of applications including statistical regression, filtering, inverse problems, and nonlinear forms of compressed sensing. The posterior density of given in the GLM model is given by
(1) 
where is a normalization constant. In the sequel, we will often omit the dependence on and simply write
(2) 
so that the dependence on in the function and the normalization constant is implicit. In this work, we consider the inference problem of estimating the posterior marginal distributions, . From these posterior marginals, one can compute the posterior means and variances
(3a)  
(3b) 
We study this inference problem in the case where the functions and are separable, in that they are of the form
(4a)  
(4b) 
for some scalar functions and . The separability assumption (4a) corresponds to the components in being a priori independent. Recalling the implicit dependence of on , the separability assumption (4b) corresponds to the observations being conditionally independent given the transform outputs .
For posterior densities of the form (2), there are several computationally efficient methods to find the maximum a posteriori (MAP) estimate, which is given by
(5) 
Under the separability assumptions (4), the MAP minimization (5) admits a factorizable dual decomposition that can be exploited by a variety of approaches, including variants of the iterative shrinkage and thresholding algorithm (ISTA) [3, 4, 5, 6, 7, 8] and the alternating direction method of multipliers (ADMM) [9, 10, 11, 12].
In contrast, the inference problem of estimating the posterior marginals and the corresponding minimum mean squared error (MMSE) estimates (3a) is often more difficult—even in the case when and are convex. As a simple example, consider the case where and each constrains to belong to some interval, so that constrains to belong to some polytope. The MAP estimate (5) is then given by any point in the polytope. Such a point can can be computed via a linear program. However, the MMSE estimate (3a) is the centroid of the polytope which is, in general, #Phard to compute [13].
GLM inference methods often use a penalized quasilikelihood method [14] or some form of Gibbs sampling [15, 16]. In recent years, Bayesian forms of approximate message passing (AMP) have been considered as a potential alternate class of methods for approximate inference in GLMs [17, 18, 19, 20, 21, 22]. AMP methods are based on Gaussian and quadratic approximations to loopy belief propagation (loopy BP) in graphical models and are both computationally simple and applicable to arbitrary separable penalty functions and . In addition, for certain large i.i.d. transforms , they have the benefit that the behavior of the algorithm can be exactly predicted by a state evolution analysis, which then provides testable conditions for Bayes optimality [23, 22, 24].
Unfortunately, for general , AMP methods may diverge [25, 26]—a situation that is not surprising given that AMP is based on loopy BP, which also may diverge. Several recent modifications have been proposed to improve the stability of AMP, including damping [25], sequential updating [27], and adaptive damping [28]. However, while these methods appear to perform well empirically, little has been proven rigorously about their convergence.
The main goal of this paper is to provide a provably convergent approach to AMP. We focus on the generalized AMP (GAMP) method of [22], which allows arbitrary separable functions for both and . Our approach to stabilizing GAMP is based on reconsidering the inference problem as a type of freeenergy minimization. Specifically, it is known that GAMP can be considered as an iterative procedure for minimizing a largesystemlimit approximation of the socalled Bethe Free Energy (BFE) [29, 30], which we abbreviate as “LSLBFE” in the sequel. The BFE also plays a central role in loopy BP [31], and we review both the BFE and LSLBFE in Section III.
In contrast to GAMP, which implicitly minimizes the LSLBFE through an approximation of the sumproduct algorithm, our proposed approach explicitly minimizes the LSLBFE. We propose a doubleloop algorithm, similar to the wellknown Convex Concave Procedure (CCCP) [32]. Specifically, the outer loop of our method successively approximates the LSLBFE by partially linearizing the LSLBFE around the current belief estimate, while the inner loop minimizes the linearized LSLBFE using ADMM [9]. Similar applications of ADMM have also been proposed for related freeenergy minimizations [33, 34]. Interestingly, our proposed doubleloop algorithm, which we dub ADMMGAMP, is similar in structure to the original GAMP method of [22], but with an additional least squares optimization. We discuss these differences in detail in Section VIII.
Our main theoretical result shows that, for strictly convex penalties, the proposed ADMMGAMP algorithm is guaranteed to converge to at least a local minimum of the LSLBFE. In this way, we obtain a variant of the GAMP method with a provable convergence guarantee for arbitrary transforms . In addition, using hardening arguments similar to [35, 36], we show that our ADMMGAMP can also be applied to the MAP estimation problem, in which case we can obtain global convergence for strictly convex, smooth penalties. Also, while our theory requires convex penalties, we present simulations that show robust behavior even in nonconvex cases.
Ii The GLM and Examples
Before describing our optimization approach, it is useful to briefly provide some examples of the model (1) to illustrate the generality of the framework. As a first simple example, consider a simple linear model
(6) 
where is a known matrix, is an unknown vector and is a noise vector. In statistics, would be the data matrix with predictors, would be the vector of regression coefficients, the vector of target or response variables and would represent the model errors. To place this model in the framework of this paper, we must impose a prior on and model the noise as a random vector independent of and . Under these assumptions, the posterior density of given will be of the form (1) if we define
(7) 
The separability assumption (4) will be valid if the components of are are independent so the prior and noise density factorizes as
If the output noise is Gaussian with independent components , the output factor in (7) has a quadratic cost,
Similarly, if has a Gaussian prior with , the input factor will be given by
Note that the estimation in this quadratic case would be given by standard least squares estimation.
However, much more general models are possible. For example, for Bayesian forms of compressed sensing problems [37], one can impose a sparse prior such as a BernoulliGaussian or a heavytailed density.
Also, for the output, any likelihood that factorizes as can be incorporated. This model would occur, for example, under any output nonlinearities as considered in [38],
where is a known, nonlinear function and is noise. The model can also include logistic regression [39] where is a binary class variable and
for some sigmoidal function . Onebit and quantized compressed sensing [40] as well as Poisson output models [41] can also be easily modeled.
Iii Bethe Free Energy Minimization
We next provide a brief review of the Bethe Free Energy (BFE) minimization approach to estimation of marginal densities in GLMs. A more complete treatment of this topic, along with related ideas in variational inference, can be found in [31, 42].
For a generic density , exact computation of the marginal densities is difficult, because it involves a potentially highdimensional integration. BFE minimization provides an approximate approach to marginal density computation for the case when the joint density admits a factorizable structure of the form
(8) 
where, for each , is a subvector of created from indices in the subset and is a potential function on that subvector. In this case, BFE minimization aims to compute the vectors of densities
where represents an estimate of the marginal density and where represents an estimate of the joint density on the subvector . These density estimates, often called “beliefs,” are computed using an optimization of the form
(9) 
where is the BFE given by
(10) 
where is the KL divergence,
(11) 
where is the entropy or differential entropy; and where (for each ) is the number of factors such that . The BFE minimization (9) is performed over the set of all whose components satisfy a particular “matching” condition: for each , the marginal density of within the belief must agree with the belief . That is, the set contains all such that
(12) 
where the integration is over the components in the subvector holding constant. Note that imposes a set of linear constraints on the belief vectors and .
The BFE minimization exactly recovers the true marginals in certain cases (e.g., when the factor graph has no cycles) and provides good estimates in many other scenarios as well; see [42] for a complete discussion. In addition, due to its separable structure, the BFE can be typically minimized “locally,” by solving a set of minimizations over the densities and . When the cardinalities of the subsets are small, these local minimizations may involve much less computation than directly calculating the marginals of the full joint density . In fact, the classic result of [31] is that loopy belief propagation can be interpreted precisely as one type of iterative local minimization of the BFE.
For the GLM in Section I, the separability assumption (4) allows us to write the density (2) in the factorized form (8) using the potentials
(13a)  
(13b) 
where is the th row of . Note that, if is a nonsparse matrix, then depends on all components in the vector . In this case, the application of traditional loopy BP—as described for example in [43]—does not generally yield a significant computational improvement.
The GAMP algorithm from [22] can be seen as an approximate BFE minimization method for GLMs with possibly dense transforms . Specifically, it was shown in [29] that the stationary points of GAMP coincide with the local minima of the constrained optimization
(14b)  
where and are product densities, i.e.,
(15) 
and the objective function is given by
(16)  
(17)  
(18)  
(19)  
(20) 
Above, and in the sequel, we use to denote the expectation of under , and we use to denote the vector whose th entry is the variance of under . Note that is not a full covariance matrix. Also, is the scale factor that makes a valid density over . Although it is not essential for this paper, we note that is an upper bound on the differential entropy of that is tight when has independent Gaussian entries with variances . It was then shown in [30] that the objective function in (16) can be interpreted as an approximation of the BFE for the GLM from Section I in a certain largesystem limit, where and has i.i.d. entries. We thus call the approximate BFE in (16) the largesystem limit Bethe Free Energy or LSLBFE.
Similar to the case of loopy BP, it has been shown in [29, 30] that the stationary points of (14) are precisely the fixed points of sumproduct GAMP. Thus, GAMP can be interpreted as an iterative procedure to find local minima of the LSLBFE, much in the same way that loopy BP can be interpreted as an iterative way to find local minima of the BFE. The trouble with GAMP, however, is that it does not always converge (see, e.g., the negative results in [25, 26, 28]). The situation is similar to the case of loopy BP. Although several modifications of GAMP have been proposed with the goal of improving convergence, such as damping [25], sequential updating [27], and adaptive damping [28], a globally convergent GAMP modification remains elusive.
Iv Minimization via Iterative Linearization
Our approach to finding a convergent algorithm for minimizing the constrained LSLBFE employs a generalization of the convexconcave procedure (CCCP) of [32] that we will refer to as Minimization via Iterative Linearization.
Iva The ConvexConcave Procedure
We first briefly review the CCCP. Observe that, in the BFE (10), the terms are convex in and the terms are concave in . Thus, the BFE (10) can be written as a sum of terms
where is convex and is concave. The CCCP finds a sequence of estimates of a BFE minimizer by iteratively linearizing the concave part of this function, i.e.,
(21a)  
(21b) 
where denotes the gradient of at . The resulting procedure is often called a “doubleloop” algorithm, since each iteration involves a minimization (21a) that is itself usually performed by an iterative procedure. Because is convex and the constraint is linear, the minimization problem (21a) is convex. Thus, the CCCP converts the nonconvex BFE minimization to a sequence of convex minimizations. In fact, it can be shown that the CCCP will monotonically decrease the BFE for arbitrary convex and concave [32].
IvB Minimization via Iterative Linearization
For the LSLBFE, it is not convenient to decompose the objective function into a convex term plus a concave term. To handle problems like LSLBFE minimization, we consider optimization problems of the form
(22) 
where now is a vector in a Hilbert space , is a closed affine subspace of , is a convex functional, is a mapping from to for some , and is an arbitrary functional. Below, we use to denote the input to . Note that the functionals and may be neither concave nor convex.
To solve (22), we propose the iterative procedure shown in Algorithm 1, which is reminiscent of the CCCP. At each iteration , an estimate of is computed by minimizing the functional
(23) 
where is a “damped” version of the gradient . In particular, when the damping parameter is set to unity, the linearization vector is exactly equal to the gradient at , i.e., , similar to CCCP. However, in Algorithm 1, we have the option of setting , which has the effect of slowing the update on . We will see that, by setting , we can guarantee convergence when and/or is nonconcave.
IvC Convergence of Minimization via Iterative Linearization
Observe that when is convex, is concave, (as when are discrete variables), is the identity map (i.e., ), and there is no damping (i.e., ), Algorithm 1 reduces to the CCCP. However, we are interested in possibly nonconcave , in which case we cannot directly apply the results of [32]. We instead consider the following alternate conditions.
Assumption 1
Consider the optimization problem (22), and suppose that the functions , , and have components that are twice differentiable with uniformly bounded second derivatives. Also, assume that there exists a convex set such that, for all :

The minimization of the linearized function,
(24) exists and is unique.

At each minimum, the linearized objective is uniformly strictly convex in the linear space in that there exists constants with such that
(25) where is the Hessian of with respect to at , i.e.,
(26) and where the constants and do not depend on .

The gradient obeys when .
Theorem 1
Proof:
See Appendix A.
The most simple case where Assumption 1 holds is the setting where is strictly convex and smooth, is linear and is smooth (but neither necessarily convex nor concave). Under these assumptions, would be strictly convex for all , thereby satisfying Assumptions (a) and (b). The assumption would also be valid for strictly convex and convex , provided we restrict to positive . In this case, to satisfy assumption (c), we would require that , i.e. is increasing in each of its component. Interestingly, in the setting we will use below, will be convex, but will be concave. Nevertheless, we will show that the assumption will be satisfied.
IvD Application to LSLBFE Minimization
To apply Algorithm 1 to the LSLBFE minimization (14), we first take to be the vector of separable density pairs satisfying the moment matching constraint
(28) 
Then, if we define the functions
(29a)  
(29b)  
(29c) 
we see that from (16) can be cast into the form in (22). Observe that, while is convex, the function is, in general, neither convex nor concave. Thus, while the CCCP does not apply, we can apply the iterative linearization method from Algorithm 1.
We will partition the linearization vector conformally with function in (29b) as
(30) 
where we use “” to denote componentwise division of two vectors and “” to denote vertical concatenation. The notation in (30) will help to clarify the connections to the original GAMP algorithm. Using the above notation, the linearized objective (23) can be written as
(31)  
Finally, we compute the gradient of the function from (29c). Similar to , we will partition the gradient into two terms,
(32) 
From (17), the derivative of with respect to is
(33) 
Similarly, using the chain rule and (33), we find
(34)  
We can then write (33) and (34) in vector form as
(35) 
Substituting the above computations into the iterative linearization algorithm, Algorithm 1, we obtain Algorithm 2. We refer to this as the outer loop, since each iteration involves a minimization of the linearized LSLBFE in line 5. We discuss this latter minimization next and show that it can itself be performed iteratively using a set of iterations that we will refer to as the inner loop.
IvE Alternative Methods
While the method proposed in this paper is based on CCCP of [32], there are other methods for direct minimization of the BFE that may apply to the LSLBFE as well. For example, for problems with binary variables and pairwise penalty functions, [44, 45] propose a clever reparametrization to convert the constrained BFE minimization to an unconstrained optimization on which gradient descent can be used. Unfortunately, it is not obvious if the LSLBFE here can admit such a reparametrization since the penalty functions are not pairwise and the variables are not binary.
V InnerLoop Minimization and ADMMGAMP
Va ADMM Principle
The outer loop algorithm, Algorithm 2, requires that in each iteration we solve a constrained optimization of the form
(36) 
We will show that this optimization can be performed by the Alternating Direction Method of Multipliers (ADMM) [9]. ADMM is a general approach to constrained optimizations of the form
(37) 
where is an objective function and is some constraint matrix. Corresponding to this optimization, let us define the augmented Lagrangian
(38) 
where is a dual vector, is a vector of positive weights and . The ADMM procedure then produces a sequence of estimates for the optimization (37) through the iterations
(39a)  
(39b) 
where creates a diagonal matrix from the vector . The algorithm thus alternately updates the primal variables and dual variables . The vector can be interpreted as a stepsize on the primal problem and an inverse stepsize on the dual problem.
The key benefit of the ADMM method is that, for any positive stepsize vector , the procedure is guaranteed to converge to a global optimum for convex functions under mild conditions on .
VB Application of ADMM to LSLBFE Optimization
The ADMM procedure can be applied to the linearized LSLBFE optimization (36) as follows. First, we replace the constraint with two constraints: and . Variable splittings of this form are commonly used in the context of monotropic programming [46]. With this splitting, the augmented Lagrangian for the LSLBFE (14) becomes
(40)  
where and represent the dual variables. Note that the vectors and that appear in the linearized LSLBFE have been used for the augmentation terms (i.e., the last two terms) in (40). This choice will be critical. From (39), the resulting ADMM recursion becomes
(41a)  
(41b)  
(41c)  
(41d) 
To compute the minimization in (41a), we first note that the second and fourth terms in (40) can be rewritten as
(42)  
where in (a) we used ; in (b) we used “” to denote componentwise multiplication between vectors; and “const” includes terms that are constant with respect to and . A similar development yields