#  Distributed Decisions on Optimal Load Balancing
in Loss Networks

Qiong Liu , Chenhao Wang* and Ce Zheng Q. Liu is with INFRE, Telecom Paris,
France. email<EMAIL_ADDRESS>C. Wang is the corresponding
author, and is with Beijing Normal University, China. email:
<EMAIL_ADDRESS>Qiong Liu1, Chenhao Wang2, Ce Zheng1  1Télécom Paris,
Institut Polytechnique de Paris, France  2Beijing Normal University, China
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>ce.zheng@telecom-
paris.fr

###### Abstract

When multiple users share a common link in direct transmission, packet loss
and link congestion may occur due to the simultaneous arrival of traffics at
the source node. To tackle this problem, users may resort to an indirect path:
the packet flows are first relayed through a sidelink to another source node,
then transmitted to the destination. This behavior brings the problems of
packet routing or load balancing: (1) how to maximize the total traffic in a
collaborative way; (2) how self-interested users choose routing strategies to
minimize their individual packet loss independently. In this work, we propose
a generalized mathematical framework to tackle the packet and load balancing
issue in loss networks. In centralized scenarios with a planner, we provide a
polynomial-time algorithm to compute the system optimum point where the total
traffic rate is maximized. Conversely, in decentralized settings with
autonomous users making distributed decisions, the system converges to an
equilibrium where no user can reduce their loss probability through unilateral
deviation. We thereby provide a full characterization of Nash equilibrium and
examine the efficiency loss stemming from selfish behaviors, both
theoretically and empirically. In general, the performance degradation caused
by selfish behaviors is not catastrophic; however, this gap is not monotonic
and can have extreme values in certain specific scenarios.

###### Index Terms:

load balancing, Nash equilibria, price of anarchy, network congestion,
sidelink

##  I Introduction

Since the seminal work of Erlang [1], loss networks have played a crucial role
in analyzing and optimizing stochastic systems involving simultaneous resource
utilization, and non-backlogging workloads (for an extensive overview, see
[2]). Meanwhile, in the post-5G era, cloud-enabled networks have emerged as a
dominant architecture, where multiple servers collect data from users and
relay it to a central hub for final processing. To guarantee network efficacy,
that is no server is either overburdened or underutilized, load balancing
strategies are well studied, e.g., [3]. In this context, loss networks provide
valuable mathematical frameworks for comprehending and enhancing load
distribution within cloud-enabled networks.

Early load balancing research for cloud-enabled networks focused on
centralized scenarios, where a centralized planner scheduled workloads to
optimize aspects like performance-energy tradeoffs [3] and algorithmic
considerations [4, 5]. However, due to the stringent latency requirement for
real-time decisions and the increasing signaling overhead caused by the large-
scale deployment of servers and massive users, distributed decisions become a
better solution. In this context, the complexity of the problem increases due
to the non-cooperative and competitive behaviors among users within the
system.

To address the challenges of load balancing in a distributed way, game theory
provides a mathematical framework that describes and analyzes scenarios with
interactive decisions [6]. Till now, some studies have demonstrated the
efficacy of game-theoretic models in addressing load balancing problems. For
instance, Mondal et al. [7] developed a game-theoretic model for load
balancing among competitive cloudlets, while Yi et al. [8] investigated a
similar problem, incorporating additional considerations of queue-aware
strategies. In [9, 10], symmetric loss models where each source has an equal
number of users are considered. However, previous studies mostly focused on
limited cases of identical user strategies, which may not reflect real-world
scenarios, i.e., different users may have different objectives and
preferences. Therefore, further research is needed to develop game-theoretic
models that can address the challenges of load balancing in a more general and
realistic manner.

In this paper, we employ game theory to address load balancing in both
distributed and centralized environments, where users have non-identical
strategies and the number of users is not evenly distributed. Specifically, we
consider the load balancing in a cloud-enabled network consisting of $m$
source nodes (servers) $\\{s_{1},\ldots,s_{m}\\}$ and one destination node
(central hub) $d$. Each source $s_{i}$ has $n_{i}$ users seeking service, and
the traffic originating from each user is assumed to follow an independent
Poisson point process with an identical rate. The nodes in the network are
connected by two types of communication links, namely sidelinks that connect
two sources, and direct links that connect a source and destination. The
sidelink has a random identical independent distribution (i.i.d) loss with a
fixed probability $q$, and the direct link has a congestion loss that depends
on the arrival rate and service rate of each server.

The user cannot split its traffic, and has to determine how to route all of
its traffic from the source node arrived at to the destination node. There are
two approaches for the traffic transmission: a direct path (DP) in which the
packet goes directly from the source arrived at to the destination, and an
indirect path (IP) in which the packet is first relayed to another source node
and then takes the direct link from that node to the destination. We treat
packet loss probability as the performance metric in load balancing, instead
of additive costs like delay or fees in classical routing games [10, 11],
resulting in a non-additive and non-convex optimization process. Each user
aims to minimize its own loss probability and engage in a game by
strategically selecting its own path. In the end, no user can reduce its loss
probability by unilateral deviation and reaches the state of Nash Equilibrium
(NE).

###  I-A Our Contributions

Our work contributes to the load balancing game in the following aspects:
First, we prove two lemmas related to the optimal solution when a centralized
planner exists. Based on these lemmas, a low-complexity algorithm that
maximizes the total traffic is proposed. Second, we study the decentralized
environment where decisions are made by autonomous and self-interested users.
The sufficient and necessary conditions on NE are derived, which depend on the
number of users on direct path and each indirect path. Moreover, since a NE
may be suboptimal, we use the price of anarchy (PoA) [12] to measure the gap
between the NE led by users’ selfish behaviors and the system optimum achieved
by the centralized planner.

The rest of the paper is structured as follows. The formal model and notations
are presented in Section II. In Section III, we provide details to compute the
optimal solution that maximizes the total traffic when a centralized planner
exists. In Section IV, we study the NE in the decentralized decision-making
scenarios, and analyzed the efficiency loss stemming from selfish behaviors.
In Section V, a fine-grained analysis is performed on the existence of NE in
various network configurations for a specific scenario involving two source
nodes. Numerical results are presented and discussed in Section VI. Finally,
Section VII concludes the paper and outlines some future work.

###  I-B Other related works

_Routing games._ As a special class of congestion games, routing games in a
network are problems of routing traffic to achieve the best possible network
performance, and have been studied within various contexts and within various
communities, for example, the mathematics community [13], the
telecommunications [14], and theoretical computer science [15, 16]. The above
references have all in common a cost framework which is additive over links,
such as delays or tolls, and is flow conserving (the amount entering a node
equals the amount leaving it). Routing games with non-additive cost in loss
networks are studied in [17].

_Braess-like paradox in distributed systems._ The Braess-like paradox is said
to occur in a network system with distributed behaviors if adding an extra
link or adding communication capacity to the system leads to a worse system
performance. It widely exists in transportation networks and queuing networks.
Bean et al. [18] show that it can occur in loss networks. Kameda et al. [19]
consider a model similar to ours in that a job (packet) can be processed
directly or indirectly; however, they do not consider the loss probability.
They identify a Braess-like paradox in which adding capacity to the channel
may degrade the system performance on the response time. Kameda and
Pourtallier [20] characterize conditions under which such paradoxical behavior
occurs, and give examples in which the degradation of performance may increase
without bound.

##  II Model and Preliminaries

We abstractly model our problem using a graph. Consider a network with $m$
source nodes $S=\\{s_{1},\ldots,s_{m}\\}$ and one destination node $d$. For
each source node $s_{i}\in S$, let $N_{i}$ be the set of users arriving at
$s_{i}$, and $n_{i}=|N_{i}|$ be the number of such users. Without loss of
generality, we assume $n_{1}\\!>\\!n_{2}\\!>\\!\dots\\!>\\!n_{m}$. Denote
$[m]=\\{1,\ldots,m\\}$. There is a total of $n=\sum_{i\in[m]}n_{i}$ users in
the system, who are self-interested players in the game. We say players and
users interchangeably throughout this paper. Each user is identified with a
flow (or traffic) of packets, which originates from the user and is assumed to
form an independent Poisson process with an identical rate $\phi$. See Fig. 1
for illustration.

Figure 1: An illustration of the game model.

Each user controls its route that all its packets should follow. For a user
associated with $s_{i}\in S$, there are only two types of routes to ship these
packets to the destination $d$: either a direct path (DP) $(s_{i},d)$, or an
indirect two-hop path (IP) $(s_{i},s_{j},d)$ for some $s_{j}\neq s_{i}$, in
which the packet is first sent to another source $s_{j}$ by the side link
$(s_{i},s_{j})$, and then passes through the direct link $(s_{j},d)$.

Strategies. For every source $s_{i}$, each user $k\in N_{i}$ decides a one-
shot strategy
$\mathbf{p}_{k}^{(i)}=\left(p_{k1}^{(i)},\ldots,p_{km}^{(i)}\right)^{T}\in[0,1]^{m}$
with $\sum_{j\in[m]}p_{kj}^{(i)}=1$, where $p_{ki}^{(i)}$ is the probability
of routing all packets through DP, and $p_{kj}^{(i)}~{}(j\in[m],j\neq i)$ is
the probability of routing all packets through IP $(s_{i},s_{j},d)$. When no
confusion arises, we simply write the strategy $\mathbf{p}_{k}^{(i)}$ as
$\mathbf{p}_{k}$. We focus on pure strategies in this paper: a strategy
$\mathbf{p}_{k}$ is pure if $\|\mathbf{p}_{k}\|_{\infty}=1$, i.e., user $k$
deterministically selects a route with probability 1 (for example,
$\mathbf{p}_{k}=(0,1,0,\ldots,0)^{T}$). Let
$\mathbf{p}=(\mathbf{p}_{1},\ldots,\mathbf{p}_{n})$ be the strategy profile of
all users.

Loss probability and loss rate. There are two types of losses: (1) Losses on
side links. We assume that a packet originating from node $s_{i}$ and relayed
to node $s_{j}$ is lost with a fixed probability $q$ for every side link
$(s_{i},s_{j})$, independently of any other loss. Denote by $\bar{q}=1-q$ the
probability that a packet is successfully relayed. (2) Congestion losses on
direct links. We assume that there is no buffer to restore the backlogged
packets, so a packet will be lost when it enters the direct path which is
occupied for the transmission of another packet. The transmission time of a
packet on a direct link $(s_{i},d)$ is a random variable $\sigma$ following a
distribution $\mathcal{X}$, which is assumed to be an identically independent
distribution (i.i.d) for all packets.

Given strategy profile $\mathbf{p}$, user $k\in N_{i}$ continuously sends
packets that follow an independent Poisson process with rate
$p_{ki}^{(i)}\cdot\phi$ to DP $(s_{i},d)$, and an independent Poisson process
of packets with rate $p_{kj}^{(i)}\cdot\phi$ to IP $(s_{i},s_{j},d)$, for any
$s_{j}\neq s_{i}$. Since there is a random loss on the side link
$(s_{i},s_{j})$, the flow of packets from user $k\in N_{i}$ that arrive at the
node $s_{j}$ is also a Poisson process with rate $\bar{q}p_{kj}^{(i)}\phi$.

Thus, for each source $s_{i}\in S$, the flow over the link $(s_{i},d)$ is
Poisson distributed with a traffic rate $T_{i}(\mathbf{p})$ given by

$\displaystyle T_{i}(\mathbf{p})=\sum_{k\in
N_{i}}p_{ki}^{(i)}\phi+\sum_{j\in[m]\backslash\\{i\\}}\sum_{k\in
N_{j}}p_{ki}^{(j)}\bar{q}\phi.$ (1)

When no confusion arises, we simply write $T_{i}(\mathbf{p})$ as $T_{i}$.

The probability of no congestion loss on the direct link $(s_{i},d)$ equals
the probability that there is no arrival during a transmission time $\sigma$,
which is given by
$\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{X}}e^{-T_{i}\sigma}$. As usual
practice, assume $\mathcal{X}$ is an exponential distribution with a rate
parameter $\mu$ (service rate) and mean $\frac{1}{\mu}$. Thus the probability
of no congestion loss on $(s_{i},d)$ is

$\displaystyle\operatorname*{\mathbb{E}}_{\sigma\sim\mathcal{X}}e^{-T_{i}\sigma}=\int_{0}^{+\infty}\mu
e^{-\mu\sigma}e^{-T_{i}\sigma}d\sigma=\frac{\mu}{T_{i}+\mu},$ (2)

and the loss probability on link $(s_{i},d)$ is $\frac{T_{i}}{T_{i}+\mu}$.

Given the strategy profile $\mathbf{p}$, for $s_{i}\in S$ and $k\in N_{i}$,
the loss rate of user $k$ is defined as

$\displaystyle
LR_{k}(\\!\mathbf{p}\\!)\\!\\!=\\!\\!\\!\left[\\!p_{ki}^{(i)}\frac{T_{i}}{T_{i}\\!+\\!\mu}\\!\\!+\\!\\!\left(\\!1\\!-\\!p_{ki}^{(i)}\\!\right)\\!q\\!+\\!(\\!1\\!-\\!q)\\!\\!\\!\\!\\!\\!\\!\\!\sum_{j\in[m]\backslash\\{i\\}}\\!\\!\\!\\!\\!\\!\\!p_{kj}^{(i)}\frac{T_{j}}{T_{j}\\!+\\!\mu}\\!\right]\\!\\!\phi,$
(3)

and the loss probability of user $k$ is $\frac{LR_{k}(\mathbf{p})}{\phi}$.

Total traffic. Regarding the system efficiency, we measure it by the total
traffic rate arriving at the destination $d$. Given the strategy profile
$\mathbf{p}$, the total traffic rate $TR(\mathbf{p})$ of the system can be
derived in two ways. The first expression is derived as the summation of
successful transmission rates on direct links:

