supp [supp-]supp.tex

# TECHNICAL REPORT:
Sparse Bayesian Lasso via a Variable-Coefficient $\ell_{1}$ Penalty

Nathan Wycoff1, Ali Arab2, Katharine M. Donato3,4, Lisa O. Singh1,5
1 The McCourt School’s Massive Data Institute, Georgetown University
2 Department of Mathematics and Statistics, Georgetown University
3 The Institute for the Study of International Migration, Georgetown
University
4 The School of Foreign Service, Georgetown University
5 Department of Computer Science, Georgetown University

###### Abstract

Modern statistical learning algorithms are capable of amazing flexibility, but
struggle with interpretability. One possible solution is sparsity: making
inference such that many of the parameters are estimated as being identically
$0$, which may be imposed through the use of nonsmooth penalties such as the
$\ell_{1}$ penalty. However, the $\ell_{1}$ penalty introduces significant
bias when high sparsity is desired. In this article, we endow the $\ell_{1}$
penalty with learnable penalty weights $\lambda_{p}$. We start by
investigating the optimization problem this poses, developing a proximal
operator associated with the $\ell_{1}$ norm. We then study the theoretical
properties of this variable-coefficient $\ell_{1}$ penalty in the context of
penalized likelihood. Next, we investigate application of this penalty to
Variational Bayes, developing a model we call the Sparse Bayesian Lasso which
allows for behavior qualitatively like Lasso regression to be applied to
arbitrary variational models. In simulation studies, this gives us the
Uncertainty Quantification and low bias properties of simulation-based
approaches with an order of magnitude less computation. Finally, we apply our
methodology to a Bayesian lagged spatiotemporal regression model of internal
displacement that occurred during the Iraqi Civil War of 2013-2017.

Keywords: variational Bayesian inference, sparsity, proximal operators,
nonsmooth optimization, gravity model, human migration

##  1 Introduction

Basic Inferential Problem: Sparsifying penalties are popular because they
allow for simultaneous variable selection and parameter estimation by
shrinking parameter values identically to $0$. Explicitly, we define the
sparse learning problem as follows:

$\begin{split}\underset{\bm{\beta},\bm{\theta}}{\min}\,\,\mathcal{L}(\bm{\beta},\bm{\theta})+\tau\|\bm{\beta}\|_{0}\,\,,\end{split}$
(1)

where $\bm{\beta}$ is a parameter vector which we wish to be sparse,
$\bm{\theta}$ is a vector of other model parameters and $\mathcal{L}$ is a
misfit term which determines the compatibility of a parameter configuration
with a dataset (such as a negative log-likelihood or log posterior density).
This is a difficult problem because the $||.||_{0}$ function, which counts the
nonzero components of its input, is nonconvex and discontinuous. The
$||.||_{1}$ norm is a useful convex approximation to $||.||_{0}$. This is the
penalty function of the celebrated Lasso (Tibshirani, 1996; Taylor et al.,
1979), which may be viewed as placing a Laplace prior on $\beta_{p}$ (Park and
Casella, 2008). The draw-back to this convexity is bias: estimates of nonzero
parameters will be shrunk towards zero, significantly so if we are to achieve
much sparsity (Zhang, 2010). One solution to this problem is to allow the
penalty coefficients to vary by parameter:

$\underset{\bm{\beta},\bm{\theta}}{\min}\,\,\mathcal{L}(\bm{\beta},\bm{\theta})+\tau\sum_{p=1}^{P}\lambda_{p}|\beta_{p}|\,\,,$
(2)

as when $\lambda_{p}\to 0$ the associated $\beta_{p}$ is unpenalized. Far from
clear is how to choose these $\lambda_{p}$, as evidenced by the assortment of
approaches discussed in the next section.

Existing Approaches and Limitations: The importance of variable selection in
nonlinear models has prompted decades of research and yielded a multitude of
approaches. In this section, we review selected penalty-based approaches to
variable selection.

Some approaches perform adaptation “offline”, that is, before or after the
optimization process. The Adaptive Lasso (Zou, 2006) specifies
$\lambda_{p}=\frac{1}{|\hat{\beta}_{p}|^{\gamma}}$, where $\hat{\beta_{p}}$ is
some initial estimate of the regression coefficient and $\gamma$ is a
hyperparameter, though an initial estimate $\hat{\beta}_{p}$ is a nontrivial
ask if $\mathcal{L}$ is complicated. Bühlmann and Meier (2008); Candes et al.
(2008) propose to iterate this procedure, updating the penalty coefficient
$\lambda_{p}$ with new $\hat{\beta_{p}}$. Alternatively, we may use Lasso for
variable selection only then proceed to an unpenalized procedure (Efron et
al., 2004; Meinshausen, 2007; Zou and Li, 2008).

In the Bayesian framework, several authors have proposed simulation-based
procedures for adaptation of $\lambda_{p}$. Kang and Guo (2009); Leng et al.
(2014); Mallick and Yi (2014) place a Gamma prior or similar for $\lambda_{p}$
and conduct inference via Gibbs sampling. Bhattacharya et al. (2012) instead
take the approach of specifying a Dirichlet prior on the regression
coefficients before conducting MCMC. Conceptually this is similar to the
Horseshoe Prior (Carvalho et al., 2010), which instead specifies a conditional
Normal prior for $\beta_{p}$. The Horseshoe prior enjoys generality as well as
fast specific implementations (Terenin et al., 2019; Makalic and Schmidt,
2015). The reader is referred to Bhadra et al. (2019) for an extensive
comparison of the Lasso and Horseshoe models. Another approach is the Spike-
Slab prior (Mitchell and Beauchamp, 1988) which uses discrete latent variables
to categorize whether $\beta_{p}$ was a priori sampled from the spike or the
slab, which would complicate gradient-based inference.

Variational Bayesian inference (Blei et al., 2017, VB) is another way to adapt
the $\lambda_{p}$. But there is a sense in which VB with a Laplace prior
differs from Lasso: the integrated cost erases the sparsifying geometry, as we
discuss in Section 3.1. Previous work has circumvented this using either a
point mass variational distribution for regression coefficients (Tung et al.,
2019), which does not directly allow for uncertainty quantification (without,
say, bootstrapping Fu and Knight (2000)), or by thresholding small
coefficients (Babacan et al., 2014), which requires setting an arbitrary
threshold (She, 2009). Kawano et al. (2015) vary the threshold while
monitoring information criteria.

Nonconvex penalties do not vary the regularization strength but are directly
constructed to impose minimal bias on nonzero coefficients, such as the
Smoothly Clipped Absolute Deviation function (Fan and Li, 2001; Hunter and Li,
2005, SCAD) or Minimax Concave Penalty (Zhang, 2010, MCP). These penalties do
not correspond to proper priors, as they place constant, positive probability
density arbitrarily far from the origin.

Other penalties, like the ridge penalty (Hoerl and Kennard, 1970) and elastic
net (Zou and Hastie, 2005), do not serve primarily to select variables, but
rather to stabilize model fit or improve predictions. These are not the
subject of this article.

Somewhat in between Bayesian and penalized likelihood approaches lies Sparse
Bayesian Learning (Tipping, 2001). Here, sparsity comes not from
nonsmooothness of the likelihood with respect to $\beta$, but rather by
shrinking prior variance terms (and hence posterior variance terms) to zero
via empirical Bayes. As originally proposed, this requires a linear model and
Gaussian error structure, though its reach can be expanded using Gaussian
mixtures (Sandhu et al., 2021). Helgøy and Li (2019) applied this framework to
the Bayesian Lasso.

Scientific Motivation: Forced human migration is at an all time high and only
increasing. In order for policymakers to plan effectively, it is important to
understand which factors affect when and where people will move during times
of crisis. Traditionally, migration researchers use gravity models, described
in Section 4.1.2. These models were developed by economist to model
international trade (the flow of goods and services) as a function of
traditional econometric indicators, typically on a coarse spatiotemporal scale
(i.e. country by year). Migration researchers use this model because it can be
effective in capturing both push and pull dynamics. In this article, we study
an approximately year-long period of the Iraqi Civil War of 2013-2017, where
the government and its international backers fought ISIL (Islamic State, also
known as ISIS) in the northeast of the country. Traditional migration
variables are not immediately available in exactly those situations where they
are needed most: dangerous, remote, and underdeveloped areas, such as our case
study. Taking advantage of the proliferation of the mobile internet, we
develop novel social media “Buzz Variables” derived from Arabic-Language
twitter (Section 4.1.1). This allows us to make predictions at finer
spatiotemporal scales than traditional country-year level administrative data,
but introduces new challenges, namely unignorable lag between cause and
effect, spatiotemporal correlation, and high dimensional and noisy predictors.
In order to meet these new challenges, we develop a Bayesian hierarchical
model which allows for zero inflation, overdispersion, spatiotemporal random
effects, and lagging and aggregation of predictors (Section 4.2) – a complex
model to perform fixed-effects selection for.

Brief Outline of Our Contributions: We begin this article with a study of the
approach of treating the Laplace inverse scale parameter
$\bm{\lambda}:=\\{\lambda_{1},\ldots,\lambda_{P}\\}$ as an additional
parameter to be optimized, one endowed with the hyperprior $p_{\lambda}$,
yielding:

$\underset{\beta,\theta,\bm{\lambda}>\mathbf{0}}{\min}\,\,\mathcal{L}(\beta,\theta)+\sum_{p=1}^{P}\big{[}\tau\lambda_{p}|\beta_{p}|-\log\lambda_{p}\big{]}+\sum_{p=1}^{P}-\log
p_{\lambda}(\lambda_{p})\,\,,$ (3)

where $\tau$ is a positive parameter controlling the sparsity level and the
$\log\lambda_{p}$ term comes from the normalization constant of the Laplace
density. The $\ell_{1}$ penalty in general may be efficiently applied to a
smooth loss function via the Iterative Shrinkage and Thresholding Algorithm
(Daubechies et al., 2004, ISTA), a proximal gradient method Parikh et al.
(2014). But ISTA assumes a known and fixed $\lambda_{p}$. Section 2 examines a
proximal operator associated with the variable-$\lambda_{p}$ optimization
problem and discusses how to deploy it to marry adaptive and nonsmooth
penalties in an analog of ISTA we term Variable ISTA (VISTA) 111As we deploy
Nesterov acceleration, the better analog is actually Fast ISTA (Beck and
Teboulle, 2009).. It also investigates the basic theoretical properties of
this procedure in penalized likelihood.

In Section 3 we pivot to a Variational Bayesian take on an adaptive Lasso
using a nonsmooth penalty we term the Sparse Bayesian Lasso (SBL). The
“Bayesian Lasso” of Park and Casella (2008) is so named as it is a Bayesian
explanation for the density used in Lasso. But the Laplace prior does not
actually encode the prior beliefs that motivate the use of Lasso, i.e., that
the coefficient vector is sparse. Rather, it encodes the belief that most of
the coefficients are near zero, which leads to poor interval coverage and bias
if they are not. As such it is not operationally a Bayesian Lasso, that is, it
is not a method that yields the sparsity properties of Lasso as part of a
Bayesian analysis. The proposed Sparse Bayesian Lasso fills this role by
allowing for both geometric sparsity and uncertainty quantification. Combined
with warm starts222A “warm start” is simply the practice of initializing an
iterative algorithm at the final step of a previous, similar algorithm., we
can efficiently calculate trajectories not just of parameter estimates, as in
the traditional Lasso, but of entire variational distributions, showing us how
uncertainty of model parameters change with sparsity. This allows us to avoid
specifying or estimating an overall model complexity parameter $\tau$, which
is difficult in practice.

In this article, we explain how to deploy our method in gradient-based
learning frameworks such as tensorflow and pytorch. After verifying the
frequentist properties of the SBL in simulation studies, we demonstrate the
potential of the SBL in complex hierarchical models using our Iraq forced
displacement case study (Section 4). We conclude by discussing research
directions suggested by the novel proximal operator in Section 5.

##  2 Coefficient-Specific $\lambda_{p}$ with Optimization

We begin this section by developing a proximal gradient method suitable for
general smooth optimization problems augmented by the variable-coefficient
$\ell_{1}$ penalty. We then describe basic theoretical properties of the
estimator in the penalized likelihood case. We simply give outlines of proofs
for theoretical results; detailed derivations are available in the
Supplementary Material.

###  2.1 Background: Nonsmooth Penalties and Proximal Operators

We are interested in minimizing a complicated but smooth loss $l$ augmented
with a simple but nonsmooth regularizer $g$:

$c(\mathbf{x})=l(\mathbf{x})+g(\mathbf{x})\,,$ (4)

where $g$ is given by the $\ell_{1}$ norm in the case of Lasso regression.
Proximal gradient descent and its relatives are the algorithms of choice in
such a situation. Given a function $g$ with domain $\mathcal{X}$ and some norm
parameterized by a positive definite matrix $\mathbf{C}$, their proximal
operator is then defined as the following mapping:

$\mathrm{prox}_{g}^{\mathbf{C}}(\mathbf{x})=\underset{\mathbf{u}\in\mathcal{X}}{\mathrm{argmin}}\,g(\mathbf{u})+\frac{1}{2}||\mathbf{x}-\mathbf{u}||_{\mathbf{C}}^{2}\,\,.$
(5)

Intuitively, the proximal operator of a function $g$ evaluated at a vector
$\mathbf{x}$ returns another vector $\mathbf{u}$ which is close to
$\mathbf{x}$ (wrt $\mathbf{C}$) but does a better job minimizing $g$.

Proximal Gradient algorithms optimize the objective function in Equation 4 via
iteration of a two step process. Given a current solution $\mathbf{x}^{k}$ and
smooth cost function $\mathcal{L}$, the next iterate is parameterized by a
step size $s$ and given by:

$\displaystyle\hat{\mathbf{x}}^{k+1}=\mathbf{x}^{k}-s\mathbf{C}^{-1}\nabla\mathcal{L}(\mathbf{x}^{k})$
(6)
$\displaystyle\mathbf{x}^{k+1}=\mathrm{prox}^{s\mathbf{C}}_{g}(\hat{\mathbf{x}}^{k+1})\,.$
(7)

Note that $\mathbf{C}$ is scaled by the step size in the proximal operator and
that we have not included a subgradient of $g$ in the gradient descent step.
The matrix $\mathbf{C}$ may in theory be any positive definite matrix and in
practice it is defined by the preconditioning strategy of the gradient descent
algorithm that the proximal operator will be deployed in.

