Max–Planck–Institut fu¨r biologische KybernetikMax Planck Institute for Biological CyberneticsTechnical Report No. TR-175Large Scale VariationalInference and ExperimentalDesign for Sparse GeneralizedLinear ModelsMatthias W. Seeger,1 Hannes Nickisch,1September 20081 Department Scho¨lkopf. Contact: seeger@tuebingen.mpg.deDagstuhl Seminar Proceedings 09181 Sampling-based Optimization in the Presence of Uncertainty http://drops.dagstuhl.de/opus/volltexte/2009/2114Large Scale Variational Inference andExperimental Design for Sparse GeneralizedLinear ModelsMatthias W. Seeger, Hannes NickischAbstract. Sparsity is a fundamental concept of modern statistics, and often the only general principle avail-able at the moment to address novel learning applications with many more variables than observations. Whilemuch progress has been made recently in the theoretical understanding and algorithmics of sparse point estima-tion, higher-order problems such as covariance estimation or optimal data acquisition are seldomly addressed forsparsity-favouring models, and there are virtually no algorithms for large scale applications of these. We providenovel approximate Bayesian inference algorithms for sparse generalized linear models, that can be used with hun-dred thousands of variables, and run orders of magnitude faster than previous algorithms in domains where eitherapply. By analyzing our methods and establishing some novel convexity results, we settle a long-standing openquestion about variational Bayesian inference for continuous variable models: the Gaussian lower bound relax-ation, which has been used previously for a range of models, is proved to be a convex optimization problem, if andonly if the posterior mode is found by convex programming. Our algorithms reduce to the same computationalprimitives than commonly used sparse estimation methods do, but require Gaussian marginal variance estimationas well. We show how the Lanczos algorithm from numerical mathematics can be employed to compute the latter.We are interested in Bayesian experimental design here (which is mainly driven by efficient approximate infer-ence), a powerful framework for optimizing measurement architectures of complex signals, such as natural images.Designs optimized by our Bayesian framework strongly outperform choices advocated by compressed sensing the-ory, and with our novel algorithms, we can scale it up to full-size images. Immediate applications of our methodlie in digital photography and medical imaging.We have applied our framework to problems of magnetic resonance imaging design and reconstruction, and partof this work appeared at a conference (Seeger et al., 2008). The present paper describes our methods in muchgreater generality, and most of the theory is novel. Experiments and evaluations will be given in a later paper.1 IntroductionGeneralized linear models are cornerstones of applied statistics, and are also very frequently used in machinelearning. In many applications from low level computer vision, bio-informatics, neuroscience, information re-trieval, adaptive filtering and control, or medical image reconstruction, a vast number of features could potentiallybe used, although the important ones for any given task are not known a priori. In such contexts, the concept ofsparsity regularization or sparsity priors is of central importance, either to select a relevant subset of features in adata-driven manner, or to improve estimation or inference by conditioning them on the assumption that a solutionfree of noise ought to be sparsely distributed.Sparsity is a fundamental concept of modern, high-dimensional statistics. For signal classes such as naturalimages, it can be motivated by direct empirical analysis: images are sparsely distributed in certain transformdomains. In general, it complements linearity as common model assumption in a way which is strictly differentfrom the traditional least squares (LS) idea. If a signal is represented as linear combination over some basis set (ordictionary), we may enforce “simplicity” of this representation in different ways. LS, or Gaussian, regularizationpenalizes the size of all linear coefficients uniformly, so that simple signals have smallish coefficients throughout.In contrast, sparsity regularization enforces most coefficients to be tiny or even zero, at the expense of allowingfor a few large ones, so that simple signals are dominantly explained by few basis functions. The rather generalsuccess of sparsity statistics at denoising signals comes from the fact that, due to the central limit theorem, noiseis typically Gaussian (coming from many small, independent error sources), while real-world signals are typically1highly structured, therefore rather sparsely distributed. More details about the statistical role of sparsity can befound in (Seeger, 2008).When sparsity is enforced in LS estimation or approximate Bayesian inference, the method is to concentrate ona small subset of explanatory variables, about which nothing is known explicitly beforehand, and one might guessthat a combinatorial problem is lurking behind the scenes. Indeed, a harsh formulation of sparse estimation, knownas L0 regularization (see Section 5.2), is NP complete. Fortunately, a surge of recent activity has established thatin many practically relevant cases, feature selection or sparse estimation can be performed by convex programs(Donoho, 2006; Cande`s et al., 2006; Donoho and Elad, 2003), which can be solved very efficiently. Modernalgorithms achieve scalability to very many variables by reducing their dominating efforts to reweighted leastsquares problems, for which efficient code is in common use.In this paper, we are interested in Bayesian inference and applications thereof, problems which are distinctlydifferent from sparse estimation. For example, if u is an unknown image, and y are linear measurements thereof,sparse estimation is useful in order to reconstruct u from y through a single point estimate. Traditional LSestimation is equivalent to maximum-likelihood estimation under the assumption of Gaussian noise. In sparseestimation, we regularize the LS estimate by minimizing the sum of the negative log likelihood − logP (y|u)and a sparsity-enforcing regularization term acting on u. From a statistical viewpoint, a single point estimateis unsatisfying, since it does not come with an assessment of uncertainty. For example, if we had an additionalmeasurement yi on top of y, would our estimate improve, and if so, how exactly? In Bayesian statistics, theposterior distribution is a complete representation of the uncertainty in a reconstruction of u from data y. Thesparsity regularization term is replaced by a sparsity prior P (u), a distribution over images u, and the posterioris obtained by Bayes formula as P (u|y) ∝ P (y|u)P (u). A sparse estimator is recovered as the posterior mode:uˆ = argmaxu P (u|y) = argminu [− logP (y|u)−logP (u)], themaximum a-posteriori (MAP) estimator. Here,− logP (u) plays the role of the sparsity regularizer. But the posterior P (u|y) is much more than just its mode,and can be used for tasks beyond point estimation. In this paper, we will demonstrate a particular application,which is driven by posterior covariances, and for which image point estimation would be of no value at all. Whilemany of the recently studied sparse estimation settings correspond to (convex) posterior maximization in sparselinear models, real Bayesian applications require us to integrate over the posterior, a computation which is lesswell understood, and is thought to be much harder to do in general.Nevertheless, our motivation is to apply Bayesian inference in domains where sparse estimation is typicallyused. This means that our algorithms have to run efficiently for very large numbers of variables, yet still produceuseful uncertainties along with point estimates. We focus on a variational relaxation of inference here, which hasbeen used before on problems of moderate size. The range of models in the scope of our results here is not re-stricted to Gaussian-linear likelihoods, but contains generalized linear models as well, such as logistic regression.We establish that a wide range of instances of these variational relaxations are convex minimization problems withunique solutions, an insight which is novel to the best of our knowledge. For example, variational Bayesian treat-ments of sparsity models with Laplace sites, of binary classification models, or of combinations thereof, are shownto be convex problems. In a nutshell, whenever potentials in the model are amenable to Fenchel duality lowerbounding (Jaakkola, 1997; Palmer et al., 2006), and the log posterior is concave (meaning that MAP estimation isconvex), the variational relaxation is proved to be convex.Beyond this insight, which strengthens a range of prior work, we present novel scalable algorithms to solvethis variational relaxation of Bayesian inference for sparse (generalized) linear models very efficiently, whetherit is convex or not.1 Just as for convex estimation, scalability means that all major computations are reduced tostandard primitives of large scale numerical mathematics and least squares estimation, which have received muchattention already in most computational application fields. This has the very considerable advantage that highlyoptimized code can be made to use, and that our methods can be imported into many large scale applicationswith minor efforts. We reduce inference to a sequence of reweighted least squares problems, as well as Gaussianmarginal variance computations. These can be reduced further (by standard numerical mathematics algorithms)to matrix-vector multiplication (MVM) primitives with model matrices, any structure exploitation of which has adirect impact on our dominating computations, with no further heuristics to be tuned.Our special interest is in Bayesian experimental design, which is a framework for improving measurement ar-chitectures automatically, with the aim of obtaining reconstructions of equivalent quality under lower cost. In theimage reconstruction example, which measurement design allows for the best sparse reconstruction of u from y?To answer this question, we compute design scores (expectations over the current posterior), whose inspection1 A local minimum is found if the relaxation is not convex.2reveals directions of improvement for the current design. The power of our approach in practice has been demon-strated for a number of applications already (Steinke et al., 2007; Seeger, 2008; Seeger and Nickisch, 2008), butthe variational approximations and algorithms used there are not scalable and cannot be used in the large scaledomains of interest here. With the novel methods presented here, these settings can be lifted to full-size images,and problems in medical image reconstruction can be addressed.On a high level, our approach can be understood as a relaxation of Bayesian inference for distinctly non-Gaussiangeneralized linear models to a (small) sequence of Gaussian linear model computations, such as computing meansand marginal variances. The considerable experience with Gaussian random fields in (say) low-level computervision can therefore be used to address inference in non-Gaussian models, which represent natural image charac-teristics such as sparsity much better than Gaussian models do, or in models for discrete observations.Beyond pure Bayesian applications, the problem of finding very good or even optimal designs for subsequentsparse image reconstruction does not have a satisfying solution yet. While much is known about good mea-surements supporting linear LS estimation, nonlinear sparse reconstruction corrects for many shortcomings of thelatter, so that the relevance of many linear design properties is most probably diminished. Recent theoretical resultsabout sparse convex estimation (Donoho, 2006; Cande`s et al., 2006; Donoho and Elad, 2003) are not helpful in thatrespect, since they focus on truly sparse signals u, while natural images are not well described by sparsity alone(Weiss et al., 2007). The inappropriateness of theoretical properties such as maximal incoherence or RIP as designprinciples for natural image reconstruction has been demonstrated in (Seeger and Nickisch, 2008). However, weobserve in practice that designs optimized by our approach support subsequent sparse MAP reconstruction suc-cessfully, indeed as well or better than other “more Bayesian” estimates linked to the posterior, such as its mean.Our method can be used for design optimization, if the objective is sparse MAP estimation. It solves the problemof “learning compressed sensing” (Weiss et al., 2007) for large scale signals. Moreover, since current scalableMAP estimation codes, or even LS estimation codes, are based on much the same underlying primitives, the addedeffort of setting up our algorithms is minor.The variational relaxation of inference employed here is not novel (Girolami, 2001; Palmer et al., 2006), butpreviously known algorithms for solving it are orders of magnitude slower than our approach on the problems con-sidered here. Moreover, their convergence behaviour is not well characterized. Our contribution basically settlesthe question under which conditions the general approach of (Palmer et al., 2006) leads to a convex optimizationproblem. Our algorithm development owes ideas to (Wipf and Nagarajan, 2008), whose interest is in aggressivelysparse estimation beyond convex MAP. Our framework moves considerably beyond their method, in generality,scope, and practical realization. Our convexity proof is novel, and our methods are applicable to more generalmodels. Our proofs are based on convexity results for certain log determinants, which are novel in machine learn-ing to our knowledge, and may have other applications there. Our interest is in estimating Bayesian uncertainties,information which is essentially destroyed in successful sparse estimation, as will be discussed in some detail. Themajor computational benefit of sparse estimation, namely many variables becoming exactly equal to zero rapidly, isresponsible for this information loss, so that successful sparse Bayesian inference has to be implemented efficientlywithout relying on exact sparsification.The structure of the paper is as follows. The sparse linear model is introduced in Section 2, where we alsodiscuss our variational relaxation of Bayesian inference. In Section 3, we prove novel convexity properties of thisrelaxation. Our scalable algorithms for solving the variational problem are introduced in Section 4. In Section 5,we present some extensions, and discuss the relationship to MAP estimation and algorithmic aspects of sparsity.Bayesian sequential experimental design is discussed in Section 6. We close with a discussion, putting our workinto context, and suggesting applications which would directly benefit from it. Experimental results for our novelmethods will be presented in a later paper.2 Sparse Bayesian Inference. Variational ApproximationsIn this section, we introduce sparse linear models, which can be used to represent image reconstructions togetherwith remaining uncertainties, as noted in Section 1. Bayesian posterior (or inference) computations, which arerequired to compute queries for optimizing the measurement designs, cannot be done exactly for these models. Weintroduce a general variational inference approximation, based on the idea of convex inequalities for the cumulantgenerating function of the posterior.Consider the problem of (natural) image reconstruction introduced in Section 1, which has real-world applica-tions in computational photography and medical imaging (for example, magnetic resonance imaging or positron-emission tomography). The image u ∈ Rn of n pixels has to be estimated from linear measurements y ∈ Rm,3wherem n in many situations of practical interest. Such measurements suggest a linear modely =Xu + ε, ε ∼ N(0, σ2I), (1)whereX ∈ Rm×n is the design matrix, and ε is Gaussian noise of variance σ2. The implied likelihood is Gaus-sian: P (y|u) = N(y|Xu, σ2I), and maximum-likelihood estimation leads to the famous normal equations oflinear least squares (LS) estimation. However, for m n, these methods can work poorly for image reconstruc-tion, mostly because they do not reflect image properties of u at all. For example, it is an established fact of naturethat the projections pTj u of a natural image under zero-mean filters pj , such as nearest neighbour differences (im-age gradient), wavelet or Fourier coefficients, are distibuted in a distictly non-Gaussian way: most coefficients areclose to zero, while a certain fraction have significant sizes (Simoncelli, 1999). If such super-Gaussian properties(also known as sparsity of natural images) are encoded in a prior distribution P (u), reconstructions are typicallyimproved. The non-Gaussianity of P (u) is important here. While a Gaussian prior leads to much simpler, analyti-cally tractable computations, the main improvement comes from non-Gaussian properties. More details about thispoint can be found in (Seeger, 2008).In this paper, we concentrate on priors of the form P (u) ∝ ∏qi=1 ti(si), where s = Bu. In the imageapplications we are interested in here, B may contain local derivative, wavelet, or Fourier filters, but many otherapplications come with this structure. For example, with time series data, B may encode temporal differences.The matricesX andB model the couplings of variables. On the lowest level, exploitable structure in these is whatrenders our algorithms scalable. An example for sparsity-favouring sites ti are Laplace (or double exponential)potentialsti(si) =τ˜i2e−τ˜i|si|, τ˜i = τi/σ > 0. (2)For this particular prior and the Gaussian likelihood (1), the MAP estimator is a (convex) quadratic program, thespecial case B = I is known as Lasso (Tibshirani, 1996). In general, if all log ti(si) are concave in si, MAPestimation reduces to a convex program. Another class of sparsity potentials are of the Student’s t type:ti(si) = (1 + (τi/νi)σ−2s2i )−(νi+1)/2, τi, νi > 0. (3)For these, log ti(si) is neither concave nor convex, and MAP estimation is in general not a convex program. Incomparison with Laplace sites, Student’s t potentials have heavier tails, and for small νi enforce sparsity morestrongly. We refer to a model with likelihood (1) and non-Gaussian sparsity prior P (u) as sparse linear model(SLM).Beyond sparsity prior sites, our framework can be used for other models, whose posterior can be written inthe form P0(u)∏i ti(si), where P0(u) has Gaussian form (but need not be normalizable), and the ti are scalarpotentials. For example, sparsity sites can be used as likelihood terms, to drive robust regression. Or the ti(si)may be binary classification likelihoods. Concrete examples will be given below.2.1 Variational Lower BoundsBayesian inference is not analytically tractable in general for sparse linear models (SLMs) with non-Gaussiansparsity potentials ti(si), and has to be approximated. As motivated in Section 1, we are interested in scalablerelaxations to standard linear optimization primitives, which is why we focus on a variational approach here, withroots in statistical physics.2 Our restriction to SLMs with priors of factorizable form facilitates the exposition.Generalizations to non-Gaussian likelihoods and coupled non-Gaussian sites will be given in later sections.The log partition function logP (y) = log∫P (y|u)P (u) du contains the gist of the posterior,3 and is easierto approximate. It cannot be computed exactly for the SLM, whose posterior is not Gaussian. At this point, weexploit a property of sparsity potentials mentioned above already: a positive even continuous function ti(si) is(strongly) super-Gaussian if gi(xi) = log ti(si), xi = s2i /σ2, is convex and nonincreasing4 for xi > 0 (Palmer2 This does not mean that other approximation techniques, such as Markov chain Monte Carlo, cannot be scalable, only thatequivalent relaxations to provably few calls of standard primitives are harder to establish.3 logP (y) is the cumulant generating function of P (u|y), with a role similar to the generating function for a convergentseries, or the characteristic function for a distribution.4 Fenchel duality works for any convex gi(xi), but if gi is not nonincreasing, the maximization would not be over pii > 0,so that at least for some xi, the closest lower bound to ti(si) would not have Gaussian form (which requires pii > 0). However,since ti(si) = egi(s2i /σ2) is in general a normalizable potential, we can constrain gi(xi) to be nonincreasing without much lossof generality.4et al., 2006). We can represent this convex function using Fenchel5 duality (Rockafellar, 1970, Sect. 12): gi(xi) =maxpii>0−xipii/2− g∗i (−pii/2), resulting inti(si) = maxpii>0e−12σ−2piis2i fi(pii), fi(pii) = e−g∗i (−pii/2).Note that log fi(pii) is concave, since the conjugate function g∗i is convex just as gi. The term “super-Gaussian”becomes clear now: ti(si) has tight Gaussian-form lower bounds of all possible widths.Palmer et al. (2006) remark several interesting facts about this “Gaussianification” step. First, there is a closerelationship to scale mixture decompositions (Gneiting, 1997; West, 1987), where a non-Gaussian density ti(si)is written as a mixture of zero mean Gaussians: ti(si) = Epii [N(si|0, pi−1i σ2)]. It is shown in (Palmer et al.,2006) that all scale mixture sites can be lower bounded as above, although fi(pii) is different from the mixturedensity. Moreover, they show that the variational approximation that arises from these lower bounds, to be detailedright below, is equivalent for scale mixture sites to a different variational principle, known as variational (meanfield) Bayes. More details about these relationships are found in (Palmer et al., 2006; Seeger, 2008), together withreferences to earlier work exploiting special cases of Fenchel duality lower bounds (Girolami, 2001; Figueiredo,2003; Jaakkola, 1997; Wipf et al., 2004). The lower bounds for the Laplace (2) aree−τ˜i|si| = maxpii>0NU (si|0, σ−2pii)e−(τ2i /2)pi−1i , τ˜i = τi/σ, (4)where we define unnormalized Gaussian functions asNU (z|b,P) := exp(−12zTPz + bTz), P 0.For Student’s t sites (3), we obtain(1 + (τi/νi)σ−2s2i )−(νi+1)/2 = maxpii>0NU (si|0, σ−2pii)e−hi(pii)/2, (5)where hi(pii) = (νi + 1)(xi − log(1 + xi))I{xi≤0} with xi = (αi/(νi + 1))pii − 1, αi = νi/τi. Finally, if ti(si)is super-Gaussian, so is a positive mixture of the form∑lαlti(βlsi), αl, βl > 0,because the logsumexp function x 7→ log 1T exp(x) is strictly convex and increasing in each argument, and thelog-mixture is the concatenation of logsumexp and (log ti(βlsi)+logαl)l, the latter being convex and nonincreas-ing for xi > 0 in each component (Boyd and Vandenberghe, 2002, Sect. 3.2.4).The variational relaxation we use here is obtained by plugging the Gaussian-form bounds into P (y), whichresults in the lower boundP (y) ≥ C˜1∫N(y|Xu, σ2I)NU (u|0, σ−2BTΠB) du, C˜1 =q∏i=1fi(pii), Π = diagpi. (6)The right hand side is a Gaussian integral and can be evaluated easily. The variational problem, to be addressedby our algorithms, is to maximize the lower bound w.r.t. the variational parameters pi 0 (i.e. , pii > 0 for alli = 1, . . . , q), with the aim of tightening the approximation to logP (y).The idea of Gaussian-form lower bounding is not restricted to symmetric sites. For example, if ti(si) =t˜i(si)eαisi , so that t˜i(si) is super-Gaussian, we can bound t˜i(si) as before, then multiply eαisi back (which isjust log-linear). In the following, commiting a slight abuse of nomenclature, we will refer to such ti(si) as super-Gaussian. For example, Bernoulli potentials used in binary classification areti(si) =(1 + e−yisi)−1 = eyisi/22 cosh(yisi/2), (7)5 Under some additional conditions on gi(xi), Fenchel duality is equivalent to Legendre duality.5and − log cosh(yisi/2) is even and convex as function of s2i (Jaakkola, 1997, Sect. 3.B). The corresponding con-jugate function is hard to compute analytically, but this is not required in our algorithms, as long as gi(xi) and itsderivative can be computed at any point (see Section 4.3).In the following, we will address two major problems for this relaxation. In the next section, we will answerthe question under which conditions the variational problem of maximizing the lower bound of (6) is a convexoptimization problem. To this end, we will derive some novel convexity results for certain log determinants. Ourfindings have impact on MAP estimation algorithms as well (see Section 5.2).However, we will see in Section 4 that the variational problem even in a convex case is more difficult to solvethan the corresponding MAP estimation. For the latter, gradient computations come at the cost of solving a singlelinear system, while computing a gradient of the variational lower bound w.r.t. pi is much more difficult than asingle system. Most previous algorithms for variational problems of this kind (Girolami, 2001; Tipping, 2001;Minka, 2001) avoid this difficulty by following a step-wise approach, optimizing w.r.t. single components of pi inturn, keeping all others fixed. While an informed scheduling of updates (Tipping and Faul, 2003; Seeger and Nick-isch, 2008) can render these algorithms feasible on problems of moderate size, their scalability is fundamentallylimited.6 In difficult image reconstruction settings, every site approximation NU (si|0, σ−2pii) has to be visited atleast once (usually several times). But for large n in the hundred thousands, every single scalar update requires theequivalent of a reweighted least squares estimation. And for problems considered here, where q can be a millionor more, O(q) LS estimations simply cannot be done.Our second contribution lies in the development of novel classes of algorithms that can cope with such large scaleproblems, by decoupling the lower bound criterion complexity in a nested double loop fashion. These algorithmscan be applied in domains where single-site updating is not an option. Also on problems of moderate size, speedupsby orders of magnitude are realized. These algorithms can be applied in general, whether the variational problemis convex or not.3 Convexity Properties of Variational InferenceWe discussed a general variational relaxation of Bayesian inference for generalized sparse linear models in theprevious section, to the maximization of the right hand side of (6). In this section, we bring this optimizationproblem into a form more amenable to algorithmic treatment. The relaxation is a special case of variational (meanfield) Bayes (Ghahramani and Beal, 2001; Attias, 2000), or of direct site bounding (Jordan et al., 1997). We answerthe question under which conditions the variational problem is a convex optimization problem. To this end, weprove a number of novel convexity properties for parts of the upper bound criterion, thereby laying groundworkfor our novel scalable algorithms introduced in Section 4 as well.Our problem is to maximize the right hand side of (6). Assume for now that BTΠB is invertible. Theend result remains valid even if this is not the case, as is easily seen by a continuity argument. Let Q(u) :=C−12 NU (u|0, σ−2BTΠB) and Q(y,u) := P (y|u)Q(u). The joint distribution is Gaussian, andQ(u|y) = N(h, σ2Σ), Σ−1 = A =XTX +BTΠB. (8)We have that∫P (y|u)Q(u) du = Q(y) = |2piσ2Σ|1/2maxuQ(y)Q(u|y) = |2piσ2Σ|1/2maxuP (y|u)Q(u), (9)because the maximum of the GaussianQ(u|y) is attained at the mean h. The crucial step here is that we can movefrom the integral over u to the maximum over u exactly, which is possible becauseQ is Gaussian. We end up withP (y) ≥ C1e−φ(pi)/2, whereφ(pi) := log |A|+ h(pi) + σ−2(minu‖y −Xu‖2 + sTΠs),h(pi) =q∑i=1hi(pii) = −2q∑i=1log fi(pii),(10)6 These comments apply to sparse Bayesian inference. For sparse estimation, forward selection schedules may well workfor very large n and q, namely because many of the pii are clamped to∞ and never once moved from there (see Section 5.2).6and C1 = (2piσ2)(n−m)/2. The variational problem now consists in maximizing the lower bound, or equivalentlyminimizing φ(pi), w.r.t. the variational parameters pi. The Gaussian posterior approximation is Q(u|y), with thefinal parameters pi plugged in. An algorithmically beneficial side effect of using a lower bound on logP (y), ratherthan just an approximation (as other variational methods, such as expectation propagation (Minka, 2001) do), isthat we can usually devise algorithms which provably converge to a local optimum. Still, φ(pi) is in general acoupled, non-convex function of non-standard form, and even ∇piφ is hard to compute.3.1 Some Convexity ResultsIn this section, we establish some convexity properties of φ in (10). A crucial term in this criterion is log |A|, whereσ2A−1 is the covariance matrix of the approximate posterior (8). At least in hindsight, much of the computationaldifficulty is caused by this term (see Section 5.2 for a detailed discussion). Let γ := pi−1,i.e. γi = 1/pii. It turnsout that log |A| has a number of convexity properties in terms of pi or γ, which are important for obtaining scalablealgorithms, or characterizations of the variational problem in the first place.Theorem 1 LetX ∈ Rm×n,B ∈ Rq×n be arbitrary matrices, so thatlog |A|, A =XTX +BTΠB, Π = diagpiexists for all pi 0 from an open set. Let γ := pi−1, Γ = diag γ.1. pi 7→ log |A| is concave.2. γ 7→ log |A| is convex.3. γ 7→ log |A|+ log |Γ| is concave.4. Let ρi(γi) be concave functions into R+. Then, pi 7→ log |XTX +BT ρ(Π)B| is concave, where ρ(Π) =diag(ρi(pii)).5. Let ρi(γi) be twice continuously differentiable functions into R+, so thatρ′′i (γi)ρi(γi) ≥ (ρ′i(γi))2for all i and γi. Then, γ 7→ log |XTX +BT ρ(Γ)B| is convex, where ρ(Γ) = diag(ρi(γi)).6. Let ρi(γi) be concave functions into R+. Then, γ 7→ log |ρ(Γ)| + log |XTX +BT ρ(Γ)−1B| is concave,where ρ(Γ) = diag(ρi(γi)).7. Let Q(u|y) be the approximate posterior of (8). Then,σ−2VarQ[si|y] = δTi BA−1BT δi ≤ γi.For the proof, (1.) is well known (Boyd and Vandenberghe, 2002, Sect. 3.1.5). The generalization (4.) followsfrom (Boyd and Vandenberghe, 2002, Sect. 3.2.4), since pi 7→ log |A| is nondecreasing in each pii. (2.) followsfrom the more general (5.), using ρi(γi) = γ−1i , and (3.) follows from (6.), using ρi(γi) = γi. To our knowledge,(5.) and (6.) are novel, at least we do not know of previous appearences in machine learning. They are proved inAppendix A.1. The proof of (7.) is a part of the proof of (5.).Note that (5.) is more general than what we require in the following. For example, it holds for all ρi(γi) = γ−βii ,βi > 0. For ρi(γi) = eγi , we obtain the convexity of γ 7→ log |XTX+BT exp(Γ)B|, generalizing the logsumexpfunction x 7→ log 1T exp(x) (Boyd and Vandenberghe, 2002, Sect. 3.1.5) to matrix values. The convexity of thelatter is behind many properties of exponential families or of maximum-likelihood estimation in log-linear models.Note also that (7.) gives a precise characterization of γi as sparsity parameter regulating the variance of si. Asimilar argument shows that the size of EQ[si|y] is also regulated by γi.What about the remaining terms in (10)?Here and in the following, we treat φ, h, and fi as functions7 of pi or γ. Based on Theorem 1, we show thatγ 7→ φ(γ)− h(γ) is a convex function.7 In general, we adopt the physics convention of treating function values as dependent variables, invariant to reparameteri-zations of the variables they depend on.7Theorem 2 Our variational relaxation of approximate inference requires the minimization of φ(γ) from (10). Thefunction γ 7→ φ(γ)− h(γ), the negative log partition function of a Gaussian, is convex for γ 0.The convexity of log |A| has been shown in Theorem 1, part (2.). ‖y−Xu‖2 is convex in u. More interestingly,(u,γ) 7→ sTΓ−1s is jointly convex, since the quadratic-over-linear function (si, γi) 7→ s2i /γi is jointly convexfor γi > 0 (Boyd and Vandenberghe, 2002, Sect. 3.1.5), and s = Bu is linear in u. Now, since minu κ(u,γ) isconvex if κ is jointly convex (Boyd and Vandenberghe, 2002, Sect. 3.2.5), we conclude that γ 7→ φ(γ) − h(γ) isa convex function.To put this result into context, note that φ − h is the negative log partition function of a Gaussian with naturalparameters being linear in pi. It is therefore well known that pi 7→ φ−h is a concave function. However, pi 7→ h isconvex (hi(pii) is the conjugate function for 2g(xi)), which does not lead to insights about pi 7→ φ. The convexityof the negative log partition function of a distribution with natural parameters linear in γ−1 seems specific to theGaussian case, and is novel to our knowledge.Given Theorem 2, if all hi(γi) are convex, the whole variational problem minγ0 φ is a convex minimization.In the next section, we establish properties for when this is the case.3.2 Convex Variational InferenceThe statement of Theorem 2 implies that whenever γ 7→ h(γ) is convex, the whole variational problemminγ0 φof interest here is a convex problem with a unique solution. In this section, we show that whenever log ti(si) areconcave functions, and ti(si) are super-Gaussian (see Section 2.1), h(γ) is in fact convex. Moreover, at least ingeneral, this property is necessary for the convexity of h(γ), and for the convexity of the variational problem. Fornotational simplicity, we do not deal with the most general case our result applies for here. For example, P (y|u)can be replaced by any factor of Gaussian form in u, and the ti(si) can be likelihood sites, or depend on multiplecomponents of s. Generalizations are given in Section 5.1.Theorem 3 Consider a model with Gaussian likelihood (1) and a prior P (u) ∝ ∏qi=1 ti(si), s = Bu. Leteach ti(si) be strongly super-Gaussian, meaning that ti(si) = eαisi t˜i(si), where t˜i(si) is an even function,and gi(xi) = log t˜i(si), xi = s2i , is strictly8 convex and nonincreasing for xi > 0. Moreover, suppose thatgi(si) = log t˜i(si) is concave9 and twice continuously differentiable for si > 0. Then, the variational problem ofminimizing (10) over γ 0 is a convex optimization problem.In general, the requirements on the ti(si) are necessary for the convexity to hold. If g′′i (si) > 0 for some si > 0,then hi(γi) is not convex at some γi > 0. For generalX,B, and y, this means that φ(γ) is not convex either.The proof is given in Appendix A.2. Our theorem provides a satisfying characterization of the variationalinference relaxation of Section 2. MAP estimation (see Section 1) is a convex problem if and only if (in general)all log ti(si) are concave. Whenever MAP estimation is convex, and the ti(si) are super-Gaussian, the variationalrelaxation is convex as well. Loosely speaking, it is the log-concavity of the posterior that renders the variationalproblem convex. This property sets the relaxation used here apart from all other approximate inference methodsfor continuous variable models we know of: most of these can be shown to be non-convex in general, even if thelog posterior is concave and has a single mode only.10With a view to Section 4, it is also interesting to ask under which conditions the hi(γi) are concave functions. Wedo not have a complete characterization for this case, but can give some examples. If ti(si) = e−τ |si|α, α ∈ (0, 2],then hi(γi) ∝ γβi with β = −α/(α − 2), which is convex iff α ≥ 1. In these cases, ti(si) is log-concave. Forα < 1, hi(γi) is concave. Moreover, if ti(si) is log-convex in si, thenhi(γi) = maxsi−σ−2s2i /γi − 2 log ti(si) = −minsi(σ−2s2i γi + 2 log ti(si)),which is concave in γi, since the argument of minsi is jointly convex in (si, γi) (Boyd and Vandenberghe,2002, Sect. 3.2.5). Therefore, given that ti(si) is super-Gaussian, its log-{concavity/convexity} implies con-vexity/concavity of hi(γi). The reverse is not true in general (for the second statement), as can be seen forti(si) = e−τ |si|αabove.8 We require a slightly stronger notion of super-Gaussianity here, in that gi(xi) has to be strictly convex.9 Here and elsewhere, we understand function values as variables dependent on their arguments, so that gi(si) = gi(xi).This convention, widely used in physics, simplifies notation and should not lead to confusions.10 An example is expectation propagation (Minka, 2001), whose log partition function approximation is non-convex (Op-per and Winther, 2005), except for trivial cases. Log-concavity of the ti(si) has important consequences for the numericalproperties of the EP algorithm (Seeger, 2008), but they do not imply convexity of the complete problem.8We close this section with some examples. For Laplace sites (2), the Fenchel duality is given by (4), wherehi(γi) = τ2i γi, a convex function as predicted by our result above. For sparse linear models with Laplace sites,MAP estimation is a convex quadratic program, with the Lasso as a special case (see Section 2). Variationalinference is a convex problem as well. While the same relaxation has been used before for these SLMs (Girolami,2001), its convexity has not been established until now.Second, binary classification Bernoulli likelihood sites (7), also known as logistic potentials, are super-Gaussian(see Section 2.1), and they are well known to be log-concave. MAP estimation for generalized linear models withthese sites is known as logistic regression, a convex problem typically solved by the iteratively reweighted leastsquares (IRLS) algorithm (also known as Fisher scoring). Variational inference for this model is a convex problem,and our algorithms introduced in Section 4 make use of IRLS as well.However, Student’s t potentials (3) are not log-concave, and hi(γi) in (5) is neither convex nor concave. NeitherMAP estimation nor variational approximate inference is a convex problem, when Student’s t sites are used.We have provided a satisfying characterization of a widely used class of variational approximate inferencemethods. For super-Gaussian sites, the variational problem is convex if and only if the search for the posteriormode is convex. This does not mean that solving the variational problem is computationally as tractable as MAPestimation (see end of Section 2.1). For example, our posthoc result that the algorithm of (Girolami, 2001) solvesa convex problem, is of little value for measurement design on full-size images, where this algorithm would notconverge in any reasonable amount of time. In the next section, we propose classes of algorithms that solve thevariational problem in a scalable way. While they are still in general more expensive than convex MAP estimationmethods, the precise relationship is clarified in Section 5.2.4 The AlgorithmsIn this section, we develop two classes of algorithms for scalable variational inference, maximizing the lowerbound of (6) for SLMs, where the prior is P (u) ∝ ∏qi=1 ti(si), s = Bu, with super-Gaussian sites. We useFenchel concave duality in order to upper bound the variational criterion by a decoupling surrogate, which can beminimized by a standard convex optimization algorithm.These algorithms are independent of our convexity analysis for the variational problem in Section 3. They canbe applied in order to find a local minimum very efficiently in general. If the problem is convex, there is a uniqueglobal optimum, which is found by our methods.4.1 The First ClassWe propose two closely related classes of scalable algorithms for minimizing φ of (10) w.r.t. pi 0 (or equiva-lently, w.r.t. γ := pi−1 0). In this section, we introduce the first class. Whether φ is convex in γ or not (seeTheorem 3), it is not obvious how to minimize φ tractably, since even the computation of ∇γφ will be seen to bea computationally expensive problem.We make use of a powerful general idea known as double loop algorithms, concave-convex algorithms, ord.c. programming (difference of convex). Special cases of such algorithms are already heavily used in machinelearning and statistics: the expectation-maximization method (Dempster et al., 1977), variational (mean field)Bayesian inference (Attias, 2000), or CCCP for approximate inference (Yuille and Rangarajan, 2003), amongmany others. The idea is to write φ as sum of a concave and a convex part. We use Fenchel-Legendre dualityonce more, in order to upper bound the concave part by a linear function, then minimize the convex upper boundto φ globally. The linear upper bounding is done iteratively, in so called outer loop steps, followed by innerloop convex minimizations. If the concave part is differentiable, the linear upper bound is a tangent plane, andunder mild conditions the resulting double loop algorithm can be shown to be globally convergent, meaning that itconverges to a local minimum point of φ, no matter from where we start.In Theorem 1, part (3.), we show that γ 7→ log |A|+log |Γ| is a concave function. Now write h(γ)− log |Γ| =h∩(γ) + h∪(γ), where h∩ is concave, and h∪ is convex. Note that h(γ) − log |Γ| decomposes as sum of scalarterms, each depending on a single γi only, and it should in general be possible to find h∩, h∪ which decomposein the same way. A concave-convex decomposition is never unique: we can always add a concave part to h∩ andsubtract it from h∪. However, since we will replace the concave part by a hyperplane, it is sensible to choose adecomposition with “minimal” h∩, as close to linear as possible. For tractability, it is also important that bothparts are simple terms composed of standard functions. For example, for Student’s t sites (3), hi(γi) is given in(5), where α = ν ◦ τ−1. A convenient decomposition is h∩(γ) = νT (log γ) and h∪,i(γi) = αi/γi + (νi +91)(log(νi + 1) − logαi − 1) for γi ≥ αi/(νi + 1), h∪,i(γi) = −(νi + 1) log γi for γi < αi/(νi + 1). Note thatthe latter is convex and twice differentiable.Now, if g(γ) := log |A|+ log |Γ|+ h∩(γ), which is concave as sum of concave functions, thenφ(γ) = g(γ) + minuh∪(γ) + σ−2(‖y −Xu‖2 + sTΓ−1s) .Using the conjugate representation g(γ) = minz0 zTγ − g∗(z), we obtainφ(γ) ≤ minuφz(u,γ), φz(u,γ) := zTγ + h∪(γ) + σ−2(‖y −Xu‖2 + sTΓ−1s)− g∗(z)We saw in the proof of Theorem 2 that φz is jointly convex. Our algorithm to minimize φ iterates between outerloop updates of z ← argminzTγ − g∗(z) and inner loop convex minimizations of φz , in order to update γ. As aby-product, we obtain minimum points u, and it is easy to see that at the end of an outer loop, u is the mean11 ofthe current posterior approximation Q(·|y).The inner loop minimization of φz can be done in any order of u and γ. If h∪(γ) =∑i h∪,i(γi), as is thecase in general, it is easiest to perform the minimization over γ first, since this decouples into q independentminimizationsminγiziγi + h∪,i(γi) + σ−2s2i /γi,leading to the equations zi + (dh∪,i)/(dγi) − σ−2s2i γ−2i = 0. For the examples given here, these equationscan be solved analytically, but in general univariate convex minimization can be used in order to solve for therequired quantities (see Section 4.3). Plugging the solutions in, we obtain convex functions h∗∪,i(si). For the caseof Student’s t sites, we have h∪,i(γi) = αi/γi, therefore h∗∪,i(si) = 2(zipi)1/2, pi = αi + s2i /σ2 (recall thatαi = νi/τi). The remaining inner loop problem isminuσ−2‖y −Xu‖2 +q∑i=1h∗∪,i(si), s = Bu. (11)This problem is of standard form, and can be solved generically by the iteratively reweighted least squares (IRLS)algorithm, a variant of the Newton-Raphson method. This method, which is discussed in detail in Section 4.3, pro-ceeds in Newton steps, each of which requires the solution of a single linear system with a matrix of the same formas A (8), only that Π is replaced by a different positive diagonal matrix, and a simple line search. Convergenceis typically attained after few iterations (less than thirty), and each step is reduced to a single (reweighted) leastsquares problem. Therefore, the inner loop minimization fulfils our criteria of scalability.The outer loop updates of z, given γ, require the minimization of zTγ − g∗(z) for the concave function g∗(z).It is not in general possible to analytically obtain g∗(z). By duality, g∗(z) can at any point be evaluated byminimizing zTγ − g(γ) over γ, which is a convex problem, but is hard in general. Fortunately, none of this isnecessary. For fixed γ, the minimizer zopt is such thatzToptγ − g(γ) = g∗(zopt) = minγ˜zToptγ˜ − g(γ˜),so that ∇γzToptγ − g(γ) = zopt − ∇γg(γ) = 0. We just have to compute the gradient of g at γ. Using the factthat d log |A| = trΣ(dA) = trBΣBT (dΠ), we arrive atzopt = ∇γg(γ) = pi +∇γh∩(γ)− pi2 ◦ diag−1(BΣBT).Since h∩(γ) =∑i h∩,i(γi), the main difficulty is the computation of the diagonal of BΣBT . A glance at (8)reveals that we need to compute Gaussian marginal variances: (VarQ[si |y])i = σ2 diag−1(BΣBT ). This isthe second primitive we require in our algorithms. As opposed to least squares estimation, this computation is notusually done in non-Bayesian applications. It is important to note that we require allmarginal variances in parallel.A single variance could be obtained as the solution of a least squares problem, but computing zopt by solving qsuch problems is certainly not feasible. Nevertheless, there are algorithms from numerical mathematics which canbe used to estimate these variances, and they reduce to the same primitives as does least squares estimation. In11 We abuse notation slightly by treating u as variable to optimize over here.10practice, the number of outer loop updates required to gain convergence is small, usually five or six iterations aresufficient. The variance estimation problem is discussed in more detail in Section 4.5.This concludes the specification of our first class of algorithms. To understand the considerable computationalbenefits of the double loop structure, note that a single gradient computation∇γφ is as expensive as an outer loopupdate in our scheme: both require all q Gaussian marginal variances. By bounding the log |A| part, which causesthe complexity (see Section 5.2), we reduce the number of expensive steps required until convergence drastically,compared to standard gradient descent algorithms. We follow the double loop strategy even if h(γ) is convex, sincethe decoupled inner loop criterion is much more efficient to minimize than φ itself. Some general characteristicsof our algorithms are given in Section 4.4.4.2 The Second ClassOur second class of algorithms is closely related to the first one, but involves another twist. Again, our aim is todevise a double loop scheme. From Theorem 1, part (1.), we know that pi 7→ log |A| is concave. We can try todecompose h = h∩(pi) + h∪(γ) (recall that γ = pi−1), h∩ concave in pi, and h∪ convex in γ. For example, forLaplace sites (4), we have that h = (τ 2)Tγ, which is convex, so that h∪(γ) = (τ 2)Tγ and h∩(pi) = 0.We can now proceed much as for the first class above, but treating g(pi) := log |A|+h∩(pi) as concave functionin pi. By Fenchel’s inequality (Rockafellar, 1970, Sect. 12): g(pi) ≤ zTpi − g∗(z) for z 0. But this upperbound is also a convex function of γ 0 (since zi > 0 for all i), which is the additional observation we make useof here. We obtain the convex inner loop problemminγminuzT (γ−1) + h∪(γ) + σ−2(‖y −Xu‖2 + sTΓ−1s)− g∗(z),which is treated much as in Section 4.1. The independent minimizations are nowminγi h∪,i(γi)+(zi+σ−2s2i )/γi.If pi = zi + σ−2s2i , the stationary equations are (dh∪,i)/(dγi) − piγ−2i = 0. Plugging the solutions in gives riseto convex functions h∗∪,i(si), and the inner loop optimization is done by the IRLS algorithm. For Laplace sites,h∪,i(γi) = τ2i γi, so that h∗∪,i(si) = 2τip1/2i . For the outer loop update,zopt = diag−1(BΣBT)+∇pih∩(pi).If h∩ decomposes, the main effort here is again to compute the Gaussian variances of Q(s|y).Characteristics of the second class, as well as a comparison with the first, are given in Section 4.4. The mainidea is the same as for the first class: difficult coupling terms in φ are bounded by simple decoupling functions, sothat the complexity of computing ∇γφ is shifted into a few outer loop updates. If R = σ−2minu [‖y −Xu‖2 +sTΓ−1s], the relationship between the two classes is given byφ[1]z ≥ log |A|+ log |Γ|+ (h− log |Γ|)︸︷︷︸zTγ−g∗(z)+h[1]∪ ≥+R = φ = log |A|+ h︸︷︷︸≤zT (γ−1)−g∗(z)+h[2]∪+R ≤ φ[2]z .The superscripts distinguishing the two classes will be dropped in the sequel, where it will always be clear fromthe context which class is used.Hybrid VariantThe coupling term log |A| is bounded in either of our classes, making use of Fenchel’s inequality for concavefunctions. In cases where the second class does not apply, the following hybrid variant can be used. Suppose thath(γ) = h∩(γ) + h∪(γ), where h∩ is concave, h∪ is convex. We can upper bound log |A| ≤ zT (γ−1) − g∗(z)and h∩(γ) ≤ zT∩γ− g∗∩(z∩). Since h∩(γ) decouples, z∩ and g∗∩(z∩) are easy to compute. In the inner loop, γi iseliminated byminγi h∪,i(γi) + pi/γi + z∩,iγi, where pi = zi + σ−2s2i . For Student’s t sites (see Section 4.1), wehave that h∗∪,i(si) = 2(z∩,i(αi+ pi))1/2, and finally γ′i = ((αi+ p′i)/z∩,i)1/2. Moreover, the outer loop update isz∩ = (ν + 1) ◦ pi.4.3 The Inner Loop OptimizationIn this section, we show how to efficiently solve the inner loop problem (11). The structure of φz(u) is the sum ofa least squares term and a decomposing penalty. Therefore, the Newton-Raphson algorithm reduces to a standard11method called iteratively reweighted least squares (IRLS). We describe a single Newton step here, starting from u.Let r := y −Xu denote the residual vector. Then,dσ−2‖r‖2 = −2σ−2rTX(du), d2σ−2‖r‖2 = 2σ−2(du)TXTX(du).If θi := (σ2/2)(dh∗∪,i)/(dsi), ρi := (σ2/2)(d2h∗∪,i)/(ds2i ), theng := ∇uφz = −2σ−2(XTr −BTθ) ,H := ∇∇uφz = 2σ−2(XTX +BT (diagρ)B).Note that ρi ≥ 0, by the convexity of h∗∪,i. The Newton search direction isd := −H−1g = (XTX +BT (diagρ)B)−1 (XTr −BTθ) .The computation of d requires to solve a system with the matrix H, which is of the same form as A. This isprecisely the computation required for least squares estimation with the likelihood P (y|u) and the Gaussian priorN(s|0, (diagρ)−1). Such systems are generally solved approximately by the linear conjugate gradients (LCG)algorithm (Golub and Van Loan, 1996). The cost per iteration of LCG is dominated by a MVM with H, whichtranslates to single MVMs withX, B, and its transposes respectively. Our scalability requirements are thereforemet.A line search along d can be run in negligible time. If f(t) := (σ2/2)φz(u + td), thenf ′(t) = −(Xd)Tr + t‖Xd‖2 + (Bd)Tθ(t).Here, θ(t) is in terms of s(t) = s + tBd. If we precompute Xd, (Xd)Tr, ‖Xd‖2, and Bd, f(t) and f ′(t)can be evaluated in O(q). No further MVMs are required during the line search. Each line search is started witht0 = 1. Note that a line search seems essential in practice. Especially in the beginning, or after X has just beenextended in a sequential design loop (see Section 6), a full Newton step (t = 1) would lead to a large increaseof the criterion, and significantly smaller t need to be taken. During later stages, the first choice t = 1 is usuallyaccepted, due to the self-scaling properties of the Newton method. Finally, once u′ = argminφz(u) is found, γ′is updated as minimizer w.r.t. γ, solving the scalar stationary equations once more.For the case of Laplace sites and the second class, recall that h∗∪,i(si) = 2τip1/2i , pi = zi + σ−2s2i . Therefore,θi = τisip−1/2i and ρi = τizip−3/2i . For Student’s t sites, h∗∪,i has the same form, but τ2i is replaced by zi, andzi by αi = νi/τi. Finally, if h∗∪,i(si) cannot be determined analytically, the procedure detailed in Appendix A.3can be used, which does not increase the computational complexity. In some cases, even hi(γi) may not be knownexplicitly, or may be cumbersome to obtain analytically (see Section 2.1). We show in Appendix A.3 how ouralgorithms can be run based on computations of gi(xi) and its derivative only.4.4 Properties of the AlgorithmsIn this section, we analyze characteristics of the two classes of algorithms, relating them to each other, and showinghow they compare with the sparse estimation method of (Wipf and Nagarajan, 2008), which inspired some of ourwork here. Properties of the variational problem they adress, are analyzed in Section 3.First, the algorithms in either class can be shown to converge globally, i.e. to find a local minimum point from anystarting point, if the outer loop updates are done exactly. This is seen just as in (Wipf and Nagarajan, 2008), whosearguments apply here as well. If the variational problem itself is convex (see Section 3), the algorithms converge tothe global optimum from any starting point. In a nutshell, the argument goes as follows. The upper bound φz hasthe same tangent plane at γ than φ. Therefore, the inner loop optimization is guaranteed to decrease φ substantiallyif ∇γφ is not equal to zero. Note that the convergence proof does not require the inner loop optimization to findthe minimum of φz . In fact, a single line search along the first Newton direction would be sufficient. On the otherhand, the outer loop updates have to be done to high accuracy to retain the guarantee, which can be problematic inlarge scale settings. This point is discussed in more detail in Section 4.5.Our first class of algorithms can be seen as generalization of the sparse estimation method in (Wipf and Nagara-jan, 2008). While they considered the special case B = I only, their method can be generalized to any B, usingthe facts proved in Theorem 1. It is obtained as limit of the Student’s t case in Section 4.1, setting νi = 0 for all12i. In this case, referred to as automatic relevance determination (ARD), the prior sites ti(si) are not normalizable.Moreover, as can be seen from Section 4.3, their inner loop criterion becomesσ−2‖y −Xu‖2 + 2σ−1q∑i=1z1/2i |si|.Solving for u is a quadratic program, whose special case for B = I is the Lasso (see Section 2). Due to thenondifferentiable criterion, special code is required (but see Section 5.2), while our inner loops can be implementedmore easily. On the other hand, for the inner loop solution u′, many components in s′ = Bu′ are exactly zero,and this fact can be exploited to find u′ very efficiently. Moreover, much effort is concentrated on designing fastalgorithms for Lasso, at least for certain B (unfortunately the most studied case, B = I , for which very efficientsoft thresholding algorithms are known, is not very useful in the context of natural images). The exact sparsityof s′ is even more useful for computing outer loop updates efficiently. Namely, since γ′i = σ−1z−1/2i |s′i|, we seethat γ is exactly sparse just as s′. If there are d non-zeros in γ, the computation12 of zopt can be done based on asystem matrix of size d×d. However, this computationally beneficial role of sparsity can be expected in the sparseestimation context only. This important point is discussed in more detail in Section 5.2. Variational approximateinference, with the aim of a useful uncertainty representation, is harder to do in practice than sparse estimation,and we show how the added complexity can be addressed computationally.How generally applicable is either class of algorithms? First, we require a Gaussian-form lower bound on ti(si).If the site is super-Gaussian, Fenchel duality provides a tight lower bound. Non-symmetric sites can also be dealtwith in our framework, as long as lower bounds are known (see Section 2.1 for binary classification Bernoullisites). If they make use of two parameters,i.e. are of the form ti(si) ≥ NU (si|bi, pii)fi(pii, bi), the extension givenin Appendix A.4 may be applicable.Second, in order to apply the first class, we require a decomposition of h(γ) − log |Γ| into a sum of concaveand convex functions. Such a decomposition can always be found in principle, although we also require that theparts have a simple tractable form, and are at least differentiable. In all cases of convex log ti(√·) we know of, thissecond requirement does not limit the applicability. We have seen above that h∗∪,i(si) may not be differentiable,which precludes the direct usage of IRLS for solving the inner loop problems. However, ARD above is the onlycase we know of where this happens, and the ARD sites ti(si) do not correspond to normalizable priors.On the other hand, in order to apply the second class of algorithms, h(γ) has to be decomposed as h∩(pi) +h∪(γ). This is simple if h(γ) is convex itself, since h∩(pi) = 0 then. Laplace sites (2) and Bernoulli potentials (7),as well as all super-Gaussian, log-concave potentials, can be treated this way. However, a simple decompositiondoes not seem to exist in many other cases.Note that whenever h(γ) is convex, the first class of algorithms can be applied as well, using h∪(γ) = h(γ)−log |Γ|. However, in this case, the second class seems to be the more direct approach. For Laplace potentials, theinner loop criterion for the second class is more sparsity-enforcing than for the first.The SLM can be configured with different sparsity prior sites. In this paper, Laplace and Student’s t sites aretreated explicitly. For the former, variational inference is convex. For the latter, the posterior is multi-modal,and our algorithms search for a local optimum. The choice of sparsity potentials should depend on the problemaddressed.How are the algorithms initialized? In practice, we found it useful to start with z = ε1 for some ε > 0 (say:1/20), and with u = 0. We also explored the alternative of setting pi ∝ 1, u = 0, and start with an update of z,but this led to less stable behaviour and required more running time.4.5 Estimation of Gaussian VariancesRecall from Section 4.1 that the outer loop steps in our algorithms require the estimation of Gaussian variances z˜ =diag−1(BΣBT ) = σ−2(VarQ[si|y])i, which would also be required for the gradient computation∇γφ. Variancecomputations are not usually required in the context of sparse MAP estimation (see Section 5.2). Certainly, thevector z˜ cannot be estimated by solving a single or few linear systems. In this section, we discuss a generic wayof estimating these variances en bulk, using the Lanczos algorithm from numerical mathematics. Further details12 The computation of zopt is slightly complicated if some γi = 0. From Section 4.1, we have that zi = pii(1−piiz˜i), wherez˜i = 0, and pii =∞. If the subscripts “0” (“1”) denote the part of γ which is zero (non-zero), andA1 =XTX+BT1 Γ−11 B1,a careful computation shows that z0 = diag−1([B0A−11 BT0 ]−1). The zi do not in general become zero.13are given in Appendix B. This section is technically more difficult than the others, and can be skipped in a firstreading.If the precision matrix A of the Gaussian posterior approximation Q(u|y) is sparse, the variances can beestimated in some cases using Gaussian belief propagation. For example, the algorithm in (Wainwright et al.,2001) can be used, which works by embedding a sequence of spanning trees, then does Gaussian propagation onthese trees. Another interesting approach for Gaussian Markov random fields is given in (Malioutov et al., 2006a).However, A is neither sparse nor has graphical model structure13 in general. In (Schneider and Willsky, 2001),the Lanczos algorithm (Golub and Van Loan, 1996) is used in order to estimate Gaussian posterior variances.The connection can be understood by noting that a single marginal variance can be estimated by running linearconjugate gradients (LCG) for the system Av = BT δi, then z˜i = δTi Bv. One way to regard the Lanczosalgorithm is that it explicitly builds up a low rank representation ofA, which allows to solve many linear systemsin parallel with the same matrix A, but different right hand sides. More specifically, k iterations of Lanczosproduce Q(k) = (q(1) . . . q(k)) ∈ Rn×k with orthonormal columns, and T (k) ∈ Rk×k tridiagonal, such thatQ(k)TAQ(k) = T (k). The extremal eigenvectors ofA can be well approximated in the column span ofQ(k) evenfor small k, and the convergence of such eigenvector estimates can be efficiently monitored within the Lanczosrecursion. A Lanczos iteration requires a single MVM withA, thus the same effort in principle than a single LCGiteration.We obtain a Lanczos estimator for z˜ by replacing Σ in the exact expression byQ(k)T (k)−1Q(k)T , resulting inz˜(k) := diag−1(BQ(k)T (k)−1Q(k)TBT ). It is shown in Appendix B that z˜(k) = z˜(k−1) + v(k) ◦ v(k), wherev(k) obeys the recursion v(k) = e−1k (Bq(k) − dk−1v(k−1)). Therefore, z˜(k) converges against z˜ monotonicallyfrom below14 in every component. We can also estimate log |A|, featuring in the criterion part g∗(z), by log |T (k)|,although the latter value is not critically required in order to run our methods. While other estimators could bederived from the Lanczos method, the monotonicity of z˜(k) would be lost (see end of this section). Moreover, thecomponents of z˜(k) correspond to the best estimate after k iterations of LCG for the systems Av = BT δi, if thesame starting vector is used for all of them.While the basic recurrence of the Lanczos algorithm is easy to describe, this apparent simplicity is treacherous.Useful Lanczos implementations involve a lot of extra technology, and in general they require O(nk) memory(details are given in Appendix B).In general, a Lanczos estimate is good if the full expression dominantly depends on the extremal eigenvalues andeigenvectors, which tend to converge rapidly. An ideal case, for which most Lanczos codes seem to be optimized,is given by a matrixA with a geometrically decaying spectrum bounded away from zero. This implies a geometricdecay of the spectrum of Σ as well, so that z˜(k) converges to z˜ rapidly. In applications with a geometric spectraldecay of theA matrices, modern Lanczos codes, such as ARPACK15, can be used. Unfortunately, for problems ofour interest here (high-quality image reconstruction from non-local measurements), the system matrices A oftenexhibit a spectral decay which is roughly linear (with the possible exception of beginning and end). In such cases,‖z˜(k) − z˜‖/‖z˜‖ converges approximately linearly, and a close approximation of z˜ requires k ≈ n iterationsof Lanczos, which is not tractable. Moreover, since the Lanczos iteration delivers optimal approximations ofeigenvalues and eigenvectors, given that k MVMs with A can be used, there may be no better generic method toapproximate z˜, if no additional knowledge about A is used. To be clear, the problem in these cases is not thatA is ill-conditioned, but that the dependence of z˜ on the eigendecomposition of A is spread across much of thespectrum. In other words, A does not have any low rank approximation which is close in the spectral norm. Innumerical mathematics, structural knowledge about the system matrix is used in order to precondition the Lanczosalgorithm, which can in principle be used to improve its spectral properties. Such preconditioning strategies arenot in the scope of this paper, but remain an important topic for further investigations.To conclude so far, we have seen that the Lanczos estimator z˜(k) for z˜ can be computed using k MVMs withA. On the other hand, for certain system matrices A arising in practice (in problems of image reconstruction),z˜(k) will not be very close to z˜ in many components. This creates a problem for the global convergence proofof Section 4.4, which relies on exact computations of z˜. In these situations, we cannot claim that our method isprovably globally convergent in practice. Moreover, since the g∗(z) term in φz cannot be computed exactly, wecannot check for sure whether an outer loop step improved the criterion. From our experience on such problems,13 To be precise,A would be the precision matrix of the Gaussian graphical modelQ(u|y), and its sparsity pattern encodesgraph structure. In cases where such an interpretation applies, this structure is normally independent of the value of pi.14 From Theorem 1, (7.), we know that z˜i ∈ [z(k)i , γi] for any i and k.15 Available at http://www.caam.rice.edu/software/ARPACK/.14the algorithm is well behaved, in that the approximate φz rapidly decreases, then jitters slightly around a value andcan be stopped. Since even∇γφ cannot be computed reliably, it is questionable whether a generic precise stoppingrule can be found. Establishing such rules in particular cases of interest remains an open problem for future work.The estimation of z˜ is the most difficult computation required in our class of methods. It is done only oncein each outer loop iteration, and typically only few iterations are needed at all. Recall from Section 4.1 thatfor a gradient descent method, z˜ would have to be estimated in every single step. In Section 5.2, we compareour variational inference methods to certain MAP estimation algorithms on the same model, showing that on apurely computational level, they mainly differ in the computation of z˜ (required for inference, but not for MAPestimation). It goes without saying that if Gaussian models with precision matrix A (under arbitrary positive pi)admit a better specific estimator for (VarQ[si|y])i, this should be preferred over our generic Lanczos solution here.Or the Lanczos method could be run with preconditioning, based on system knowledge. For example, if A hasa sparse graphical model structure, Gaussian belief propagation on modified graphs is used to precondition LCG(Malioutov et al., 2006b), and these ideas could be used to precondition the Lanczos algorithm as well.Surprisingly, in our experiments of main interest here, the bad relative accuracy of our estimator of z˜ does notseem to have a major impact at all. Although the spectral decay of A is linear, and the variance estimates havesignificant errors, the L2 reconstruction errors are often slightly better for smaller numbers of Lanczos iterationsk than for exact computations. Here is an idea why this might be the case. Importantly, z˜(k) approaches z˜componentwise from below, so we generally use underestimates. Moreover, the inner loop sparsity penalty forLaplace sites is 2∑i τi(zi + s2i /σ2)1/2 (see Section 4.2; recall that zi = z˜i in this case), which is stronger forsmaller zi. An underestimate of zi leads to a stonger sparsity penalty on si in the subsequent inner loop, andthis amplification happens mostly on those si, whose true zi values would be moderately small. More sparsity isimplied for s than what is specified by the prior. While this effect could be self-amplifying, it seems to be benignor even slightly beneficial in our applications, maybe because the Laplace sparsity potentials are not strong enoughin the first place.One may not be this lucky in other applications, or with other models. Indeed, it would be preferable to workwith proper estimators, and adjust the model and its potentials for the right degree of sparsity. Certainly, findingbetter Gaussian variance estimators is an important point for future work.165 ExtensionsIn this section, we collect a number of extensions of the approximate inference algorithms described in Section 4.Moreover, we discuss the precise relationship to closely related sparse MAP estimation methods, pinning downthe added complexity of variational inference. In this context, we also comment on compressed sensing for naturalimages, and on the misguided advice to measure natural or medical images with randomly drawn designs.5.1 Direct GeneralizationsIt should be obvious that our framework is not limited to sparse linear models with factorizing Gaussian likelihood.It can be applied whenever the posterior can be written as the normalized product of univariate sites, whose argu-ments are linear combinations of u, and if each of the sites has Gaussian-form lower bounds. For example, Laplaceor Student’s t sites can be used as likelihood terms for robust regression (Tipping and Lawrence, 2005). Binaryclassification with Bernoulli sites is discussed in Section 2.1. We can also deal with multivariate non-Gaussiansites, as long as corresponding multivariate lower bounds are given.17 Π and Γ become block-diagonal in sucha setting. All convexity properties discussed in Section 3.1, as well as Theorem 3, remain valid if positivity re-quirements on scalars are replaced by requirements of positive definiteness on square blocks. The sites can haveoverlapping supports, since our methods can certainly deal with matricesB that do not have full rank.We can also in principle accommodate fully coupled sites, as long as they are Gaussian. For example, thelikelihood P (y|u) may come with a general covariance matrixΠ(0)−1. This case is accommodated by replacing16 However, in the absence of such, or of a different approximate inference method as scalable as ours, the only alternativeseems to not do Bayesian experimental design on a large scale at all. Fortunately, this narrow “fundamentalist” view onanything vaguely Bayesian (often ignoring real-world aspects such as running time or user-friendlyness) does not hamper workin machine learning.17 For example, Fenchel bounds may be generalized to multivariate scale mixtures. The latter are useful to specify correla-tions in non-Gaussian priors. In the context of natural images, certain filter responses in s are known to be typically correlated,and multivariate scale mixture priors have been used in this context (Portilla et al., 2003).15σ−2‖y−Xu‖2 above by (y−Xu)TΠ(0)(y−Xu), and σ−2XTX byXTΠ(0)X in the system matrices. In thiscase, MVMs withΠ(0) (the inverse covariance matrix) have to be computed inside LCG and Lanczos iterations.In this paper, our interest in scalable approximate Bayesian inference is mainly driven by experimental designapplications to improve linear measurement architectures for images (see Section 6). Bayesian inference has manyother applications, and while we do not pursue any of them in detail here, we close this section with some remarksabout how our algorithms could be employed. First, the partition function P (y) of Section 2.1 is an importantconcept on its own, besides being a favourable target for variational relaxations. It is the marginal likelihood ofthe data, where the unknown parameters u have been integrated out. Bayes factors (Kaas and Raftery, 1995)are differences of log partition functions for different models of the same data, and they are routinely used formodel selection (they are the Bayesian equivalent to likelihood ratio statistics). The variational relaxation weemploy results in a bound on the log partition function, and our algorithms can be used to evaluate these boundsrapidly.18 Moreover, a powerful method for adjusting free hyperparameters, such as the τi scale parameters in (2),consists of maximizing the marginal likelihood P (y) of the data. This is implemented easily within our variationalframework, by just maximizing the lower bound to logP (y) instead. Technically, it is equivalent to maximumlikelihood learning of parameters in an undirected graphical model, and the derivatives of the lower bound w.r.t.hyperparameters have the usual form of simple posterior expectations. The maximization of our bound w.r.t. thenoise variance σ2 is a convex problem if the relaxation is convex, as will be shown elsewhere. However, in general,hyperparameter optimization by marginal likelihood techniques is a non-convex problem.5.2 Relationship to MAP Estimation. Added Complexity of Variational InferenceOur variational approximate inference methods are, from a purely computational viewpoint, closely related toMAP estimation algorithms for the same underlying posterior. A related point is made in (Wipf and Nagarajan,2008), where they compare convex MAP with non-convex sparsity estimators. In this section, we will see that itis precisely the step from MAP estimation to variational inference which makes the outer loop updates hard. Theadded complexity of variational inference versus MAP estimation can be quantified precisely in this case. Withthis analogy in mind, we can give weight to the core messages of this paper, towards the end of this section.The problem of maximum a posteriori (MAP) estimation for the sparse linear model ismaxuN(y|Xu, σ2I)∏iti(si).It is convex if the ti(si) are log-concave functions. Given that the ti(si) are super-Gaussian, as we assume in therest of this section, Theorem 3 states that the variational inference approximation we employ, is convex as well. Infact, it is closely related to particular MAP estimation techniques. A key step in Section 3 is (9), where∫(. . . ) duis replaced by |2piσ2Σ|1/2maxu(. . . ), an equality for Gaussian integrals. This implies that the MAP problem canbe written in much the same form as the approximate inference problem, only that the log |A| term vanishes. Forthe general case, the MAP estimation problem is equivalent tominuminγh(γ) + σ−2‖y −Xu‖2 + σ−2sTΓ−1s, h(γ) = −2q∑i=1log fi(1/γi).By Theorem 3, h(γ) is convex, therefore the whole criterion is jointly convex in (u,γ), which allows us tointerchange the ordering of minu and minγ , and we can solve the MAP problem by iteratively updating u and γ.This algorithm has been proposed in (Figueiredo, 2003) for the case of Laplace sites. Our Theorem 3 implies thatthe same algorithmic recipe can be applied to other models as well, for example (sparse) logistic regression.19In general, for sites that may not be log-concave, h(γ) = h∩(γ) + h∪(γ). The concave part h∩(γ) is upperbounded by Fenchel’s inequality, whence we repeatedly need to solve inner loop convex problems of the formminuminγzTγ + h∪(γ) + σ−2‖y −Xu‖2 + σ−2sTΠs.18 Some caution is necessary here, due to the problems noted in Section 4.5.19 These MAP estimation algorithms may not be globally convergent for sparsity-enforcing sites, at least no proof has beengiven in (Figueiredo, 2003). The problem is that many γi have to become exactly zero eventually, since the MAP estimator forsparse linear models is provably sparse. However, if the algorithm attends some u such that si = 0 for s = Bu, then γi andsi remain clamped at zero ever after. A global convergence proof for these MAP estimation methods is complete only if it isshown that only such si become zero, which are in fact zero at the true solution (assuming, of course, that all si 6= 0 initially).This problem does not occur for variational inference applications, since no γi can ever become zero there.16These inner loop steps are exactly the same as for the approximate inference method discussed in Section 4.1.However, the outer loop updates of z are different: zopt = ∇γh∩(γ). Since h∩ typically decouples, these updatesare simple. The log |A| term, coupling the sites, is not present in the MAP estimation context. For example, forStudent’s t sites (3), we have that zopt = (ν + 1) ◦ pi.Our observations here are related to (Wipf and Nagarajan, 2008, Sect. 3). Recall from Section 4.4 that theirmethod is useful for sparse estimation, yet does not compute the MAP estimator for the underlying model withStudent’s t sites, but also features a log |A| coupling term. In their context, the coupling term helps to be veryaggressively sparse. Our goal is not sparsity feature selection, but rather a faithful approximation of posteriorcovariances. As shown below, explicit sparsity destroys such covariance information. Since the MAP estimatorfor the sparse linear model exhibits exact sparsity in general, as does the method of (Wipf and Nagarajan, 2008),neither is useful for our purposes. The added complexities of obtaining meaningful uncertainty estimates requiredfor Bayesian experimental design (see Section 6) are precisely the outer loop computations, which consist of bulkmarginal variance estimations. As shown in Section 4.5, this added complexity is considerable, at least in terms oftechnology required. While the algorithm of (Wipf and Nagarajan, 2008) features a log |A| term as well, its roleis to enhance the degree of sparsity (i.e. , the number of γi = 0) beyond convex MAP estimation. The high degreeof exact sparsity in γ helps them in turn to compute the marginal variances more efficiently (see Section 4.4).We have seen that if γ becomes exactly sparse to a high degree, it is easy to implement algorithms such as oursefficiently, because the coupling matricesA are effectively only as large as the number of non-zeros in γ. In sparseestimation, γ has many zeros at the final solution, and in most algorithms, γ is sparse from the beginning. In thissetting, single-site updating algorithms (see Section 2.1) can be applicable to large problems, since single updatesare cheap due to the sparsity, and the true sparse optimum may be found in a moderate number of steps. However,when the goal is approximate inference for estimating uncertainties such as posterior variances or covariances,exact sparsity of γ cannot be expected. By Theorem 1, (7.), the posterior variance estimate VarQ[si|y] is upperbounded by σ2γi. If γi = 0, the method asserts that there is no posterior variance in si at all. By setting γi = 0, siis fixed to zero exactly, with absolute posterior certainty. This means that also the estimated correlation betweensi and any other sj is zero. Since computational savings through exact sparsity can only be expected if most γi areset to zero, this means that in the corresponding posterior Q(u|y), there is no uncertainty about the large majorityof the si, and no correlations between these and the few coefficients that survive. Basically, most coefficients arejust eliminated. Q(u|y) exists on the hyperplane corresponding to the few surviving coefficients only. The verysensible and important question about how sure the method is in switching any of the si off, cannot be answered,not even a ranking among the eliminated coefficients can be extracted. In our opinion, it makes little sense toapproximate Bayesian inference with such drastic side conditions, which do not come from the model or the data,but are nothing but artefacts due to the overly sparse approximation technique.One of the main messages of our paper is that sparsity is an important statistical principle, not only for pointestimation, but also if data analysis with meaningful uncertainty estimates is the goal. However, in this case, thecommon practice of clamping many variables to zero exactly, which is successfully used in sparse estimation, to allof our knowledge cannot be used. The goal of faithful posterior approximations has to be reached without relyingon exact sparsity, and other large scale numerical techniques have to be used. The algorithms we use here, suchas Lanczos or LCG, do not require exact sparsity in the variables to scale to very large problems. They are fast,because structure in the model matricesX andB is exploited for efficient MVMs. This structure may be sparsity(in the matrices, not in γ or s), but does not have to be. For example, if an image is to be reconstructed from Fouriersamples, MVMs withX can make use of the Fast Fourier Transform and related signal processing code. At themoment, the dominating interest in sparsity in machine learning and statistics, both in theory and in practice, seemsto be to obtain ever sparser estimators, as close as possible to the most extreme case of all: L0 regularization.20In some settings, such as some applications of sparse estimation, this may be the real goal (although by far themost experiments in sparse estimation papers we have seen, work with artificially generated data). But in others,certainly of at least equal importance in practice, it is not. Uncertainty estimates are basically destroyed by overlyaggressive sparsification. Beyond experimental design, uncertainty estimates are acknowledged to be important inoptimal decision theory. Decision-making systems based on aggressive sparse estimation may run fast, but alsorun the risk of getting it wrong with absolute confidence.20 The NP-hardness of this problem helps to understand some part of the apparent attractiveness of aggressive sparsification.But coming ever closer to a computationally hard problem with efficient relaxations does not automatically mean that real-worldproblems are solved better.17Compressed Sensing of Natural ImagesOur experimental design method is about compressed sensing of real-world signals, and we close this sectionwith a comment about compressed sensing and sparse estimation, which is in line with the core message of our pa-per. Natural images are approximately sparse in transform domains, such as wavelet or Fourier: the filter responsesfollow a power law, with a concentration close to zero (see Section 2). But for non-synthetic images, they are neverexactly zero, not even in single coefficients, and certainly the dominating coefficients are not distributed uniformlyat random. To see this, just transform an image, permute the wavelet coefficients randomly, and transform back:you will never retrieve anything like an image. In (Cande`s and Romberg, 2006, Sect. 2.2), a natural image is arti-ficially sparsified by setting small wavelet coefficients to zero, and this latter signal (which is not a natural image)is then reconstructed from random measurements. It goes without saying that such a pre-sparsifying oracle is notsomething you can buy for your camera or your MR scanner: it is not realizable, and examples like this cannot tellus much about how MAP reconstruction from random measurements performs on real natural images.Comprehensive studies on natural images can do that. The results in (Seeger and Nickisch, 2008) indicate thatrandom measurements perform poorly on natural images, and some theoretical arguments for why this might thecase, have been given in (Weiss et al., 2007). They show that for signals with a spectral decay (power as a functionof spatial frequency) exhibited by natural images, the signal-to-noise ratio (SNR) of random measurements tendsto zero rapidly. But the full story is likely to be even more interesting. In results to be reported elsewhere, wecompare different ways of sampling 2d Fourier coefficients, which have the same density of samples as a functionof distance from the Fourier space origin, so should give rise to the same SNR, according to the arguments of(Weiss et al., 2007). However, there is a large spread in reconstruction errors from these designs, and again,randomized designs work worst.The problem of optimizing designs for image measurement devices is of high importance in practice, in compu-tational photography, and even more so in medical imaging, or with cameras operating beyond the visible wave-lenghts. It is a fascinating one to study, owing to the complexity of natural images, and the constraints and errorsources coming with the devices. But to all of our present knowledge, it is not solved by uniformly randomizingthe measurement design. We hope that our work here, which allows to optimize measurement designs for fullimages, will help along the recognition that compressed sensing is a problem about real-world signals, not abouttruly sparse, unstructered random vectors, and that more work than uniform random sampling is needed in orderto solve it adequately.6 Bayesian Experimental DesignIn this paper, the main motivation for our scalable approximate inference method, which can maintain posteriorsover entire images, is that this allows us to optimize the image measurement designX through Bayesian sequentialdesign. This method is developed in the present section, where we also show how many large candidates for adesign extension can be scored efficiently.The framework we use here has been described in detail in (Seeger et al., 2007; Seeger, 2008), and its useful-ness for optimizing image measurement architectures has been demonstrated in (Seeger and Nickisch, 2008). Aclear outcome from the latter study was that while significant reductions in reconstruction error are realized byswitching from linear to sparse MAP reconstruction, it is the optimization of the measurement design specificallyfor sparse MAP estimation that allows for much larger gains. In fact, once good designs are used, the differencesin reconstruction errors between MAP and least squares reconstruction can be minor. In previous work, imagesof moderate size (such as 64 × 64) were dealt with, but the inference methods used there are not scalable. Ournovel variational algorithms can be used for significantly larger images, and for the model used in these references(SLM with Laplace prior sites) solves a convex problem. In this section, we demonstrate how the design scorecomputations can be scaled up accordingly, making use of the Lanczos algorithm once more.6.1 Sequential Design Score ComputationIn sequential experimental design21 (Chaudhuri and Mykland, 1993; Fedorov, 1972),X is built up through severalrounds. In each round, a set of candidates X∗ ∈ Rd×n of equal size is scored, and the winner (maximizing the21 We use the term “experimental design” in a narrow sense, compared to what readers from statistics may be familiar with.We wish to quantify amounts of information in parts of experiments, and to exploit this inference to optimize the measurementsautomatically, where the sole aim is to obtain faithful reconstructions faster or at a lower cost. Classical ED (in our sense)concentrates on the Fisher information matrix of an estimator. For sparse MAP estimation, this matrix is not well defined, dueto the shrinkage-to-zero properties (see also Section 5.2). Bayesian ED is not plagued by these problems.18score) is appended to X. A candidate scores highly if its measurements are deemed to reveal as much novelinformation about u as possible, given what is already known at the beginning of the round. In the following, wediscuss the scoring for a single round, starting from (X,y) and P (u|y). We employ the entropy difference scorehere,H[P (u|y)]− EP (y∗|y) [H[P (u|y,y∗)]] ,which quantifies the reduction in posterior uncertainty (Seeger, 2008). Note that reduction in uncertainty is scoredglobally over all of u. Importantly, the posterior correlations are fully contained in this score, setting it apart fromscores based purely on marginals of u. As discussed in (Seeger, 2008), its maximizer depends on the full posteriorcovariance matrix.If Q(u|y) = N(h, σ2Σ) is the approximation (8) to the current posterior P (u|y), we approximate this scoreby∆(X∗) := log |Σ|+ log∣∣Σ−1 +XT∗ X∗∣∣ = log ∣∣I +X∗ΣXT∗ ∣∣ , (12)the entropy difference betweenQ(u|y) and an updated Gaussian, whereX∗ is appended toX, and y∗ ∼ Q(y∗|y)to y. This is an approximation, because we do not adapt the variational parameters pi of Q(u|y) to the new data,but keep them at their old values. This approximation is mainly done for efficiency reasons, since our aim is toscore many candidates in each round. Note that no integration over y∗ is required for the ∆(X∗) computation.Details are given in (Seeger, 2008).Suppose there are N candidates of d rows each. The computation of (12) for these requires the solution ofN d linear systems with A, but different right hand sides, which is not feasible to do with LCG in applicationsof our interest. We came across a related problem in Section 4.5 already, and our approach once more involvesthe Lanczos algorithm. Recalling the notation there, if T (k) = L(k)L(k)T (Cholesky decomposition; L(k) isbidiagonal) and Q˜(k):= Q(k)L(k)−T (this computation is O(nk) only), thenI +X∗ΣXT∗ ≈ I + V T∗ V∗, V∗ := Q˜(k)TXT∗ ,ifΣ is replaced by its Lanczos low rank approximation. Finally, the score is computed using a Cholesky decompo-sition of this d× d matrix, or of I + V∗V T∗ (which has the same determinant) if k < d. In the applications we areinterested in, MVMs with largeX matrices are efficient, so the V∗ for the different candidates are best computeden bulk, given sufficient memory. On the other hand, these score computations can directly be parallelized, givenQ˜(k).The approximation described here manifests the role of the Lanczos matrices Q(k), T (k) of A as principalrepresentation of the approximate posterior itself. Once these are computed, they can be used to answer ratherarbitrary posterior queries, such as marginal variances of s = Bu (in Section 4.5) or entropy difference designscores.6.2 Other Design Setups. Relation to Classical DesignIn the previous section, we showed how to compute the design score (12) for many unrelated candidates. In thissection, we consider some other optimization settings for∆(X∗).First, consider d = 1, so ∆(x∗) = log(1 + xT∗Σx∗). As noted in (Seeger and Nickisch, 2008), the globalmaximizer of ∆(x∗) among all unit-norm vectors x∗ is the eigenvector corresponding to the largest eigenvalue ofΣ (Horn and Johnson, 1985, Sect. 4.2). This makes sense intuitively, since the uncertainty in xT∗ u is largest for thisdirection.22 It should be clear now that finding optimal extensions ofX with Bayesian experimental design needsan estimate of the whole posterior covariance matrix: its single node marginals are not enough. We can use theLanczos algorithm (see Section 4.5) in order to find the optimal direction. In fact, finding extremal eigenvectors istypically the main application of this algorithm. Recalling the discussion of convergence properties in Section 4.5,we note that if Lanczos is applied to A with a geometrically decaying spectrum, its minimal eigenvector (themaximal eigenvector of Σ = A−1) comes out last. In such a situation, packages like ARPACK run Lanczoson A−1 instead, which requires LCG in each iteration. However, for applications of our interest here, the linearspectral decay of A means that eigenvectors from both ends of the spectrum converge rapidly, so we can just runLanczos onA, until the minimal eigenvector converges.22 However, our simplifying assumption leading to (12) is also witnessed here, in that the outcome of a measurement x∗is assumed to mainly reduce the uncertainty in xT∗ u alone, while any further non-Gaussian “spread” of uncertainty reduction(which could happen, if pi was updated for the score computation) is ignored.19It is also possible to compute ∂∆(X∗)/(∂α) if ∂X∗/(∂α) is known. If M = I + X∗ΣXT∗ and W∗ =Q˜(k)T(∂X∗/(∂α))T , then∂∆(X∗)/(∂α) = trM−1(W T∗ V∗ + VT∗ W∗).If k < d, we can work with I + V∗V T∗ instead.How difficult is the optimization of∆(X∗) in general, over infinitely many candidatesX∗? For the case d = 1,we have to maximize the convex quadratic xT∗Σx∗, which is easy if the feasible set of x∗ is a Euclidean ball orellipse: the solution is a generalized eigenproblem. However, in general, convex quadratic maximization is a hardproblem, even for linear constraints on x∗. This is not surprising, since finding the optimal design in a Gaussianlinear model is already a hard problem in principle, at least for large n and q.From a purely computational viewpoint, our optimization of∆(X∗) is related to classical sequential D-optimaldesign, or its Bayesian analogue for Gaussian linear models (Chaloner and Verdinelli, 1995), with the importantdifference that the Fisher information matrix (which would be XTX for the Gaussian linear model) is replacedby the posterior precision matrix A here. Therefore, the work referenced in (Chaudhuri and Mykland, 1993) canbe used in our framework as well. In the nomenclature of (Chaloner and Verdinelli, 1995), our sparse Bayesiandesign framework is nonlinear. While they note that P (u|y) is often approximated by a Gaussian, none of themethods they refer to uses a modern variational approximation. With smooth non-Gaussian sites, the Laplaceapproximation is typically employed in statistics to obtain a Gaussian posterior approximation. However, in thecase of sparse linear models, the log posterior is strongly singular at its mode, so that the Laplace approximation isnot well defined. Moreover, the fact m n invalidates the typical justification for this approximation, as well asall asymptotic results about it we know of. Classical D-optimal design, as well as its Bayesian variant for Gaussianlinear models, are linear techniques, meaning that the design score to optimize does not depend on the observationsy. The design optimization is done without ever looking at any real data. In the applications of our interest, whilethe linearity of the measurements and the sparsity properties of the signal u can be motivated well, the modelsetup certainly does not perfectly represent the true data-generating process. The dependence of design decisionson data gathered along the process is an important feature of the our method, rendering it robust against modelmismatch, which is surely present. Nonlinear sequential design based on maximum likelihood estimation has beenanalyzed in (Chaudhuri and Mykland, 1993), and an interesting point for future research would be to extend theiranalysis to settings such as ours, where modern variational approximations of the posterior covariance matrix areused in place of the inverse Fisher information matrix, and MAP or posterior mean estimators replace maximumlikelihood techniques.7 DiscussionMany modern applications of statistical inference and estimation come with a large number of latent variables, of-ten many more than the number of independent datapoints. While a recent surge of activity has established efficientconvex methods for sparse point estimation from such data, little work has been done on higher-order problems forsparsity-favouring models, such as estimating confidences and dependencies, or optimizing measurement architec-tures. These problems can be addressed by Bayesian inference, but the commonly used standard approaches, suchas Laplace approximation or Markov chain Monte Carlo, either do not apply (the Hessian is strongly singular atthe mode; see also Section 5.2) or do not run fast enough for many real-world setups. Variational approximationsto Bayesian inference have been applied to sparse linear models (Tipping, 2001; Girolami, 2001; Figueiredo, 2003;Wipf et al., 2004), but the algorithms known so far scale up to large problems only if the main objective is sparseestimation once more. Moreover, the variational relaxations as well as the algorithms for solving them, have notbeen given satisfying characterizations.In this paper, we concentrate on a widely used general variational relaxation based on Gaussian-form lowerbounds (Jaakkola, 1997; Palmer et al., 2006; Girolami, 2001), which for a range of models is equivalent to thevariational (mean field) Bayes technique. We settle a long-standing question for this approximation, by showingthat the variational problem is a convex minimization if and only if the search for the posterior mode is. Moreover,and probably of more importance in practice, we provide the first truly scalable algorithms for solving this min-imization (whether convex or not) on large scale models. Scalability is achieved, just as in many convex sparseestimation algorithms, by reducing the dominating efforts to problems of standard form, with very well developedsolutions in numerical mathematics: least squares estimation, and variance estimation in Gaussian random fields.If code for these primitives is in place, the implementation of our methods is straightforward. Our setup is genericand can be configured with little effort to models featuring sparsity potentials (Laplace, Student’s t; see (Palmer20et al., 2006) for further super-Gaussian potentials), binary classification likelihoods, or many other exponentialfamily sites. Since the dominating computations are spent in standard primitives, no additional heuristics have tobe tuned. Moreover, structure in the model coupling matrices (measurement design, prior filters) can be exploitedvery effectively.Our main interest in this paper is Bayesian experimental design, with the aim of optimizing measurement ar-chitectures for natural images. Uniform random sampling is not enough to find useful measurement designs forreal-world signals (this point is discussed in more detail at the end of Section 5.2). The ability of estimating poste-rior covariances is crucial for finding good designs, and methods for aggressive sparse estimation cannot in generalbe used towards this end (they are important for reconstruction, once a good design has been found). Moreover,the queries required to improve a design can be approximated using the same primitives our inference algorithmsrely on.Our algorithms are of the double loop, or difference of convex type, which holds much promise for manyproblems in machine learning. The idea is to decouple a criterion φ, for which even the gradient is hard to compute,by upper bounding a critical part, so that the resulting bound φz can be minimized efficiently. This decouplingis successful if much of the criterion structure is still contained in the bound, so that only a few outer loop steps(re-fits of the bound) are required until convergence. This may not happen for all instances of our algorithms. Ifstrong couplings are eliminated by the bounding step, the inner loop optimization tends to converge rapidly withoutmuch progress in φ, and many outer loop steps are required. Since double loop algorithms are increasingly usedin machine learning and statistics (EM algorithm, CCCP, variational Bayes), it is an important point for furtherresearch to understand under which conditions they work well, and what can be done if they do not.23The algorithms proposed here have been shown to work well on challenging large scale medical imaging prob-lems, but are certainly useful in many other applications as well. In time series, filtering, or tracking problems,large state spaces could be endowed with sparsity priors, involving potentials between subsequent time points.Systems biology applications create datasets with very many latent variables, and experimental design can be usedthere to save on expensive experiments (Steinke et al., 2007). Sparsity models can also be used to analyze datafrom neural cell recordings, allowing for sharper predictions than traditional second order correlation analysis(Gerwinn et al., 2008). Our methods should be especially interesting for low-level computer vision applications,since sparsity-enforcing models represent natural image statistics much better than purely Gaussian ones, yet ourmethods draw exclusively on primitives such as mean and variance estimation in GaussianMarkov random fields,for which recently very efficient algorithms have been proposed (Malioutov et al., 2006a,b). In this sense, ouralgorithms reduce inference for certain non-Gaussian MRFs to (repeated) computations in Gaussian MRFs. TheSLMs we employ in this paper, are already used for problems such as image denoising or super-resolution, byway of sparse estimation. Beyond these applications close to machine learning, our techniques can be importedrather easily into application fields, where least squares estimation or lower-level signal processing techniques aredominatingly used, simply because our methods are configured solely in terms of such well-studied computationalprimitives.Finally, we hope that our work here contributes to the wider recognition of scalable Bayesian approximationsin computational application fields. It should become clear from our results that maximum a posteriori (MAP)estimation is not Bayesian inference (see Section 5.2 for a discussion), although even in machine learning or com-puter vision, Bayesian techniques are often equated with MAP throughout. Bayesian inference is about integration,while MAP is about optimization towards a single point estimate. While MAP estimates of the unknowns them-selves can often be computed by convex programming, covariances between variables or other uncertainty queriescannot properly be estimated this way. Variational approximations reduce Bayesian integration to optimizationproblems, and although by far the most work in this field is done on discrete graphical models, similar principlesapply to continuous variable models as well. In fact, our work shows that resulting algorithms can fit in seemlessly,running on the same computational primitives than MAP or least squares estimation, and thereby deal with long-range couplings and strongly non-local measurements in a way that is not currently possible for discrete randomfields methods.23 Some hints may come from loopy belief propagation for discrete graphical models, where double loop methods like CCCP(Yuille and Rangarajan, 2003) have been proposed, yet are much too slow to be practical.21AcknowledgmentsWe thank Rolf Pohmann, Bernhard Scho¨lkopf, Florian Steinke, and David Wipf for helpful discussions andinput. Supported in part by the IST Programme of the European Community, under the PASCAL Network ofExcellence, IST-2002-506778.A Further DetailsIn this section, we collect details of arguments which have been ommitted in the text.A.1 Proof of Theorem 1In this section, we provide the proof of Theorem 1. We begin with (5.). Define A = XTX + BT ρ(Γ)B fornow, where ρi(γi) > 0 for all i, and ψ1 = log |A|. We have that dψ1 = trSD with S = BA−1BT andD = ρ′(Γ)(dΓ). Now, since dA−1 = −A−1(dA)A−1, we have thatd2ψ1 = − trSDSD + trSρ′′(Γ)(dΓ)2 = trDSD(g(Γ)− S),where gi(γi) = ρ′′i (γi)/(ρ′i(γi))2. Since S 0 (positive semi-definite) we can write S = VV T with some matrixV, and d2ψ1 = tr(DV)T (g(Γ) − S)DV. If we can show that g(Γ) − S 0, then for γ(t) = γ + t(∆γ), wehave that ψ′′1 (0) = trMT (g(Γ)−S)M ≥ 0 for all small∆γ, whereM = ρ′(Γ)(∆Γ)V. This implies convexityof ψ1.Next, we show that ρ(Γ)−1 − S 0. Our proof employs the identityrTM−1r = maxx2rTx − xTMx, (13)which holds wheneverM 0 (positive definite). For any vector r ∈ Rq, we have thatrTBA−1BTr = maxx2rTBx − xT (XTX +BT ρ(Γ)B)x ≤ maxk=Bx2rTk − kT ρ(Γ)k,using (13) and xTXTXx = ‖Xx‖2 ≥ 0. Now, if the maximum is taken over all k ∈ Rq, the expression cannotbecome smaller, sorTSr ≤ maxk2rTk − kT ρ(Γ)k = rT ρ(Γ)−1r,using (13) once more. Therefore, ρ(Γ)−1 −S 0. Note that this argument, applied to ρi(γi) = γ−1i and r = δi,proves (7.).Collecting all parts, the convexity of γ 7→ ψ1 is implied by g(Γ) − ρ(Γ)−1 0. An elementary computationshows that the latter is implied by ρi(γi)ρ′′i (γi) ≥ (ρ′i(γi))2 for all γi. This completes the proof of (5.).We continue with (6.). Define A = XTX +BT ρ(Γ)−1B, and ψ2 = log |A| + log |ρ(Γ)|. The concavity ofγ 7→ ψ2 is shown by induction on q, the number of rows of B. Assume for now thatXTX is nonsingular. First,let q = 1, and b = BT . Then,log ρ1(γ1) + log∣∣XTX + ρ1(γ1)−1bbT ∣∣ = log ∣∣XTX∣∣+ log (ρ1(γ1) + bT (XTX)−1b) .Now, log(·) is concave and nondecreasing, so the concavity for q = 1 follows from (Boyd and Vandenberghe,2002, Sect. 3.2.4). If q > 1, letB = (BT 0,κ(ρ +∆δi) = bT(A 0. Forany γ 0 s.t. ψ2(γ) > −∞, ψε2 converges uniformly to ψ2 on a closed environment of γ (ψ2 and all ψε2 arecontinuous), so that ψ2 is concave around γ. This completes the proof of (6.).A.2 Proof of Theorem 3In this section, we prove Theorem 3. In fact, by Theorem 2, we only need to establish the convexity of h(γ) =−2∑i log fi(1/γi). Recall super-Gaussianity from Section 2.1. To simplify notation, we ignore some positivescaling, they do not alter convexity properties. We also ignore additional linear terms αisi in log ti(si), since theycan be dealt with as noted at the end of Section 2.1, leading to an additional linear term in φ(γ). In the following,we pick an index i ∈ {1, . . . , q}, and drop the corresponding subscript.Let x = s2 and g(s) = g(x) = log t(s). g(s) is odd, and we only deal with s ≥ 0 in the following. s 7→ g(s)is concave and twice continuously differentiable for s > 0, and x 7→ g(x) is strictly convex and nonincreasing forx > 0. By Fenchel duality,h(γ) = g∗(−1/γ) = sups≥0f = supx≥0f, f = −s2/γ − g(s) = −x/γ − g(x).We start with a simple, general observation. Let 0 < γ < γ′, so that −g(0) < h(γ), h(γ′) < ∞. If s∗ :=argmaxs f(s, γ) is a maximum point,24 then s∗ > 0, and h(γ′) ≥ −s2∗/γ′ − g(s∗) > −s2∗/γ − g(s∗) = h(γ).Therefore, if γ0 = sup{γ | f(s, γ) ≤ −g(0) ∀s} (γ0 = 0 if this set is empty), then s∗ = 0, h(γ) = −g(0) for0 < γ ≤ γ0, and for γ > γ0, h is strictly increasing, and s∗ > 0 (note that s∗ = ∞ is allowed). Therefore, itsuffices to show that h is convex at all γ where s∗ ∈ (0,∞).Here and in the following, g′ = dg/(ds), etc. We use the notation fs = ∂f/(∂s), functions are evaluated at(s∗, γ) if nothing else is said. Now, fs = −2s∗/γ − g′(s∗) = 0, so thatg′(s∗) = −2s∗/γ. (14)Next, x 7→ g(x) is twice continuously differentiable as well. We have x∗ = s2∗ at γ. fx is continuously dif-ferentiable, and g′′(x) > 0 by the strict convexity of g(x). By the implicit function theorem, x∗(·) is contin-uously differentiable at γ, and since h(γ) = f(x∗(γ), γ), h′(γ) exists. Moreover, 0 = (d/dγ)fx(x∗(γ), γ) =fx,γ + fx,x(dx∗)/(dγ), so that (dx∗)/(dγ) = γ−2/g′′(x∗) > 0, and x∗ is increasing in γ.From fs = 0, we have that h′(γ) = fγ = s2∗/γ2 = (g′(s∗))2/4 by (14). Now, g′(s) is nonincreasing by theconcavity of g(s), and g′(s∗) < 0 by (14), which means that s∗ 7→ h′(γ) is nondecreasing. Since s2∗ is increasingin γ, so is s∗. Therefore, γ 7→ h′(γ) is nondecreasing, thus h(γ) is convex for s∗ ∈ (0,∞).It remains to show necessity of the concavity of g(s). Suppose that g′′(s˜) > 0 for some s˜ > 0. From (14)we see that g′(s∗) < 0 whenever s∗ > 0, and by the same equation, there exists a γ˜ > 0 so that s∗(γ˜) = s˜.But if g′′(s∗) > 0 at γ˜, then s∗ 7→ h′(γ) is decreasing at s∗ = s˜, and just as above h′(γ) is decreasing at γ˜,so that h(γ) is not convex. Reverting to subscripts, suppose that hi is not convex at γ˜i > 0. Setting y = 0 in(10) leaves us with φ = h(γi) + log |A|, viewed as a function of γi. We can easily construct a setup so that∂2 log |A|/(∂γ2i ) < −h′′i (γ˜i)/2 at γi = γ˜i, showing that φ is not in general convex.A.3 Generic Inner Loop CriterionRecall the inner loop criterion (11), and its minimization by IRLS (see Section 4.3). In this section, we show howIRLS can be run if h∗∪,i(γi) cannot be determined analytically. In this case, its value and θi, ρi have to be computedon demand for each si encountered. In the following, we fix i ∈ {1, . . . , q} and drop the subscript.We only treat first class inner loops, the procedure for the second class is very similar. Let k(s, γ) = zγ +h∪(γ)+(s2/σ2)γ−1. Then h∗∪(γ) = k(s, γ∗), where the minimizer γ∗ is found by univariate convex minimization.Since kγ(s, γ∗) = 0, we have that θ = (σ2/2)(2/σ2)s/γ∗ = s/γ∗. Next, because kγ(s, γ∗) = 0 for all s, then0 = (d/ds)kγ(s, γ∗) = ks,γ + kγ,γ · (dγ∗)/(ds) (always evaluated at s, γ∗), leading todγ∗ds=sγ∗s2 + γ∗κ, κ =σ22γ2∗h′′∪(γ∗). (15)24 We do not require that s∗ is unique (f(s∗, γ) is unique, since x 7→ f(x, γ) is concave), but only that f(s∗, γ) > f(0, γ). Ifthe supremum is not attained at any point, then s∗ =∞ (by the continuity of f , this happens only ifmaxs∈[0,s0] f(s, γ) < h(γ)for all s0).23Finally, ρ = (dθ)/(ds) = (γ∗)−1[1− θ(dγ∗)/(ds)] = κ/(s2 + γ∗κ).Moreover, not even h∪(γ) needs to be known analytically. We show how to deal with the case h∪(γ) = h(γ),other cases being similar. All we need is code in order to compute g(x) and its first and second derivative. Recallthat g(x) is convex, and h(γ) = −minx≥0 l(x, γ), with l(x, γ) = x/γ + 2g(x). The procedure requires thecomputation of h(γ), h′(γ), and h′′(γ) on demand. First, h(γ) = −l(x∗, γ), where x∗ is found by univariateconvex minimization. Since lx = 0 (all functions evaluated at (x∗, γ)), we have that h′(γ) = −lγ = x∗/γ2.Moreover, (d/dγ)lx(x∗(γ), γ) = 0, so that (dx∗)/(dγ) = γ−2/(2g′′(x∗)) (recall that g′′(x∗) > 0). Therefore,h′′(γ) = γ−4((2g′′(x∗))−1 − 2x∗γ). The second derivative is needed in (15) only, we have thatκ = γ−1∗ σ2(14γ∗g′′(x∗)− x∗),where x∗ = x∗(γ∗).A.4 Two-Parameter Gaussian Site BoundsIf ti(si) is not even, the Gaussian-form lower bound may not be centered at zero, as happens in the case of Bernoullipotentials (7). More generally, the Gaussian form may come with a second parameter b˜i controlling its position.In this case, Q(u) = C−12 NU (s|σ−2b˜, σ−2Π). Our criterion φ remains as in (10), but sTΠs is replaced bysTΠs − 2b˜Ts. Let b := Π−1b˜ and s˜ := Bu − b. Then,φ(γ, b) = log |A| − σ−2bTΠb + h(γ, b) + σ−2minu(‖y −Xu‖2 + s˜TΠs˜) .Here, (u,γ, b) 7→ s˜TΠs˜ is jointly convex, giving rise to an inner loop much as in Section 4.3. Moreover,(γ, b) 7→ −bTΠb is jointly concave, so that either of our classes of algorithms can be generalized (using boundslinear in b and γ).B The Lanczos AlgorithmThe Lanczos algorithm (Golub and Van Loan, 1996) is a standard tool of numerical mathematics. For a linearsystem Ax = c, where A ∈ Rn×n is positive definite, the linear conjugate gradients (LCG) algorithm (Goluband Van Loan, 1996) is a method to approximate x∗ = A−1c, by minimizing the convex quadratic q(x) :=2cTx−xTAx. In each iteration, a single matrix-vector multiplication (MVM) withA is required. The sequenceof Krylov subspaces is defined as Kk := span{c,Ac, . . . ,Ak−1c}. The outcome of LCG after k iterations is theminimizer of q(x) within Kk. The Lanczos algorithm is an extension of LCG. In iteration k, a unit norm vectorq(k+1) is generated, which is orthogonal to all previous q(j), j ≤ k, so that Kk = span{q(1), . . . , q(k)}. In theLanczos sequence, q(k+1) is generated by a recurrence involving q(k) and q(k−1) only, and still only a single MVMwith A per iteration is required. On finite-precision computers, a number of additional techniques are required toobtain a useful algorithm (see Appendix B.1). In general, the quantities of interest to be estimated by the Lanczosmethod do not depend on the right hand side c. We follow the custom of using randomly drawn unit norm vectorsc.The Lanczos method constructs an orthonormal Q(k) := (q(1) . . . q(k)) and a tridiagonal T (k) with main diag-onal α ∈ Rk and subdiagonal β ∈ Rk−1, so thatQ(k)TAQ(k) = T (k). Here, leading eigenvalues ofA are closeto leading eigenvalues of T (k) rapidly, and this convergence of eigenvalues and eigenvectors can be tested withinthe algorithm itself (Parlett and Scott, 1979). The rough idea behind Lanczos estimates of linear functions ofA orΣ = A−1 is to plug in the following low rank approximationsA 7→ Q(k)T (k)Q(k)T , Σ 7→ Q(k)T (k)−1Q(k)T .If the function of interest depends on the spectrum of A only, the spectrum of T (k) is used instead. For example,the Lanczos estimate of z˜ = diag−1(BΣBT ) in Section 4.5 is given byz˜(k) := diag−1(BQ(k)T (k)−1Q(k)TBT).Since the T (k) are tridiagonal and nested, it is easy to derive a recurrence for these estimates. T (k) is positivedefinite just as A. Let T (k) = L(k)L(k)T be its Cholesky decomposition, where L(k) is lower triangular. The24factors are bidiagonal, say e := diag−1(L(k)), and d ∈ Rk−1 the subdiagonal. The following recurrence takes usfrom L(k−1) to L(k):dk−1 = βk−1/ek−1, ek =√αk − d2k−1.Here, d0 = 0. Let V (k) := BQ(k)L(k)−T . It is easy to see that the sequence V (k) is nested, so that V (k) =(v(1) . . .v(k)). Moreover, we have the recurrencev(k) = e−1k(Bq(k) − dk−1v(k−1)),which depends on the last recent v(k−1) only. Finally,z˜(k) = diag−1(V (k)V (k)T ) = z(k−1) + v(k) ◦ v(k).Since z(k)i = z(k−1)i + (v(k)i )2, the estimates z˜(k) are nondecreasing in each component, and converge towards z˜from below. Recall from Section 4.5 that this monotonicity property has important implications in practice. In ourimplementation, we also estimate log |A|, required in the computation of g∗(z), by log |T (k)| (the recurrence forL(k) leads to a recurrence for these estimates), although this is not required in order to run our algorithms.B.1 Lanczos ImplementationsA word of warning: the Lanczos algorithm is a powerful method to approximate spectral information of verylarge structured matrices, but the simple three-point recurrence it is based on, hides much of the complicationsin practice. Without a number of additional mechanisms, significantly more difficult than the recurrence, an im-plementation on a finite-precision computer fails almost surely. Orthogonality is rapidly lost among the Lanczosvectors, which leads to the method getting stuck in a subspace. Ironically, this degradation is a byproduct of therecurrence being so successful: it is caused by Ritz vectors (estimates of eigenvectors of A in Krylov subspacesKk) converging. The full story, one of the fascinating ones in numerical mathematics, is given in (Parlett and Scott,1979). But even the three-term recurrence itself is easily done wrong, as pointed out by Paige. We use the variantapproved in (Paige, 1976).We ignore preconditioning for now. Let A be the system matrix, and c be a unit norm starting vector, corre-sponding to the right hand side in linear conjugate gradients. The algorithm is given in Algorithm 1.Algorithm 1 Lanczos algorithmRequire: OperatorA. Initial c, ‖c‖ = 1q(1) = c. u(1) = Aq(1)for k = 1, 2, . . . doαk = q(k)Tu(k)Update estimatesr(k) = u(k) − αkq(k)Re-orthogonalize r(k)βk = ‖r(k)‖. Stop if too smallq(k+1) = r(k)/βku(k+1) = Aq(k+1) − βkq(k)end forFrom a practical viewpoint, the most important point in Algorithm 1 is the re-orthogonalization step: r(k) (tobecome the new q(k+1)) has to be orthogonalized (or deflated) against all previous q(j), j ≤ k. For a moderatenumber of iterations, this can be done naively, say by Gram-Schmidt, but we incur a cost of O(nk2) for thisstrategy, which for large k dominates the running time of the whole method. More advanced techniques have beenproposed, where r(k) is orthogonalized only against previously converged Ritz vectors (Parlett and Scott, 1979).However, they require more than the full matrix Q(k) to be stored in memory, and eigendecompositions of someof the T (k) have to be done. In most Lanczos applications, the spectrum ofA decays geometrically, and queries ofinterest depend on few leading eigenvectors only. Generic codes such as ARPACK (see Section 4.5; Matlab eigscalls this code) seem to be tailored towards such spectral behaviour. If the system matrices A in a variationalinference application are of this sort, these standard codes can be used. Unfortunately, as discussed in Section 4.5,25the spectral behaviour ofA in natural image reconstruction problems is different, and standard codes may be slowor even fail. Our strategy in such cases is described at the end of this section.Both LCG and Lanczos can typically be improved by preconditioning. The idea is to design some invertible C,so that systems withC andCT can be solved very rapidly (say, in O(n)), and at the same time, A˜ = C−1AC−Thas better properties as system matrix than A itself. For example, even the simple choice C = (diagA)1/2can lead to A˜ being better conditioned than A. If A is sparse, C can be computed by an incomplete Choleskydecomposition with small amount of fill-in. The Lanczos algorithm of Algorithm 1 can be run with a preconditionerC, by simply replacing A by A˜ and c by C−1c everywhere. The system matrices A in the application of ourinterest here are sufficiently well conditioned, and we do not use preconditioning at present. As noted in (SchneiderandWillsky, 2001), preconditioning could in principle be used to improve spectral properties of A˜ in comparison toA. For specialA coming from structured graphical models, such strategies have been proposed (see Section 4.5).However, in the image reconstruction applications of main interest here, such structure is not present. Research inpreconditioners for these setups is an important topic for future work.Lanczos with Lazy Selective OrthogonalizationThe dilemma of linear spectral decay ofA in applications of our interest here has been discussed in Section 4.5.The comments at the end of that section show that the resulting inaccuracy of Lanczos estimates do not invali-date reconstruction estimates or Bayesian design scores. However, other applications, or similar applications ondifferent data, may rely more strongly on accurate posterior variance estimates.We do not know of numerical mathematics work discussing this problem. Standard codes such as ARPACK justfail in these cases, or are very slow. Once more, the problem is not that leading eigenvectors of A are not found,but rather (ironically) that they converge rapidly. Time and memory requirements of modern Lanczos codes scalebadly with the number of converged Ritz vectors. But the variance estimates we require depend significantly on alarge part of the eigenspectrum ofA: they will be accurate only once many Ritz vectors have converged. Withoutadditional structural knowledge about A, there is probably little hope for principled improvement. But even forthe specific structure of main interest here (image reconstruction from non-local Fourier measurements and localfinite difference gradient potentials), which is heavily used in medical image reconstruction, we do not know ofany prior analysis of LCG or Lanczos that would be helpful here.Our present solution is to use Lanczos with complete (naive) re-orthogonalization, whenever the total run-ning time is not dominating by the deflations. This variant is easiest to code and requires the least amount ofmemory. We use Gram-Schmidt deflation.25 Our implementation also contains a modification of selective re-orthogonalization (Parlett and Scott, 1979), which is faster than the naive approach, but not drastically so (in theimage reconstruction applications). In the language of Parlett and Scott (1979), a direct implementation of theirscheme pauses in almost every iteration, and the running time is (perversely!) dominated by the eigendecompo-sitions of T (k) (required to test for converged Ritz vectors), and the re-computation of these vectors. Moreover,the strategies to avoid frequent pauses given there basically do not work in our case. A closer look reveals thattheir bounds are approximate themselves, relying on a favourable spectral decay ofA as well. Our most importantmodification is that at a pause, previously converged Ritz vectors are assigned to present ones, without having tore-compute them. This assignment strategy is based on the assumption that a converged Ritz vector will appear atall subsequent pauses just as well. We still have to recompute the eigendecomposition of T (k) at every pause.26Second, we do not use the heuristics of (Parlett and Scott, 1979), since they are not reliable in our case, but ratherallow for pauses only at every fourth (or so) iteration.ReferencesH. Attias. A variational Bayesian framework for graphical models. In S. Solla, T. Leen, and K.-R. Mu¨ller, editors,Advances in Neural Information Processing Systems 12, pages 209–215. MIT Press, 2000.S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2002.E. Cande`s and J. Romberg. Sparsity and incoherence in compressive sampling. Inverse Problems, 23:969–985,2006.25 Some experiments showed that Householder orthogonalization, advocated in (Golub and Van Loan, 1996), is not moreaccurate in our case, but runs significantly slower, presumably because fast BLAS primitives can be exploited to a lesser degree.26 Even though T (k) is tridiagonal, this needs an iterative computation. We have a very good initial guess about thisdecomposition from the last recent pause, which could be used (Gu and Eisenstat, 1994).26E. Cande`s, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal reconstruction from highly incom-plete frequency information. IEEE Transactions on Information Theory, 52(2):489–509, 2006.K. Chaloner and I. Verdinelli. Bayesian experimental design: A review. Statistical Science, 10(3):273–304, 1995.P. Chaudhuri and P. Mykland. Nonlinear experiments: Optimal design and inference based on likelihood. Journalof the American Statistical Association, 88(422):538–546, 1993.A. Dempster, N. Laird, and D. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journalof Roy. Stat. Soc. B, 39:1–38, 1977.D. Donoho. Compressed sensing. IEEE Transactions on Information Theory, 52(4):1289–1306, 2006.D. Donoho and M. Elad. Optimally sparse representation in general (nonorthogonal) dictionaries via l1 minimiza-tion. Proc. Natl. Acad. Sci. USA, 100:2197–2202, 2003.V. Fedorov. Theory of Optimal Experiments. Academic Press, 1972.M. Figueiredo. Adaptive sparseness for supervised learning. IEEE Transactions on Pattern Analysis and MachineIntelligence, 25(9):1050–1059, 2003.S. Gerwinn, J. Macke, M. Seeger, and M. Bethge. Bayesian inference for spiking neuron models with a sparsityprior. In Platt et al. (2008).Z. Ghahramani and M. Beal. Propagation algorithms for variational Bayesian learning. In Leen et al. (2001), pages507–513.M. Girolami. A variational method for learning sparse and overcomplete representations. Neural Computation,13:2517–2532, 2001.T. Gneiting. Normal scale mixtures and dual probability densities. J. Statist. Comput. Simul., 59:375–384, 1997.G. Golub and C. Van Loan. Matrix Computations. Johns Hopkins University Press, 3rd edition, 1996.M. Gu and S. Eisenstat. A stable and efficient algorithm for the rank-one modification of the symmetric eigen-problem. SIAM. J. Matrix Anal., 15(4):1266–1276, 1994.R. Horn and C. Johnson. Matrix Analysis. Cambridge University Press, 1st edition, 1985.T. Jaakkola. Variational Methods for Inference and Estimation in Graphical Models. PhD thesis, MassachusettsInstitute of Technology, 1997.M. Jordan, Z. Ghahramani, T. Jaakkola, and L. Saul. An introduction to variational methods in graphical models.In M. I. Jordan, editor, Learning in Graphical Models. Kluwer, 1997.R. E. Kaas and A. E. Raftery. Bayes factors and model uncertainty. Journal of the American Statistical Association,90:773–795, 1995.T. Leen, T. Dietterich, and V. Tresp, editors. Advances in Neural Information Processing Systems 13, 2001. MITPress.D. Malioutov, J. Johnson, and A. Willsky. Low-rank variance estimation in large-scale GMRF models. In Pro-ceedings of ICASSP, 2006a.D.Malioutov, J. Johnson, and A.Willsky. Walk-sums and belief propagation in Gaussian graphical models. Journalof Machine Learning Research, 7:2031–2064, 2006b.T. Minka. Expectation propagation for approximate Bayesian inference. In J. Breese and D. Koller, editors,Uncertainty in Artificial Intelligence 17. Morgan Kaufmann, 2001.M. Opper and O. Winther. Expectation consistent approximate inference. Journal of Machine Learning Research,6:2177–2204, 2005.C. Paige. Error analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix. IMA Journal of AppliedMathematics, 18(3):341–349, 1976.A. Palmer, D. Wipf, K. Kreutz-Delgado, and B. Rao. Variational EM algorithms for non-Gaussian latent variablemodels. In Y. Weiss, B. Scho¨lkopf, and J. Platt, editors, Advances in Neural Information Processing Systems 18.MIT Press, 2006.B. Parlett and D. Scott. The Lanczos algorithm with selective orthogonalization. Mathematics of Computation, 33(145):217–238, 1979.27J. Platt, D. Koller, Y. Singer, and S. Roweis, editors. Advances in Neural Information Processing Systems 20, 2008.MIT Press.J. Portilla, V. Strela, M. Wainwright, and E. Simoncelli. Image denoising using Gaussian scale mixtures in thewavelet domain. IEEE Transactions on Image Processing, 12:1338–1351, 2003.R. Rockafellar. Convex Analysis. Princeton University Press, 1970.M. Schneider and A. Willsky. Krylov subspace estimation. SIAM Journal on Scientific Computing, 22(5):1840–1864, 2001.M. Seeger. Bayesian inference and optimal design for the sparse linear model. Journal of Machine LearningResearch, 9:759–813, 2008.M. Seeger and H. Nickisch. Compressed sensing and Bayesian experimental design. In A. McCallum, S. Roweis,and R. Silva, editors, International Conference on Machine Learning 25. Omni Press, 2008.M. Seeger, F. Steinke, and K. Tsuda. Bayesian inference and optimal design in the sparse linear model. In M. Meilaand X. Shen, editors,Workshop on Artificial Intelligence and Statistics 11, 2007.M. Seeger, H. Nickisch, R. Pohmann, and B. Scho¨lkopf. Bayesian experimental design of magnetic resonanceimaging sequences. To appear at Neural Information Processing Systems 21, 2008.E. Simoncelli. Modeling the joint statistics of images in the Wavelet domain. In Proceedings 44th SPIE, pages188–195, 1999.F. Steinke, M. Seeger, and K. Tsuda. Experimental design for efficient identification of gene regulatory networksusing sparse Bayesian models. BMC Systems Biology, 1(51), 2007.R. Tibshirani. Regression shrinkage and selection via the Lasso. Journal of Roy. Stat. Soc. B, 58:267–288, 1996.M. Tipping. Sparse Bayesian learning and the relevance vector machine. Journal of Machine Learning Research,1:211–244, 2001.M. Tipping and A. Faul. Fast marginal likelihood maximisation for sparse Bayesian models. In C. Bishop andB. Frey, editors, Workshop on Artificial Intelligence and Statistics 9, 2003. Electronic Proceedings (ISBN 0-9727358-0-1).M. Tipping and N. Lawrence. Variational inference for Student-t models: Robust Bayesian interpolation andgeneralised component analysis. NeuroComputing, 69:123–141, 2005.M. Wainwright, E. Sudderth, and A. Willsky. Tree-based modeling and estimation of Gaussian processes on graphswith cycles. In Leen et al. (2001), pages 661–667.Y. Weiss, H. Chang, and W. Freeman. Learning compressed sensing. Snowbird Learning Workshop, Allerton, CA,2007.M. West. On scale mixtures of normal distributions. Biometrika, 74(3):646–648, 1987.D. Wipf and S. Nagarajan. A new view of automatic relevance determination. In Platt et al. (2008).D. Wipf, J. Palmer, and B. Rao. Perspectives on sparse Bayesian learning. In S. Thrun, L. Saul, and B. Scho¨lkopf,editors, Advances in Neural Information Processing Systems 16. MIT Press, 2004.A. Yuille and A. Rangarajan. The concave-convex procedure. Neural Computation, 15(4):915–936, 2003.28

MLPerf Inference BERT sparse model on SQuAD v1.1 dataset

Abstract

Extracted data

MLPerf Inference BERT sparse model on SQuAD v1.1 dataset

Abstract

Extracted data

Related items

Related items