$\displaystyle TR(\mathbf{p})=\sum_{i\in[m]}T_{i}\cdot\frac{\mu}{T_{i}+\mu}$
$\displaystyle=\mu\\!\left[\\!m-\sum_{i\in[m]}\frac{\mu}{\sum\limits_{k\in
N_{i}}p_{ki}^{(i)}\phi+\\!\\!\\!\sum\limits_{j\in[m]\backslash\\{i\\}}\sum\limits_{k\in
N_{j}}p_{ki}^{(j)}\bar{q}\phi+\mu}\right]$ (4)

where $T_{i}$ is the traffic rate over link $(s_{i},d)$, and
$\frac{\mu}{T_{i}+\mu}$ is the probability of no congestion loss on
$(s_{i},d)$.

The second expression is from users’ perspective:

$\displaystyle TR(\mathbf{p})$ $\displaystyle:=\sum_{i\in[m]}\sum_{k\in
N_{i}}(\phi-LR_{k}),$ (5)

where $\phi-LR_{k}(\mathbf{p})$ is the traffic rate of user $k\in N_{i}$ that
successfully arrive at $d$. It is not hard to see that (II) and (5) are
equivalent.

Nash equilibria. A Nash equilibrium (NE) is a strategy profile where no player
can decrease its loss probability by unilaterally deviating to any other
strategy. Formally, we give a definition.

######  Definition 1.

A strategy profile $\mathbf{p}$ is a Nash equilibrium, if for any source
$s_{i}\in S$ and any player $k\in N_{i}$, we have

$LR_{k}(\mathbf{p}_{k},\mathbf{p}_{-k})\leq
LR_{k}(\mathbf{p}_{k}^{\prime},\mathbf{p}_{-k}),$

where $\mathbf{p}_{k}^{\prime}$ can be any feasible strategy of player $k$,
and $\mathbf{p}_{-k}$ is the strategy profile of all other players.

We measure the efficiency of NEs by the price of anarchy (PoA) [12], which is
defined as the ratio between social efficiencies in an optimal solution and in
the worst NE. Formally, given an instance $\Gamma$ of this game, we define

$PoA(\Gamma)=\frac{TR(opt)}{\min_{\mathbf{p}\in\mathbb{NE}}TR(\mathbf{p})}.$

where $opt$ is an optimal solution of $\Gamma$, and $\mathbb{NE}$ is the set
of all NEs. The PoA of the whole game is defined as the maximum over all
instances, that is, $PoA=\max_{\Gamma}PoA(\Gamma)$.

##  III Centralized Analysis

The main technical results of the paper are presented now. We show how to
compute an optimal solution that maximizes the total traffic. Note that the
total traffic rate depends on the number of users working on each source by DP
or IP, but not the users’ identity. Given a strategy profile $\mathbf{p}$, let
$u_{i}=|\\{k\in N_{i}~{}|~{}p_{ki}^{(i)}=1\\}|$ be the number of users working
with DP $(s_{i},d)$, and let $v_{i}=|\\{k\in
N_{j},j\in[m]\backslash{i}|~{}p_{ki}^{(j)}=1\\}|$ be the users working with IP
through link $(s_{i},d)$. Define $y_{i}=u_{i}+v_{i}$ as the number of users
who choose source $s_{i}$ (including both DP and IP).

######  Lemma 1.

In any optimal solution, for any source $s_{i}$, either $u_{i}=n_{i}$ or
$v_{i}=0$ or both hold.

###### Proof.

Let $\mathbf{p}$ be an optimal solution. Suppose for contradiction that
$u_{i}<n_{i},v_{i}>0$ for some source $s_{i}$. Then there exists a user (say,
$k$) in $N_{i}$ who chooses IP (say, $(s_{i},s_{i^{\prime}},d)$ for some
$i^{\prime}\neq i$). Also, since $v_{i}>0$, there exist a source $s_{j}\neq
s_{i}$ and a user $l\in N_{j}$ who chooses IP $(s_{j},s_{i},d)$. The total
traffic rate is $TR(\mathbf{p})=\frac{\mu
T_{i}}{T_{i}+\mu}+\sum_{w\in[m]\backslash\\{i\\}}\left(\frac{\mu
T_{w}}{T_{w}+\mu}\right)$.

Now we show that the total traffic rate can be improved by revising
$\mathbf{p}$. Let user $k\in N_{i}$ choose DP, and let user $l\in N_{j}$
choose IP $(s_{j},s_{i^{\prime}},d)$. Fixing all others’ strategies, denote
the new strategy profile by $\mathbf{p}^{\prime}$, and define
$u_{i}^{\prime},v_{i}^{\prime}$ accordingly. Note that
$u_{i}^{\prime}=u_{i}+1,v_{i}^{\prime}=v_{i}-1$, and
$T_{w}(\mathbf{p}^{\prime})=T_{w}(\mathbf{p})$ for all source $s_{w}\neq
s_{i}$. Since $q>0$, we have

$T_{i}(\mathbf{p}^{\prime})=(u_{i}+1)\phi+(v_{i}-1)\bar{q}\phi>u_{i}\phi+v_{i}\bar{q}\phi=T_{i}(\mathbf{p}).$

So $TR(\mathbf{p}^{\prime})>TR(\mathbf{p})$, contradicting to the optimality.
∎

Lemma 1 indicates that if a source (say $s_{i}$) provides service to the users
of other sources, then all users of $s_{i}$ choose DP.

######  Lemma 2.

In any optimal solution, there must exist $\tilde{i}\in[m]$, such that
$v_{l}=0$ for all $l\leq\tilde{i}$, and $u_{j}=n_{j}$ for all $j>\tilde{i}$.

###### Proof.

Given an optimal solution $\mathbf{p}$, suppose for contradiction that there
exist $i,j\in[m]$ ($i<j$) such that $v_{i}>0$ and $u_{j}<n_{j}$. By Lemma 1,
we have $u_{i}=n_{i}$ and $v_{j}=0$. There exists a source $s_{i^{\prime}}$
and a user $k\in N_{i^{\prime}}$ selecting IP $(s_{i^{\prime}},s_{i},d)$.
There exists a source $s_{j^{\prime}}$ ($j^{\prime}\neq j$) and a user
$k^{\prime}\in N_{j}$ selecting IP $(s_{j},s_{j^{\prime}},d)$. Note that when
$i^{\prime}=j$ and $j^{\prime}=i$, users $k$ and $k^{\prime}$ may coincide.
The total traffic rate is

$\displaystyle TR(\mathbf{p})$ $\displaystyle=\frac{\mu
T_{i}}{T_{i}+\mu}+\frac{\mu
T_{j}}{T_{j}+\mu}+\sum_{w\in[m]\backslash\\{i,j\\}}\left(\frac{\mu
T_{w}}{T_{w}+\mu}\right)$
$\displaystyle=\mu(2-\frac{1}{T_{i}+\mu}-\frac{1}{T_{j}+\mu})+\sum_{w\in[m]\backslash\\{i,j\\}}\left(\frac{\mu
T_{w}}{T_{w}+\mu}\right).$

Now we show that the total traffic rate can be improved by revising
$\mathbf{p}$. Let user $k$ choose IP $(s_{i^{\prime}},s_{j^{\prime}},d)$ if
$i^{\prime}\neq j^{\prime}$ and choose DP $(s_{i^{\prime}},d)$ if
$i^{\prime}=j^{\prime}$. Let user $k^{\prime}\in N_{j}$ choose DP. Fixing all
others’ strategies, denote the new strategy profile by $\mathbf{p}^{\prime}$,
and define $u_{i}^{\prime},v_{i}^{\prime}$ accordingly. Note that
$v_{i}^{\prime}\\!=\\!v_{i}-1$, $u_{j}^{\prime}\\!=\\!u_{j}+1$, and
$T_{w}(\mathbf{p}^{\prime})=T_{w}(\mathbf{p})$ for all other sources
$s_{w}\\!\neq\\!s_{i},s_{j}$. Since $i\\!<\\!j$, it follows that
$n_{i}\\!\geq\\!n_{j}\\!>\\!u_{j}$. Therefore, we have

$\displaystyle\frac{1}{T_{i}+\mu}+\frac{1}{T_{j}+\mu}=\frac{1}{n_{i}\phi+v_{i}\bar{q}\phi+\mu}+\frac{1}{u_{j}\phi+\mu}$
$\displaystyle>~{}$
$\displaystyle\frac{1}{n_{i}\phi+v_{i}^{\prime}\bar{q}\phi+\mu}+\frac{1}{u_{j}^{\prime}\phi+\mu}=\frac{1}{T_{i}^{\prime}+\mu}+\frac{1}{T_{j}^{\prime}+\mu},$

which indicates that $TR(\mathbf{p})<TR(\mathbf{p}^{\prime})$, a
contradiction. ∎

Lemma 2 shows that there exists a threshold $\tilde{i}$: 1) if $i>\tilde{i}$,
all users from $s_{i}$ chose DP; 2) if $i\leq\tilde{i}$, portion users chose
DP, and portion users chose IP.

Now we are ready to present Algorithm 1. The main idea is searching for
$\tilde{i}$ in Lemma 2. For each candidate of $\tilde{i}$, let $B$ be the
number of users selecting IP, all of whom come from
$L\\!=\\!\\{s_{l}~{}|~{}l\leq\tilde{i}\\}$, and go to
$R\\!=\\!\\{s_{j}~{}|~{}j>\tilde{i}\\}$. For every possible value of $B$, we
compute the best possible way for extracting the $B$ users from $L$ and
distributing them over $R$.

0: $m$ source with $n_{1}\geq n_{2}\geq\cdots\geq n_{m}$, $\phi,\mu,q$

0: $(u^{*}_{i},v^{*}_{i})_{i\in N}$

1: Initialize $u_{i}=v_{i}=0$ for all $i\in[m]$. $TR^{*}=0$

2: for $\tilde{i}=1,\ldots,m$ do

3: Let $v_{l}=0$ for all $l=1,2,\ldots,\tilde{i}$

4: Let $u_{j}=n_{j}$ for all $j=\tilde{i}+1,\ldots,m$

5: for $B=1,2,\ldots,\sum_{l=1}^{\tilde{i}}n_{l}$ do

6: (a) Compute $(u_{l})_{l\in[\tilde{i}]}$ s.t.
$\sum_{l\\!=\\!1}^{\tilde{i}}(n_{l}\\!-\\!u_{l})\\!=\\!B$

and the values of $u_{l}$ are as equal as possible.

7: (b) Compute $(v_{j})_{j>\tilde{i}}$ s.t.
$\sum_{j=\tilde{i}+1}^{m}v_{j}\\!=\\!B$ and the values of
$n_{j}\phi\\!+\\!v_{j}\bar{q}\phi\\!+\\!\mu$ are as equal as possible.

8: (c) Compute $TR$ with respect to $(u_{i},v_{i})_{i\in N}$

9: if $TR>TR^{*}$ then

10: $TR^{*}\leftarrow TR$, and $(u^{*}_{i},v^{*}_{i})_{i\in
N}\leftarrow(u_{i},v_{i})_{i\in N}$

11: end if

12: end for

13: end for

ALGORITHM 1 Computing an optimal solution.

In Algorithm 1, step (a) is to make $T_{l}$ (and thus no congestion
probability $\frac{\mu}{T_{l}+\mu}$) as equal as possible for
$l\\!\in\\![\tilde{i}]$. This can be realized by initializing
$u_{l}\\!=\\!n_{l}$, and then removing players one by one from the highest
$u_{l}$ and updating until $B$ players have been removed. The goal of step (b)
is to make $T_{j}$ (and thus $\frac{\mu}{T_{j}+\mu}$) as even as possible.
This can be realized by initializing $v_{j}=0$, then adding users one by one
to $v_{j^{\prime}}=\arg\min_{j>\tilde{i}}n_{j}\phi+v_{j}\bar{q}\phi+\mu$ and
updating, until $B$ players have been added. These two steps guarantee that
the $B$ loads are distributed in an optimal way to maximize the traffic rate.

######  Theorem 1.

Algorithm 1 returns an optimal solution for maximizing the total traffic, and
runs in $O(mn^{2})$ time.

###### Proof.

In the first loop, we traverse all indexes in $[m]$ to find the $\tilde{i}$ in
Lemma 2. In the second loop, we traverse all possible numbers of users who
select IP, and given any such a number $B$, we extract the $B$ users from
$\\{s_{l}~{}|~{}l\leq\tilde{i}\\}$ and distribute them over
$\\{s_{j}~{}|~{}j>\tilde{i}\\}$ in an optimal way to maximize the traffic
rate. So all possible optimal solutions have been searched by the algorithm,
giving the optimality. For the time complexity, we have $m$ iterations in the
first loop, at most $n$ iterations in the second loop, and the time for each
iteration is $O(n)$. ∎

Intuitively, when the transmission loss probability is sufficiently large, all
packets should go through DP; when there is no transmission loss, the load of
packets should be distributed evenly over all sources. We verify the intuition
as follows.

######  Corollary 1.

If $q=1$, the unique optimal solution is that all users choose DP (i.e.,
$u_{i}=n_{i},v_{i}=0,\forall i\in M$). If $q=0$, a strategy profile
$\mathbf{p}$ is optimal if and only if $|y_{i}-y_{j}|\leq 1$ for all
$i,j\in[m]$.

###### Proof.