The proximal operator is most useful when it can be computed efficiently. In
this article, we will assume a diagonal preconditioner
$\mathbf{C}=\mathrm{diag}(c_{1},\ldots,c_{P})$, which is conducive to breaking
the proximal problem into subproblems defined along each axis. For example,
when $g:\mathbb{R}\to\mathbb{R}$ is given by $g(x)=\lambda|x|$, the proximal
operator is given by elementwise application of the Soft Thresholding Operator
(STO):

$[\mathrm{prox}^{\mathrm{sdiag(c_{1},\ldots,c_{P})}}_{\lambda|.|}(\mathbf{x})]_{p}=(|x_{p}|-sc_{p}\lambda)^{+}\mathrm{sgn}(x_{p})\,,$
(8)

where $(a)^{+}$ gives $\max(0,a)$ and $\mathrm{sgn}(a)$ gives the sign of $a$.

###  2.2 The Variable-Coefficient $\ell_{1}$ Proximal Operator

Figure 1: The Proximal Cost Left: The function $g(x,\lambda)=\lambda|x|$. The
proximal cost Center: of $(x,\lambda)$ Right: marginal for $\lambda$ with
$\lambda_{0}=x_{0}=1;s_{\lambda}=s_{x}=2$ yielding two optima.

Because the $\ell_{1}$ regularization coefficient $\lambda$ is now being
optimized over, the STO is no longer the pertinent proximal operator. Indeed,
the $\ell_{1}$ regularization function, considered formally as a function of
both $\lambda$ and $\mathbf{x}$, is nonconvex (see Figure 1, Left), somewhat
complicating proximal operator computation. In fact, many authors such as
Parikh et al. (2014) define the $\mathrm{prox}$ operator as one that acts on
convex functions. However, there has been work on extending proximality to
larger classes of functions (see e.g. (Hare and Sagastizábal, 2009)) and, of
particular interest to the statistical community, development of proximal
operators for the nonconvex “Bridge Penalties” $|\beta_{p}|^{q}$ for
$q\in(0,1)$ (Marjanovic and Solo, 2013); see Polson et al. (2015) for more on
proximal methods in statistics. Perhaps because of this focus on convexity
(and despite the popularity of the $\ell_{1}$ norm and adaptive penalty
methods), the proximal operator of $\lambda|x|$ as a
$\mathbb{R}^{+}\times\mathbb{R}\to\mathbb{R}^{+}$ function has not to our
knowledge been previously examined in the literature. It turns out that the
action $\mathrm{prox}(\mathbf{x},\bm{\lambda})$ of this proximal operator is
available in closed form and is single-valued for almost all inputs and always
for sufficiently small step sizes $s_{x}$ and $s_{\lambda}$ such that
$s_{x}s_{\lambda}<1$. We will assume in this section that $\tau=1$, since a
different $\tau$ simply scales the step sizes.

Consider the proximal operator of the variable-coefficient $\ell_{1}$ norm
function $g(\mathbf{x},\bm{\lambda})=\sum_{p=1}^{P}\lambda_{p}|x_{p}|$. Since
this function decomposes into additive functions of each $(\lambda_{p}$,
$x_{p})$ individually, its proximal operator acts on each block independently
of the others. Therefore, for the remainder of this section, we consider a
single block $(\lambda,x)$, dropping the index $p$, and consider the proximal
operator of the 2-dimensional function $g(x,\lambda)=\lambda|x|$:

$\displaystyle\mathrm{prox}^{s_{\lambda},s_{x}}_{g}(x_{0},\lambda_{0})=\underset{x\in\mathbb{R},\lambda>0}{\mathrm{argmin}}\,\,\lambda|x|+\frac{(x-x_{0})^{2}}{2s_{x}}+\frac{(\lambda-\lambda_{0})^{2}}{2s_{\lambda}}\,\,.$
(P1)

######  Lemma 1.

The marginal cost of P1 with respect to $\lambda$ (i.e. with $x$ profiled out)
is the following piecewise quadratic expression:

$\underset{\lambda>0}{\mathrm{argmin}}\begin{cases}\frac{1}{2}(\frac{1}{s_{\lambda}}-s_{x})\lambda^{2}+(|x_{0}|-\frac{\lambda_{0}}{s_{\lambda}})\lambda+\frac{\lambda_{0}^{2}}{2s_{\lambda}}&\lambda<\frac{|x_{0}|}{s_{x}}\\\
\frac{(\lambda-\lambda_{0})^{2}}{2s_{\lambda}}+\frac{x_{0}^{2}}{2s_{x}}&\lambda\geq\frac{|x_{0}|}{s_{x}}\,\,,\\\
\end{cases}$ (9)

where the changepoint $\lambda=\frac{|x_{0}|}{s_{x}}$ is the point where
$\lambda$ is just large enough to push $x$ to zero.

###### Proof.

Convert to nested optimization and exploit the known solution for fixed
$\lambda$ given by the soft thresholding operator. ∎

The quadratic polynomial in the interval $[\frac{|x_{0}|}{s_{x}},\infty)$ is
always convex. When $s_{\lambda}s_{x}<1$, the quadratic polynomial in the
other interval is convex as is the overall expression. But when
$s_{\lambda}s_{x}>1$, the coefficient of the quadratic term is negative, and
that polynomial is concave, yielding a nonconvex piecewise function (see
Figure 1, center and right).

We are now prepared to develop the proximal operator.

######  Theorem 1.

The optimizing $\lambda$ for the proximal program P1 is given by, when
$s_{x}s_{\lambda}<1$:

$\lambda^{*}=\begin{cases}\lambda_{0}&\lambda_{0}\geq\frac{|x_{0}|}{s_{x}}\\\
\frac{(\lambda_{0}-s_{\lambda}|x_{0}|)^{+}}{1-s_{\lambda}s_{x}}&o.w.\,\,\,\,,\end{cases}$
(10)

and, when $s_{x}s_{\lambda}\geq 1$, by
$\lambda^{*}=\mathbbm{1}_{\big{[}\frac{\lambda_{0}}{\sqrt{s_{\lambda}}}>\frac{|x_{0}|}{\sqrt{s_{x}}}\big{]}}\lambda_{0}$
(here $\mathbbm{1}$ denotes the indicator function). In either case
$x^{*}=(|x_{0}|-s_{x}\lambda^{*})^{+}\mathrm{sgn}(x_{0})$.

###### Proof.

We need only compare the optima of the quadratic functions of Lemma 1. ∎

Due to the nonconvexity of the proximal cost, this proximal program may have
two global optima. Thus the proximal operator is discontinuous and multi-
valued at the discontinuity, as visualized in the top right of Figure 2.

######  Remark 1.

When $s_{x}s_{\lambda}<0$ and $\lambda_{0}<s_{\lambda}|x_{0}|$ or when
$s_{x}s_{\lambda}>1$ and
$\frac{\lambda_{0}}{\sqrt{s_{\lambda}}}>\frac{|x_{0}|}{\sqrt{s_{x}}}$, the
solution to the proximal problem gives $\lambda=0$, which would lead to no
shrinkage on $\beta$.

Remark 1 is interesting, as it implies that it is possible to develop a
procedure with “dual sparsity”: on the regression coefficient, when
appropriate, or on the penalty coefficient. However, in our application of
this operator to the Laplace penalty, the $\lambda$ normalization term gives
this a density of zero, precluding that point being a penalized maximizer. We
look forward to examining other models which allow for zero penalties.

Fixed $s_{x}s_{\lambda}$:

Fixed $\frac{|x_{0}|}{s_{x}}$:

Fixed $\lambda_{0}$:

Figure 2:  The Action of the Proximal Operator: Plots of the reduced proximal
operator (Eq 11) considering two variables at a time and fixing the third.
Top: for various fixed $b:=s_{x}s_{\lambda}<1$ and with
$\lambda_{0},\frac{|x_{0}|}{s_{x}}\in(0,2)$. Values
$b=s_{x}s_{\lambda}\in\\{0.1,0.35,0.65,0.99\\}$ are shown left to right. Mid:
$a=\frac{|x_{0}|}{s_{x}}\in\\{0,0.5,1,1.8\\}$. Bottom:
$\lambda_{0}\in\\{0.2,0.75,1.25,1.75\\}$.

This proximal operator has been conceptualized as a mapping of
$(\lambda_{0},x_{0})\to(\lambda^{*},x^{*})$ parameterized by $s_{x}$ and
$s_{\lambda}$, but to visualize it we will briefly study it as a function of
these four quantities mapping to an optimizing $\lambda^{*}$.

######  Remark 2.

Since when $s_{x}s_{\lambda}>1$ the $\lambda^{*}$ is either $0$ or
$\lambda_{0}$, we will focus on the case where $s_{x}s_{\lambda}<1$. Then, let
$a:=\frac{|x_{0}|}{s_{x}}$ and $b:=s_{x}s_{\lambda}$ yielding a function of
just three variables:

$\lambda(\lambda_{0},a,b)=\begin{cases}\lambda_{0}&\lambda_{0}\geq a\\\
\frac{(\lambda_{0}-ab)^{+}}{1-b}&o.w.\,\,\,\,\,.\end{cases}$ (11)

This function is visualized in Figure 2. For given step size product $b$ the
function is stepwise linear, and converges to the identity mapping with
respect to its $\lambda_{0}$ input as $b\to 0$. As $b\to 1$, the mapping
becomes more and more steep for $\lambda_{0}\in(ab,a)$, gradually converging
to the discontinuous mapping
$\lambda(\lambda_{0},a,b)\to\lambda_{0}\mathbbm{1}_{\lambda_{0}>a}$ . The
proximal operator is multi-valued at $\lambda_{0}=a$ when $b=1$.

Comparison to Other Penalties in Statistics: Existing nonconvex alternatives
to our proposed biconvex penalty in the form of the MCP, SCAD and Bridge
penalties all come with an extra hyperparameter controlling how close the
penalty gets to approximating $||.||_{0}$ directly. The proposed approach
retains only the single global penalty strength parameter, which in practice
means it is possible to simply try all pertinent parameter values via a
sequence of warm starts. However, it does require specification of a prior on
$\lambda$, though we have found the standard Half-Cauchy to be generally
sufficient. Candes et al. (2008) and Fan et al. (2014) propose to set
$\lambda_{p}$ in order to locally approximate a prespecified nonconvex
penalty, and so also updates $\lambda_{p}$ each iteration.

###  2.3 Deploying the $\mathrm{prox}$ operator with VISTA

We can deploy the proximal operator of Equation 10 as part of a proximal
gradient method:

$\displaystyle\tilde{\bm{\beta}}^{t+1}\leftarrow\bm{\beta}^{t}-s_{\bm{\beta}}\nabla_{\bm{\beta}}\mathcal{L}(\bm{\beta}^{t},\bm{\lambda}^{t})$
(12)
$\displaystyle\tilde{\bm{\lambda}}^{t+1}\leftarrow\bm{\lambda}^{t}-s_{\bm{\lambda}}\nabla_{\bm{\lambda}}[\mathcal{L}(\bm{\beta}^{t},\bm{\lambda}^{t})-\log\bm{\lambda}^{t}-\log
P_{\bm{\lambda}}(\bm{\lambda}^{t})]$ (13)
$\displaystyle\bm{\beta}^{t+1},\bm{\lambda}^{t+1}\leftarrow\mathrm{prox}_{g}^{s_{\bm{\lambda}},s_{\bm{\beta}}}(\tilde{\bm{\beta}}^{t+1},\tilde{\bm{\lambda}}^{t+1})$
(14)

In practice, success with any gradient descent method relies on choosing good
step sizes, preconditioners, and acceleration, as we detail in Supplementary
2.

###  2.4 Basic Theoretical Properties for Penalized Likelihood

In this section, we will assume that $\bm{\lambda}$ is excluded from the
misfit term $\mathcal{L}(\bm{\beta})$ such that it only appears in the penalty
and, eventually, that this misfit term $\mathcal{L}(\bm{\beta})$ is given by a
negative likelihood $-L(\mathbf{y};\bm{\beta})$. In this case, we can rewrite
our penalty as such (we will assume all parameters are penalized for ease of
discussion):

$\displaystyle\underset{\bm{\beta},\bm{\lambda}>\mathbf{0}}{\min}\,\,-L(\bm{\beta})+\sum_{p=1}^{P}\big{[}\tau\lambda_{p}|\beta_{p}|-\log\lambda_{p}\big{]}+\sum_{p=1}^{P}-\log
p_{\lambda}(\lambda_{p})$ (15)
$\displaystyle\iff\underset{\bm{\beta}}{\min}\,\,-L(\bm{\beta})+\sum_{p=1}^{P}\underset{\lambda_{p}>0}{\min}\big{[}\tau\lambda_{p}|\beta_{p}|-\log\lambda_{p}-\log
p_{\lambda}(\lambda_{p})\big{]}$ (16)

Therefore, we may profile over $\lambda$ to develop a penalty
$g_{\tau}(|\beta|)=\underset{\lambda>\mathbf{0}}{\min}\,\tau\lambda|\beta|+\log\lambda+\rho(\lambda)$,
where $\rho(\lambda)=-\log p_{\lambda}(\lambda)$. We begin with some basic
properties of this penalty.

######  Lemma 2.

The following hold, where $\lambda^{*}$ denotes the optimizing $\lambda$, and
is formally a function of $\tau$ and $|\beta|$:

  1. 1.

$\lambda^{*}=\frac{1}{\tau|\beta|+\rho^{\prime}(\lambda^{*})}$.

  2. 2.

$\frac{\partial\lambda^{*}}{\partial|\beta|}=-\frac{\tau}{\frac{1}{\lambda^{*2}}+\rho^{\prime\prime}(\lambda)}$.

  3. 3.

$g_{\tau_{n}}^{\prime}(|\beta|)=\tau\lambda^{*}$.

  4. 4.

$g_{\tau_{n}}^{\prime\prime}(|\beta|)=-\frac{\tau_{n}^{2}}{(\tau_{n}+\rho^{\prime}(\lambda^{*}))^{2}+\rho^{\prime\prime}(\lambda^{*})}$.

###### Proof.

These follow from implicit differentiation on first order optimality
conditions. ∎

This allows us to quantify the behavior of this penalty as follows:

######  Theorem 2.

Assume that the logarithmic derivative of the hyperprior density on $\lambda$
is bounded ($|\rho^{\prime}(\lambda)|<M_{1}\,\,\forall\lambda\geq 0$) and that
the density is decreasing on $(0,\infty)$. Then:

  1. 1.

$g^{\prime}_{\tau}(|\beta|)\approx\frac{1}{|\beta|}$ for large $\beta$.

  2. 2.

The minimum of $|\beta|+g^{\prime}_{\tau}(|\beta|)$ is achieved at $\beta=0$
with value $\lambda_{a}\tau$.

###### Proof.