If $q=1$, $TR=\sum_{i\in[m]}\frac{\mu u_{i}\phi}{u_{i}\phi+\mu}$ is increasing
with respect to every $u_{i}$. By the monotonicity, the optimum is achieved
when $u_{i}=n_{i}$. If $q=0$, suppose for contradiction that there exist
$i,j\in[m]$ in an optimal solution $\mathbf{p}$ such that $y_{i}-y_{j}\geq 2$.
The total traffic rate is $TR=\mu(m-\sum_{k\in[m]}\frac{\mu}{y_{k}\phi+\mu})$.
Consider a new strategy profile $\mathbf{p}^{\prime}$ with
$y_{i}^{\prime}=y_{i}-1,y_{j}^{\prime}=y_{j}+1$, i.e., a user who chooses
source $s_{i}$ deviates to $s_{j}$. Then the total traffic rate becomes
$TR^{\prime}=\mu(m-\sum_{k\in[m]\backslash\\{i,j\\}}\frac{\mu}{y_{k}\phi+\mu}-\frac{\mu}{y_{i}^{\prime}\phi+\mu}-\frac{\mu}{y_{j}^{\prime}\phi+\mu})>TR$,
a contradiction. ∎

##  IV Decentralized Analysis

In this section, we study the Nash equilibria in the decentralized decision-
making scenario where each user makes a decision on the choice of DP or IP.

###  IV-A Characterization of NEs

A NE should satisfy that: for a user selecting DP, its loss rate will not
decrease if it deviates to any IP; for a user selecting IP, its loss rate will
not decrease if it deviates to DP or another IP. We formalize it as the
following characterization.

######  Theorem 2.

Given an arbitrary strategy profile $\mathbf{p}$ with
$(u_{i},v_{i})_{i\in[m]}$, let
$i^{*}\in\arg\min_{i\in[m]}\\{u_{i}+v_{i}\bar{q}\\}$, and let
$x_{ij}\in\\{0,1\\}$ be an indicator where $x_{ij}=1$ if there exists at least
one user selecting IP $(s_{i},s_{j},d)$. Then, $\mathbf{p}$ is a NE, if and
only if the following conditions are satisfied:

  * (i)

for all $i\in[m]$ with $u_{i}>0$, we have

$\bar{q}(u_{i}+v_{i}\bar{q})\leq
u_{i^{*}}+v_{i^{*}}\bar{q}+\bar{q}+\frac{q\mu}{\phi};$ (6)

  * (ii)

for all $i,l\in[m]$ with $x_{il}=1$, we have

$\displaystyle
u_{l}\\!+\\!v_{l}\bar{q}\leq\min\\!\left\\{\\!\bar{q}(\\!u_{i}\\!+\\!1\\!+\\!v_{i}\bar{q})\\!-\\!\frac{q\mu}{\phi},u_{i^{*}}\\!+\\!v_{i^{*}}\bar{q}\\!+\\!\bar{q}\\!\right\\}$
(7)

###### Proof.

Suppose $\mathbf{p}$ is a NE. Consider an arbitrary source $s_{i}\in S$ and
arbitrary user $k\in N_{i}$.

Case 1. User $k$ selects DP in $\mathbf{p}$ (denoted as $\mathbf{p}_{k}^{(i)}$
where $p_{ki}^{(i)}=1$). If it deviates to IP $(s_{i},s_{j},d)$ where $j\neq
i$ (denoted as $\mathbf{p^{\prime}}_{k}^{(i)}$ where $p_{kj}^{i}=1$), by
Definition 1, we have $LR_{k}(\mathbf{p}_{k}^{(i)},\mathbf{p}_{-k})\leq
LR_{k}(\mathbf{p^{\prime}}_{k}^{(i)},\mathbf{p}_{-k})$. It is equivalent to

$\displaystyle
1\\!-\\!\frac{\mu}{u_{i}\phi+v_{i}\bar{q}\phi+\mu}\\!\leq\\!q+\bar{q}\cdot\left(1\\!-\\!\frac{\mu}{u_{j}\phi+(v_{j}+1)\bar{q}\phi+\mu}\right)$
$\displaystyle\Leftrightarrow$
$\displaystyle~{}\frac{\bar{q}}{u_{j}\phi+(v_{j}+1)\bar{q}\phi+\mu}\leq\frac{1}{u_{i}\phi+v_{i}\bar{q}\phi+\mu}$
$\displaystyle\Leftrightarrow$
$\displaystyle~{}\bar{q}(u_{i}+v_{i}\bar{q})-\frac{q\mu}{\phi}\leq
u_{j}+(v_{j}+1)\bar{q},$

The above inequality should hold for all $j\neq i$, and thus is equivalent to
Equation (6).

Case 2. User $k$ selects IP $(s_{i},s_{l},d)$ in $\mathbf{p}$. If it deviates
to DP $(s_{i},d)$, by Definition 1, we should have
$LR_{k}(\mathbf{p}_{k}^{(i)},\mathbf{p}_{-k})\leq
LR_{k}(\mathbf{p^{\prime}}_{k}^{(i)},\mathbf{p}_{-k})$. It is equivalent to

$\displaystyle~{}q\\!+\\!\bar{q}\cdot\left(1\\!-\\!\frac{\mu}{u_{l}\phi+v_{l}\bar{q}\phi+\mu}\right)\leq
1\\!-\\!\frac{\mu}{(u_{i}+1)\phi+v_{i}\bar{q}\phi+\mu}$
$\displaystyle\Leftrightarrow$
$\displaystyle~{}\frac{1}{(u_{i}+1)\phi+v_{i}\bar{q}\phi+\mu}\leq\frac{\bar{q}}{u_{l}\phi+v_{l}\bar{q}\phi+\mu}$
$\displaystyle\Leftrightarrow$
$\displaystyle~{}u_{l}+v_{l}\bar{q}\leq\bar{q}(u_{i}+1+v_{i}\bar{q})-\frac{q\mu}{\phi}.$

Moreover, a NE must guarantee that user $k$ will not deviate to another IP
$(s_{i},s_{j},0)$, and thus we should have

$\displaystyle~{}1-\frac{\mu}{u_{l}\phi+v_{l}\bar{q}\phi+\mu}\leq
1-\frac{\mu}{u_{j}\phi+(v_{j}+1)\bar{q}\phi+\mu}$
$\displaystyle\Leftrightarrow$ $\displaystyle~{}u_{l}+v_{l}\bar{q}\leq
u_{j}+(v_{j}+1)\bar{q}.$

Note that the above inequality should hold for all $j\neq i,l$. Therefore, we
obtain Equation (7). ∎

###  IV-B Price of Anarchy

We investigate the price of anarchy in this section, which measures the
efficiency of NE. We give an upper bound on the optimal total traffic rate,
and a lower bound on the total traffic rate of any NE.

######  Lemma 3.

In an optimal solution $\mathbf{p}$, the total traffic rate is
$TR_{\text{}}(\mathbf{p})\leq\mu m\left(1-\frac{\mu}{n_{1}\phi+\mu}\right)$.

###### Proof.

Let $i$ be the index stated in Lemma 2. It suffices to show with the proof by
contradiction that in the optimal solution $\mathbf{p}$,
$u_{i}+v_{i}\bar{q}\leq n_{i}$. First, for $i=1$, we have $v_{1}\\!=\\!0$, and
thus it satisfies $u_{1}+v_{1}\bar{q}\\!=\\!u_{1}\leq n_{1}$. For any $i>1$,
suppose for contradiction that $u_{i}+v_{i}\bar{q}>n_{i}$. Then $v_{i}>0$, and
there exists a source $s_{j}$ and a user $k\in N_{j}$ that chooses the IP
$(s_{j},s_{i},d)$, i.e., $u_{j}<n_{j}$. By Lemma 1, it must be $v_{j}=0$, and
thus $T_{j}=u_{j}\phi$. Denote the strategy as $\mathbf{p}$ with
$p_{ki}^{(j)}=1$. The total traffic rate is

$\displaystyle
TR(\mathbf{p})=\mu\left(m-\frac{\mu}{T_{j}+\mu}-\frac{\mu}{T_{i}+\mu}-\sum_{w\in[m]\backslash\\{j\\}}\frac{\mu}{T_{w}+\mu}\right)$
(8)

We show that the total traffic rate can be improved with user $k\in N_{j}$
deviating from IP $(s_{j},s_{i},d)$ to DP $(s_{j},d)$. Fixing the strategies
of all others, denote by $\mathbf{p}^{\prime}$ the new strategy profile, and
define $(u^{\prime}_{w},v^{\prime}_{w},T^{\prime}_{w})_{w\in[m]}$ accordingly.
Note that $u_{j}^{\prime}=u_{j}+1$, $v_{i}^{\prime}=v_{i}-1$,
$u_{i}^{\prime}=u_{i}$, and $T_{w}^{\prime}=T_{w}$ for any
$w\in[m]\backslash\\{j\\}$. Since $u_{i}+v_{i}\bar{q}>n_{1}\geq n_{j}\geq
u_{j}$, we have

$\displaystyle\frac{1}{T_{j}+\mu}+\frac{1}{T_{i}+\mu}=\frac{1}{u_{j}\phi+\mu}+\frac{1}{u_{i}\phi+v_{i}\bar{q}\phi+\mu}$
$\displaystyle>$
$\displaystyle~{}\frac{1}{u_{j}^{\prime}\phi+\mu}+\frac{1}{u_{i}^{\prime}\phi+v_{i}^{\prime}\bar{q}\phi+\mu}=\frac{1}{T_{j}^{\prime}+\mu}+\frac{1}{T_{i}^{\prime}+\mu}.$

It indicates that $TR(\mathbf{p}^{\prime})\\!>\\!TR(\mathbf{p})$ is a
contradiction. Consequently,
$u_{i}+v_{i}\bar{q}\\!\leq\\!n_{i}\\!\leq\\!n_{1}$. According (II), we have
$TR_{\text{}}(\mathbf{p})\leq\mu m(1-\frac{\mu}{n_{1}\phi+\mu})$. ∎

######  Lemma 4.

Let $z=\min\\{n_{m},\frac{n}{4m}-\bar{q}-\frac{q\mu}{\phi}\\}$. For every NE
$\mathbf{p}$, the total traffic rate satisfies
$TR(\mathbf{p})\geq\mu(m-\frac{m\mu}{z\phi+\mu})$.

###### Proof.

Let $i^{*}=\arg\min_{i\in[m]}\\{u_{i}+v_{i}\bar{q}\\}$. Since
$TR(\mathbf{p})\geq\mu m(1-\frac{\mu}{(u_{i^{*}}+v_{i^{*}}\bar{q})\phi+\mu})$,
it suffices to prove that $u_{i^{*}}+v_{i^{*}}\bar{q}\geq z$. If
$u_{i^{*}}+v_{i^{*}}\bar{q}\geq n_{m}$, it is done. We only need to consider
the case when $u_{i^{*}}+v_{i^{*}}\bar{q}<n_{m}\leq n_{i^{*}}$. There exists
some users in $N_{i^{*}}$ selecting IP. By Equation (7), we have

$u_{i^{*}}+v_{i^{*}}\bar{q}\leq\bar{q}(u_{i^{*}}+1+v_{i^{*}}\bar{q})-\frac{q\mu}{\phi}.$

By Theorem 2, for each $i\in[m]$, if $u_{i}>0$, then
$\bar{q}(u_{i}+v_{i}\bar{q})\leq
u_{i^{*}}+v_{i^{*}}\bar{q}+\bar{q}+\frac{q\mu}{\phi}$; if $v_{i}>0$, then
$u_{i}+v_{i}\bar{q}\leq u_{i^{*}}+v_{i^{*}}\bar{q}+\bar{q}$. In both cases, we
obtain $u_{i}+v_{i}/2\leq
2(u_{i^{*}}+v_{i^{*}}\bar{q}+\bar{q}+\frac{q\mu}{\phi})$. Summing up over all
$i\in[m]$, we have

$n/2\leq\sum_{i\in[m]}(u_{i}+v_{i}/2)\leq
2m(u_{i^{*}}+v_{i^{*}}\bar{q}+\bar{q}+\frac{q\mu}{\phi}),$

which implies that
$u_{i^{*}}+v_{i^{*}}\bar{q}\geq\frac{n}{4m}-\bar{q}-\frac{q\mu}{\phi}$. ∎

######  Theorem 3.

For any instance with $m$ sources, the price of anarchy is $PoA\leq
1+\frac{n_{1}\mu}{n_{1}z\phi+z\mu}$, where
$z=min\\{n_{m},\frac{n}{4m}-\bar{q}-\frac{q\mu}{\phi}\\}$.

###### Proof.

Combining the upper bound in Lemma 3 and the lower bound on $TR(\mathbf{p})$
for any NE $\mathbf{p}$ in Lemma 4, it follows

$PoA\leq\frac{m-\frac{m\mu}{n_{1}\phi+\mu}}{m-\frac{m\mu}{z\phi+\mu}}=\frac{n_{1}(z\phi+\mu)}{z(n_{1}\phi+\mu)}=1+\frac{n_{1}\mu}{n_{1}z\phi+z\mu}.$

∎

##  V A Particular Case: Two Sources

In this section, we focus on the special case of $m=2$. That is, there are
only two sources $s_{1}$ and $s_{2}$. Assume w.l.o.g. that $n_{1}\geq n_{2}$.
For each user $k\in N_{i}$, there is only one IP. Accordingly, its strategy
becomes $\mathbf{p}_{k}^{(i)}=\left(p_{k1}^{(i)},p_{k2}^{(i)}\right),i=1,2$.
And we have

$\displaystyle n_{1}=u_{1}+v_{2};$ (9) $\displaystyle n_{2}=u_{2}+v_{1}.$ (10)

The traffic rate $T_{i}(\mathbf{p})$ in (1) is rephrased as

$\displaystyle T_{i}(\mathbf{p})=\sum_{k\in N_{i}}p_{ki}^{(i)}\phi+\sum_{k\in
N_{j},j\neq i}p_{ki}^{(j)}\bar{q}\phi=u_{i}\phi+v_{i}\bar{q}\phi.$ (11)

Given strategy profile $\mathbf{p}$, the set $N$ is further partitioned into 4
subsets $(V_{1},V_{2},V_{3},V_{4})$ where $V_{1}=\\{k\in
N_{1}~{}|~{}\mathbf{p}_{k}^{(1)}=(1,0)\\}$, $V_{2}=\\{k\in
N_{1}~{}|~{}\mathbf{p}_{k}^{(1)}=(0,1)\\}$, $V_{3}=\\{k\in
N_{2}~{}|~{}\mathbf{p}_{k}^{(2)}=(0,1)\\}$ and $V_{4}=\\{k\in
N_{2}~{}|~{}\mathbf{p}_{k}^{(2)}=(1,0)\\}$. Clearly, users in $V_{1}$ and
$V_{3}$ choose DP, and users in $V_{2}$ and $V_{4}$ choose IP.

Suppose $\mathbf{p}$ is a NE. We study the deviation of users in
$V_{1},V_{2},V_{3},V_{4}$, respectively. For user $k\in V_{1}$, the strategy
is $\mathbf{p}_{k}^{(1)}=\left(p_{k1}^{(1)},p_{k2}^{(1)}\right)=(1,0)$, and
the loss rate in (3) is

$\displaystyle LR_{k}(\mathbf{p})$ $\displaystyle\\!=\\!\frac{\phi
T_{1}(\mathbf{p})}{T_{1}(\mathbf{p})+\mu}=\left[1-\frac{\mu/\phi}{u_{1}+v_{1}\bar{q}+\mu/\phi}\right]\phi.$

When user $k\in V_{1}$ deviates to IP, the strategy profile becomes
$\mathbf{p^{\prime}}=\left(\mathbf{p^{\prime}}_{k}^{(1)},\mathbf{p}_{-k}\right)$
where $\mathbf{p}_{k}^{\prime(1)}=(0,1)$. The loss rate of user $k$ becomes

$\displaystyle
LR_{i}(\mathbf{p}^{\prime})\\!=\\!q\phi\\!+\\!\bar{q}\phi\frac{T_{2}(\mathbf{p}^{\prime})}{T_{2}(\mathbf{p}^{\prime})\\!+\\!\mu}\\!=\\!\left[1\\!-\\!\frac{\bar{q}\mu/\phi}{u_{2}\\!+\\!(v_{2}\\!+\\!1)\bar{q}\\!+\\!\mu/\phi}\right]\phi.$

Since $\mathbf{p}$ is NE, $k$ has no incentive to deviate, and thus
$LP_{k}(\mathbf{p})\leq LR_{k}(\mathbf{p}^{\prime})$, which is equivalent to

$\displaystyle
t_{1}(u_{2}):=\frac{q\mu/\phi\\!+\\!u_{2}(1\\!+\\!\bar{q}^{2})\\!+\\!(n_{1}\\!+\\!1)\bar{q}\\!-\\!n_{2}\bar{q}^{2}}{2\bar{q}}\geq
u_{1},$ (12)

where $t_{1}(u_{2})$ is a function with respect to variable $u_{2}$.

For user $k\in V_{2}$ with strategy $\mathbf{p}_{k}^{(1)}=(0,1)$, the loss
rate is

$\displaystyle
LR_{k}(\mathbf{p})=q\phi+\bar{q}\phi\frac{T_{2}(\mathbf{p})}{T_{2}(\mathbf{p})+\mu}=\left[1-\frac{\bar{q}\mu/\phi}{u_{2}+v_{2}\bar{q}+\mu/\phi}\right]\phi.$

When user $k\in V_{2}$ deviates to DP, the strategy profile becomes
$\mathbf{p^{\prime}}=\left(\mathbf{p^{\prime}}_{k}^{(1)},\mathbf{p}_{-k}\right)$
where $\mathbf{p}_{k}^{\prime(1)}=(1,0)$. The loss rate of $k$ becomes

$\displaystyle LR_{k}(\mathbf{p}^{\prime})=\frac{\phi
T_{1}(\mathbf{p}^{\prime})}{T_{1}(\mathbf{p}^{\prime})+\mu}=\left[1-\frac{\mu/\phi}{u_{1}+1+v_{1}\bar{q}+\mu/\phi}\right]\phi.$

Since $\mathbf{p}$ is NE, we have $LR_{k}(\mathbf{p})\leq
LR_{k}(\mathbf{p}^{\prime})$, that is,

$u_{1}\geq\frac{q\mu/\phi\\!+\\!u_{2}(1\\!+\\!\bar{q}^{2})\\!+\\!(n_{1}\\!-\\!1)\bar{q}\\!-\\!n_{2}\bar{q}^{2}}{2\bar{q}}=t_{1}(u_{2})\\!-\\!1.$
(13)

Symmetrically, for each user $k\in V_{3}$ and $k\in V_{4}$, since $\mathbf{p}$
is NE, we have

$\displaystyle t_{2}(u_{1})-1\leq u_{2}\leq t_{2}(u_{1}),$ (14)

where

$\displaystyle
t_{2}(u_{1}):=\frac{q\mu/\phi+u_{1}(1+\bar{q}^{2})+(n_{2}+1)\bar{q}-n_{1}\bar{q}^{2}}{2\bar{q}}.$

Note that Eq. (12) - (14) are the sufficient and necessary conditions for an
arbitrary strategy $\mathbf{p}$ to achieve NE. Now we are ready to give a
characterization of NEs.

######  Theorem 4.

Let $\mathbf{p}$ be an arbitrary strategy profile for the game with two
sources. Let $u_{1}$ and $u_{2}$ be the number of users in $N_{1}$ and $N_{2}$
who choose DP under $\mathbf{p}$, respectively. We have

  * [1]

when (a) $u_{1}=n_{1},u_{2}<n_{2}$, or (b) $u_{1}=0,u_{2}>0$, $\mathbf{p}$
cannot be a NE;

  * [2]

when $u_{1}\in[0,n_{1}),u_{2}\in[0,n_{2})$, $\mathbf{p}$ is NE if and only if
$u_{1}\geq t_{1}(u_{2})-1$ and $u_{2}\geq t_{2}(u_{1})-1$;

  * [3]

when $u_{1}\in(0,n_{1}),u_{2}=n_{2}$, $\mathbf{p}$ is NE if and only if
$u_{1}\in[t_{1}(u_{2})-1,t_{1}(u_{2})]$;

  * [4]

when $u_{1}=n_{1},u_{2}=n_{2}$, $\mathbf{p}$ is NE if and only if
$n_{1}\bar{q}\leq q\mu/\phi+n_{2}+\bar{q}$.

###### Proof.

Given $\mathbf{p}$, let $(V_{1},V_{2},V_{3},V_{4})$ be a partition of $N$ as
defined above. We discuss the four cases.

Case 1. When (a) $u_{1}=n_{1}$ and $u_{2}<n_{2}$, $V_{4}$ is nonempty. If
$\mathbf{p}$ is a NE, it must satisfy $t_{2}(u_{1})-1\leq u_{2}$. However,
because $\bar{q}u_{2}\leq n_{1},\bar{q}u_{2}\leq(n_{2}-1)\bar{q}$ and
$q\mu/\phi>0$, it cannot hold.

When (b) $u_{1}=0$ and $u_{2}>0$, $V_{2}$ is nonempty. If $\mathbf{p}$ is a
NE, it must satisfy Eq. (13), that is, $u_{1}\geq t_{1}(u_{2})-1$. It follows
that $0=2\bar{q}u_{1}\geq
q\mu/\phi+u_{2}(1+\bar{q}^{2})+(n_{1}-1)\bar{q}-n_{2}\bar{q}^{2}\geq
q\mu/\phi+1+(n_{1}-1)\bar{q}-(n_{2}-1)\bar{q}^{2}\geq q\mu/\phi+1>0$, a
contradiction.

Case 2. When $u_{1}\in[0,n_{1}),u_{2}\in[0,n_{2})$, $V_{2}$ and $V_{4}$ are
nonempty. It is easy to see that $\mathbf{p}$ is NE if and only if $u_{1}\geq
t_{1}(u_{2})-1$ and $u_{2}\geq t_{2}(u_{1})-1$ are satisfied simultaneously.

Case 3. When $u_{2}=n_{2},u_{1}\in(0,n_{1})$, $V_{1},V_{2},V_{3}$ are
nonempty, and $V_{4}$ is empty. $\mathbf{p}$ is NE if and only if
$t_{1}(u_{2})-1\leq u_{1}\leq t_{1}(u_{2})$ and $u_{2}\leq t_{2}(u_{1})$ hold
simultaneously. Moreover, note that $u_{2}\leq t_{2}(u_{1})$ is implied by
$u_{1}\geq t_{1}(u_{2})-1$. Therefore, the sufficient and necessary condition
for NE is $u_{1}\in[t_{1}(u_{2})-1,t_{1}(u_{2})]$.

Case 4. When $u_{1}=n_{1},u_{2}=n_{2}$, $V_{1},V_{3}$ are nonempty, and
$V_{2},V_{4}$ are empty. $\mathbf{p}$ is NE if and only if $u_{1}\leq
t_{1}(u_{2})$ and $u_{2}\leq t_{2}(u_{1})$. It is easy to see that, it is
equivalent to $n_{1}\bar{q}\leq q\mu/\phi+n_{2}+\bar{q}$.

Note that every situation of $u_{1},u_{2}$ is included in the above four
cases. So we complete a characterization. ∎

Case 4 can be intuitively explained by considering the sidelink loss
probability $q$ over link $(s_{1},s_{2})$. If $q$ is sufficiently high, no
user would prefer the indirect path, and selecting the direct path would be a
NE for all users. Conversely, when there is no transmission loss over sidelink
$(s_{1},s_{2})$ (i.e., $q=0$), every user would prefer to use the source with
fewer users. Therefore, the profile of all users selecting DP is a NE only if
the user distribution between the two sources is as even as possible, with
$n_{1}\leq n_{2}+1$. Based on Theorem 4, we give some interesting conclusions.

######  Corollary 2.

If a strategy profile with $u_{1}=n_{1},u_{2}=n_{2}$ is optimal, then it is
also a NE.

######  Corollary 3.

A strategy profile with $u_{1}=0,u_{2}=0$ is a NE, if and only if (a)
$n_{1}=n_{2}+1,\bar{q}=1$, or (b)
$n_{1}=n_{2},n_{1}(1-\bar{q}^{2})\leq\bar{q}-\frac{q\mu}{\phi}$.

Note that $u_{1}\geq t_{1}(u_{2})-1$ and $u_{2}\geq t_{2}(u_{1})-1$ cannot
hold simultaneously when $q>\frac{2}{n}$, and $u_{1}\geq t_{1}(n_{2})-1$
cannot hold when $n_{1}\bar{q}<q\mu/\phi+n_{2}+\bar{q}$.

######  Corollary 4.

When $n_{1}\bar{q}<q\mu/\phi+n_{2}+\bar{q}$ and $q>\frac{2}{n}$, the unique NE
is that all users choose DP, i.e., $u_{1}=n_{1},u_{2}=n_{2}$.

We end this section by proving the existence of NE.

######  Theorem 5.

For any game instance with two sources, there exists a NE with $u_{1}>0$ and
$u_{2}=n_{2}$.

###### Proof.

By Theorem 4 (4), if $n_{1}\bar{q}\leq q\mu/\phi+n_{2}+\bar{q}$, then the
strategy profile that all users choose DP (i.e., $u_{1}=n_{1},u_{2}=n_{2}$) is
a NE. Otherwise, $n_{1}\bar{q}>q\mu/\phi+n_{2}+\bar{q}$. Let $\tilde{m}$ be an
integer in interval
$[\frac{q\mu/\phi+n_{2}+n_{1}\bar{q}-\bar{q}}{2\bar{q}},\frac{q\mu/\phi+n_{2}+n_{1}\bar{q}+\bar{q}}{2\bar{q}}]=[t_{1}(n_{2})-1,t_{1}(n_{2})]$,
which always admits at least one integer. Note that
$n_{1}>\frac{q\mu/\phi+n_{2}+n_{1}\bar{q}+\bar{q}}{2\bar{q}}\geq\tilde{m}>0$.
By Theorem 4, a strategy profile with $u_{1}=\tilde{m}$ and $u_{2}=n_{2}$ is a
NE. ∎

##  VI Numerical experiments

(a) $q\\!\in\\![0,1]$ with $n_{1}\\!=\\!1000$, $n_{2}\\!=\\!100$,
$\phi\\!=\\!1$, $\mu=300$.

(b) $\mu\\!\in\\![1,6000]$ with $n_{1}\\!=\\!1000$, $n_{2}\\!=\\!100$,
$\phi\\!=\\!1$, $q\\!=\\!0.3$.

(c) $n_{1}\\!\in\\![500,8000]$, $n_{1}$ is an integer with $n_{2}\\!=\\!100$,
$\phi\\!=\\!1$, $\phi=1$, $q\\!=\\!0.7$.

Figure 2:  The PoA and total traffics in $opt$ and NE in two-source network

Through numerical simulations, we explore the impact of traffic condition on
network performance, i.e., the total traffic rate and PoA. Recall that the
traffic flow originating from each user is Poisson with rate $\phi$, the
service rate of each direct link is $\mu$, and the loss probability over each
side link is $q$. Assume $\phi=1$ for normalization.