Follows from Lemma 2. ∎

######  Remark 3.

$g^{\prime}_{\tau}(|\beta|)\approx\frac{1}{|\beta|}$ approximates the gradient
of the adaptive lasso (Zou, 2006) procedure with optimal weights with
hyperparameter $\gamma=1$ as well as the iteration described by Candes et al.
(2008).

######  Remark 4.

Fan and Li (2001) describe three desirable properties of nonconcave penalties:
first, that they have bias decreasing quickly in nonzero parameter size,
second, that they induce sparsity, and third, that they be continuous in the
data. The shrinking gradient size with parameter norm is sufficient to ensure
unbiasedness for large parameters. The fact that the minimum of
$|\beta|+g^{\prime}_{\tau}(|\beta|)$ is strictly positive ensures sparsity,
while the fact that the minimum occurs at zero ensures continuity. This latter
condition is not satisfied by, for example, the bridge penalty.

######  Remark 5.

Bounded logarithmic derivatives are satisfied by the densities of, for
example, the Cauchy and Exponential distributions, but not the Gaussian
distribution.

We next consider the asymptotic distribution of the penalized likelihood
estimator. In particular, we demonstrate that there exists a local minimizer
of the penalized loss which satisfies the oracle property of Fan and Li
(2001). We assume that the first $r$ entries of the true parameter vector
$\bm{\beta}_{0}$ are nonzero, and the rest 0, such that
$\bm{\beta}_{0}={\tiny\begin{pmatrix}\bm{\beta}_{10}\\\
\bm{\beta}_{20}=\mathbf{0}\end{pmatrix}}$.

######  Theorem 3.

Let $\tau_{n}=n\tau_{0}$ for $\tau_{0}>0$, and further assume that
$|\rho^{\prime\prime}(|\lambda|)|<M_{2}$ (bounded second logarithmic
derivative). Then, under the standard regularity conditions on the likelihood
enumerated in the supplementary material, there is a local minimum of 16 that
satisfies the following:

  1. 1.

$\hat{\bm{\beta}}_{2}=\mathbf{0}$ with probability approaching 1 as
$n\to\infty$.

  2. 2.

$\hat{\bm{\beta}}_{1}$ is asymptotically normal with covariance given
approximately by $\frac{1}{n}I^{-1}(\bm{\beta}_{10})$, the Fisher information
matrix considering only active variables.

###### Proof.

See Supplementary 1\. ∎

The oracle property tells us that this estimator has the same asymptotic
distribution as that estimator with truly zero $\beta_{j}$ clamped to zero.

##  3 The Sparse Bayesian Lasso

Having developed a general-purpose optimization algorithm and developed basic
properties in the penalized likelihood case, we pivot to develping a Bayesian
Lasso procedure which provides full uncertainty quantification and penalty
coefficient adaptation. We achieve this by deploying the VISTA procedure of
the previous section on a novel nonsmooth Variational Bayesian methodology
which we call the Sparse Bayesian Lasso.

###  3.1 Variational Inference and Nonsmooth Penalties

Variational Bayes searches for a distribution over an unknown parameter vector
$\bm{\theta}$ that is 1) analytically tractable, and 2) sufficiently close to
the posterior such that what each has to say about the quantities of interest
are approximately the same. It does so by defining an optimization problem
$\underset{Q_{\theta}\in\mathcal{Q}}{\mathrm{argmin}}\,d(Q_{\theta},P_{\theta|\mathbf{y}})$
where $P_{\theta|\mathbf{y}}$ is our true posterior probability with density
$p$ and $Q_{\theta}$ is a candidate variational distribution with density $q$
from a space of possible distributions $\mathcal{Q}$ parameterized by
variational parameters, and $d$ is some measure of dissimilarity between
distributions. In this article, for each parameter $\theta$, the variational
distribution is defined by a location parameter $\eta_{\theta}$ and a scale
parameter $\nu_{\theta}$. We will impose variational independence as is
common:
$q_{\bm{\eta},\bm{\nu}}(\bm{\theta})=\prod_{m=1}^{M}q_{\eta_{m},\nu_{m}}(\theta_{m})$.

The most common $d$ is the KL divergence from variational to posterior:

$d(Q,P)=\mathrm{KL}(Q_{\theta}||P_{\theta|\mathbf{y}})=\mathbb{E}_{\bm{\theta}\sim
Q}[\log\bigg{(}\frac{q(\bm{\theta})}{p(\bm{\theta}|\mathbf{y})}\bigg{)}]\,=-\mathbb{E}_{\theta\sim
Q}[\log\mathcal{L}(\mathbf{y}|\bm{\theta})]+\mathrm{KL}(Q_{\theta}||P_{\theta})\,.$
(17)

We see this $d(Q_{\theta},P_{\theta|\mathbf{y}})$ may be viewed as the
negative expected log likelihood (which we can view as a model misfit term)
penalized by the KL divergence between the variational and prior distributions
(which we can view as a model complexity term). Since $Q$ was chosen to be
simple, we can use a Monte Carlo sample to estimate the expected likelihood.

For example, for regression coefficients $\beta_{p}$, we might specify Laplace
distributions both for variational and prior distributions:
$\beta_{p}\overset{Q}{\sim}\mathrm{L}(\eta_{\beta_{p}},\nu_{\beta_{p}})$ and
$\beta_{p}\overset{P}{\sim}\mathrm{L}(0,\frac{1}{\lambda_{p}\tau})$, leading
to the following penalty function on our variational parameters
$\eta_{\beta_{p}}$ and $\nu_{\beta_{p}}$:

$g_{\mathrm{KL}}(\eta_{\beta_{p}},\nu_{\beta_{p}},\lambda_{p})=\mathrm{KL}(\mathrm{L}(\eta_{\beta_{p}},\nu_{\beta_{p}})||\mathrm{L}(0,\frac{1}{\lambda\tau}))=\tau\lambda(\nu_{\beta_{p}}e^{-\frac{|\eta_{\beta_{p}}|}{\nu_{\beta_{p}}}}+|\eta_{\beta_{p}}|)-\log(\nu_{\beta_{p}})-\log\lambda\,.$
(18)

One might be forgiven for thinking that Variational Bayes, a procedure built
on optimization, would be the approach that allows us to combine the
generality afforded by Bayesian inference with the sparsity of penalized
likelihood procedures. Unfortunately, this is not the case. The proof of this
is straightforward, but we have not seen it explicitly mentioned in the
academic literature.

######  Theorem 4.

Penalty functions defined by KL-cost Variational Bayesian procedures with
continuous prior and variational distributions do not induce sparsity.

###### Proof.

As Fan and Li (2001, Section 2) discuss in the context of penalized
likelihood, any penalty function which induces sparsity and is continuous in
the data must have a singularity at the origin. Since Variational Bayesian
penalties are defined by a prior-variational KL divergence, and this in turn
is defined by an integral, an application of the fundamental theorem of
calculus is sufficient to reach the desired conclusion. ∎

The penalty given in Equation 18 is smooth despite the two absolute values
because they cancel at zero as visualized in Figure 3 and demonstrated
algebraically by Meyer (2021), who studied this function because of its
similarity to the Huber loss function.

###  3.2 A Nonsmooth Approximation

Owen (2007), in the context of penalized M estimation, advocates to
reintroduce nonsmoothness by exchanging the quadratic and linear parts of the
penalty. We propose to approximate
$\nu_{\beta_{p}}e^{-\frac{|\eta_{\beta_{p}}|}{\nu_{\beta_{p}}}}+|\eta_{\beta_{p}}|$
by $\nu_{\beta_{p}}+|\eta_{\beta_{p}}|$ since it is a good approximation for
small $\nu_{\beta_{p}}$:

$g_{\mathrm{NS}}(\eta_{\beta_{p}},\nu_{\beta_{p}},\lambda_{p}):=\tau\lambda_{p}(\nu_{\beta_{p}}+|\eta_{\beta_{p}}|)-\log(\nu_{\beta_{p}})-\log\lambda_{p}-1\,.$
(19)

Figure 3 show that $g_{\mathrm{NS}}$ retains the qualitative properties of the
KL penalty, and $g_{\mathrm{NS}}(\eta_{\beta_{p}},\nu_{\beta_{p}})-1$ bounds
$g_{\mathrm{KL}}(\eta_{\beta_{p}},\nu_{\beta_{p}})$ from below. Notice that
the term $\log\eta_{\beta_{p}}$ serves as a logarithmic barrier preventing
$\eta_{\beta_{p}}\to 0$. Therefore, the sparsity induced by this penalty is
qualitatively different from that induced by the Horseshoe, or Spike-Slab,
which compress the entire posterior distribution to zero. In contrast, this
proposed penalty sets only the variational mode exactly to zero while allowing
for nonzero, and sometimes significant, variational variance.

Figure 3: Nonsmooth Approximation Left: Bivariate KL penalty (blue) and
nonsmooth approximation (orange). The marginal penalties, ($g$ gives KL) for
Center: $\eta$, Right: $\nu$.

In the case of a linear regression with Gaussian error, the cost is closed
form whenever the first two moments of $Q_{\bm{\beta}}$ and the inverse and
log expectation of $Q_{\sigma^{2}}$ are tractable (Supplementary 3). We cannot
expect this to be the case in general, and must estimate the expected log
likelihood term $-\mathbb{E}_{\bm{\theta}\sim
Q_{\bm{\theta}}}[\log\mathcal{L}(\mathbf{y}|\bm{\theta})]$. Most common in
“blackbox” VB (Ranganath et al., 2014) is to optimize this with stochastic
optimization over a Monte Carlo sample. We instead use the Sample Average
Approximation (Robinson, 1996, SAA), which converts the optimization to a
deterministic one by employing the same Monte Carlo draw throughout the
optimization procedure, and use antithetic sampling for variance reduction
(Owen, 2013, Chapter 8.2). Collecting those parameters which are not desired
to be sparse into the vector $\bm{\theta}$ with associated variational
parameters $\eta_{\theta},\nu_{\theta}$, our cost is given as:

$\underset{\bm{\eta}}{\min}\,\,\frac{1}{B}\underset{\bm{\theta}_{b}\sim
Q_{\bm{\theta}},\bm{\beta}_{b}\sim
Q_{\bm{\beta}}}{\sum}\mathcal{L}(\bm{\theta}_{b},\bm{\beta}_{b})+\sum_{p=1}^{P}g_{\mathrm{NS}}(\eta_{\beta_{b,p}},\nu_{\beta_{b,p}},\lambda_{p})+\mathrm{KL}(Q_{\bm{\beta},\bm{\theta}}||P_{\bm{\beta},\bm{\theta}})$
(20)

Note that to compute gradients we must differentiate through the sampling
procedure; this is generally implemented in automatic differentiation packages
like tensorflow and pytorch, though quality of the gradients vary. When the
distribution is a location-scale family, this differentiation is particularly
straightforward (Kingma and Welling, 2014). With the gradient of the smooth
part computed, we can apply the VISTA algorithm to yield sparse learning with
any smooth likelihood. Since $\lambda_{p}$ is of secondary interest, we
estimate it via maximization rather than specifying a distribution for it.

###  3.3 Building Coefficient Trajectories by Varying $\tau$

Figure 4: Qualitative Behavior: Solid lines indicate coefficient estimates
(y-axis) as function of regularization strength for Lasso (left) SBL (right).
Dotted blue lines give true nonzero coefficient values (-2,2), and solid blue
lines give estimates. Orange lines give estimates of coefficients which are 0.
Dotted green line (second y-axis) gives percent of estimated coefficients
estimated as 0 (sparsity).

A practically useful aspect of Lasso regression is the ability to efficiently
form trajectories of parameter estimates while varying the penalty strength
$\tau$ by warm-starting optimizations. We also deploy VISTA estimating the SBL
with Gaussian error structure and a Half-Cauchy prior, illustrated in Figure
4, evaluate a trajectory for varying $\tau$, and then compare to the Lasso on
a simple simulated dataset with a true Gaussian-linear model in 50 dimensions.
Though it successfully thresholds the active parameters last, there is no
regularization strength which allows Lasso to identify the correct model
without bias. In the SBL, by contrast, the variational mean of the active
parameters are near their true values for the whole range of penalty values
until they are thresholded out, closer to the behavior of nonconvex penalties
such as MCP and SCAD. Unlike penalized likelihood procedures, however, we can
build trajectories of entire variational distributions, with credible
intervals pictured here.

###  3.4 Empirical Frequentist Properties

Likelihood | Method |  $\beta$-FNR |  $\beta$-FPR |  $\beta$-coverage |  $\sigma^{2}$-coverage | Time (s)
---|---|---|---|---|---|---
bernoulli | VISTA-MAP | 0.03 | 0.00 | - | - | 4.43
| VISTA-SBL | 0.00 | 0.00 | 0.91 | 1.00 | 31.58
| Horseshoe | 0.00 | 0.04 | 0.93 | 1.00 | 400.09
cauchy | SSLasso | 0.88 | 0.00 | - | - | 0.02
| VISTA-MAP | 0.30 | 0.00 | - | - | 6.43
| VISTA-SBL | 0.00 | 0.03 | 0.94 | 0.97 | 40.86
| Horseshoe | 0.00 | 0.00 | 0.94 | 0.87 | 458.67
nb | VISTA-MAP | 0.00 | 0.02 | - | - | 6.31
| VISTA-SBL | 0.00 | 0.00 | 0.95 | 1.00 | 28.77
| Horseshoe | 0.00 | 0.00 | 0.94 | 1.00 | 962.86
normal | SSLasso | 0.00 | 0.00 | - | - | 0.02
| VISTA-MAP | 0.00 | 0.02 | - | - | 7.12
| VISTA-SBL | 0.00 | 0.00 | 0.94 | 1.00 | 34.93
| Horseshoe | 0.00 | 0.00 | 0.95 | 1.00 | 39.10
poisson | VISTA-MAP | 0.00 | 0.00 | - | - | 18.01
| VISTA-SBL | 0.00 | 0.00 | 0.80 | 1.00 | 46.77
| Horseshoe | 0.00 | 0.00 | 0.98 | 1.00 | 394.34
Table 1: FNR is False Negative Rate, FPR is False Positive Rater, COVR is 95%
empirical interval coverage, time is elapsed real time in seconds for given
$\tau$.

In this section, we examine the empirical performance of our proposed method
on selected GLMs. We use the mean function described in our working example.
We deploy the VISTA algorithm, implemented in tensorflow_probability, on for
both MAP and VB inference with Bernoulli, Negative Binomial, Normal, and
Poisson likelihoods using the canonical link functions and with the Cauchy
likelihood with an identity link. We tuned sparsity level $\tau$ by hand which
is what we advocate doing in practice by examining the entire trajectory
produced by varying $\tau$ (Supplementary 5), fixing it across Monte Carlo
replicates.