We first present the simulation results in two-source networks. In Fig. 2, the
PoA and the total traffics are plotted under different $q$, $\mu$ and $n_{1}$,
showing a PoA of less than 1.08. Such a little gap between the optimal
solution and the worst NE suggests that the gain of centralized-decision
making over decentralized-decision making is trivial most of the time. As
shown in Fig. 2(a), the total traffic decreases with the increase of $q$,
i.e., the increased loss rate on sidelink. On the other hand, the PoA is first
increasing from $1$ at $q=0$, implying that the NE and optimal solution are
the same with $u_{1}=n_{2}+v_{1}$. That is, we have the equal number of users
on $(s_{1},d)$ and $(s_{2},d)$ in terms of both IP and DP. With the increase
of $q$, the benefit of centralized-decision making is gradually unveiled.
However, when $p$ reaches a certain value, the PoA goes down to $1$ quickly.
An intuitive explanation is that, when $q$ becomes larger than the loss
probability on DP, no users will choose IP in NE. And this strategy is optimal
as well.

In Fig. 2(b), the traffic rates grow with the increase of $\mu$ due to the
increased probability of no congestion in (2). This is, a high service rate
help clear the collision and relieve congestion on both DP and IP. The PoA
curve indicates that either in overloaded or less congested scenarios, there
is little improvement of centralized-decision making. In Fig. 2(c), the
increased number of users leads to an increase of traffic rate in spite of the
rise in loss rate. What is more, PoA tends towards $1$ for small and large
$n_{1}$. As the strategies in $opt$ and NE are much similar for users at
source $s_{1}$, i.e., DP in less biased scenario and IP severely biased
scenario.

Figure 3: The impact of $n_{1}$ on total traffics in multi-source network.
Figure 4: Comparison of NE strategy and optimal strategy in terms of traffic
rate.

In the multi-source network, while the optimal solution can be easily computed
by Algorithm 1, it is difficult to find all NEs even given Theorem 2. Hence,
we merely consider a small value of $m$ and $n$ (i.e., $m=3$). The service
rate and traffic arrival rate are fixed as $\mu=1,\phi=1$. Results are given
in Fig. 3, which shows similar results in Fig, 2(c). It is obvious that the
growth of the total traffic slows down gradually, because given the service
rate, an increase of $n_{1}$ aggravates the network congestion. Second, the
increase of loss rate on sidelink, leads to the increase of loss rate on IP.
As a result, more users choose DP instead, which in turn worsens the network
congestion.

Figure 4 plots the performances for a range of $q$. When $q=0$ and $q=1$, the
PoA is exactly 1. The PoA converges to 1 when $q$ goes to 1, because when the
service rate is large enough compared with arrival rate, there is a
sufficiently small congestion loss and all users like to choose DP.

##  VII Conclusion

In this work, we give a theoretical analysis of a load balancing game in
cloud-enabled networks, in which the users want to minimize the loss
probability of their packets with suitable routing strategies. In the
centralized analysis, an efficient algorithm for maximizing the total traffic
rate is proposed, according to Lemma 1 and Lemma 2. In the decentralized
analysis, a characterization of Nash equilibrium is given, and the PoA is
investigated. Numerical experiments show that the efficiency loss due to
selfish behaviors is relatively small in most cases.

There are many future directions that are worth exploring. First, we only
focus on pure strategies of players in this work, and an immediate and natural
question is how the users act when mixed strategies are allowed. Second, it
would be interesting to investigate heterogeneous servers (source nodes) where
each $s_{i}$ serves a different purpose or has a different service rate
$\mu_{i}$. Moreover, while we only consider direct path and one-hop indirect
paths, a more general scenario where players can choose multi-hop indirect
paths to the destination can be taken into consideration.

## References

  * [1] A. K. Erlang, “Solution of some problems in the theory of probabilities of significance in automatic telephone exchanges,” _Post Office Electrical Engineer’s Journal_ , vol. 10, pp. 189–197, 1917.
  * [2] K. Jung, Y. Lu, D. Shah, M. Sharma, and M. S. Squillante, “Revisiting stochastic loss networks: Structures and approximations,” _Mathematics of Operations Research_ , vol. 44, no. 3, pp. 890–918, 2019.
  * [3] J. Cao, K. Li, and I. Stojmenovic, “Optimal power allocation and load distribution for multiple heterogeneous multicore server processors across clouds and data centers,” _IEEE Transactions on Computers_ , vol. 63, no. 1, pp. 45–58, 2013.
  * [4] E. J. Ghomi, A. M. Rahmani, and N. N. Qader, “Load-balancing algorithms in cloud computing: A survey,” _Journal of Network and Computer Applications_ , vol. 88, pp. 50–71, 2017.
  * [5] J. Zhao, K. Yang, X. Wei, Y. Ding, L. Hu, and G. Xu, “A heuristic clustering-based task deployment approach for load balancing using bayes theorem in cloud environment,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 27, no. 2, pp. 305–316, 2015.
  * [6] G. Scutari, D. P. Palomar, F. Facchinei, and J.-S. Pang, “Convex optimization, game theory, and variational inequality theory,” _IEEE Signal Processing Magazine_ , vol. 27, no. 3, pp. 35–49, 2010.
  * [7] S. Mondal, G. Das, and E. Wong, “A game-theoretic approach for non-cooperative load balancing among competing cloudlets,” _IEEE Open Journal of the Communications Society_ , vol. 1, pp. 226–241, 2020.
  * [8] C. Yi, J. Cai, K. Zhu, and R. Wang, “A queueing game based management framework for fog computing with strategic computing speed control,” _IEEE Transactions on Mobile Computing_ , vol. 21, no. 5, pp. 1537–1551, 2020\.
  * [9] E. Altman, J. Kuri, and R. El-Azouzi, “A routing game in networks with lossy links,” in _the 7th International Conference on Network Games, Control and Optimization_. IEEE, 2014, pp. 91–94.
  * [10] B. Toure, S. Paturel, and E. Altman, “Congestion load balancing game with losses,” in _the 8th International Conference on Wireless Networks and Mobile Communications (WINCOM)_. IEEE, 2020, pp. 1–5.
  * [11] M. Patriksson, _The traffic assignment problem: models and methods_. Courier Dover Publications, 2015.
  * [12] E. Koutsoupias and C. Papadimitriou, “Worst-case equilibria,” in _Annual Symposium on Theoretical Aspects of Computer Science_. Springer, 1999, pp. 404–413.
  * [13] R. W. Rosenthal, “A class of games possessing pure-strategy nash equilibria,” _International Journal of Game Theory_ , vol. 2, no. 1, pp. 65–67, 1973.
  * [14] D. Acemoglu, A. Makhdoumi, A. Malekian, and A. Ozdaglar, “Informational braess’ paradox: The effect of information on traffic congestion,” _Operations Research_ , vol. 66, no. 4, pp. 893–917, 2018.
  * [15] X. Chen, X. Hu, and C. Wang, “The equilibrium existence of a robust routing game under interval uncertainty,” in _International Symposium on Algorithmic Game Theory_. Springer, 2018, pp. 189–200.
  * [16] X. Chen, X. Hu, C. Wang, and X. Wu, “The efficiency of nash equilibria in the load balancing game with a randomizing scheduler,” _Theoretical Computer Science_ , vol. 838, pp. 180–194, 2020.
  * [17] S. Chowdhury and C. Giri, “Non-cooperative game theory based congestion control in lossy wsn,” in _2018 IEEE Global Communications Conference (GLOBECOM)_. IEEE, 2018, pp. 1–7.
  * [18] N. Bean, F. Kelly, and P. Taylor, “Braess’s paradox in a loss network,” _Journal of Applied Probability_ , pp. 155–159, 1997.
  * [19] H. Kameda, E. Altman, T. Kozawa, and Y. Hosokawa, “Braess-like paradoxes in distributed computer systems,” _IEEE Transactions on Automatic Control_ , vol. 45, no. 9, pp. 1687–1691, 2000.
  * [20] H. Kameda and O. Pourtallier, “Paradoxes in distributed decisions on optimal load balancing for networks of homogeneous computers,” _Journal of the ACM (JACM)_ , vol. 49, no. 3, pp. 407–433, 2002.


# A Sparse and Adaptive Prior for Time-Dependent Model Parameters

Dani Yogatama Bryan R. Routledge Noah A. Smith
_draft in review; do not cite or circulate_ Dani Yogatama
Language Technologies Institute
Carnegie Mellon University
Pittsburgh, PA 15213, USA
<EMAIL_ADDRESS>&Bryan R. Routledge
Tepper School of Business
Carnegie Mellon University
Pittsburgh, PA 15213, USA
<EMAIL_ADDRESS>&Noah A. Smith
Language Technologies Institute
Carnegie Mellon University
Pittsburgh, PA 15213, USA
<EMAIL_ADDRESS>
###### Abstract

We consider the scenario where the parameters of a probabilistic model are
expected to vary over time. We construct a novel prior distribution that
promotes sparsity and adapts the strength of correlation between parameters at
successive timesteps, based on the data. We derive approximate variational
inference procedures for learning and prediction with this prior. We test the
approach on two tasks: forecasting financial quantities from relevant text,
and modeling language contingent on time-varying financial measurements.

##  1 Introduction

When learning from streams of data to make predictions in the future, how
should we handle the timestamp associated with each instance? Ignoring
timestamps and assuming data are i.i.d. is scalable but risks distracting a
model with irrelevant “ancient history.” On the other hand, using only the
most recent portion of the data risks overfitting to current trends and
missing important time-insensitive effects. In this paper, we seek a general
approach to learning model parameters that are overall sparse, but that adapt
to variation in how different effects change over time.

Our approach is a prior over parameters of an exponential family (e.g.,
coefficients in linear or logistic regression). We assume that parameter
values shift at each timestep, with correlation between adjacent timesteps
captured using a multivariate normal distribution whose precision matrix is
restricted to a tridiagonal structure. We (approximately) marginalize the
(co)variance parameters of this normal distribution using a Jeffreys prior,
resulting in a model that allows smooth variation over time while encouraging
overall sparsity in the parameters. (The parameters themselves are not given a
fully Bayesian treatment.)

We demonstrate the usefulness of our model on two tasks, showing gains over
alternative approaches. The first is a text regression problem in which an
economic variable (volatility of returns) is forecast from financial reports
(Kogan et al., 2009). The second forecasts text by constructing a language
model that conditions on highly time-dependent economic variables.

Notation is given in §2. Our prior distribution is presented in §3. We draw
connections to related work in §4. §5 presents our inference algorithm and §6
our experimental results.

##  2 Notation

We assume data of the form $\\{(x_{n},y_{n})\\}_{n=1}^{N}$, where each $x_{n}$
includes a timestamp denoted $t\in\\{1,\ldots,T\\}$.111In this work we assume
timestamps are discretized. The aim is to learn a predictor that maps input
$x_{N+1}$, assumed to be at timestep $T$, to output $y_{N+1}$. In the
probabilistic setting we adopt here, the prediction is MAP inference over r.v.
$Y$ given $X=x$ and a model parameterized by
$\boldsymbol{\beta}\in\mathbb{R}^{I}$. Learning is parameter estimation to
solve:

$\operatorname*{argmax}_{\boldsymbol{\beta}}\log
p(\boldsymbol{\beta})+\overbrace{\sum_{n=1}^{N}\log\underbrace{p(y_{n}\mid
x_{n},\boldsymbol{\beta})}_{\mathrm{link}^{-1}(\boldsymbol{f}(x)^{\top}\boldsymbol{\beta})}}^{L(\boldsymbol{\beta})}$
(1)

The focus of the paper is on the prior distribution $p(\boldsymbol{\beta})$.
Throughout, we will denote the task-specific log-likelihood (second term) by
$L(\boldsymbol{\beta})$ and assume a generalized linear model such that a
feature vector function $\boldsymbol{f}$ maps inputs $x$ into $\mathbb{R}^{I}$
and $\boldsymbol{f}(x)^{\top}\boldsymbol{\beta}$ is “linked” to the
distribution over $Y$ using, e.g., a logit or identity. We will refer to
elements of $\boldsymbol{f}$ as “features” and to $\boldsymbol{\beta}$ as
“coefficients.” We assume $T$ discrete timesteps.

##  3 Time-Series Prior

Our time-series prior draws inspiration from the probabilistic interpretation
of the sparsity-inducing lasso (Tibshirani, 1996) and group lasso (Yuan & Lin,
2007). In non-overlapping group lasso, features are divided into groups, and
the coefficients within each group $m$ are drawn according to:

  1. 1.

Variance $\sigma^{2}_{m}\sim$ an exponential distribution.222The exponential
distribution can be replaced by the (improper) Jeffreys prior, although then
the familiar Laplace distribution interpretation no longer holds (Figueiredo,
2002).

  2. 2.

$\boldsymbol{\beta}_{m}\sim\mathrm{Normal}(\boldsymbol{0},\sigma^{2}_{m}\mathbf{I})$.

We seek a prior that lets each coefficient vary smoothly over time. A high-
level intuition of our prior is that we create copies of $\boldsymbol{\beta}$,
one at each timestep:
$\langle\boldsymbol{\beta}^{(1)},\boldsymbol{\beta}^{(2)},\ldots,\boldsymbol{\beta}^{(T)}\rangle$.
For each feature $i$, let the sequence
$\langle\beta^{(1)}_{i},\beta^{(2)}_{i},\ldots,\beta^{(T)}_{i}\rangle$ form a
group, denoted $\boldsymbol{\beta}_{i}$. Group lasso does not view
coefficients in a group as explicitly correlated; they are independent given
the variance parameter. Given the sequential structure of
$\boldsymbol{\beta}_{i}$, we replace the covariance matrix
$\sigma^{2}\mathbf{I}$ to capture autocorrelation. Specifically, we assume the
vector $\boldsymbol{\beta}_{i}$ is drawn from a multivariate normal
distribution with mean zero and a $T\times T$ precision matrix
$\mathbf{\Lambda}$ with the following tridiagonal form:333We suppress the
subscript $i$ for this discussion; each feature $i$ has its own
$\mathbf{\Lambda}_{i}$.

$\displaystyle\mathbf{\Lambda}$ $\displaystyle=\frac{1}{\lambda}\mathbf{A}$
$\displaystyle=\frac{1}{\lambda}\left[\begin{array}[]{cccccc}1&\alpha&0&0&\dots\\\
\alpha&1&\alpha&0&\dots\\\ 0&\alpha&1&\alpha&\dots\\\ 0&0&\alpha&1&\dots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right]$ (7)

$\lambda\geq 0$ is a scalar multiplier whose role is to control sparsity in
the coefficients, while $\alpha$ dictates the degree of correlation between
coefficients in adjacent timesteps (autocorrelation). Importantly, $\alpha$
and $\lambda$ (and hence $\mathbf{A}$ and $\mathbf{\Lambda}$) are allowed to
be different for each group $i$.

We need to ensure that $\mathbf{A}$ is positive definite. Fortunately, it is
easy to show that for $\alpha\in(-0.5,0.5)$, the resulting $\mathbf{A}$ is
positive definite.

###### Proof sketch.

To show this, since $\mathbf{A}$ is a symmetric matrix, we verify that each of
its principal minors have strictly positive determinants. The principal minors
of $\mathbf{A}$ are uniform tridiagonal symmetric matrices, and the
determinant of a uniform tridiagonal $N\times N$ matrix can be written as
$\prod_{n=1}^{N}\left\\{1+2\alpha\cos\left(\frac{(n+1)\pi}{N+1}\right)\right\\}$
(see, e.g., Volpi (2003) for the proof). Since $\cos(x)\in[-1,1]$, if
$\alpha\in(-0.5,0.5)$, the determinant is always positive. Therefore,
$\mathbf{A}$ is always p.d. for $\alpha\in(-0.5,0.5)$. ∎

###  3.1 Generative Model

Our generative model for the group of coefficients
$\boldsymbol{\beta}_{i}=\langle\beta^{(1)}_{i},\beta^{(2)}_{i},\ldots,\beta^{(T)}_{i}\rangle$
is given by:

  1. 1.

$\lambda_{i}\sim$ an improper Jeffreys prior
($p(\lambda)\propto\lambda^{-1}$).

  2. 2.

$\alpha_{i}\sim$ a truncated exponential prior with parameter $\tau$. This
distribution forces $\alpha_{i}$ to fall in $(-C,0]$, so that $\mathbf{A}_{i}$
is p.d. and autocorrelations are always positive:

$p(\alpha\mid\tau)=\frac{\tau\exp(-\tau(\alpha+C))\boldsymbol{1}\\{-C<\alpha\leq
0\\}}{(1-\exp(-\tau C))}.$ (8)

We fix $C=\frac{1}{2}-10^{-5}$.

  3. 3.

$\boldsymbol{\beta}_{i}\sim\mathrm{Normal}(\boldsymbol{0},\mathbf{\Lambda}_{i}^{-1})$,
with the precision matrix $\mathbf{\Lambda}_{i}$ as defined in Eq. 7.

During estimation of $\boldsymbol{\beta}$, each $\lambda_{i}$ and $\alpha_{i}$
are marginalized, giving a sparse and adaptive estimate for
$\boldsymbol{\beta}$.

###  3.2 Scalability

Our design choice of the precision matrix $\mathbf{\Lambda}_{i}$ is driven by
scalability concerns. Instead of using, e.g., a random draw from a Wishart
distribution, we specify the precision matrix to have a tridiagonal structure.
This induces dependencies between coefficients in adjacent timesteps (first-
order dependencies) and allows the prior to scale to fine-grained timesteps
more efficiently. Let $N$ denote the number of training instances, $I$ the
number of base features, and $T$ the number of timesteps. A single pass of our
variational algorithm (discussed in §5) has runtime $\mathcal{O}(I(N+T))$ and
space requirement $\mathcal{O}(I(N+T))$, instead of $\mathcal{O}(I(N+T^{2}))$
for both if each $\mathbf{\Lambda}_{i}$ is drawn from a Wishart distribution.
This can make a big difference for applications with large numbers of features
($I$). Additionally, we choose the off-diagonal entries to be uniform, so we
only need one $\alpha_{i}$ for each base feature. This design choice restricts
the expressive power of the prior but still permits flexibility in adapting to
trends for different coefficients, as we will see. The prior encourages
sparsity at the group level, essentially performing feature selection: some
feature coefficients $\boldsymbol{\beta}_{i}$ may be driven to zero across all
timesteps, while others will be allowed to vary over time, with an expectation
of smooth changes.

Note that this model introduces only one hyperparameter, $\tau$, since we
marginalize $\boldsymbol{\alpha}=\langle\alpha_{1},\ldots,\alpha_{I}\rangle$
and $\boldsymbol{\lambda}=\langle\lambda_{1},\ldots,\lambda_{I}\rangle$.

##  4 Related Work

Our model is related to autoregressive integrated moving average approaches to
time-series data (Box et al., 2008), but we never have access to _direct_
observations of the time-series. Instead, we observe data ($x$ and $y$)
assumed to have been sampled using time-series-generated variables as
_coefficients_ ($\boldsymbol{\beta}$). During learning, we therefore use
probabilistic inference to reason about the variables at all timesteps
together. In §5, we describe a scalable variational inference algorithm for
inferring coefficients at all timesteps, enabling prediction of future data
and inspection of trends.

We follow Yogatama et al. (2011) in creating time-specific copies of the base
coefficients, so that
$\boldsymbol{\beta}=\langle\boldsymbol{\beta}^{(1)},\boldsymbol{\beta}^{(2)},\ldots,\boldsymbol{\beta}^{(T)}\rangle$.
As a prior over $\boldsymbol{\beta}$, they used a multivariate Gaussian
imposing non-zero covariance between each $\beta_{i}^{(t)}$ and its time-
adjacent copies $\beta_{i}^{(t-1)}$ and $\beta_{i}^{(t+1)}$. The strength of
that covariance was set for each base feature by a global hyperparameter,
which was tuned on held-out development data along with the global variance
hyperparameter. Yogatama et al.’s model can be obtained from ours by fixing
the same $\alpha$ and $\lambda$ for all features $i$. Our approach differs in
that (i) we marginalize the hyperparameters, (ii) we allow each coefficient
its own autocorrelation, and (iii) we encourage sparsity.

There are many related Bayesian approaches for time-varying model parameters
(Belmonte et al., 2012; Nakajima & West, 2012; Caron et al., 2012), as well as
work on time-varying signal estimation (Angelosante & Giannakis, 2009;
Angelosante et al., 2009; Charles & Rozell, 2012). Each provides a different
probabilistic interpretation of parameter generation. Our model has a
distinctive generative story in that correlations between parameters of
successive timesteps are encoded in a precision matrix. Additionally, unlike
these fully Bayesian approaches that infer full posterior distributions, we
only obtain posterior mode estimates of coefficients, which has computational
advantages at prediction time (e.g., straightforward MAP inference and
sparsity) and interpretability of $\boldsymbol{\beta}$.

As noted, our grouping together of each feature’s instantiations at all
timesteps,
$\langle\beta_{i}^{(1)},\beta_{i}^{(2)},\ldots,\beta_{i}^{(T)}\rangle$ and
seeking sparsity, bears clear similarity to _group lasso_ (Yuan & Lin, 2007),
which encourages whole groups of coefficients to collectively go to zero. A
probabilistic interpretation for lasso as a two level exponential-normal
distribution that generalizes to (non-overlapping) group lasso was introduced
by Figueiredo (2002). He also showed that the exponential distribution prior
can be replaced with an improper Jeffreys prior for a parameter-free model, a
step we follow as well.

Our model is also related to the fused lasso (Tibshirani et al., 2005), which
penalizes a loss function by the $\ell_{1}$-norm of the coefficients and their
differences. Our prior has a more clear probabilistic interpretation and
adapts the degree of autocorrelation for each coefficient, based on the data.

Zhang & Yeung (2010) proposed a regularization method using a matrix-variate
normal distribution prior to model task relationships in multitask learning.
If we consider timesteps as tasks, the technique resembles our regularizer.
Their method jointly optimizes the covariance matrix with the feature
coefficients; we choose a Bayesian treatment and encode our prior belief to
the (inverse) covariance matrix, while still allowing the learned feature
coefficients to modify the matrix by posterior inference. As a result, our
method allows different base features to have different matrices.

##  5 Learning and Inference

We marginalize $\boldsymbol{\lambda}$ and $\boldsymbol{\alpha}$ and obtain a
maximum _a posteriori_ estimate for $\boldsymbol{\beta}$, which includes a
coefficient for each base feature $i$ at each timestep $t$. Specifically, we
seek to maximize:

$\displaystyle L(\boldsymbol{\beta})+\sum_{i=1}^{I}\log\int d\alpha_{i}\int
d\lambda_{i}p(\boldsymbol{\beta}_{i}\mid\alpha_{i},\lambda_{i})p(\alpha_{i}\mid\tau)p(\lambda_{i})$
(9)

Exact inference in this model is intractable. We use mean-field variational
inference to derive a lower bound on the above log-likelihood function. We
then apply a standard optimization technique to jointly optimize the
variational parameters and the coefficients $\boldsymbol{\beta}$.

We introduce fully factored variational distributions for each $\lambda_{i}$
and $\alpha_{i}$. For $\lambda_{i}$, we use a Gamma distribution with
parameters $a_{i},b_{i}$ as our variational distribution:

$\displaystyle q_{i}(\lambda_{i}\mid
a_{i},b_{i})=\frac{\lambda_{i}^{a_{i}-1}\exp(-\lambda_{i}/b_{i})}{b_{i}^{a_{i}}\Gamma(a_{i})}$

Therefore, we have $\mathbb{E}_{q_{i}}[\lambda_{i}]=a_{i}b_{i}$,
$\mathbb{E}_{q_{i}}[\lambda_{i}^{-1}]=((a_{i}-1)b_{i})^{-1}$, and
$\mathbb{E}_{q_{i}}[\log\lambda_{i}]=\Psi(a_{i})+\log b_{i}$ ($\Psi$ is the
digamma function).

For $\alpha_{i}$, we choose the form of our variational distribution to be the
same truncated exponential distribution as its prior, with parameter
$\kappa_{i}$, denoting this distribution $q_{i}(\alpha_{i}\mid\kappa_{i})$. We
have

$\displaystyle\mathbb{E}_{q_{i}}[\alpha_{i}]$
$\displaystyle=\int^{0}_{-C}\alpha_{i}\frac{\kappa_{i}\exp(-\kappa_{i}(\alpha_{i}+C))}{1-\exp(-\kappa_{i}C)}d\alpha_{i}$
$\displaystyle=\frac{1}{\kappa_{i}}-\frac{C}{1-\exp(-\kappa_{i}C)}$ (10)

We let $q$ denote the set of all variational distributions over
$\boldsymbol{\lambda}$ and $\boldsymbol{\alpha}$.

The variational bound $B$ that we seek to maximize is given in Figure 1. Our
learning algorithm involves optimizing with respect to variational parameters
$\boldsymbol{a}$, $\boldsymbol{b}$, and $\boldsymbol{\kappa}$, and the
coefficients $\boldsymbol{\beta}$. We employ the L-BFGS quasi-Newton method
(Liu & Nocedal, 1989), for which we need to compute the gradient of $B$. We
turn next to each part of this gradient.

$\displaystyle
B(\boldsymbol{a},\boldsymbol{b},\boldsymbol{\kappa},\boldsymbol{\beta})$
$\displaystyle\propto$ $\displaystyle
L(\boldsymbol{\beta})+\sum_{i=1}^{I}\left\\{\frac{1}{2}(-T\mathbb{E}_{q}[\log\lambda_{i}]\framebox{$-\mathbb{E}_{q}[\log\det{\mathbf{A}_{i}^{-1}}]$})-\mathbb{E}_{q}[\lambda_{i}^{-1}]\frac{1}{2}\boldsymbol{\beta}_{i}^{\top}\mathbb{E}_{q}[\mathbf{A}_{i}]\boldsymbol{\beta}_{i}\right\\}$
$\displaystyle+\sum_{i=1}^{I}\left\\{-(\mathbb{E}_{q}[\alpha_{i}]+C)\tau-\mathbb{E}_{q}[\log\lambda_{i}]\right\\}-\sum_{i=1}^{I}\left\\{(a_{i}-1)\mathbb{E}_{q}[\log\lambda_{i}]-\frac{\mathbb{E}_{q}[\lambda_{i}]}{b_{i}}-\log\Gamma({a_{i}})-a_{i}\log
b_{i}\right\\}$
$\displaystyle-\sum_{i=1}^{I}\left\\{\log\kappa_{i}-\kappa_{i}(\mathbb{E}_{q}[\alpha_{i}]+C)-\log(1-\exp(-\kappa_{i}C))\right\\}$
Figure 1:  The variational bound on Equation 1 that is maximized to learn
$\boldsymbol{\beta}$. The boxed expression is further bounded by
$-\log\det\mathbb{E}_{q}[\mathbf{A}_{i}]$ using Jensen’s inequality, giving a
new lower bound we denote by $B^{\prime}$.

###  5.1 Coefficients $\boldsymbol{\beta}$

For $1<t<T$, the first derivative with respect to time-specific coefficient
$\beta_{i}^{(t)}$ is:

$\frac{\partial B}{\partial\beta_{i}^{(t)}}=\frac{\partial
L}{\partial\beta_{i}^{(t)}}-\frac{1}{2}\mathbb{E}[\lambda_{i}^{-1}]\left(\mathbb{E}[\alpha_{i}](\beta_{i}^{(t-1)}+\beta_{i}^{(t+1)})+2\beta_{i}^{(t)}\right)$
(11)