We use a Monte Carlo sample size of 40. As a comparator, we show a JAGS
implementation of the Horseshoe prior, run for 1,000 iterations after
adaptation on two chains, which should represent a lower bound on the amount
of computation required for inference in practice. Additionally, we compare to
the Spike-Slab Lasso using the R package SSLASSO when the link function is the
identity, namely on the Gaussian and Cauchy cases. The Spike-Slab Lasso
penalty could of course be combined with a non-Gaussian likelihood, but
solving the marginal problem for coordinate descent will no longer in general
be in closed form and is not implemented in SSLASSO. In contrast, the proximal
operator is the same for any differentiable likelihood.

As Table 1 shows, SBL is able to maintain high interval coverage using
generally about an order of magnitude the elapsed real time of the Horseshoe
prior. The performance gap is smaller on the Normal likelihood where the Gibbs
sampler enjoys conjugacy. On the Negative Binomial likelihood, which is the
likelihood we use in the case study of Section 4.3, we observe about a
$30\times$ speedup. The SSLasso is fast but not yet applicable to non-Gaussian
likelihoods and unable provide UQ.

##  4 Modeling Migration in Iraq

We begin this section with an overview of the classical Gravity Model. We then
pivot to discussing the data set we collected, and finally present the results
of our analysis.

###  4.1 ISIL and the 2013-2017 Iraqi Civil War

In late 2013, tensions between Sunni and Shia Muslims in the Anbar province of
Iraq boiled over into violence between the Shia government and its allies
against the Islamic State. The violence led to significant internal
displacement which is the subject of our case study. In the absence of
traditional migration variables, we consider the use of social media for
capturing indirect indicators of migrations. Specifically, we focus on novel
conversation buzz and insecurity predictors based on Twitter data.

####  4.1.1 Data Collection and Processing

Figure 5: Iraq Dataset: Top Two Rows: Buzz Variables broken up by category and
normalized to have max value 1. Bottom Left Quadrant: The cumulative buzz by
each category is shown above the UN IOM data giving the families fleeing
during each period. Bottom Right Quadrant: Spatial totals for Buzz (upper row)
and flow (lower row); gray districts not observed.

Social Media Indicators We collected Arabic language tweets between March 2016
and June 2017 from Twitter’s Streaming Application Program Interface (API).
Buzz variables were formed by monitoring counts of keywords selected by Arabic
language experts related to various migration-relevant topics, such as ethnic
angst and discussion about the economy. We also formed sentiment indicators.
Variables such as Social buzz for different demographic groups, (e.g. Social-
shia_buzz), Insecurity, (e.g. Insecurity-emotions), Environmental (e.g.
Environmental-weather) are examples of the 68 variables identified.

Movement Data We use the flow data provided by the International Organization
for Migration’s Iraq Displacement Tracking Matrix, which gives stocks of
displaced families in one of 97 Districts, as well as which of the 18
Governorates they originated from. We differenced these data with respect to
time to estimate flows. The counts are aggregated over variable length
periods, with the shortest period consisting of 11 days and longest of 31
days. There are 48 total time periods in our study, leading to a total of
57,024 observations. 96.4% of response values are zero, with a mean of 6.71
families displaced for all time across origin-destination pairs.

####  4.1.2 New Challenges for the Gravity Model

Gravity models of migration (Poprawe, 2015; Karemera et al., 2000; Ramos and
Suriñach, 2017) measure the exchange between the origin $o$ and destination
$d$ for a given time period $t$, denoted as $y_{t,o,d}$; in our case study
$y_{t,o,d}$ is measured in number of families. The gravity model specifies a
generalized linear model relating $\mathbf{y}\in\mathbb{R}^{N_{o}N_{d}T}$ to
predictors $\mathbf{X}\in\mathbb{R}^{N_{o}N_{d}T\times P}$ as well as
resistance to exchange terms $\omega_{t}$, $\omega_{o}$, $\omega_{d}$ which
represent the effects associated with each time periods, origin and
destination. Poisson gravity models have been fit using MCMC for hospital
patient flows (Congdon, 2000) and trade (Ranjan and Tobias, 2007). Chen et al.
(2018) perform explicit Bayesian Model Selection on a gravity model of trade
via Normal approximations to their likelihood. The fine timescale of our data
introduces temporal misalignment between cause and effect: if an event which
causes migration is recorded by our Twitter Buzz variables on a given day, we
do not expect migrants to register at a camp the same day. At the yearly scale
which gravity models are typically deployed, this is not a concern, but we
need to develop a mechanism for estimating lags in this specific case study.

###  4.2 Spatiotemporal ZINB Gravity Model with Gaussian Lags

In this section, we develop the following hierarchical Bayesian model to
relate the organic variables to the observed migration:

$\displaystyle\lambda_{p}\overset{\mathrm{iid}}{\sim}\mathrm{C}(0,1)^{+}\,\,\mathrm{for}\,\,p\in\\{1,\ldots,P\\}$
(21)
$\displaystyle\beta_{p}|\lambda_{p}\overset{\mathrm{indep}}{\sim}\mathrm{L}(0,\frac{1}{\tau\lambda_{p}})^{+}\,\,\mathrm{for}\,\,p\in\\{1,\ldots,P\\}$
(22)
$\displaystyle\phi_{\mathcal{I}}\sim\mathrm{logitN}(a_{\phi},b_{\phi})\,\,\mathrm{for}\,\,\mathcal{I}\in\\{\mathcal{T},\mathcal{O},\mathcal{D}\\}$
(23)
$\displaystyle\rho^{2}_{\mathcal{I}}\sim\mathrm{logN}(a_{\rho},b_{\rho})\,\,\mathrm{for}\,\,\mathcal{I}\in\\{\mathcal{T},\mathcal{O},\mathcal{D}\\}$
(24)
$\displaystyle\bm{\omega}_{\mathcal{I}}|\phi_{\mathcal{I}},\rho^{2}_{\mathcal{I}}\sim\mathrm{N}\big{(}\mathbf{0},\rho^{2}_{\mathcal{I}}\Sigma_{\mathcal{I}}(\phi_{\mathcal{I}})\big{)}\,\,\mathrm{for}\,\,\mathcal{I}\in\\{\mathcal{T},\mathcal{O},\mathcal{D}\\}$
(25) $\displaystyle\pi\sim\mathrm{logitN}(a_{\pi},b_{\pi})$ (26)
$\displaystyle y_{t,o,d}\sim\begin{cases}\delta_{0}&\mathrm{w.p.}\,\,\,\pi\\\
\mathrm{N.B.}(g(\mathbf{z}_{t,o,d}(\mu_{l},\sigma_{l})^{\top}\beta+\omega_{t}+\omega_{o}+\omega_{d}),\frac{l_{t}}{\sigma^{2}})&\mathrm{o.w.}\end{cases}$
(27)

where $\phi_{\mathcal{T}}\in[0,1]$ is the correlation between successive
temporal effects $\bm{\omega}_{t}$ giving the AR1 correlation matrix
$\Sigma_{t}(\phi_{t})$, while $\phi_{o},\phi_{d}$ are the correlations for the
origin and destination spatial effects
$\bm{\omega}_{\mathcal{O}},\bm{\omega}_{\mathcal{D}}$, respectively, which lie
between the inverse maximum and minimum eigenvalues of the adjacency matrix
(De Oliveira, 2012). These matrices are scaled by variance terms
$\rho_{\mathcal{I}}$ endowed with LogNormal priors. The number of families
moving between two regions on a given day is given by a negative binomial
distribution parameterized by the canonical link transformation $g$ of our
linear predictions. The positive “number of failures” parameter controlling
overdispersion is estimated by $\frac{1}{\sigma^{2}}$ and is scaled by the
duration of the time period $l_{t}$. Due to the many of zeros, we include the
possibility of zero inflation via a parameter $\pi$ with logitNormal prior. We
account for unknown and variable lag between cause and effect through a
“flattened” Gaussian filter with parameter $\mu_{l}$ determining how far back
to shift the time series and the parameter $\sigma_{l}$ determining how far
outside the bounds of the period length to integrate and producing the
aggregated effects $\mathbf{z}$ (Supplementary 5). We use a neighborhood
structure that has weight 1 if two regions touch and 0 otherwise.

####  4.2.1 Variational Distribution and Inference

We optimize the divergence via the VISTA algorithm. We place logitNormal
variational distributions on $\phi_{t},\phi_{o},\phi_{d},\pi$ and logNormal
variationals on $\rho_{t},\rho_{o},\rho_{d},\sigma^{2}$. The KL divergence
between prior and variational in all of these cases is given by the KL
divergence of the underlying normal distributions and thus is available in
closed form. The $\omega$’s are given a Normal variational distribution. While
the KL divergence cannot completely be determined in closed form due to the
spatiotemporal correlation parameters $\bm{\phi}$, the $\omega$ and $\rho$
integrals may be computed analytically, and this expression may be evaluated
for Monte Carlo samples of $\phi$ without any matrix decompositions beyond an
initial decomposition of the neighborhood structure matrix (Supplementary 5).
As is done in the context of Gibbs sampling for CAR models, we “center on the
fly” (Ferreira, 2019) by projecting the variational mean of
$\bm{\omega}_{\mathcal{I}}$ to have mean 0 after each iteration:
$\bm{\eta}_{\omega_{\mathcal{I}}}\leftarrow\bm{\eta}_{\omega_{\mathcal{I}}}-\mathrm{mean}(\bm{\eta}_{\omega_{\mathcal{I}}})\mathbf{1}$.

###  4.3 Results

Figure 6: Parameter Trajectories vs $\tau$. Top Left: Regression Coefficient
Variational Trajectory; dotted lines give 95% credible intervals; intervals
are faded when the posterior mean is $0$. Top Right: Spatiotemporal
Correlations fairly constant. Bottom Left: Gaussian convolution parameters
vary with fixed effects. Bottom Left: Zero inflation parameter is constant.

In this section, we apply our methodology to our case study. We begin by
setting the regularization coefficient $\tau=250$, where the model chooses to
include only the destination-level variable insecurity-emotions-buz. We then
decreased this coefficient along a logarithmically spaced grid to $\tau=10$
with 100 points in between. Figure 6 shows the trajectory borne out by the
regression coefficients, random effect variances and lag and zero inflation
parameters. Begin with the top-left three panels, which show the variational
distributions of the regression coefficients. Qualitatively, we notice that
variances of zeroed out coefficients decrease smoothly as the regularization
strength increases, reflecting the increased prior distribution concentration.
The first variable to enter the model is destination insecurity-emotions-buzz,
followed by road-dist, the dyadic distance between origin and destination,
which actually seems to indicate that migration occurs over longer distances
(see Supplementary 5.4 for more discussion). Next relief-aid enters with a
positive coefficient, indicating either that good relief aid attracts or is
located where people are attracted, followed by environment-negative sentiment
at the destination level and environment-negative and environment-agriculture
at the origin level.

The lag-aggregation parameters are shown in the bottom left of Figure 6. As
the number of variables in the model increases, so too does the lookback
window size $\sigma_{l}$, from about $9$ days to about $15$, while the number
of days to look back travels from around $18$ to around $22$. This varies from
a model with sharper predictions from fewer variables to one with smoother
predictions incorporating more variables and time. The spatiotemporal
correlation parameters’ variational distributions are shown in the top 3 right
panels of Figure 6. The origin parameter is not well pinned down by our only
18 governorates, the Destination is also not strongly identified though it
does rule out substantial negative spatial correlation, and the temporal
correlation seems minimal. Finally, the bottom right panel indicates that the
zero-inflation parameter is also fairly constant with respect to $\tau$.

Figure 7: Spatiotemporal Random Effects: At the largest $\tau$ value.

The most positive origin random effects belongs to Ninewa Governorate with
capital Mosul, which experienced urban combat against ISIL, followed by Maysan
and Kirkuk governorates. For destinations, the Baghdad area has some of the
highest effects, with its constituent Al-Sadr and Kut districts well separated
from the rest. On the lower end we see a different story: many districts which
observed no inflow of migration have random effects which are experiencing
separation and similar, negative values. The time period with the largest
random effect was from 2016-05-12 to 2016-05-26, a period which saw a series
of ISIL bomb attacks which left over 100 dead in Baghdad (May 12) as well as
combat between Iraqi and ISIL forces (Wilson Center, 2019), while lowest
random effect was 2016-09-01 to 2016-09-15, during which ISIL was more active
abroad than in Iraq (ibid.).

##  5 Discussion and Future Work

We began this article by examining the density of the Bayesian Lasso when
$\lambda_{p}$ is given a hyperprior and treated as a variable to optimize.
This lead to a biconvex proximal operator that is well defined for
sufficiently small step sizes. We then showed how to plug this proximal
operator into ISTA, producing the VISTA algorithm which allows for minimal
bias even in the presence of high sparsity. Next, we examined the asymptotic
behavior of penalized likelihood estimates, finding that there exists a local
penalized likelihood maximizer with the oracle property. We then changed
gears, showing how to deploy the penalty in the context of Variational Bayes
and developing the Sparse Bayesian Lasso, which builds trajectories of entire
variational distributions. Simulation studies showed approximately nominal
coverage of its credible intervals. Finally, we deployed the VISTA algorithm
and Sparse Bayesian Lasso on a sophisticated hierarchical model and were able
to generate trajectories of variational distributions of all model parameters,
varying the penalty strength.

We look forward to deploying the VISTA algorithm on other applications such as
the inverse and imaging problems where ISTA made its name. And we are even
more excited to see what other possibilities viewing penalty coefficients as
optimization quantities can open up. In particular, Group Lasso (Yuan and Lin,
2006) and Fused Lasso (Tibshirani et al., 2005), among other Lasso analogs,
enjoy fast inference via proximal gradient methods as well; perhaps the
variable-penalty versions of these proximal operators are available in closed
form. Variable-coefficient low-rank penalties for matrices (e.g. Koltchinskii
et al. (2011)) and their associated proximal operators are also of future
interest. Recently, Quaini and Trojani (2022) developed basic inferential
properties of estimators defined via a proximal operator; it would be
interesting to see which of these properties apply to the biconvex proximal
operator we have developed here. We are also curious about different ways to
motivate nonsmoothness in VB, perhaps via different divergences or
discontinuous priors.

## Acknowledgements

The authors gratefully acknowledge funding from the Massive Data Institute and
McCourt Institute. We also acknowledge the Georgetown Data Lab and MDI
Technical Team for support creating migration indices, and Douglas Post for
his preliminary analysis of the Iraq data. We would like to thank Stephen
Becker, Warren Hare, Nicholas Polson, Mahlet Tadesse, and Stefan Wild for
valuable conversations and input. Any errors are our own.