We can interpret the first derivative as including a penalty scaled by
$\mathbb{E}[\lambda_{i}^{-1}]$. We rewrite this penalty as:

$\displaystyle\mathbb{E}[\lambda_{i}^{-1}]\left(\vphantom{1-\mathbb{E}[\alpha_{i}])2\beta_{i}^{(t)}}\right.$
$\displaystyle(1-\mathbb{E}[\alpha_{i}])$ $\displaystyle\cdot
2\beta_{i}^{(t)}$ $\displaystyle+\mathbb{E}[\alpha_{i}]$
$\displaystyle\cdot(\beta_{i}^{(t)}-\beta^{(t-1)}_{i})$
$\displaystyle+\mathbb{E}[\alpha_{i}]$
$\displaystyle\left.\cdot(\beta_{i}^{(t)}-\beta^{(t+1)}_{i})\right)$

This form makes it clear that the penalty depends on $\beta^{(t-1)}_{i}$ and
$\beta^{(t+1)}_{i}$, penalizing the difference between $\beta^{(t)}_{i}$ and
these time-adjacent coefficients proportional to $\mathbb{E}[\alpha_{i}]$.

The form bears strong similarity to the first derivative of the time-series
(log-)prior introduced in Yogatama et al. (2011), which depends on fixed,
global hyperparameters analogous to our $\alpha$ and $\lambda$. Because our
approach does not require us to specify scalars playing the roles of
“$\mathbb{E}[\lambda_{i}^{-1}]$” and “$\mathbb{E}[\alpha_{i}]$” in advance, it
is possible for each feature to have its own autocorrelation. Obtaining the
same effect in their model would require careful tuning of $\mathcal{O}(I)$
hyperparameters, which is not practical.

It also has some similarities to the fused lasso penalty (Tibshirani et al.,
2005), which is intended to encourage sparsity in the differences between
features coefficients across timesteps. Our prior, on the other hand,
encourages smoothness in the differences, with additional sparsity at the
feature level.

###  5.2 Variational Parameters for $\boldsymbol{\alpha}$ and
$\boldsymbol{\lambda}$

Recall that the variational distribution for $\lambda_{i}$ is a Gamma
distribution with parameters $a_{i}$ and $b_{i}$.

##### Precision matrix scalar $\boldsymbol{\lambda}$.

The first derivative for variational parameters $\boldsymbol{a}$ is easy to
compute:

$\frac{\partial B}{\partial
a_{i}}=\left(-\frac{T}{2}-a_{i}\right)\Psi_{1}(a_{i})+\frac{\boldsymbol{\beta}_{i}^{\top}\mathbb{E}[\mathbf{A}_{i}]\boldsymbol{\beta}_{i}}{2b_{i}(a_{i}-1)^{2}}+1$
(12)

where $\Psi_{1}$ is the trigamma function. We can solve for $\boldsymbol{b}$
in closed form given the other free variables:

$b_{i}=\frac{\boldsymbol{\beta}_{i}^{\top}\mathbb{E}[\mathbf{A}_{i}]\boldsymbol{\beta}_{i}}{(a_{i}-1)T}$
(13)

We therefore treat $\boldsymbol{b}$ as a function of $\boldsymbol{a}$,
$\boldsymbol{\kappa}$, and $\boldsymbol{\beta}$ in optimization.

##### Off-diagonal entries $\boldsymbol{\alpha}$.

First, notice that using Jensen’s inequality:
$\mathbb{E}[\log\det{\mathbf{A}_{i}^{-1}}]=\mathbb{E}[-\log\det{\mathbf{A}_{i}}]\geq-\log\det\mathbb{E}[\mathbf{A}_{i}]$
due to the fact that $-\log\det\mathbf{A}_{i}$ is a convex function.
Furthermore, for a uniform symmetric tridiagonal matrix like $\mathbf{A}_{i}$,
the log determinant can be computed in closed form as follows (Volpi, 2003):

$\displaystyle\log\det\mathbb{E}[\mathbf{A}_{i}]=$
$\displaystyle\log\left(\prod_{t=1}^{T}1+2\mathbb{E}[\alpha_{i}]\cos\left(\frac{(t+1)\pi}{T+1}\right)\right)$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\log\left(1+2\mathbb{E}[\alpha_{i}]\cos\left(\frac{(t+1)\pi}{T+1}\right)\right)$

We therefore maximize a lower bound on $B$, making use of the above to
calculate first derivatives with respect to $\kappa_{i}$:

$\displaystyle\frac{\partial B^{\prime}}{\partial\kappa_{i}}=$
$\displaystyle-\tau\frac{\partial\mathbb{E}[\alpha_{i}]}{\partial\kappa_{i}}-\frac{1}{\kappa_{i}}+C+\mathbb{E}[\alpha_{i}]+\frac{\partial\mathbb{E}[\alpha_{i}]}{\partial\kappa_{i}}\kappa_{i}$
$\displaystyle+\frac{C\exp(-C\kappa_{i})}{1-\exp(-C\kappa_{i})}+\frac{1}{2}\frac{\partial\log\det\mathbb{E}[\mathbf{A}_{i}]}{\partial\kappa_{i}}$
$\displaystyle-\frac{1}{2}\mathbb{E}[\lambda_{i}^{-1}]\frac{\partial\boldsymbol{\beta}_{i}^{\top}\mathbb{E}[\mathbf{A}_{i}]\boldsymbol{\beta}_{i}}{\partial\kappa_{i}}$

The partial derivatives
$\frac{\partial\mathbb{E}[\alpha_{i}]}{\partial\kappa_{i}}$,
$\frac{\partial\log\det\mathbb{E}[\mathbf{A}_{i}]}{\partial\kappa_{i}}$, and
$\frac{\partial\boldsymbol{\beta}_{i}^{\top}\mathbb{E}[\mathbf{A}_{i}]\boldsymbol{\beta}_{i}}{\partial\kappa_{i}}$
are easy to compute. We omit them for space.

###  5.3 Implementation Details

A well-known property of numerical optimizers like the one we use (L-BFGS; Liu
& Nocedal (1989)) is the failure to reach optimal values exactly at zero.
Although theoretically strongly sparse, our prior only produces weak sparsity
in practice. Future work might consider a more principled proximal-gradient
algorithm to obtain strong sparsity (Bach et al., 2011; Liu & Ye, 2010; Duchi
& Singer, 2009).

If we expect feature coefficients at specific timesteps to be sparse as well,
it is straightforward to incorporate additional terms in the objective
function that encode this prior belief (analogous to an extension from group
lasso to _sparse_ group lasso). For the tasks we consider in our experiments,
we found that it does not substantially improve the overall performance.
Therefore, we keep the simpler bound given in Figure 1.

##  6 Experiments

We report two sets of experiments, one with a continuous $y$, the other a
language modeling application where $y$ is text. Each timestep in our
experiments is one year.

###  6.1 Baselines

On both tasks, we compare our approach to a range of baselines. Since this is
a forecasting task, at each test year, we only used training examples that
come from earlier years. Our baselines vary in how they use this earlier data
and in how they regularize.

  * •

ridge-one: ridge regression (Hoerl & Kennard, 1970), trained on only examples
from the year prior to the test data (e.g., for the 2002 task, train on
examples from 2001)

  * •

ridge-all: ridge regression trained on the full set of past examples (e.g.,
for the 2002 task, train on examples from 1996–2001)

  * •

ridge-ts: the non-adaptive time-series ridge model of Yogatama et al. (2011)

  * •

lasso-one: lasso regression (Tibshirani, 1996), trained on only examples from
the year prior to the test data444Brendan O’Connor (personal communication)
has established the superiority of the lasso to the support vector regression
method of Kogan et al. (2009) on this dataset; lasso is a strong baseline for
this problem.

  * •

lasso-all: lasso regression trained on the full set of past examples

In all cases, we tuned hyperparameters on a development data. Note that, of
the above baselines, only ridge-ts replicates the coefficients at different
timesteps (i.e., $IT$ parameters); the others have only $I$ time-insensitive
coefficients.

The model with our prior always uses all training examples that are available
up to the test year (this is equivalent to a sliding window of size infinity).
Like ridge-ts, our model trusts more recent data more, allowing coefficients
farther in the past to drift farther away from those most relevant for
prediction at time $T+1$. Our model, however, adapts the “drift” of each
coefficient separately rather than setting a global hyperparameter.

###  6.2 Forecasting Risk from Text

In the first experiment, we apply our prior to a forecasting task. We consider
the task of predicting volatility of stock returns from financial reports of
publicly-traded companies, similar to Kogan et al. (2009).

Table 1:  MSE on the 10-K dataset (various test sets). The first test year (2002) was used as our development data. Our model uses the sparse adaptive prior described in §3. The overall differences between our model and all competing models are statistically significant (Wilcoxon signed-rank test, $p<0.01$). year | # examples | ridge-one | ridge-all | ridge-ts | lasso-one | lasso-all | our model
---|---|---|---|---|---|---|---
2002(dev) | 2,845 | 0.182 | 0.176 | 0.171 | 0.165 | 0.156 | 0.158
2003 | 3,611 | 0.185 | 0.173 | 0.171 | 0.164 | 0.176 | 0.164
2004 | 3,558 | 0.125 | 0.137 | 0.129 | 0.116 | 0.119 | 0.113
2005 | 3,474 | 0.135 | 0.133 | 0.136 | 0.124 | 0.124 | 0.122
overall | 13,488 | 0.155 | 0.154 | 0.151 | 0.141 | 0.143 | 0.139

In finance, _volatility_ refers to a measure of variation in a quantity over
time; for stock returns, it is measured using the standard deviation during a
fixed period (here, one year). Volatility is used as a measure of financial
risk. Consider a linear regression model for predicting the log
volatility555Similar to Kogan et al. (2009) and as also the standard practice
in finance, we perform a log transformation, since log volatilities are
typically close to normally distributed. of a stock from a set of features
(see §6.2.1 for a complete description of our features). We can interpret a
linear regression model probabilistically as drawing $y\in\mathbb{R}$ from a
normal distribution with $\boldsymbol{\beta}^{\top}\boldsymbol{f}(x)$ as the
mean of the normal. Therefore, in this experiment:
$L(\boldsymbol{\beta})=-\sum_{t=1}^{T}\sum_{i=1}^{N_{t}}(y_{i}^{(t)}-\boldsymbol{\beta}^{(t)\top}\boldsymbol{f}(x_{i}^{(t)}))^{2}$.

We apply the time-series prior to the feature coefficients
$\boldsymbol{\beta}$. When making a prediction for the test data, we use
$\boldsymbol{\beta}^{(T)}$, the set of feature coefficients for the last
timestep in the training data.

####  6.2.1 Dataset

We used a collection of Securities Exchange Commission-mandated annual reports
from 10,492 publicly traded companies in the U.S. There are 27,159 reports
over a period of ten years from 1996–2005 in the corpus. These reports are
known as “Form 10-K.”666See Kogan et al. (2009) for a complete description of
the dataset; it is available at http://www.ark.cs.cmu.edu/10K. For the feature
set, we downcased and tokenized the texts and selected the 101st–10,101st most
frequent words as binary features. The feature set was kept the same across
experiments for all models. It is widely known in the financial community that
the past history of volatility of stock returns is a good indicator of the
future volatility. Therefore, we also included the log volatility of the
stocks twelve months prior to the report as a feature. Our response variable
$y$ is the log volatility of stock returns over a period of twelve months
after the report is published.

####  6.2.2 Results

The first test year (i.e., 2002) was used as our development data for
hyperparameter tuning ($\tau$ was selected to be $1.0$). We initialized all
the feature coefficients by the coefficients from training a lasso regression
on the last year of the training data (lasso-one). Table 1 provides a summary
of experimental results. We report the results in mean squared error on the
test set: $\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}$, where $y_{i}$ is
the true response for instance $i$ and $\hat{y}_{i}$ is the predicted
response.

Our model consistently outperformed ridge variants, including the one with a
time-series penalty (Yogatama et al., 2011). It also outperformed the lasso
variants without any time-series penalty, on average and in three out of four
test sets apiece.

One of the major challenges in working with time-series data is to choose the
right window size, in which the data is still relevant to current predictions.
Our model automates this process with a Bayesian treatment of the strength of
each feature coefficient’s autocorrelation. The results indicate that our
model was able to learn when to trust a longer history of training data, and
when to trust a shorter history of training data, demonstrating the
adaptiveness of our prior. Figure 2 shows the distribution of the expected
values of the autocorrelation paramaters under the variational distributions
$\mathbb{E}_{q_{i}}[\alpha_{i}]$ for 10,002 features, learned by our model
from the last run (test year 2005).

In future work, an empirical Bayesian treatment of the hyperprior $\tau$,
fitting it to improve the variational bound, might lead to further
improvements.

Figure 2:  The distribution of expected values of the autocorrelation
paramaters under the variational distributions
$\mathbb{E}_{q_{i}}[\alpha_{i}]$ for 10,002 features used in our experiments
(10,000 unigram features, the previous year log volatility feature, and a bias
feature).

###  6.3 Text Modeling in Context

In the second experiment, we consider a hard task of modeling a collection of
texts over time conditioned on economic measurements. The goal is to predict
the probability of words appearing in a document, based on the “state of the
world” at the time the document was authored. Given a set of macroeconomic
variables in the U.S. (e.g., unemployment rate, inflation rate, average
housing prices, etc.), we want to predict what kind of texts will be produced
at a specific timestep. These documents can be written by either the
government or publicly-traded companies directly or indirectly affected by the
current economic situation.

####  6.3.1 Model

Our text model is a sparse additive generative model (SAGE; Eisenstein et al.
(2011)). In SAGE, there is a background lexical distribution that is perturbed
additively in the log-space. When the effects are due to a (sole) feature
$f(x)$, the probability of a word is:

$\displaystyle
p(w\mid\boldsymbol{\theta},\boldsymbol{\beta},x)=\frac{\exp(\theta_{w}+\beta_{w}f(x))}{\sum_{w^{\prime}\in
V}\exp(\theta_{w^{\prime}}+\beta_{w^{\prime}}f(x))}$

where $V$ is the vocabulary, $\boldsymbol{\theta}$ (always observed) is the
vector of background log-frequencies of words in the corpus, $f(x)$ (observed)
is the feature derived from the context $x$, and $\boldsymbol{\beta}$ is the
feature-specific deviation.

Notice that the formulation above is easily extended to multiple effects with
coefficients $\boldsymbol{\beta}$. In our experiment, we have 117 effects
(features), each with its own $\boldsymbol{\beta}_{i}$. The first 50
correspond to U.S. states, plus an additional feature for the entire U.S., and
they are observed for each text since each text is associated with a known set
of states (discussed below). We assume that texts that are generated in
different states have distinct characteristics; for each state, we have a
binary indicator feature. The other 66 features depend on observed
macroeconomics variables at each timestep (e.g., unemployment rate, inflation
rate, house price index, etc.). Given an economic state of the world, we
hypothesize that there are certain words that are more likely to be used, and
each economic variable has its own (sparse) deviation from the background word
frequencies. The generative story for a word at timestep $t$ associated with
(observed) features $\boldsymbol{f}(x^{(t)})$ is:

  * •

Given observed real-world observed variables $x^{(t)}$, draw word $w$ from a
multinomial distribution
$p(w\mid\boldsymbol{\theta}^{(t)},\boldsymbol{\beta}^{(t)},x^{(t)})\propto\exp(\theta_{w}^{(t)}+\boldsymbol{\beta}_{w}^{(t)\top}\boldsymbol{f}(x^{(t)}))$.

Our $L(\boldsymbol{\beta})$ is simply the negative log-loss function commonly
used in multiclass logistic regression:
$L(\boldsymbol{\beta})=\sum_{t=1}^{T}\sum_{i=1}^{N_{t}}\log
p(\boldsymbol{w}_{i}^{(t)}\mid\boldsymbol{\theta}^{(t)},\boldsymbol{\beta}^{(t)},x^{(t)}_{i})$.
We apply our time-series prior from §3 to $\boldsymbol{\beta}$.
$\boldsymbol{\theta}^{(t)}$ is fixed to be the log frequencies of words at
timestep $t$. For a single feature, coefficients over time for different
classes (words) are assumed to be generated from the same prior.

####  6.3.2 Dataset

There is a great deal of text that is produced to describe current
macroeconomic events. We conjecture that the connection between the economy
and the text will have temporal dependencies (e.g., the amount of discussion
about housing or oil prices might vary over time). We use three sources of
text commentary on the economy. The first is a subset of the 10-K reports we
used in §6.2. We selected the 10-K reports of 200 companies chosen randomly
from the top quintile of size (measured by beginning-of-sample market
capitalization). This gives us a manageable sample of the largest U.S.
companies. Each report is associated with the state in which the company’s
head office is located. Our next two data sources come from the Federal
Reserve System, the primary body responsible for monetary policy in the
U.S.777 For an overview of the Federal Reserve System, see the Federal
Reserve’s “Purpose and Functions” document at
http://www.federalreserve.gov/pf/pf.htm. The Federal Open Market Committee
(FOMC) meets roughly eight times per year to discuss economic conditions and
set monetary policy. Prior to each meeting, each of the twelve regional banks
write an informal “anecdotal” description of economic activity in their region
as well as a national summary. This “Beige Book” is akin to a blog of economic
activity released prior to each meeting. Each FOMC meeting also produces a
transcript of the discussion. For our experiments here, we focus on text from
1996–2006.888All the text is freely available at
http://www.federalreserve.gov. The Beige Book is released to the public prior
to each meeting. The transcripts are released five years after the meetings.
As a result, we have 2,075 documents in the final corpus, consisting of 842
documents of the 10-K reports, 89 documents of the FOMC meeting transcripts,
and 1,144 documents of the Beige Book summaries.

We use the 501st–5,501st most frequent words in the dataset. We associated the
FOMC meeting transcripts with all states. The “Beige Book” texts were produced
by the Federal Reserve Banks. There are twelve Federal Reserve Banks in the
United States, each serving a collection of states. We associated texts from a
Federal Reserve Bank with the states that it serves.

Table 2:  Negative log-likelihood of the documents on various test sets (lower is better). The first test year (2003) was used as our development data. Our model uses the sparse adaptive prior in §3. | # tokens | ridge-one | ridge-all | ridge-ts | lasso-one | lasso-all | our model
---|---|---|---|---|---|---|---
year | ($\times 10^{6}$) | ($\times 10^{3}$) | ($\times 10^{3}$) | ($\times 10^{3}$) | ($\times 10^{3}$) | ($\times 10^{3}$) | ($\times 10^{3}$)
2003(dev) | 1.1 | 2,736 | 2,765 | 2,735 | 2,736 | 2,765 | 2,735
2004 | 1.5 | 2,975 | 3,004 | 2,975 | 2,975 | 3,004 | 2,974
2005 | 1.9 | 2,999 | 3,027 | 2,997 | 2,998 | 3,027 | 2,997
2006 | 2.3 | 2,916 | 2,922 | 2,913 | 2,912 | 2,922 | 2,912
overall | 6.8 | 11,626 | 11,718 | 11,619 | 11,620 | 11,718 | 11,618

Quantitative U.S. macroeconomic data was obtained from the Federal Reserve
Bank of St. Louis data repository (“FRED”). We used standard measures of
economic activity focusing on output (GDP), employment, and specific markets
(e.g., housing).999For growing output series, like GDP, we calculate growth
rates as log differences. We use equity market returns for the U.S. market as
a whole and various industry and characteristic portfolios.101010Returns are
monthly, excess of the risk-free rate, and continuously compounded. The data
are from CRSP and are available for these portfolios at
http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html. They
are used as $\boldsymbol{f}(x)$ in our model; in addition to state indicator
variables, there are 66 macroeconomic variables in total.

We compare our model to the baselines in §6.1. The lasso variants are
analogous to the original formulation of SAGE (Eisenstein et al., 2011),
except that our model directly conditions on macroeconomic variables instead
of a Dirichlet-multinomial compound.

####  6.3.3 Results

We score models by computing the negative log-likelihood on the test
dataset:111111Out-of-vocabulary items are ignored. $-\sum_{i=1}^{N}\log
p(\boldsymbol{w}_{i}^{(T+1)}\mid\boldsymbol{\theta}^{(T)},\boldsymbol{\beta}^{(T)},x^{(T+1)}_{i})$.
We initialized all the feature coefficients by the coefficients by training a
lasso regression on the last year of the training data (lasso-one). The first
test year (i.e., 2003) was used as our development data for hyperparameter
tuning ($\tau$ was selected to be $.001$). Table 2 shows the results for the
six models we compared. Similar to the forecasting experiments, at each test
year, we trained only on documents from earlier years. When we collapsed all
the training data and ignored the temporal dimension (ridge-all and lasso-
all), the background log-frequencies $\boldsymbol{\theta}^{(t)}$ are computed
using the entire training data, which is different compared to the background
log-frequencies for only the last timestep of the training data. Our model
outperformed all ridge and lasso variants, including the one with a time-
series penalty (Yogatama et al., 2011), in terms of negative log-likelihood on
unseen dataset.

In addition to improving predictive accuracy, the prior also allows us to
discover trends in the feature coefficients and gain insight. We manually
examined the model from the last run (test year 2006). Examples of temporal
trends learned by our model are shown in Figure 3. The plot illustrates
feature coefficients for words that contain the string employ. For comparison,
we also included the percentage of unemployment rate in the U.S. (which was
used as one of the features $f(x)$), scaled to fit into the plot. We can see
that there is a correlation between feature coefficients for the word
unemployment and the actual unemployment rate. On the other hand, the
correlations are less evident for other words.

Figure 3:  Temporal trends learned by our model for the words that contains
employ in our dataset, as well as the actual unemployment rate (scaled by
$10^{-16}$ for ease of reading). The $y$-axis denotes coefficients and the
$x$-axis is years. See the text for explanation.

##  7 Conclusions

We presented a time-series prior for the parameters of probabilistic models;
it produces sparse models and adapts the strength of temporal effects on each
coefficient separately, based on the data, without an explosion in the number
of hyperparameters. We showed how to do inference under this prior using
variational approximations. We evaluated the prior for the task of forecasting
volatility of stock returns from financial reports, and demonstrated that it
outperforms other competing models. We also evaluated the prior for the task
of modeling a collection of texts over time, i.e., predicting the probability
of words given some observed real-world variables. We showed that the prior
achieved state-of-the-art results as well.

## Acknowledgments

The authors thank several anonymous reviewers for helpful feedback on earlier
drafts of this paper. This research was supported in part by a Google research
award to the second and third authors. This research was supported in part by
the Intelligence Advanced Research Projects Activity via Department of
Interior National Business Center contract number D12PC00347. The U.S.
Government is authorized to reproduce and distribute reprints for Governmental
purposes notwithstanding any copyright annotation thereon. The views and
conclusions contained herein are those of the authors and should not be
interpreted as necessarily representing the official policies or endorsements,
either expressed or implied, of IARPA, DoI/NBC, or the U.S. Government.

## References

  * Angelosante & Giannakis (2009) Angelosante, Daniele and Giannakis, Georgios B.  Rls-weighted lasso for adaptive estimation of sparse signals.  In _Proc. of ICASSP_ , 2009.
  * Angelosante et al. (2009) Angelosante, Daniele, Giannakis, Georgios B., and Grossi, Emanuele.  Compressed sensing of time-varying signals.  In _Proc. of ICDSP_ , 2009.
  * Bach et al. (2011) Bach, Francis, Jenatton, Rodolphe, Mairal, Julien, and Obozinski, Guillaume.  _Convex optimization with sparsity inducing norms_.  MIT Press, 2011.
  * Belmonte et al. (2012) Belmonte, Miguel A. G., Koop, Gary, and Korobilis, Dimitris.  Hierarchical shrinkage in time-varying parameter models, 2012.  Working paper.
  * Box et al. (2008) Box, George E. P., Jenkins, Gwilym M., and Reinsel, Gregory C.  _Time Series Analysis: Forecasting and Control_.  Wiley Series in Probability and Statistics, 2008.
  * Caron et al. (2012) Caron, François, Bornn, Luke, and Doucet, Arnaud.  Sparsity-promoting bayesian dynamic linear models, 2012.  arXiv 1203.0106.
  * Charles & Rozell (2012) Charles, Adam S. and Rozell, Christopher J.  Re-weighted $\ell_{1}$ dynamic filtering for time-varying sparse signal estimation, 2012.  arXiv 1208.0325.
  * Duchi & Singer (2009) Duchi, John and Singer, Yoram.  Efficient online and batch learning using forward backward splitting.  _Journal of Machine Learning Research_ , 10:2873–2908, 2009\.
  * Eisenstein et al. (2011) Eisenstein, Jacob, Ahmed, Amr, and Xing, Eric P.  Sparse additive generative models of text.  In _Proc. of ICML_ , 2011.
  * Figueiredo (2002) Figueiredo, Mario A. T.  Adaptive sparseness using jeffreys’ prior.  In _Proc. of NIPS_ , 2002.
  * Hoerl & Kennard (1970) Hoerl, Arthur E. and Kennard, Robert W.  Ridge regression: Biased estimation for nonorthogonal problems.  _Technometrics_ , 12(1):55–67, 1970.
  * Kogan et al. (2009) Kogan, Shimon, Levin, Dimitry, Routledge, Bryan R., Sagi, Jacob S., and Smith, Noah A.  Predicting risk from financial reports with regression.  In _Proc. of NAACL_ , 2009.
  * Liu & Nocedal (1989) Liu, Dong C. and Nocedal, Jorge.  On the limited memory BFGS method for large scale optimization.  _Mathematical Programming B_ , 45(3):503–528, 1989.
  * Liu & Ye (2010) Liu, Jun and Ye, Jieping.  Moreau-yosida regularization for grouped tree structure learning.  In _Proc. of NIPS_ , 2010.
  * Nakajima & West (2012) Nakajima, Jouchi and West, Mike.  Bayesian analysis of latent threshold dynamic models.  _Journal of Business and Economic Statistics_ , 2012.
  * Tibshirani (1996) Tibshirani, Robert.  Regression shrinkage and selection via the lasso.  _Journal of Royal Statistical Society B_ , 58(1):267–288, 1996.
  * Tibshirani et al. (2005) Tibshirani, Robert, Saunders, Michael, Rosset, Saharon, Zhu, Ji, and Knight, Keith.  Sparsity and smoothness via the fused lasso.  _Journal of the Royal Statistical Society B_ , 67(1):91–108, 2005.
  * Volpi (2003) Volpi, Leonardo.  Eigenvalues and eigenvectors of tridiagonal uniform matrices.  2003\.
  * Yogatama et al. (2011) Yogatama, Dani, Heilman, Michael, O’Connor, Brendan, Dyer, Chris, Routledge, Bryan R., and Smith, Noah A.  Predicting a scientific community’s response to an article.  In _Proc. of EMNLP_ , 2011.
  * Yuan & Lin (2007) Yuan, Ming and Lin, Yi.  Model selection and estimation in regression with grouped variables.  _Journal of the Royal Statistical Society, Series B_ , 68(1):49–67, 2007.
  * Zhang & Yeung (2010) Zhang, Yu and Yeung, Dit-Yan.  A convex formulation for learning task relationships in multi-task learning.  In _Proc. of UAI_ , 2010.