SUPPLEMENTARY MATERIAL

Python Package:

PyPI-package sbl containing code to perform the the methods outlined in this
article.

Supplementary Discussion and Results

PDF Document containing additional numerical results, mathematical details,
and background.

Supplementary Material: Sparse Bayesian Lasso via a Variable-Coefficient
$\ell_{1}$ Penalty

##  1 Extended Proofs for Section 2

######  Lemma 1.

The marginal cost of P1 with respect to $\lambda$ (i.e. with $x$ profiled out)
is the following piecewise quadratic expression:

$\underset{\lambda>0}{\mathrm{argmin}}\begin{cases}\frac{1}{2}(\frac{1}{s_{\lambda}}-s_{x})\lambda^{2}+(|x_{0}|-\frac{\lambda_{0}}{s_{\lambda}})\lambda+\frac{\lambda_{0}^{2}}{2s_{\lambda}}&\lambda<\frac{|x_{0}|}{s_{x}}\\\
\frac{(\lambda-\lambda_{0})^{2}}{2s_{\lambda}}+\frac{x_{0}^{2}}{2s_{x}}&\lambda\geq\frac{|x_{0}|}{s_{x}}\,\,,\\\
\end{cases}$ (28)

where the changepoint $\lambda=\frac{|x_{0}|}{s_{x}}$ is the point where
$\lambda$ is just large enough to push $x$ to zero.

###### Proof.

Convert to nested optimization and exploit the fact that the solution for
known $\lambda$ is given by the soft thresholding operator:

$\displaystyle\underset{x\in\mathbb{R},\lambda>0}{\mathrm{argmin}}\,\,\lambda|x|+\frac{(x-x_{0})^{2}}{2s_{x}}+\frac{(\lambda-\lambda_{0})^{2}}{2s_{\lambda}}$
(29)
$\displaystyle\iff\underset{\lambda>0}{\mathrm{argmin}}\,\,\frac{(\lambda-\lambda_{0})^{2}}{2s_{\lambda}}+\underset{x\in\mathbb{R}}{\mathrm{argmin}}\,\,\lambda|x|+\frac{(x-x_{0})^{2}}{2s_{x}}$
(30)
$\displaystyle\iff\underset{\lambda>0}{\mathrm{argmin}}\,\,\frac{(\lambda-\lambda_{0})^{2}}{2s_{\lambda}}+\lambda(|x_{0}|-s_{x}\lambda)^{+}+\frac{(|x_{0}|-s_{x}\lambda)^{+,2}-2|x_{0}|(|x_{0}|-s_{x}\lambda)^{+}}{2s_{x}}\,\,.$
(31)

∎

######  Theorem 1.

The optimizing $\lambda$ for the proximal program P1 is given by, when
$s_{x}s_{\lambda}<1$:

$\lambda^{*}=\begin{cases}\lambda_{0}&\lambda_{0}\geq\frac{|x_{0}|}{s_{x}}\\\
\frac{(\lambda_{0}-s_{\lambda}|x_{0}|)^{+}}{1-s_{\lambda}s_{x}}&o.w.\,\,\,\,,\end{cases}$
(32)

and by
$\lambda^{*}=\mathbbm{1}_{\big{[}\frac{\lambda_{0}}{\sqrt{s_{l}}}>\frac{|x_{0}|}{\sqrt{s_{x}}}\big{]}}\lambda_{0}$
(here $\mathbbm{1}$ denotes the indicator function) otherwise. In either case
$x^{*}=(|x_{0}|-s_{x}\lambda^{*})^{+}\mathrm{sgn}(x_{0})$.

###### Proof.

We need only find the optimum of each interval. When
$\lambda\geq\frac{|x_{0}|}{s_{x}}$, the optimum is simply as close as we can
get to $\lambda_{0}$, namely
$\lambda\leftarrow\max[\lambda_{0},\frac{x_{0}}{s_{x}}]$. On the other hand,
when $\lambda\leq\frac{|x_{0}|}{s_{x}}$, if $s_{x}s_{\lambda}<1$, the optimum
is as close as we can get to the stationary point
$\frac{(\lambda_{0}-s_{\lambda}|x_{0}|)}{1-s_{\lambda}s_{x}}$, explicitly
$\lambda\leftarrow\min[\frac{(\lambda_{0}-s_{\lambda}|x_{0}|)^{+}}{1-s_{\lambda}s_{x}},\frac{|x_{0}|}{s_{x}}]$.
When $s_{x}s_{\lambda}\geq 1$, however, the solution is at one of the interval
boundaries $[0,\frac{|x_{0}|}{s_{x}}]$; the boundaries have costs of
$\frac{\lambda_{0}^{2}}{2s_{\lambda}}$ and
$\frac{(\frac{|x_{0}|}{s_{x}}-\lambda_{0})^{2}}{2s_{\lambda}}+\frac{x_{0}^{2}}{2s_{x}}$,
respectively, and so we choose $\lambda\leftarrow 0$ if
$\frac{\lambda_{0}^{2}}{2s_{\lambda}}<\frac{(\frac{|x_{0}|}{s_{x}}-\lambda_{0})^{2}}{2s_{\lambda}}+\frac{x_{0}^{2}}{2s_{x}}$
and $\lambda\leftarrow\frac{|x_{0}|}{s_{x}}$ otherwise. But the cost at
$\lambda=\lambda_{0}$ is only $\frac{x_{0}^{2}}{2s_{x}}$, so the choice is
between $0$ and $\lambda_{0}$ with costs
$\frac{\lambda_{0}^{2}}{2s_{\lambda}}$ and $\frac{x_{0}^{2}}{2s_{x}}$. ∎

######  Lemma 2.

The following hold, where $\lambda^{*}$ denotes the optimizing $\lambda$, and
is formally a function of $\tau$ and $|\beta|$:

  1. 1.

$\lambda^{*}=\frac{1}{\tau|\beta|+\rho^{\prime}(\lambda^{*})}$.

  2. 2.

$\frac{\partial\lambda^{*}}{\partial|\beta|}=-\frac{\tau}{\frac{1}{\lambda^{*2}}+\rho^{\prime\prime}(\lambda)}$.

  3. 3.

$g_{\tau}^{\prime}(|\beta|)=\tau\lambda^{*}$.

  4. 4.

$g_{\tau}^{\prime\prime}(|\beta|)=-\frac{\tau_{n}^{2}}{(\tau_{n}+\rho^{\prime}(\lambda^{*}))^{2}+\rho^{\prime\prime}(\lambda^{*})}$.

###### Proof.

It will be convenient to develop notation for the cost function inside our
penalty:
$g_{\tau}(|\beta|)=\underset{\lambda>0}{\min}\,\,\big{[}\tau\lambda|\beta|-\log\lambda+\rho(\lambda)\big{]}:=\underset{\lambda>0}{\min}\,\,c^{p}(|\beta|,\lambda)$.
For 1, since $\lambda^{*}$ is the optimizing $\lambda$, and due to the
$-\log\lambda$ constraining the optimum to be an interior point, we know that
$0=\frac{\partial}{\partial\lambda}\big{[}\tau\lambda|\beta|-\log\lambda+\rho(\lambda)\big{]}=\tau|\beta|-\frac{1}{\lambda}+\rho^{\prime}(\lambda)$.
For 2, we can use implicit differentiation on this same equation. For 3, we
simply note that
$\frac{\partial}{\partial|\beta|}g_{\tau}(|\beta|)=\frac{\partial}{\partial|\beta|}\big{[}\tau\lambda^{*}|\beta|-\log\lambda^{*}+\rho(\lambda^{*})\big{]}=\tau\lambda^{*}+\frac{\partial\lambda^{*}}{\partial|\beta|}\frac{\partial
c^{p}(|\beta|,\lambda)}{\partial\lambda}\Bigr{|}_{\lambda^{*}}=\tau\lambda^{*}$.
4 proceeds by differentiating 3 and plugging in 1 and 2. ∎

######  Theorem 2.

Assume that the logarithmic derivative of the hyperprior density on $\lambda$
is bounded ($|\rho^{\prime}(\lambda)|\leq M_{1}\,\,\forall\lambda\geq 0$) and
that the density is decreasing on $(0,\infty)$. Then:

  1. 1.

$g^{\prime}_{\tau}(|\beta|)\approx\frac{1}{|\beta|}$ for large $\beta$.

  2. 2.

The minimum of $|\beta|+g^{\prime}_{\tau}(|\beta|)$ is achieved at $\beta=0$
with value $\lambda_{a}\tau$.

###### Proof.

For 1, $g^{\prime}(|\beta|)=\frac{\tau}{\tau|\beta|+\rho(\lambda^{*})}$, and
since $\rho(\lambda)$ is bounded,
$\underset{|\beta|\to\infty}{\lim}\frac{\tau}{\tau|\beta|+\rho(\lambda^{*})}=\frac{1}{|\beta|}$.
For 2, let $\lambda_{a}$ be the $\lambda$ such that
$\frac{1}{\lambda_{a}}=\rho^{\prime}(\lambda_{a})$ (which is unique by the
assumption that $\rho$ is increasing). Note that $\lambda_{a}\leq\lambda^{*}$
and $\lambda_{a}=\lambda^{*}(0)$. So each term of $|\beta|+\lambda^{*}\tau$ is
decreasing in $|\beta|$ individually, and so the minimum of their sum must
occur at $0$, yielding value $\lambda_{a}\tau$. ∎

######  Theorem 3.

Let $\tau_{n}=n\tau_{0}$ for $\tau_{0}>0$, and further assume that
$|\rho^{\prime\prime}(|\lambda|)|<M_{2}$ (bounded second logarithmic
derivative). Then, under the following standard regularity conditions on the
likelihood:

  1. 1.

The data $\mathbf{y}_{i}$ are i.i.d. with density function
$f(\mathbf{y};\bm{\beta})$ providing for common support and model
identifiability. We assume it has a score function with expectation zero
$\mathbb{E}_{\bm{\beta}}\Big{[}\nabla_{\bm{\beta}}\log
f(\mathbf{y};\bm{\beta})\Big{]}=\mathbf{0}$ and a Fisher information
expressible in terms of second derivatives:
$I(\bm{\beta})=\mathbb{E}_{\bm{\beta}}\Big{[}\nabla_{\bm{\beta}}^{2}\log
f(\mathbf{y};\bm{\beta})\Big{]}$.

  2. 2.

The information matrix is finite and positive definite when
$\bm{\beta}=\bm{\beta}_{0}$, with $\bm{\beta}_{0}$ the true parameter vector.

  3. 3.

For some open subset $\mathcal{B}$ containing $\bm{\beta}_{0}$, for almost all
$\mathbf{y}$, the density is thrice differentiable
$\forall\bm{\beta}\in\mathcal{B}$ and that $\Bigr{|}\frac{\partial^{3}\log
f(\mathbf{y};\bm{\beta})}{\partial\beta_{i}\partial\beta_{j}\partial\beta_{k}}\Bigr{|}\leq
M_{i,j,k}(\mathbf{y})$, also over $\mathcal{B}$, where the functions $M$ are
such that $\mathbb{E}_{\bm{\beta}_{0}}[M_{i,j,k}(\mathbf{y})]<\infty$.

there is a local minimum of
$Q(\bm{\beta})=-L(\bm{\beta})+\sum_{p=1}^{P}g_{\tau_{n}}(|\beta_{p}|)$ that
satisfies the following:

  1. 1.

$\hat{\bm{\beta}}_{2}=\mathbf{0}$ with probability approaching 1 as
$n\to\infty$.

  2. 2.

$\hat{\bm{\beta}}_{1}$ is asymptotically normal with covariance given
approximately by $\frac{1}{n}I(\bm{\beta}_{1})$, the Fisher information matrix
considering only active variables.

###### Proof.

We begin by establishing the existence of sparse local minima. Let
$\mathcal{A}_{n}=\\{\bm{\beta}_{1}^{\prime}:||\bm{\beta}_{1}^{\prime}-\bm{\beta}_{10}||<\frac{C_{1}}{\sqrt{n}}\\}$.
We want to show that, asymptotically, for $\bm{\beta}_{1}\in\mathcal{A}_{n}$:

$\underset{\bm{\beta}_{2}}{\min}\,\,Q\Big{(}{\begin{bmatrix}\bm{\beta}_{1}\\\
\bm{\beta}_{2}\end{bmatrix}}\Big{)}=Q\Big{(}{\begin{bmatrix}\bm{\beta}_{1}\\\
\mathbf{0}\end{bmatrix}}\Big{)}$ (33)

We can do this by showing that $\textrm{sgn}(\frac{\partial
Q(\bm{\beta})}{\partial\beta_{j}})=\textrm{sgn}(\beta_{j})$ for $j>r$ and for
$\bm{\beta}_{2}$ sufficiently small. Starting with the expression for the
gradient and applying a series expansion on the likelihood, we see that:

$\displaystyle\frac{\partial Q(\bm{\beta})}{\partial\beta_{j}}=-\frac{\partial
L(\bm{\beta})}{\partial\beta_{j}}+g_{\tau_{n}}^{\prime}(|\beta_{j}|)\textrm{sgn}(\beta_{j})=-\frac{\partial
L(\bm{\beta}_{0})}{\partial\beta_{j}}-\sum_{p=1}^{P}\frac{\partial^{2}L(\bm{\beta}_{0})}{\partial\beta_{j}\partial\beta_{p}}(\beta_{l}-\beta_{l,0})$
(34)
$\displaystyle-\sum_{p_{1},p_{2}=1}^{P}\frac{\partial^{3}L(\bm{\beta}^{*})}{\partial\beta_{j}\partial\beta_{p_{1}}\partial\beta_{p_{2}}}(\beta_{p_{1}}-\beta_{p_{1},0})(\beta_{p_{2}}-\beta_{p_{2},0})+g^{\prime}_{\tau_{n}}(|\beta_{j}|)\textrm{sgn}(\beta_{j})$
(35)

where $\bm{\beta}^{*}$ is associated with the Cauchy form for the remainder.
Since $||\bm{\beta}_{1}-\bm{\beta}_{1,0}||=O_{p}(\frac{1}{\sqrt{n}})$ by
assumption and plugging in the expression for
$g_{\tau_{n}}^{\prime}(|\beta_{j}|)$ from Lemma 2, when
$\|\bm{\beta}_{2}\|<\frac{1}{n}$, we have that:

$\displaystyle\frac{\partial
Q(\bm{\beta})}{\partial\beta_{j}}=O_{p}(\sqrt{n})+\frac{n\tau_{0}}{n\tau_{0}|\beta_{j}|+\rho^{\prime}(\lambda^{*})}\textrm{sgn}(\beta_{j})$
(36)

but when $\beta_{j}<\frac{1}{n}$, we have that
$g_{\tau_{n}}^{\prime}(|\beta_{j}|)\geq
n\frac{\tau_{0}}{\tau_{0}+\rho^{\prime}(\lambda^{*})}$, which dominates the
$O_{p}(\sqrt{n})$ likelihood term and means that $\textrm{sgn}(\beta_{j})$
determines the sign of the partial derivative, and establishes the desired
existence of sparse local minima.

Next, we want to show that within this set $\mathcal{A}_{n}$ which induces
sparse minima in $\bm{\beta}_{2}$ is a local minimizer of the likelihood with
respect to $\bm{\beta}_{1}$. To this end, let vector
$\mathbf{u}\in\mathbb{R}^{r}$ be such that $\|\mathbf{u}\|_{2}=C_{1}$. We
consider the difference in penalized likelihood between $\bm{\beta}_{10}$ and
$\bm{\beta}_{10}+\frac{\mathbf{u}}{\sqrt{N}}$, and again take the series
expansion about $\bm{\beta}_{0}$:

$\displaystyle
Q(\bm{\beta}_{0}+\frac{\mathbf{u}}{\sqrt{n}})-Q(\bm{\beta}_{0})=-L(\begin{bmatrix}\bm{\beta}_{10}+\frac{\mathbf{u}}{\sqrt{n}}\\\
\mathbf{0}\end{bmatrix})+L(\begin{bmatrix}\bm{\beta}_{10}\\\
\mathbf{0}\end{bmatrix})+\sum_{i=1}^{r}g_{\tau_{n}}(|\beta_{j}|)-g_{\tau_{n}}(|\beta_{j0}|)$
(37) $\displaystyle=-\frac{1}{\sqrt{n}}\frac{\partial
L(\bm{\beta}_{0})}{\partial\bm{\beta}_{1}}^{\top}\mathbf{u}+\frac{1}{2}\mathbf{u}^{\top}I_{1}(\bm{\beta}_{1})\mathbf{u}\big{[}1+o_{p}(1)\big{]}+\sum_{p=1}^{r}\Big{[}\frac{1}{\sqrt{n}}g_{\tau_{n}}^{\prime}(|\beta_{p}|)\textrm{sgn}(\beta_{p})u_{p}+\frac{1}{n}g_{\tau_{n}}^{\prime\prime}(|\beta_{p}|)u_{p}^{2}\big{[}1+o_{p}(1)\big{]}\Big{]}$
(38)

Notice that
$g_{\tau_{n}}^{\prime}(|\beta_{p}|)=\frac{n\tau_{0}}{n\tau_{0}|\beta_{p}|+\rho^{\prime}(\lambda^{*})}\overset{n\to\infty}{\to}\frac{1}{|\beta_{p}|}$
so
$\frac{1}{\sqrt{n}}g_{\tau_{n}}^{\prime}(|\beta_{p}|)\overset{n\to\infty}{\to}0$
and that
$g_{\tau_{n}}^{\prime\prime}(|\beta_{p}|)=-\frac{\tau_{n}^{2}}{(\tau_{n}+\rho^{\prime}(\lambda^{*}))^{2}+\rho^{\prime\prime}(\lambda^{*})}$
also converges to a finite constant as $n\to\infty$, so
$\frac{1}{n}g_{\tau_{n}}^{\prime\prime}(|\beta_{p}|)\to 0$ as well. For
sufficiently large $C_{1}$, the quadratic likelihood term thus dominates all
other terms and the difference is positive, and there is therefore a local
minimizer with respect to $\bm{\beta}_{1}$ within $\mathcal{A}_{n}$.

To determine its asymptotic distribution and conclude the proof, we examine
the stationarity conditions for $j\in\\{1,\ldots,r\\}$ of this local optimum:

$\displaystyle 0=\frac{\partial Q({\tiny\begin{bmatrix}\bm{\beta}_{1}\\\
\mathbf{0}\end{bmatrix}})}{\partial\beta_{j}}\Bigr{|}_{\beta=\big{(}\begin{subarray}{c}\hat{\bm{\beta}_{1}}\\\
\mathbf{0}\end{subarray}\big{)}}=-\frac{\partial
L(\bm{\beta})}{\partial\beta_{j}}\Bigr{|}_{\beta=\big{(}\begin{subarray}{c}\hat{\bm{\beta}_{1}}\\\
\mathbf{0}\end{subarray}\big{)}}+g_{\tau_{n}}^{\prime}(|\hat{\beta}_{j}|)\textrm{sgn}(\beta_{j})=\frac{\partial
L(\bm{\beta}_{0})}{\partial\beta_{j}}$ (39)
$\displaystyle+\sum_{p=1}^{r}\Big{[}\frac{\partial
L(\bm{\beta}_{0})}{\partial\beta_{j}\partial\beta_{p}}+o_{p}(1)\Big{]}(\hat{\beta}_{p}-\beta_{p0})+g^{\prime}_{\tau_{n}}(|\beta_{j}|)\textrm{sgn}(\beta_{j})+\Big{[}g_{\tau_{n}}^{\prime\prime}(|\beta_{j}|)+o_{p}(1)\Big{]}(\hat{\beta_{j}}-\beta_{j,0})\,\,.$
(40)

Slutsky’s together with the Central Limit theorems thus give us that:

$\sqrt{n}\big{(}I_{1}(\bm{\beta}_{1,0})+\frac{1}{n}\mathbf{D}\big{)}\Big{(}\hat{\bm{\beta}}_{1}-\bm{\beta}_{1,0}+\frac{1}{n}\big{(}I_{1}(\bm{\beta}_{1,0})+\frac{1}{n}\mathbf{D}\big{)}^{-1}\mathbf{g}\Big{)}\overset{d}{\to}\mathcal{N}\big{(}\mathbf{0},I_{1}(\bm{\beta}_{1,0})\big{)}\,\,,$
(41)

denoting:

$\mathbf{D}=\mathrm{diag}([g^{\prime\prime}(|\beta_{1}|),\ldots,g^{\prime\prime}(|\beta_{r}|)]);\,\,\,\,\mathbf{g}=[g^{\prime}_{\tau_{n}}(|\beta_{j}|)\textrm{sgn}(\beta_{1}),\ldots,g^{\prime}_{\tau_{n}}(|\beta_{j}|)\textrm{sgn}(\beta_{1})]$
(42)

∎

##  2 Optimization Details

Step Size Selection with Trust Regions: We deploy an adaptive step size based
on Trust Regions. The trust region framework involves comparing the observed
function reduction with that implied, typically, by the second order model
itself implied by our gradient descent:

$\hat{f}(\mathbf{x})\approx f(\mathbf{x}^{t})+\nabla
f(\mathbf{x}^{t})^{\top}(\mathbf{x}-\mathbf{x}^{t})+||\mathbf{x}^{t}-\mathbf{x}||_{\mathbf{C}}^{2}$
(43)

with the expected reduction in cost given by
$\rho_{tr}:=\frac{f(\mathbf{x}^{t+1})-f(\mathbf{x}^{t})}{\hat{f}(\mathbf{x}^{t+1})-\hat{f}(\mathbf{x}^{t})}$,
and where the norm $\mathbf{C}$ is given in our case by our diagonal step
sizes:
$||\mathbf{x}||_{\mathbf{C}}^{2}=\sum_{p=1}^{P}\frac{x_{p}^{2}}{\eta_{p}}$.

In our case, our proximal approach means we have to include our regularization
term in the trust region surrogate, and burden our notation by splitting the
optimization variables into unpenalized variables $\theta$, penalized
variables $\beta$, and penalizing variables $\lambda$, so
$\mathbf{x}=[\theta,\beta,\lambda]$:

$\displaystyle\hat{f}([\theta,\beta,\lambda])\approx
f([\theta^{t},\beta^{t},\lambda^{t}])+\nabla
f([\theta^{t},\beta^{t},\lambda^{t}])^{\top}([\theta^{t},\beta^{t},\lambda^{t}]-[\theta,\beta,\lambda])$
$\displaystyle+||[\theta^{t},\beta^{t},\lambda^{t}]-[\theta,\beta,\lambda]||_{\mathbf{C}}^{2}+\sum_{p=1}^{P}\lambda_{p}|\beta_{p}|,$

but is otherwise identical. When the $\rho_{rt}$ is too small, it means the
expected reduction was much less than the actual reduction, and the trust
region shrinks. Conversely, a ratio closer to $1$ is indicative of an overly
cautious stepsize. A trust region procedure proceeds by shrinking the step
size by some factor $m_{s}$ whenever $\rho_{tr}<r_{l}$ is too small, expanding
it by $m_{e}$ if $\rho_{tr}>r_{u}$ is too big, and leaving it unchanged when
$\rho_{tr}$ is within a prespecified interval $[l_{\rho},u_{\rho}]$. We use
the default parameters recommended by the McCormick School’s optimization Wiki
444https://optimization.mccormick.northwestern.edu/index.php/Trust-
region_methods, which are $m_{s}=0.25,m_{e}=2,r_{l}=0.25,r_{u}=0.75$. Steps
which worsen the cost function are rejected and the step size is again shrunk
by $m_{s}$.

Nesterov Acceleration: Along with careful step-size selection, we implement
Nesterov acceleration (Nesterov, 1983), which the unfamiliar reader should
think of as a careful kind of momentum. Nesterov acceleration has been shown
to make a big difference in proximal gradient methods in the form of the ISTA
algorithm, which is called Fast ISTA (Beck and Teboulle, 2009) when endowed
with Nesterov acceleration.

One point of implementation we should mention is that Nesterov’s search
direction does not guarantee descent, unlike a gradient or proximal gradient
direction. This can be problematic since we are using a trust region to select
step sizes, which expects our model to provide descent directions for
sufficiently small step sizes. Various strategies exist to deal with this; we
take the approach of simply resetting the Nesterov variables each time the
trust region shrinks 3 times in a row (and resetting the step size to its
previous value).

Preconditioning: Nesterov Acceleration gets us a faster asymptotic convergence
rate. The other half of speedy optimization is preconditioning, which
essentially involves improving the constants in the convergence rate. By
preconditioner, we simply mean a positive definite linear transformation of
the gradient: $\tilde{\mathbf{g}}\leftarrow\mathbf{C}\mathbf{g}$. We studied
the Natural Gradient Method (Amari, 1998), which uses the Fisher information
associated with the variational distribution as a preconditioner, as well as
“bespoke” preconditioners for various likelihood terms. But ultimately, we
found maximum generalizability in simply using an exponential moving average
of gradient norms as our preconditioner, á la Adam (Kingma and Ba, 2015). It
is very likely that significant improvements in execution time on general
problems can be attained by better preconditioners.

Ablation Study: We demonstrate the importance of Nesterov acceleration and
proper preconditioning on our working problem. As Figure 8 shows, the gradient
search direction with as large a step size as we observed to be stable is much
slower than that with a trust region-selected step size. In turn, Nesterov
acceleration is much faster than the trust region alone. And finally,
preconditioning improves optimization speed even more.

Figure 8: Accelerating Gradient Descent: Using adaptive step sizes via a Trust
Region, better conditioning via a preconditioner, and acceleration via
Nesterov’s method yield much improved optimization performance on this
problem.

##  3 SBL in Linear Gaussian Models

Figure 9: Nonconvex Behavior. Top Left and Mid: MAP inference with uniform
prior leads to reasonable local optima, but likelihood is unbounded. Lower and
Mid Left: Penalty $\tau$ varies: gray region is correct sparsity level, Lasso
is biased within; SBL is debiased (like nonconvex penalties for likelihoods)
and provides credible intervals. Top Right: Instability in RNG for selection
of correlated parameters. Bottom Right: Cost vs Initialization, true parameter
value in red.

Denote by $\mathbf{P}_{X}$ the orthogonal projector onto the range of
$\mathbf{X}$. If we assume
$\mathbf{y}\sim\mathrm{N}(\mathbf{X}\bm{\beta},\sigma^{2}\mathbf{I})$ with
$\bm{\beta}\sim Q_{\bm{\beta}}$ and $\sigma^{2}\sim Q_{\sigma^{2}}$, then we
have the following expected log likelihood:

$\displaystyle-\underset{\bm{\beta}\sim Q_{\bm{\beta}},\sigma^{2}\sim
Q_{\sigma^{2}}}{\mathbb{E}}[-\mathcal{L}(\bm{\beta},\sigma^{2}|\mathbf{y})]=\underset{\bm{\beta}\sim
Q_{\bm{\beta}},\sigma^{2}\sim
Q_{\sigma^{2}}}{\mathbb{E}}[\frac{N}{2}\log\sigma^{2}+\frac{||\mathbf{y}-\mathbf{X}\bm{\beta}||_{2}^{2}}{2\sigma^{2}}]+C$
(44) $\displaystyle=\underset{\sigma^{2}\sim
Q_{\sigma^{2}}}{\mathbb{E}}[\mathrm{tr}[\frac{\mathbf{X}^{\top}\mathbf{X}}{2\sigma^{2}}\mathbb{V}_{Q}[\bm{\beta}]]+(\mathbb{E}_{Q}[\bm{\beta}]-\hat{\bm{\beta}})^{\top}\frac{\mathbf{X}^{\top}\mathbf{X}}{2\sigma^{2}}(\mathbb{E}_{Q}[\bm{\beta}]-\hat{\bm{\beta}})+\frac{\mathbf{y}^{\top}(\mathbf{I}-\mathbf{P}_{X})\mathbf{y}}{2\sigma^{2}}+\frac{N}{2}\log\sigma^{2}]$
(45) $\displaystyle=\underset{\sigma^{2}\sim
Q_{\sigma^{2}}}{\mathbb{E}}[\frac{\mathrm{tr}[\mathbf{X}^{\top}\mathbf{X}\mathrm{diag}(2\nu^{2})]+(\eta-\hat{\bm{\beta}})^{\top}\mathbf{X}^{\top}\mathbf{X}(\eta-\hat{\bm{\beta}})+\mathbf{y}^{\top}(\mathbf{I}-\mathbf{P}_{X})\mathbf{y}}{2\sigma^{2}}+\frac{N}{2}\log\sigma^{2}]\,.$
(46)

To evaluate the expectation wrt $\sigma^{2}$, we’ll define the numerator of
the first term as $\xi$ which yields the following:

$=\underset{\sigma^{2}\sim
Q_{\sigma^{2}}}{\mathbb{E}}[\frac{\xi}{2\sigma^{2}}+\frac{N}{2}\log\sigma^{2}]=\frac{\xi}{2}\underset{\sigma^{2}\sim
Q_{\sigma^{2}}}{\mathbb{E}}[\frac{1}{\sigma^{2}}]+\frac{N}{2}\underset{\sigma^{2}\sim
Q_{\sigma^{2}}}{\mathbb{E}}[\log\sigma^{2}]\,.$ (47)

We see that if we want a closed form estimate of KL we need to choose
$Q_{\bm{\beta}}$ to have closed form mean and variance and $Q_{\sigma^{2}}$ to
have a closed form reciprocal and logarithmic expectation. The Gamma and
Lognormal distribution satisfy this. However, we observed difficulties
differentiating through the Gamma simulation in tensorflow_probability, so we
settled on the Lognormal which had better numerical behavior.

Normal and lognormal variational-prior models have the same well known KL
divergence expression.

##  4 $\lambda$-Prior Specification in Orthogonal Designs

Working example $\mathbf{X}\in\mathbb{R}^{N\times P}=\mathbb{R}^{250\times
10}$ with entries sampled from $i.i.d$ standard Gaussians,
$\bm{\beta}\in\mathbb{R}^{50}$ with only 6 nonzero entries, of
$[-2.5,-2,-1.5,1.5,2,2.5]$ (similar $\beta$ as Ročková and George (2018)), and
$\mathbf{y}=\mathbf{X}\bm{\beta}+\bm{\epsilon}$ with again iid standard normal
$\epsilon_{n}$. Figure 9, a shows a more readable problem with only $[-2,2]$
nonzero in $\beta$ and $N=100$.

We now consider the orthogonal linear case where the problem decomposes axis-
by-axis, starting again by converting to a nested optimization marginalizing
over $x$, yielding:

$\displaystyle\underset{\beta\in\mathbb{R},\lambda\in\mathbb{R}^{+}}{\min}\frac{\sum_{i=1}^{N}(y_{i}-x_{i}\beta)^{2}}{2\sigma^{2}}+\tau\lambda|\beta|-\log\lambda-\log
p_{\lambda}(\lambda)$ (48)
$\displaystyle\iff\underset{\beta\in\mathbb{R},\lambda\in\mathbb{R}^{+}}{\min}\frac{(\beta-\hat{\beta})^{2}}{2\frac{\sigma^{2}}{\sum_{i=1}^{N}x_{i}^{2}}}+\tau\lambda|\beta|-\tau\log\lambda-\log
p_{\lambda}(\lambda)\,\,,$ (49)

where $\hat{\beta}$ gives the least squares estimate. Again, the STO gives us
profile cost for $\lambda$ with
$\beta^{*}(\lambda)=(|\hat{\beta}|-\frac{\sigma^{2}\tau\lambda}{\sum_{i=1}^{N}x_{i}^{2}})^{+}\mathrm{sgn}(\hat{\beta})$.
Using the fact that
$(\beta^{*}(\lambda)-\hat{\beta})^{2}=\min(\frac{\sigma^{2}\tau\lambda}{\sum_{i=1}^{N}x_{i}^{2}},|\hat{\beta}|)^{2}$
, we can express the profile cost as a piecewise log-quadratic (plus the prior
term):

$\displaystyle\underset{\lambda\in\mathbb{R}^{+}}{\min}\,\,\frac{\tau^{2}\min(\lambda,\frac{|\hat{\beta}|\sum_{i=1}^{N}x_{i}^{2}}{\tau\sigma^{2}})^{2}}{2}+\tau\lambda\Big{(}|\hat{\beta}|-\frac{\sigma^{2}\tau\lambda}{\sum_{i=1}^{N}x_{i}^{2}}\Big{)}^{+}-\log\lambda-\log
p_{\lambda}(\lambda)$ (50)
$\displaystyle\iff\underset{\lambda\in\mathbb{R}^{+}}{\min}\begin{cases}\tau^{2}\bigg{(}\frac{1}{2}-\frac{\sigma^{2}}{\sum_{i=1}^{N}x_{i}^{2}}\bigg{)}\lambda^{2}+\tau|\hat{\beta}|\lambda-\log(\lambda)-\log
p_{\lambda}(\lambda);&\lambda\leq\frac{|\hat{\beta}|\sum_{i=1}^{N}x_{i}^{2}}{\tau\sigma^{2}}\\\
\frac{\hat{\beta}^{2}[\sum_{i=1}^{N}x_{i}^{2}]^{2}}{2\sigma^{4}}-\log\lambda-\log
p_{\lambda}(\lambda)\end{cases}$ (51)

Uniform $p_{\lambda}\propto 1$: Without a great idea of what prior we ought to
put on $\lambda$, we might be tempted to place a uniform prior on it. However,
this leads to an unbounded density as
$(\bm{\beta},\bm{\lambda})=(\mathbf{0},\infty)$: the term
$\frac{||\mathbf{y}-\mathbf{X}\bm{\beta}||_{2}^{2}}{2\sigma^{2}}$ becomes
$\frac{||\mathbf{y}||_{2}^{2}}{2\sigma^{2}}$, while the $-\log\lambda_{p}$
terms are unbounded.

Undeterred, we performed MAP inference with a uniform prior on our working
problem using the Adam optimizer (Kingma and Ba (2015)) initialized with
$\bm{\beta}=\mathbf{0}$, $\bm{\lambda}=\mathbf{1}$ and subgradient descent,
where it is seemingly able to perform debiased variable selection (see Figure
9a). The “zero $\bm{\beta}$ big $\bm{\lambda}$” strategy only surpasses the
posterior density of local estimates on this particular example at
$\lambda\approx 10^{200}$; since the $\log\lambda$ term is unbounded, we
expect this behavior for any $\mathbf{X},\mathbf{y},\sigma^{2}$. Evidently,
local optimizations may converge to a reasonable local minimum despite this.
Though the biconvex penalty proposed here has its advantages over convex
penalties, nonconvexity can be a treacherous property.

Power Inverse $p_{\lambda}\propto\frac{1}{\lambda^{a}}$: In this case, setting
$\bm{\beta}=0$ leaves us with a cost of $\sum_{p=1}^{P}(a-1)\log\lambda_{p}$.
We see that for $a<1$, we get similar behavior to the uniform case, where cost
is unbounded below as $\lambda\to 0$. $a>1$ leads to a positive coefficient
for the $\log\lambda$ term, in turn leading to the opposite behavior, where
cost is unbounded below as $\lambda\to 0$. When $a=1$, the log terms cancel,
and we are left with a cost of
$\frac{||\mathbf{y}-\mathbf{X}\bm{\beta}||_{2}^{2}}{2\sigma^{2}}+\sum_{p=1}^{P}\lambda_{p}|\beta_{p}|$.
Without any motivation to stay positive, the $\lambda_{p}$ term will vanish,
yielding the unpenalized problem.

Half-Gaussian $p_{\lambda}\propto e^{-\frac{(\lambda-
m_{\lambda})^{2}}{b_{\lambda}^{2}}}$: Plugging in yields:

$\underset{\lambda\in\mathbb{R}^{+}}{\min}\begin{cases}\bigg{(}\tau^{2}(\frac{1}{2}-\frac{\sigma^{2}}{\sum_{i=1}^{N}x_{i}^{2}})+\frac{1}{b_{\lambda}^{2}}\bigg{)}\lambda^{2}+(\tau|\hat{\beta}|-2\frac{m_{\lambda}}{b_{\lambda}^{2}})\lambda-\log(\lambda)+\frac{m_{\lambda}^{2}}{b_{\lambda}^{2}};&\lambda\leq\frac{|\hat{\beta}|\sum_{i=1}^{N}x_{i}^{2}}{\tau\sigma^{2}}\\\
\frac{\hat{\beta}^{2}[\sum_{i=1}^{N}x_{i}^{2}]^{2}}{2\sigma^{4}}-\log\lambda+\frac{(\lambda-
m_{\lambda})^{2}}{b_{\lambda}^{2}}&o.w.\end{cases}$ (52)

This function, like the proximal cost, may be nonconvex when
$\frac{\sigma^{2}}{\sum_{i=1}^{N}x_{i}^{2}}>\frac{1}{\tau^{2}b_{\lambda}^{2}}+\frac{1}{2}$,
as this is when the quadratic term is concave. But the $-\log\lambda$ term
contributes to convexity near the origin such that the overall function is
convex until
$\frac{1}{\lambda^{2}}>\tau\big{(}\frac{2\tau\sigma^{2}}{\sum_{i=1}^{N}x_{i}^{2}}-\tau\big{)}-\frac{1}{b^{2}_{\lambda}}$.
The stationary point of the upper convex region is
$\lambda^{*}=\frac{m_{\lambda}+\sqrt{m_{\lambda}^{2}+2b_{\lambda}^{2}}}{2}$
(if this lies in the region), and is

$\frac{2b_{\lambda}^{2}s_{x}}{\beta b_{\lambda}^{2}\tau
s_{x}+\sqrt{s_{x}\left(b_{\lambda}^{4}\tau^{2}\left(\left(\beta^{2}+4\right)s_{x}-8\sigma\right)-4b_{\lambda}^{2}s_{x}(\beta
m_{\lambda}\tau-2)+4m_{\lambda}^{2}s_{x}\right)}-2m_{\lambda}s_{x}}$ (53)

in the lower region. The $\log\lambda$ term prevents optima at $\lambda=0$, so
we need only check the cost at these two points to determine which is optimal.

Half-Cauchy $p_{\lambda}\propto\frac{1}{1+\frac{\lambda^{2}}{a}}$: We next
turn our attention to the case $p_{\lambda}=\mathrm{C}(0,a_{\lambda})^{+}$, a
Half-Cauchy prior, in the style of the Horseshoe Prior (but with a Laplace
rather than Normal conditional prior for $\beta$). This leads us to the
following cost function:

$\displaystyle\underset{\lambda\geq
0}{\min}\begin{cases}\tau^{2}\bigg{(}\frac{1}{2}-\frac{\sigma^{2}}{\sum_{i=1}^{N}x_{i}^{2}}\bigg{)}\lambda^{2}+\tau\lambda|\hat{\beta}|+\log\Big{[}\frac{1+\frac{\lambda^{2}}{a_{\lambda}}}{\lambda}\Big{]};&\lambda\in\lambda\leq\frac{|\hat{\beta}|\sum_{i=1}^{N}x_{i}^{2}}{\tau\sigma^{2}}\\\
\frac{\hat{\beta}^{2}}{2\sigma^{4}}+\log\Big{[}\frac{1+\frac{\lambda^{2}}{a_{\lambda}}}{\lambda}\Big{]};&o.w.\end{cases}$
(54)

Finding the critical points of the function in
$[0,\frac{|\hat{\beta}|\sum_{i=1}^{N}x_{i}^{2}}{\tau\sigma^{2}}]$ may be cast
as solving a quartic polynomial equation. While we are fortunate that such an
equation can be solved in closed form, we have not been able to extract much
intuition from it. Nevertheless, the variable penalty $\ell_{1}$ proximal
operator enables iterative numerical solution of this problem.

##  5 Additional Experimental Results and Details

###  5.1 Simulation Study

$\tau$ values were selected manually which gave good performance in simulation
studies. In practice, we recommend evaluating the SBL variational distribution
for a whole spectrum of $\tau$ rather than choosing a particular one.

####  5.1.1 Nonconvexity and Instability

Since our penalty is nonconvex, the final parameters of an optimization are
dependent on their initializations. Figure 9c shows the final costs associated
with a grid of initializations on a 2D toy problem, where $\beta=[2,0]$. The
initial values are used on a grid of 51 points between [-4,4] and with $\tau$
set to 100. We perform MAP inference to focus on the nonconvexity in the
coefficients (rather than on variational standard deviations). We see that in
this case, the model lands on the true solution when $\beta_{1}$ is
initialized at large values. For this simulation and throughout the article,
$\lambda_{p}$ were all initialized at 1 and $\beta$ at 0. The classical Lasso
is famously unstable insofar as it essentially randomly selects from groups of
correlated parameters. We conduct a numerical experiment to confirm the
inheritance of this property, using the same toy problem, but now with
$\mathbf{X}$ data simulated from a correlated normal distribution with high
correlation, and this time using the SBL. Figure 9d shows how, starting at a
correlation of $0.999$, as we approach a correlation of $1$, the model begins
changing from nearly always selecting $\beta_{1}$ as nonzero to doing so only
$50\%$ of the time.

###  5.2 Additional Case Study Implementation Details

Symbol | Generating Dist | Variational Dist | Name
---|---|---|---
$\beta$ | $\mathrm{L}(0,\frac{1}{\lambda_{\ell_{1}}})$ | $\mathrm{L}(\eta_{\beta},\nu_{\beta})$ | Regression Coefficients
$\sigma^{2}$ | $\log\mathrm{N}(0,b_{\sigma})$ | $\log\mathrm{N}(\eta_{\sigma},\nu_{\sigma})$ | Overdispersion Parameter
$\omega$ | $\mathrm{N}(0,\rho^{2})$ | $\mathrm{N}(\eta_{\omega},\nu_{\omega})$ | Spatiotemporal Random Effects
$\rho^{2}$ | $\log\mathrm{N}(0,b_{\rho})$ | $\log\mathrm{N}(\eta_{\rho},\nu_{\rho})$ | Random Effects Variance
$\phi$ | $\mathrm{U}(0,1)$ | $\mathrm{logit}\mathrm{N}(\eta_{\phi},\nu_{\phi})$ | Spatial/Temporal Correlation
$\pi$ | $\mathrm{U}(0,1)$ | $\mathrm{logit}\mathrm{N}(\eta_{\pi},\nu_{\pi})$ | Zero Inflation Probability
Table 2: Model Estimands along with their assumed prior and variational
distributions.

####  5.2.1 Lagging and Aggregation

In order to account for lag between tweet measurement and recording of a
migration event, we use a Gaussian filter. Recall that we have periods $r$ of
variable length, with beginning $l_{r}$ and end $u_{r}$, and duration
$s_{r}:=u_{r}-l_{r}$. The filter, parameterized by a mean lookback time
$\mu_{l}$ and smoothing window size $\sigma_{l}$, is defined as:

$z_{r,o,d,p}\propto\sum_{t^{\prime}=t-w}^{t+w}\frac{f((|\frac{u_{t}+l_{t}}{2}-t^{\prime}|-\frac{s_{t}}{2})^{+},\mu_{l},\sigma_{l}^{2})x_{t^{\prime},o,d,p}}{\sum_{t^{\prime}=t-w}^{t+w}f((|\frac{u_{t}+l_{t}}{2}-t^{\prime}|-\frac{s_{t}}{2})^{+}),\mu_{l},\sigma_{l}^{2})}$
(55)

where $w$ is the maximum window size, $f(.,\mu,\sigma^{2})$ gives the Gaussian
density function. In other words, all tweets within the period $r$ when
shifted by $\mu_{l}$ is assigned maximum density. Events outside of this
window less density, with speed of decay give by $\sigma_{l}^{2}$

####  5.2.2 Spatiotemporal KL Penalty Estimation

In this section we describe how to form the KL terms related to the spatial
random effects $\omega$, the spatial correlation parameter $\phi$ and scaling
parameter $\rho^{2}$. We assume the CAR model
$\omega\overset{P}{\sim}N(\mathbf{0},\rho(\mathbf{I}-\phi\mathbf{M})^{-1})$,
where $\mathbf{M}$ is a spatial connectivity matrix.

We begin with the normal KL, denoting by $B$ the dimension of $\bm{\omega}$:

$\displaystyle\frac{1}{2}(\log|\Sigma_{P}|-\log|\Sigma_{Q}|)-B+B\mathbb{E}[\log\rho^{2}]+\big{[}\mathrm{tr}[\Sigma_{Q}]+\mathbb{E}[\bm{\omega}]^{\top}\Sigma_{P}^{-1}\mathbb{E}[\bm{\omega}]\big{]}\mathbb{E}[\frac{1}{\sigma^{2}}]\,\,,$
(56)

and plugging in the following:

$\displaystyle\Sigma_{P}^{-1}=\mathbf{I}-\phi\mathbf{M}$ (57)
$\displaystyle\log|\Sigma_{P}|=-\sum_{i}\log(1-\phi\lambda_{i}(\mathbf{M}))$
(58) $\displaystyle\log|\Sigma_{Q}|=-\sum_{i}\nu_{\omega_{i}}\,.$ (59)

In the lognormal variational case, the two expectations of $\rho^{2}$ are
available in closed form as $\mathbb{E}[\log\rho^{2}]=\eta_{\rho}$,
$\mathbb{E}[\frac{1}{\rho^{2}}]=e^{-\eta_{\rho}+\frac{\nu_{\rho}}{2}}$. It
remains to integrate out $\phi$, which is not analytical. This we accomplish
via SAA, just as we estimate the negative log likelihood.

###  5.3 Case Study Random Effects

The variational estimates of the random effects are given in Figure 10. The
right set of panels, which give us the random effect variance and thus average
magnitude, indicate that there is not much change in the magnitude of effects
overall as $\tau$ is varied. On the other hand, some individual random effects
do undergo nontrivial changes as $\tau$ varies, and the four most changing
random effects are highlighted in each figure.

Figure 10: Trajectory of Random Effects and Random Effects Variance parameters
in Iraq case study. Figure 11: Antigravitation: Left: Scatterplots do not
reveal a clear relationship between pairwise distance and migrant exchange.
Mid: Boxplots don’t either. Right: y-axis gives the mean distance between
those regions with a migrant exchange of at least the number given by the
x-axis. For some exchange levels, this is above the overall overage, and for
others, below.

###  5.4 The Antigravity Model

The gravity model is so named due to its consideration of the distances
between exchanging regions, supposing that greater distance means less
exchange. But this might not be the case in subnational regions: we observe a
negative coefficient, meaning exchange is, ceteris paribus, increasingly
likely with distance. In a small country where one side of it is at war, it’s
not too surprising that migrants would want to go to the “other side”. Without
random effects and without penalties, distance still assumes a negative,
though smaller, value.

## References

  * Amari (1998) Amari, S.-I. (1998).  Natural gradient works efficiently in learning.  Neural computation 10(2), 251–276.
  * Babacan et al. (2014) Babacan, S. D., S. Nakajima, and M. N. Do (2014).  Bayesian group-sparse modeling and variational inference.  IEEE transactions on signal processing 62(11), 2906–2921.
  * Beck and Teboulle (2009) Beck, A. and M. Teboulle (2009).  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM journal on imaging sciences 2(1), 183–202.
  * Bhadra et al. (2019) Bhadra, A., J. Datta, N. G. Polson, and B. Willard (2019).  Lasso meets horseshoe: A survey.  Statistical Science 34(3), 405–427.
  * Bhattacharya et al. (2012) Bhattacharya, A., D. Pati, N. S. Pillai, and D. B. Dunson (2012).  Bayesian shrinkage.  arXiv preprint arXiv:1212.6088.
  * Blei et al. (2017) Blei, D. M., A. Kucukelbir, and J. D. McAuliffe (2017).  Variational inference: A review for statisticians.  Journal of the American statistical Association 112(518), 859–877.
  * Bühlmann and Meier (2008) Bühlmann, P. and L. Meier (2008).  Discussion: One-step sparse estimates in nonconcave penalized likelihood models.  The Annals of Statistics 36(4), 1534 – 1541.
  * Candes et al. (2008) Candes, E. J., M. B. Wakin, and S. P. Boyd (2008).  Enhancing sparsity by reweighted l 1 minimization.  Journal of Fourier analysis and applications 14(5), 877–905.
  * Carvalho et al. (2010) Carvalho, C. M., N. G. Polson, and J. G. Scott (2010).  The horseshoe estimator for sparse signals.  Biometrika 97(2), 465–480.
  * Chen et al. (2018) Chen, H., A. Mirestean, and C. G. Tsangarides (2018).  Bayesian model averaging for dynamic panels with an application to a trade gravity model.  Econometric Reviews 37(7), 777–805.
  * Congdon (2000) Congdon, P. (2000).  A Bayesian approach to prediction using the gravity model, with an application to patient flow modeling.  Geographical analysis 32(3), 205–224.
  * Daubechies et al. (2004) Daubechies, I., M. Defrise, and C. De Mol (2004).  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  Communications on Pure and Applied Mathematics 57(11), 1413–1457.
  * De Oliveira (2012) De Oliveira, V. (2012).  Bayesian analysis of conditional autoregressive models.  Annals of the Institute of Statistical Mathematics 64(1), 107–133.
  * Efron et al. (2004) Efron, B., T. Hastie, I. Johnstone, and R. Tibshirani (2004).  Least angle regression.  The Annals of Statistics 32(2), 407 – 499.
  * Fan and Li (2001) Fan, J. and R. Li (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association 96(456), 1348–1360.
  * Fan et al. (2014) Fan, J., L. Xue, and H. Zou (2014).  Strong oracle optimality of folded concave penalized estimation.  Annals of statistics 42(3), 819.
  * Ferreira (2019) Ferreira, M. A. (2019).  The limiting distribution of the Gibbs sampler for the intrinsic conditional autoregressive model.  Brazilian Journal of Probability and Statistics 33(4), 734–744.
  * Fu and Knight (2000) Fu, W. and K. Knight (2000).  Asymptotics for lasso-type estimators.  The Annals of statistics 28(5), 1356–1378.
  * Hare and Sagastizábal (2009) Hare, W. and C. Sagastizábal (2009).  Computing proximal points of nonconvex functions.  Mathematical Programming 116(1), 221–258.
  * Helgøy and Li (2019) Helgøy, I. M. and Y. Li (2019).  A Bayesian lasso based sparse learning model.  arXiv preprint arXiv:1908.07220.
  * Hoerl and Kennard (1970) Hoerl, A. E. and R. W. Kennard (1970).  Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics 12(1), 55–67.
  * Hunter and Li (2005) Hunter, D. R. and R. Li (2005).  Variable selection using MM algorithms.  The Annals of Statistics 33(4), 1617 – 1642.
  * Kang and Guo (2009) Kang, J. and J. Guo (2009).  Self-adaptive lasso and its Bayesian estimation.  Technical report, Working Paper.
  * Karemera et al. (2000) Karemera, D., V. I. Oguledo, and B. Davis (2000).  A gravity model analysis of international migration to north america.  Applied Economics 32(13), 1745–1755.
  * Kawano et al. (2015) Kawano, S., I. Hoshina, K. Shimamura, and S. Konishi (2015).  Predictive model selection criteria for Bayesian lasso regression.  Journal of the Japanese Society of Computational Statistics 28(1), 67–82.
  * Kingma and Ba (2015) Kingma, D. P. and J. Ba (2015).  Adam: A method for stochastic optimization.  In ICLR (Poster).
  * Kingma and Welling (2014) Kingma, D. P. and M. Welling (2014).  Auto-Encoding Variational Bayes.  In 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings.
  * Koltchinskii et al. (2011) Koltchinskii, V., K. Lounici, and A. B. Tsybakov (2011).  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion.  The Annals of Statistics 39(5), 2302–2329.
  * Leng et al. (2014) Leng, C., M.-N. Tran, and D. Nott (2014).  Bayesian adaptive lasso.  Annals of the Institute of Statistical Mathematics 66(2), 221–244.
  * Makalic and Schmidt (2015) Makalic, E. and D. F. Schmidt (2015).  A simple sampler for the horseshoe estimator.  IEEE Signal Processing Letters 23(1), 179–182.
  * Mallick and Yi (2014) Mallick, H. and N. Yi (2014).  A new Bayesian lasso.  Statistics and its interface 7(4), 571.
  * Marjanovic and Solo (2013) Marjanovic, G. and V. Solo (2013).  On exact l q denoising.  In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 6068–6072. IEEE.
  * Meinshausen (2007) Meinshausen, N. (2007).  Relaxed lasso.  Computational Statistics & Data Analysis 52(1), 374–393.
  * Meyer (2021) Meyer, G. P. (2021, jun).  An alternative probabilistic interpretation of the huber loss.  In 2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), Los Alamitos, CA, USA, pp. 5257–5265. IEEE Computer Society.
  * Mitchell and Beauchamp (1988) Mitchell, T. J. and J. J. Beauchamp (1988).  Bayesian variable selection in linear regression.  Journal of the American Statistical Association 83(404), 1023–1032.
  * Nesterov (1983) Nesterov, Y. E. (1983).  A method for solving the convex programming problem with convergence rate o (1/k^ 2).  In Dokl. akad. nauk Sssr, Volume 269, pp. 543–547.
  * Owen (2007) Owen, A. B. (2007).  A robust hybrid of lasso and ridge regression.  Contemporary Mathematics 443(7), 59–72.
  * Owen (2013) Owen, A. B. (2013).  Monte Carlo theory, methods and examples.
  * Parikh et al. (2014) Parikh, N., S. Boyd, et al. (2014).  Proximal algorithms.  Foundations and Trends in Optimization 1(3), 127–239.
  * Park and Casella (2008) Park, T. and G. Casella (2008).  The Bayesian lasso.  Journal of the American Statistical Association 103(482), 681–686.
  * Polson et al. (2015) Polson, N. G., J. G. Scott, and B. T. Willard (2015).  Proximal algorithms in statistics and machine learning.  Statistical Science 30(4), 559–581.
  * Poprawe (2015) Poprawe, M. (2015, Jun).  On the relationship between corruption and migration: empirical evidence from a gravity model of migration.  Public Choice 163(3), 337–354.
  * Quaini and Trojani (2022) Quaini, A. and F. Trojani (2022).  Proximal estimation and inference.  arXiv preprint arXiv:2205.13469.
  * Ramos and Suriñach (2017) Ramos, R. and J. Suriñach (2017).  A gravity model of migration between the ENC and the EU.  Tijdschrift voor Economische en Sociale Geografie 108(1), 21–35.
  * Ranganath et al. (2014) Ranganath, R., S. Gerrish, and D. Blei (2014, 22–25 Apr).  Black Box Variational Inference.  In S. Kaski and J. Corander (Eds.), Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, Volume 33 of Proceedings of Machine Learning Research, Reykjavik, Iceland, pp. 814–822. PMLR.
  * Ranjan and Tobias (2007) Ranjan, P. and J. L. Tobias (2007).  Bayesian inference for the gravity model.  Journal of Applied Econometrics 22(4), 817–838.
  * Robinson (1996) Robinson, S. M. (1996).  Analysis of sample-path optimization.  Mathematics of Operations Research 21(3), 513–528.
  * Ročková and George (2018) Ročková, V. and E. I. George (2018).  The spike-and-slab lasso.  Journal of the American Statistical Association 113(521), 431–444.
  * Sandhu et al. (2021) Sandhu, R., M. Khalil, C. Pettit, D. Poirel, and A. Sarkar (2021).  Nonlinear sparse Bayesian learning for physics-based models.  Journal of Computational Physics 426, 109728.
  * She (2009) She, Y. (2009).  Thresholding-based iterative selection procedures for model selection and shrinkage.  Electronic Journal of statistics 3, 384–415.
  * Taylor et al. (1979) Taylor, H. L., S. C. Banks, and J. F. McCoy (1979).  Deconvolution with the $\ell_{1}$ norm.  Geophysics 44(1), 39–52.
  * Terenin et al. (2019) Terenin, A., S. Dong, and D. Draper (2019).  GPU-accelerated Gibbs sampling: a case study of the horseshoe probit model.  Statistics and Computing 29(2), 301–310.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological) 58(1), 267–288.
  * Tibshirani et al. (2005) Tibshirani, R., M. Saunders, S. Rosset, J. Zhu, and K. Knight (2005).  Sparsity and smoothness via the fused lasso.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67(1), 91–108.
  * Tipping (2001) Tipping, M. E. (2001).  Sparse Bayesian learning and the relevance vector machine.  Journal of machine learning research 1(Jun), 211–244.
  * Tung et al. (2019) Tung, D. T., M.-N. Tran, and T. M. Cuong (2019).  Bayesian adaptive lasso with variational Bayes for variable selection in high-dimensional generalized linear mixed models.  Communications in Statistics-Simulation and Computation 48(2), 530–543.
  * Wilson Center (2019) Wilson Center (2019).  Timeline: the rise, spread, and fall of the islamic state.
  * Yuan and Lin (2006) Yuan, M. and Y. Lin (2006).  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), 49–67.
  * Zhang (2010) Zhang, C.-H. (2010).  Nearly unbiased variable selection under minimax concave penalty.  The Annals of Statistics 38(2), 894 – 942.
  * Zou (2006) Zou, H. (2006).  The adaptive lasso and its oracle properties.  Journal of the American statistical association 101(476), 1418–1429.
  * Zou and Hastie (2005) Zou, H. and T. Hastie (2005).  Regularization and variable selection via the elastic net.  Journal of the royal statistical society: series B (statistical methodology) 67(2), 301–320.
  * Zou and Li (2008) Zou, H. and R. Li (2008).  One-step sparse estimates in nonconcave penalized likelihood models.  The Annals of Statistics 36(4), 1509 – 1533.