# Best Arm Identification in Graphical Bilinear Bandits

Geovani Rizk  Albert Thomas  Igor Colin  Rida Laraki  Yann Chevaleyre

###### Abstract

We introduce a new graphical bilinear bandit problem where a learner (or a
_central entity_) allocates arms to the nodes of a graph and observes for each
edge a noisy bilinear reward representing the interaction between the two end
nodes. We study the best arm identification problem in which the learner wants
to find the graph allocation maximizing the sum of the bilinear rewards. By
efficiently exploiting the geometry of this bandit problem, we propose a
_decentralized_ allocation strategy based on random sampling with theoretical
guarantees. In particular, we characterize the influence of the graph
structure (e.g. star, complete or circle) on the convergence rate and propose
empirical experiments that confirm this dependency.

Machine Learning, ICML

##  1 Introduction

In many multi-agent systems the contribution of an agent to a common team
objective is impacted by the behavior of the other agents. The agents must
coordinate (or be coordinated) to achieve the best team performance. Consider,
for instance, the problem of configuring antennas of a wireless cellular
network to obtain the best signal quality over the whole network (Siomina et
al.,, 2006). The signal quality of the region covered by a given antenna might
be degraded by the behavior of its neighboring antennas due to an increase of
interferences or bad user handovers. Another example is the adjustment of the
turbine blades of a wind farm where the best adjustment for one turbine may
generate turbulence for its neighboring turbines and thus be suboptimal for
the global wind farm objective (Bargiacchi et al.,, 2018).

These real-life problems can be viewed as instances of a _stochastic multi-
agent multi-armed bandit_ problem (Robbins,, 1952; Bargiacchi et al.,, 2018)
where a learner (or a central entity) sequentially pulls a joint arm, one arm
for each agent (e.g., all the configuration parameters of the antennas), and
receives an associated global noisy reward (e.g., the signal quality over the
whole network). The goal of the learner can either be to maximize the
accumulated reward, implying a trade-off between exploration and exploitation,
or to find the joint arm maximizing the reward, known as _pure exploration_ or
_best arm identification_ (Bubeck et al.,, 2009; Audibert and Bubeck,, 2010).

In this paper we focus on the best arm identification problem in a multi-agent
system for which we assume the knowledge of a coordination graph
$\mathcal{G}=(V,E)$ representing the agent interactions (Guestrin et al.,,
2002).

At each round $t$, a learner

  1. 1.

chooses for each node $i\in V$ an arm $x_{t}^{(i)}$ in an finite arm set
$\mathcal{X}\subset\mathbb{R}^{d}$,

  2. 2.

observes for each edge $(i,j)\in E$ a bilinear reward
$r_{t}^{(i,j)}=x_{t}^{(i)\top}{\mathbf{M}_{\star}}x_{t}^{(j)}+\eta_{t}^{(i,j)}$.

Here, we denote by $\mathbf{M}_{\star}\in\mathbb{R}^{d\times d}$ the unknown
parameter matrix, and $\eta^{(i,j)}_{t}$ a zero-mean $\sigma$-sub-Gaussian
random variable for all edges $(i,j)\in E$ and round $t$.

The goal of the central entity is to find, within a minimum number of rounds,
the joint arm $(x^{(1)}_{\star},\dots,x^{(|V|)}_{\star})$ such that the
expected global reward $\sum_{(i,j)\in
E}x_{\star}^{(i)\top}{\mathbf{M}_{\star}}x_{\star}^{(j)}$ is maximized.

The reward $r^{(i,j)}_{t}$ reflects the quality of the interaction between the
neighboring nodes $i$ and $j$ when pulling respectively the arm $x^{(i)}_{t}$
and $x^{(j)}_{t}$ at time $t$. For instance, when configuring handover
parameters of a wireless network, $r^{(i,j)}_{t}$ can be any criterion
assessing the handover quality between antenna $i$ and antenna $j$, the
parameters selected by each antenna both impacting this quantity. The bilinear
setting appears as a natural extension of the commonly studied linear setting
to model the interaction between two agents. Furthermore, instead of a global
reward being a sum of independent linear agent rewards, the global reward is
now the result of the interactions between neighboring agents.

As exposed in Jun et al., (2019), the bilinear reward can be written as a
linear reward in a higher dimensional space:

$\displaystyle
r^{(i,j)}_{t}=\operatorname*{vec}{\left(x^{(i)}_{t}x^{(j)\top}_{t}\right)}^{\top}\operatorname*{vec}{\left(\mathbf{M}_{\star}\right)}+\eta^{(i,j)}_{t}\enspace,$
(1)

where for any matrix $\mathbf{A}\in\mathbb{R}^{d\times d}$,
$\operatorname*{vec}{\left(\mathbf{A}\right)}$ denotes the vector in
$\mathbb{R}^{d^{2}}$ which is the concatenation of all the columns of
$\mathbf{A}$.

Since the unknown parameter $\mathbf{M}_{\star}$ is common to all the edges
$(i,j)$ of the graph, the expected global reward at time $t$ can also be
written as the scalar product $\left\langle\sum_{(i,j)\in
E}\operatorname*{vec}{\left(x^{(i)}_{t}x_{t}^{(j)\top}\right)},\operatorname*{vec}{(\mathbf{M}_{\star})}\right\rangle$.
Hence, solving the best arm identification problem in the described graphical
bilinear bandit boils down to solving the same problem in a global linear
bandit. Although this trick allows to use classical algorithms in linear
bandits, the number of joint arms is growing exponentially with the number of
nodes, making such methods impractical.

Another possible way to address this problem based on equation (1) is to
consider one linear bandit per edge, with constraints between edges. For more
clarity, let us define the arm set
$\mathcal{Z}=\\{\operatorname*{vec}{(xx^{\prime\top})}\ |\
(x,x^{\prime})\in\mathcal{X}^{2}\\}$, and let us refer to any
$z\in\mathcal{Z}$ as an _edge-arm_ and to any $x\in\mathcal{X}$ as an _node-
arm_. At each round $t$, the learner chooses for each edge $(i,j)$ an edge-arm
in $\mathcal{Z}$ with the constraint that for any pair of edges $(i,j)$ and
$(i,k)$, if the edge-arm $\operatorname*{vec}{(xx^{\prime\top})}$ is assigned
to the edge $(i,j)$ and the edge-arm
$\operatorname*{vec}{(x^{\prime\prime}x^{\prime\prime\prime\top})}$ is
assigned to the edge $(i,k)$, then it must be that $x=x^{\prime\prime}$.

Given this constraint, how do we choose the appropriate sequence of edge-arms
in order to build a good estimate of
$\operatorname*{vec}{(\mathbf{M}_{\star})}$? Moreover, assuming we have built
such a good estimator, is there a tractable algorithm to identify the best
joint arm, or at least to find a joint arm yielding a high expected reward? In
this paper, we answer these questions and provide algorithms and theoretical
guarantees.

We show that even with a perfect estimator
$\operatorname*{vec}{(\hat{\mathbf{M}})}=\operatorname*{vec}{(\mathbf{M}_{\star})}$,
identifying the best joint arm is NP-Hard. To address this issue, we design a
polynomial time twofold algorithm. Given
$\operatorname*{vec}{(\hat{\mathbf{M}})}$, it first identifies the best edge-
arm $z_{\star}\in\mathcal{Z}$ maximizing $\langle
z_{\star},\operatorname*{vec}{(\hat{\mathbf{M}})}\rangle$. Then, it allocates
$z_{\star}$ to a carefully chosen subset of edges. We show that this yields a
good approximation ratio in Section 4.

To build our estimator $\operatorname*{vec}{(\hat{\mathbf{M}})}$, we rely on
the G-Allocation strategy, as in Soare et al., (2014). We show that there
exists a sampling procedure over the node-arms such that the associated edge-
arms follow the optimal G-allocation strategy developed in the linear bandit
literature. This procedure allows us to avoid the difficulty of having to
satisfy the edge-arm constraints explicitly. Furthermore, we analyze the
sample complexity of this method. This is detailed in Section 5.

In addition, we highlight the impact of the graph structure in Section 6 and
provide the explicit repercussion on the convergence rate of the algorithm for
different types: star, complete, circle and matching graphs. In particular, we
show that for favorable graph structures (e.g. circles), our convergence rate
matches that of standard linear bandits. Finally, Section 7 evidences the
theoretical findings on numerical experiments.

##  2 Related Work

Best arm identification in linear bandits.  There exists a vast literature on
the problem of best arm identification in linear bandits (Soare et al.,, 2014;
Xu et al.,, 2018; Degenne et al.,, 2020; Kazerouni and Wein,, 2019; Zaki et
al.,, 2020; Jedra and Proutiere,, 2020), would it be by using greedy
strategies (Soare et al.,, 2014), rounding procedures (Fiez et al.,, 2019) or
random sampling (Tao et al.,, 2018). Although our problem can be formulated as
a linear bandit problem, none of the existing methods would scale-up with the
number of agents. Nevertheless, we will be relying on classical techniques,
and more specifically those developed in Soare et al., (2014).

Bilinear bandits.  Bandits with bilinear rewards have been studied in Jun et
al., (2019). The authors derived a no-regret algorithm based on Optimism in
the Face of Uncertainty Linear bandit (OFUL) (Abbasi-Yadkori et al.,, 2011),
using the fact that a bilinear reward can be expressed as a linear reward in
higher dimension. Our work extends their setting by considering a set of
dependent bilinear bandits. Besides, the goal here is to find the best arm
rather than minimizing the regret.

Bandits and graphs.  Graphs are often used to bring structure to a bandit
problem. In Valko et al., (2014) and Mannor and Shamir, (2011), the arms are
the nodes of a graph and pulling an arm gives information on the rewards of
the neighboring arms. The reader can also refer to Valko, (2020) for an
account on such problems. In Cesa-Bianchi et al., (2013) each node is an
instance of a linear bandit and the neighboring nodes are assumed to have
similar unknown regression coefficients. The main difference with our setting
is that the rewards of the nodes are independent.

Combinatorial and multi-agent bandits.  Allocating arms to each node of a
graph to then observe a global reward is a combinatorial bandit problem (Cesa-
Bianchi and Lugosi,, 2012), the number of joint arms scaling exponentially
with the number of nodes. This has been extensively studied both in the
regret-based (Chen et al.,, 2013; Perrault et al.,, 2020) and the pure
exploration context (Chen et al.,, 2014; Cao and Krishnamurthy,, 2019; Jourdan
et al.,, 2021; Du et al.,, 2020). Our problem is closer to the one presented
in Amin et al., (2011) and Bargiacchi et al., (2018), where several agents
want to maximize a global team reward that can be decomposed into a sum of
observable local rewards as in a _semi-bandit game_ (Audibert et al.,, 2011;
Chen et al.,, 2013). However, we study a more structured context as we assume
observable bilinear rewards for each edge of the graph. Furthermore, note that
our problem can be solved by the algorithm presented in Du et al., (2020) with
a sample complexity increasing in the number of nodes. On the contrary, we
propose in this paper an algorithm with a sample complexity decreasing in the
number of nodes exploiting the structure of the bilinear reward and the graph.
Finally, most of the algorithms developed for combinatorial bandits assume the
availability of an oracle to solve the combinatorial optimization problem
returning the arm to play or the final best arm recommendation. We make no
such assumption.

##  3 Preliminaries and Notations

Let $\mathcal{G}=(V,E)$ be a directed graph with $V$ the set of nodes, $E$ the
set of edges where we assume that if $(i,j)\in E$ then $(j,i)\in E$, and
$\mathcal{N}(i)$ the set containing the neighbors of a node $i\in V$. We
denote by $n=|V|$ the number of nodes and $m=|E|$ the number of edges. We
define the _graphical bilinear bandit_ on the graph $\mathcal{G}$ as the
setting where a learner sequentially pulls at each round $t$ a joint arm
$\left(x^{(1)}_{t},\dots,x^{(n)}_{t}\right)\in\mathcal{X}^{n}$, also called
graph allocation or simply allocation when it is clear from the context, and
then receives a bilinear reward $r^{(i,j)}_{t}$ for each edge $(i,j)\in E$. At
each round, the joint arm can be constructed simultaneously or sequentially,
however all the bilinear rewards are only revealed after the joint arm has
been pulled.

We denote $K=|\mathcal{X}|$ the number of node-arms and it is assumed that
$\mathcal{X}$ spans $\mathbb{R}^{d}$. For each round $t$ of the learning
procedure and each node $i\in V$, $x^{(i)}_{t}\in\mathcal{X}$ represents the
node-arm allocated to the node $i\in V$. For each edge $(i,j)\in E$, we denote
$z^{(i,j)}_{t}=\operatorname*{vec}{(x^{(i)}_{t}x^{(j)\top}_{t})}\in\mathcal{Z}$
the associated chosen edge-arm.

The goal is to derive an algorithm that minimizes the number of pulled joint
arms required to find the one maximizing the sum of the associated expected
bilinear rewards, for a given confidence level. For the sake of simplicity, we
assume that the unknown parameter matrix $\mathbf{M}_{\star}$ in the bilinear
reward is symmetric. We provide an analysis of the non-symmetric case in
Appendix E.

For any finite set $X$,
$\mathcal{S}_{X}\triangleq\\{\lambda\in[0,1]^{|X|}\text{, }\sum_{x\in
X}\lambda_{x}=1\\}$ denotes the simplex in $\mathbb{R}^{|X|}$. For any vector
$x\in\mathbb{R}^{d}$, $\|x\|$ will denote the $\ell_{2}$-norm of $x$. For any
square matrix $\mathbf{A}\in\mathbb{R}^{d\times d}$, we denote by
$\|\mathbf{A}\|\triangleq\sup_{x:\|x\|=1}\|\mathbf{A}x\|$ the spectral norm of
$\mathbf{A}$. Finally, for any vector $x\in\mathbb{R}^{d}$ and a symmetric
positive-definite matrix $\mathbf{A}\in\mathbb{R}^{d\times d}$, we define
$\|x\|_{\mathbf{A}}\triangleq\sqrt{x^{\top}\mathbf{A}x}$.

##  4 An NP-Hard Problem

In this section, we address the problem of finding the best joint arm given
$\mathbf{M}_{\star}$ or a good estimator $\hat{\mathbf{M}}$. If the best edge-
arm $z_{\star}$ is composed of a single node-arm $x_{\star}$, that is
$z_{\star}=\operatorname*{vec}{(x_{\star}x_{\star}^{\top})}$, then finding the
best joint arm is trivial and the solution is to assign $x_{\star}$ to all
nodes. Conversely, if $z_{\star}$ is composed of two distinct node-arms
$(x_{\star},x_{\star}^{\prime})$, the problem is harder.

The following theorem states that, even with the knowledge of the true
parameter $\mathbf{M}_{\star}$, identifying the best join-arm is NP-Hard with
respect to the number of nodes $n$.

######  Theorem 4.1.

Consider a given matrix $\mathbf{M}_{\star}\in\mathbb{R}^{d\times d}$ and a
finite arm set $\mathcal{X}\subset\mathbb{R}^{d}$. Unless P=NP, there is no
polynomial time algorithm guaranteed to find the optimal solution of

$\displaystyle\max_{\left(x^{(1)},\ldots,x^{(n)}\right)\in\mathcal{X}^{n}}\sum_{(i,j)\in
E}x^{(i)\top}\mathbf{M}_{\star}\ x^{(j)}\enspace.$

The proof of this theorem is in Appendix A and relies on a reduction to the
Max-Cut problem. Hence, no matter which estimate $\hat{\mathbf{M}}$ of
$\mathbf{M}_{\star}$ one can build, the learner is not guaranteed to find in
polynomial time the joint arm
$\left(x^{(1)}_{\star},\dots,x^{(n)}_{\star}\right)$ maximizing the expected
global reward. However, one can notice that, given the matrix
$\mathbf{M}_{\star}$ or even a good enough estimate $\hat{\mathbf{M}}$,
identifying the edge-arm
$z^{\star}=\operatorname*{vec}{(x_{\star}x_{\star}^{\prime\top})}\in\mathcal{Z}$
that maximizes the reward
$z_{\star}^{\top}\operatorname*{vec}{(\mathbf{M}_{\star})}$ requires only
$K^{2}$ reward estimations (we simply estimate all the linear reward
associated to each edge-arm in $\mathcal{Z}$). Thus, instead of looking for
the best joint arm explicitely, we will first identify the best edge-arm
$z^{\star}$, and then allocate $z^{\star}$ to the largest number of edges in
the graph. We will also show that this approach gives a guarantee on its
associated global reward.

Let us consider the graph allocation that places the maximum number of edge-
arms $z_{\star}$ in $\mathcal{G}$. It is easy to show that the subgraph
containing only the edges where $z_{\star}$ has been pulled is the largest
bipartite subgraph included in $\mathcal{G}$. Recall that a graph
$\mathcal{G}^{\prime}=(V^{\prime},E^{\prime})$ is a bipartite if and only if
one can partition the node set $V^{\prime}=(V^{\prime}_{1},V^{\prime}_{2})$
such that

$\displaystyle(i,j)\in E^{\prime}\Rightarrow(i,j)\in V^{\prime}_{1}\times
V^{\prime}_{2}\text{ or }(j,i)\in V^{\prime}_{1}\times
V^{\prime}_{2}\enspace.$

Notice, that if $\mathcal{G}^{\prime}$ is the largest bipartite subgraph in
$\mathcal{G}$, the number of edges in $E^{\prime}$ is the maximal number of
edge-arms $z_{\star}$ that can be allocated with a single graph allocation.

Hence, finding the joint arm with the largest number of edge-arms $z_{\star}$
allocated in the graph is equivalent to finding the largest bipartite subgraph
$\mathcal{G}^{\prime}=(V^{\prime},E^{\prime})$ in $\mathcal{G}$. Once that
subgraph is determined, we just need to allocate to all the nodes in
$V^{\prime}_{1}$ the node-arm $x_{\star}$ and to all the nodes of
$V^{\prime}_{2}$ the node-arm $x_{\star}^{\prime}$ (which is equivalent to
allocating to all the edges in $E^{\prime}$ the edge-arm $z_{\star}$).

Furthermore, we know that every $m$-edge graph contains a bipartite subgraph
of at least $m/2$ edges (Erdos,, 1975). Therefore, we propose Algorithm 1
which iteratively constructs a bipartite subgraph and allocates the nodes
accordingly to create at least $m/2$ edge-arms $z_{\star}$.

Input :

$\mathcal{G}=(V,E)$, $\mathcal{X}$, $\mathbf{M}$

Find
$(x_{\star},x^{\prime}_{\star})\in\mathop{\mathrm{arg\,max}}_{(x,x^{\prime})\in\mathcal{X}^{2}}x^{\top}\mathbf{M}\
x^{\prime}$

Set $V_{1}=\emptyset$, $V_{2}=\emptyset$

for  _$i$ in $V$_ do

Set $n_{1}$ the number of neighbors of $i$ in $V_{1}$ Set $n_{2}$ the number
of neighbors of $i$ in $V_{2}$ if  _$n_{1} >n_{2}$_ then

$x^{(i)}=x^{\prime}_{\star}$ $V_{2}\leftarrow V_{2}\cup\\{i\\}$

else

$x^{(i)}=x_{\star}$ $V_{1}\leftarrow V_{1}\cup\\{i\\}$

end if

end for

return $\mathbf{x}=\left(x^{(1)},\dots,x^{(n)}\right)$

Algorithm 1 Bipartite graph algorithm for Best Arm Identification in Graphical
Bilinear Bandits

The following result gives the guarantee on the global reward associated to
the joint arm returned by Algorithm 1. We refer the reader to Appendix A for
the proof.

######  Theorem 4.2.

Let us consider the graph $\mathcal{G}=(V,E)$, a finite arm set
$\mathcal{X}\subset\mathbb{R}^{d}$ and the matrix $\mathbf{M}_{\star}$ given
as input to Algorithm 1. Then, the expected global reward $r=\sum_{(i,j)\in
E}x^{(i)\top}\mathbf{M}_{\star}x^{(j)}$ associated to the returned allocation
$\mathbf{x}=\left(x^{(1)},\dots,x^{(n)}\right)\in\mathcal{X}^{n}$ verifies:

$\displaystyle\frac{r-r_{\mathrm{min}}}{r_{\star}-r_{\mathrm{min}}}\geq\frac{1}{2}\enspace.$

where $r_{\star}$ and $r_{\mathrm{min}}$ are respectively the highest and
lowest global reward one can obtain with the appropriate joint arm. Finally,
the complexity of the algorithm is in $\mathcal{O}(K^{2}+n)$.

This type of approximation result is sometimes referred to as _differential
approximation_ or _$z$ -approximation_, and is often viewed as a more subtle
analysis than standard approximation ratio. We emphasize that finding a better
ratio than $\frac{1}{2}$ is a very hard task: such a finding would immediately
yield an improved differential approximation ratio for the Max-Cut problem,
which is an opened problem since 2001 (Hassin and Khuller,, 2001).

##  5 Construction of the Estimate $\hat{\mathbf{M}}$

In the previous section, we designed a polynomial time method that computes a
$1/2$-approximation to the NP-Hard problem of finding the best joint arm given
$\mathbf{M}_{\star}$. Notice that $\mathbf{M}_{\star}$ is only used to
identify the best edge-arm $z_{\star}$. Thus, using an estimate
$\hat{\mathbf{M}}$ of $\mathbf{M}_{\star}$ having the following property:

$\displaystyle\mathop{\mathrm{arg\,max}}_{z\in\mathcal{Z}}z^{\top}\operatorname*{vec}{(\hat{\mathbf{M}})}=\mathop{\mathrm{arg\,max}}_{z\in\mathcal{Z}}z^{\top}\operatorname*{vec}{(\mathbf{M}_{\star})}\enspace,$
(2)

would still allow us to identify $z_{\star}$, and would thus give us the same
guarantees.

In this section we tackle the problem of pulling the edge-arms during the
learning procedure such that the estimated unknown parameter verifies (2) in
as few iterations as possible. To do so, we first formalize the objective
related to the linearized version of the problem. Then, we propose an
algorithm reaching the given objective with high probability while satisfying
the edge-arms constraints.

We denote by
$\theta_{\star}=\operatorname*{vec}{\left(\mathbf{M}_{\star}\right)}$ the
parameter of the linearized problem and $\hat{\theta}_{t}$ the Ordinary Least
Squares (OLS) estimate of $\theta_{\star}$ computed with all the data
collected up to round $t$. The empirical gap between two edge-arms $z$ and
$z^{\prime}$ in $\mathcal{Z}$ is denoted
$\hat{\Delta}_{t}\left(z,z^{\prime}\right)\triangleq\left(z-z^{\prime}\right)^{\top}\hat{\theta}_{t}$.

###  5.1 A Constrained G-Allocation

The goal here is to define the optimal sequence
$\left(z_{1},\dots,z_{mt}\right)\in\mathcal{Z}^{mt}$ that should be pulled in
the first $t$ rounds so that (2) is reached as soon as possible. A natural
approach is to rely on classical strategies developed for best arm
identification in linear bandits. Most of the known strategies (see e.g.,
Soare et al., (2014); Xu et al., (2018); Fiez et al., (2019)) are based on a
bound of the gap error
$|(\theta_{\star}-\hat{\theta}_{t})^{\top}(z-z^{\prime})|$ for all
$z,z^{\prime}\in\mathcal{Z}$. This bound is then used to derive a stopping
condition, indicating a sufficient number of rounds $t$ after which the OLS
estimate $\hat{\theta}_{t}$ is precise enough to ensure the identification of
the best edge-arm, with high probability.

Let $\delta\in(0,1)$ and let $\mathbf{A}_{t}=\sum_{i=1}^{mt}z_{i}z_{i}^{\top}$
be the matrix computed with the $mt$ edge-arms constructed during $t$ rounds.
Following the steps of Soare et al., (2014), we can show that if there exists
$z\in\mathcal{Z}$ such that for all $z^{\prime}\in\mathcal{Z}$ the following
holds:

$\displaystyle\|z-z^{\prime}\|_{\mathbf{A}_{t}^{-1}}\sqrt{8\sigma^{2}\log{\left(\frac{6m^{2}t^{2}K^{4}}{\delta\pi^{2}}\right)}}\leq\hat{\Delta}_{t}\left(z,z^{\prime}\right)\enspace,$
(3)

then with probability at least $1-\delta$, the OLS estimate $\hat{\theta}_{t}$
leads to the best edge-arm. Details of the derivation are given in Appendix B.

As mentioned in Soare et al., (2014), by noticing that
$\max_{(z,z^{\prime})\in\mathcal{Z}^{2}}\left\|z-z^{\prime}\right\|_{\mathbf{A}_{t}^{-1}}\leq
2\max_{z^{\in}\mathcal{Z}}\left\|z\right\|_{\mathbf{A}_{t}^{-1}}$, an
admissible strategy is to pull edge-arms minimizing
$\max_{z\in\mathcal{Z}}\left\|z\right\|_{\mathbf{A}_{t}^{-1}}$ in order to
satisfy the stopping condition as soon as possible. More formally, one wants
to find the sequence of edge-arms
$\mathbf{z}_{mt}^{\star}=\left(z_{1}^{\star},\dots,z_{mt}^{\star}\right)$ such
that:

$\displaystyle\mathbf{z}_{mt}^{\star}\in\mathop{\mathrm{arg\,min}}_{\left(z_{1},\dots,z_{mt}\right)}\max_{z^{\prime}\in\mathcal{Z}}\;{z^{\prime}}^{\top}\left(\sum_{i=1}^{mt}z_{i}z_{i}^{\top}\right)^{-1}{z^{\prime}}\enspace.$
(G-opt-$\mathcal{Z}$)

This is known as _G-allocation_ (see e.g., Pukelsheim, (2006); Soare et al.,
(2014)) and is NP-hard to compute (Çivril and Magdon-Ismail,, 2009; Welch,,
1982). One way to find an approximate solution is to rely on a convex
relaxation of the optimization problem (G-opt-$\mathcal{Z}$) and first compute
a real-valued allocation $\lambda^{\star}\in\mathcal{S}_{\mathcal{Z}}$ such
that

$\displaystyle\lambda^{\star}\in\mathop{\mathrm{arg\,min}}_{\lambda\in\mathcal{S}_{\mathcal{Z}}}\max_{z^{\prime}\in\mathcal{Z}}{z^{\prime}}^{\top}\left(\sum_{z\in\mathcal{Z}}\lambda_{z}zz^{\top}\right)^{-1}z^{\prime}\enspace.$
(G-relaxed-$\mathcal{Z}$)

One could either use random sampling to draw edge-arms as i.i.d. samples from
the $\lambda^{\star}$ distribution or rounding procedures to efficiently
convert each $\lambda^{\star}_{z}$ into an integer. However, these methods do
not take into account the graphical structure of the problem, and at a given
round, the $m$ chosen edge-arms may result in two different assignments for
the same node. Therefore, random sampling or rounding procedures cannot be
straightforwardly used to select edge-arms in $\mathcal{Z}$. Nevertheless,
(G-relaxed-$\mathcal{Z}$) still gives a valuable information on the number of
times, in proportion, each edge-arm $z\in\mathcal{Z}$ must be allocated to the
graph. In the next section, we present an algorithm satisfying both the
proportion requirements and the graphical constraints.

###  5.2 Random Allocation over the Nodes

Our algorithm is based on a randomized method directly allocating node-arms to
the nodes and thus avoiding the difficult task of choosing edge-arms and
trying to allocate them to the graph while ensuring that every node has an
unique assignment. The validity of this random allocation is based on Theorem
5.1 below showing that one can draw node-arms in $\mathcal{X}$ and allocate
them to the graph such that the associated edge-arms follow the probability
distribution $\lambda^{\star}$ solution of (G-relaxed-$\mathcal{Z}$).

######  Theorem 5.1.

Let $\mu^{\star}$ be a solution of the following optimization problem:

$\displaystyle\min_{\mu\in\mathcal{S}_{\mathcal{X}}}\max_{x^{\prime}\in\mathcal{X}}{x^{\prime}}^{\top}\left(\sum_{x\in\mathcal{X}}\mu_{x}xx^{\top}\right)^{-1}x^{\prime}\enspace.$
(G-relaxed-$\mathcal{X}$)

Let $\lambda^{\star}\in\mathcal{S}_{\mathcal{Z}}$ be defined for all
$z=\operatorname*{vec}{\left(xx^{\prime\top}\right)}\in\mathcal{Z}$ by
$\lambda_{z}^{\star}=\mu^{\star}_{x}\mu^{\star}_{x^{\prime}}$. Then,
$\lambda^{\star}$ is a solution of (G-relaxed-$\mathcal{Z}$).

#### Sketch of proof.

The objective at the optimum in (G-relaxed-$\mathcal{X}$) and
(G-relaxed-$\mathcal{Z}$) are respectively equal to $d$ and $d^{2}$ which is
the dimension of their respective problem, a result known as the Equivalence
Theorem (Kiefer and Wolfowitz,, 1960). Thus, by multiplying the optimum value
of $\eqref{eq:relaxed_G_allocation_X}$ by itself, we can show that for all
$z\in\mathcal{Z}$ where $z=\operatorname*{vec}{\left(xx^{\prime\top}\right)}$
with $(x,x^{\prime})\in\mathcal{X}^{2}$, $\lambda^{\star}_{z}$ can be written
as the product $\mu_{x}^{\star}\mu_{x^{\prime}}^{\star}$. We refer to the
Appendix C for the detailed proof.

This theorem implies that, at each round $t>0$ and each node $i\in V$, if
$x_{t}^{(i)}$ is drawn from $\mu^{\star}$, then for all pairs of neighbors
$(i,j)\in E$ the probability distribution of the associated edge-arms
$z^{(i,j)}_{t}$ follows $\lambda^{\star}$. Moreover, as $\mu^{\star}$ is a
distribution over the node-arm set $\mathcal{X}$, $\lambda^{\star}$ is a joint
(product) probability distribution on $\mathcal{X}^{2}$ with marginal
$\mu^{\star}$.

We apply the Frank-Wolfe algorithm (Frank et al.,, 1956) to compute the
solution $\mu^{\star}$ of (G-relaxed-$\mathcal{X}$), as it is more suited to
optimization tasks on the simplex than projected gradient descent. Although we
face a min-max optimization problem, we notice that the function
$h(\mu)=\max_{x^{\prime}\in\mathcal{X}}x^{\prime\top}\left(\sum_{x\in\mathcal{X}}\mu_{x}xx^{\top}\right)^{-1}x^{\prime}$
is convex. We refer the reader to Appendix F and references therein for a
proof on the convexity of $h$ and a discussion about using Frank-Wolfe for
solving (G-relaxed-$\mathcal{X}$).

Given the characterization in Theorem 5.1 and our objective to verify the
stopping condition in (3), we present our sampling procedure in Algorithm 2.
We also note that at each round the sampling of the node-arms can be done in
parallel.

Input :

graph $\mathcal{G}=(V,E)$, arm set $\mathcal{X}$

Set $A_{0}=I$ ; $b_{0}=0$ ; $t=1$;

Apply the Frank-Wolfe algorithm to find $\mu^{\star}$ solution of
(G-relaxed-$\mathcal{X}$).

while  _stopping condition ( 3) is not verified_ do

// Sampling the node-arms  Draw
$x_{t}^{(1)},\ldots,x_{t}^{(n)}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\mu^{\star}$
and obtain for all $(i,j)$ in $E$ the rewards $r_{t}^{(i,j)}$;  // Estimating
$\hat{\theta}_{t}$ with the associated edge-arms
$\mathbf{A}_{t}=\mathbf{A}_{t-1}+\sum_{(i,j)\in
E}z_{t}^{(i,j)}z_{t}^{(i,j)\top}$; $b_{t}=b_{t-1}+\sum_{(i,j)\in
E}z_{t}^{(i,j)}r_{t}^{(i,j)}$; $\hat{\theta}_{t}=\mathbf{A}_{t}^{-1}b_{t}$
$t\leftarrow t+1$;

end while

return $\hat{\theta}_{t}$

Algorithm 2 Randomized G-Allocation strategy for Graphical Bilinear Bandits

This sampling procedure implies that each edge-arm follows the optimal
distribution $\lambda^{\star}$. However, if we take the number of times each
$z\in\mathcal{Z}$ appears in the $m$ pulled edge-arms of a given round, we
might notice that the observed proportion is not close to
$\lambda^{\star}_{z}$, regardless of the size of $m$. This is due to the fact
that the $m$ edge-arms are not independent because of the graph structure (cf.
Section 6). Conversely, since each group of $m$ edge-arms are independent from
one round to another, the proportion of each $z\in\mathcal{Z}$ observed among
the $mt$ pulled edge-arms throughout $t$ rounds is close to
$\lambda^{\star}_{z}$.

One may wonder if deterministic rounding procedures could be used instead of
random sampling on $\mu^{\star}$, as it is done in many standard linear bandit
algorithms (Soare et al.,, 2014; Fiez et al.,, 2019). Applying rounding
procedure on $\mu^{\star}$ gives the number of times each node-arm
$x\in\mathcal{X}$ should be allocated to the graph. However, it does not
provide the actual allocations that the learner must choose over the $t$
rounds to optimally pull the associated edge-arms (_i.e.,_ pull edge-arms
following $\lambda^{\star}$). Thus, although rounding procedures give a more
precise number of times each node-arm should be pulled, the problem of
allocating them to the graph remains open, whereas by concentration of the
measure, randomized sampling methods imply that the associated edge-arms
follow the optimal probability distribution $\lambda^{\star}$. In this paper,
we present a simple and standard randomized G-allocation strategy, but other
more elaborated methods could be considered, as long as they include the
necessary randomness.

On the choice of the G-allocation problem.  We have considered the
G-allocation optimization problem (G-opt-$\mathcal{Z}$), however, one could
want to directly minimize
$\max_{(z,z^{\prime})\in\mathcal{Z}^{2}}\left\|z-z^{\prime}\right\|_{\mathbf{A}_{t}^{-1}}$,
known as the XY-allocation (Soare et al.,, 2014; Fiez et al.,, 2019). Hence,
one may want to construct edge-arms that follow the distribution
$\lambda^{\star}_{\mathrm{XY}}$ solution of the relaxed XY-allocation problem:

$\displaystyle\min_{\lambda}\max_{z^{\prime},z^{\prime\prime}}\;\left(z^{\prime}-z^{\prime\prime}\right)^{\top}\left(\sum_{z\in\mathcal{Z}}\lambda_{z}zz^{\top}\right)^{-1}\left(z^{\prime}-z^{\prime\prime}\right)\enspace.$

Although efficient in the linear case, this approach outputs a distribution
$\lambda^{\star}_{\mathrm{XY}}$ which is not a joint probability distribution
of two independent random variables, and so cannot be decomposed as the
product of its marginals. Hence, there is no algorithm that allocates
_identically_ and _independently_ the nodes of the graph to create edge-arms
following $\lambda^{\star}_{\mathrm{XY}}$. Thus, we will rather deal with the
upper bound given by the G-allocation as it allows sampling over the nodes.

Static design versus adaptive design. Adaptive designs as proposed for example
in Soare et al., (2014) and Fiez et al., (2019) provide a strong improvement
over static designs in the case of linear bandits. In our particular setting
however, it is crucial to be able to adapt the edge-arms sampling rule to the
node-arms, which is possible thanks to Theorem 5.1. This result requires a set
of edge-arms $\mathcal{Z}$ expressed as a product of node-arms set
$\mathcal{X}$. Extending the adaptive design of Fiez et al., (2019) to our
setting would eliminate edge-arms from $\mathcal{Z}$ at each phase, without
trivial guarantees that the newly obtained edge-arms set
$\mathcal{Z}^{\prime}\subset\mathcal{Z}$ could still be derived from another
node-arms set $\mathcal{X}^{\prime}\subset\mathcal{X}$. An adaptive approach
is definitely a natural and promising extension of our method, and is left for
future work.

###  5.3 Convergence Analysis

We now prove the validity of the random sampling procedure detailed in
Algorithm 2 by controlling the quality of the approximation
$\max_{z\in\mathcal{Z}}z^{\top}\mathbf{A}_{t}^{-1}z$ with respect to the
optimum of the G-allocation optimization problem
$\max_{z^{\prime}\in\mathcal{Z}}{z^{\prime}}^{\top}\left(\sum_{i=1}^{mt}z_{i}^{\star}z_{i}^{\star\top}\right)^{-1}z^{\prime}$
described in (G-opt-$\mathcal{Z}$). As is usually done in the optimal design
literature (see e.g., Pukelsheim, (2006); Soare et al., (2014); Sagnol,
(2010)) we bound the relative error $\alpha$:

$\max_{z\in\mathcal{Z}}z^{\top}\mathbf{A}_{t}^{-1}z\leq(1+\alpha)\max_{z^{\prime}\in\mathcal{Z}}{z^{\prime}}^{\top}\left(\sum_{i=1}^{mt}z_{i}^{\star}z_{i}^{\star\top}\right)^{-1}z^{\prime}\enspace.$

Our analysis relies on several results from matrix concentration theory. One
may refer for instance to Tropp et al., (2015) and references therein for an
extended introduction on that matter. We first introduce a few additional
notations.

Let $f_{\mathcal{Z}}$ be the function such that for any non-singular matrix
$\mathbf{Q}\in\mathbb{R}^{d^{2}\times d^{2}}$,
$f_{\mathcal{Z}}(\mathbf{Q})=\max_{z\in\mathcal{Z}}z^{\top}\mathbf{Q}^{-1}z$
and for any distribution $\lambda\in\mathcal{S}_{\mathcal{Z}}$ let
$\Sigma_{\mathcal{Z}}(\lambda)\triangleq\sum_{z\in\mathcal{Z}}\lambda_{z}zz^{\top}$
be the associated covariance matrix. Finally let
$\mathbf{A}^{\star}_{t}=\sum_{i=1}^{mt}z_{i}^{\star}z_{i}^{\star\top}$ be the
G-optimal design matrix constructed during $t$ rounds.

For $i\in\\{1,\ldots,n\\}$ and $s\in\\{1,\ldots,t\\}$, let $X^{(i)}_{s}$ be
i.i.d. random vectors in $\mathcal{X}$ such that for all $x\in\mathcal{X}$,

$\mathbb{P}\left(X_{1}^{(1)}=x\right)=\mu^{\star}_{x}\enspace.$

Each $X^{(i)}_{s}$ is to be viewed as the random arm pulled at round $s$ for
the node $i$. Using this notation, the random design matrix $\mathbf{A}_{t}$
can be defined as

$\displaystyle\mathbf{A}_{t}=\sum_{s=1}^{t}\sum_{(i,j)\in
E}\operatorname*{vec}{\left(X^{(i)}_{s}X^{(j)\top}_{s}\right)}\operatorname*{vec}{\left(X^{(i)}_{s}X^{(j)\top}_{s}\right)}^{\top}\enspace.$

One can first observe that $f_{\mathcal{Z}}(\mathbf{A}_{t})$ can be bounded by
the following quantity:

$\displaystyle f_{\mathcal{Z}}(\mathbf{A}_{t})$
$\displaystyle=\max_{z\in\mathcal{Z}}\;z^{\top}\left(\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}+\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right)z$
$\displaystyle\leq\max_{z\in\mathcal{Z}}\;z^{\top}\left(\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right)z$
$\displaystyle\phantom{\leq}+f_{\mathcal{Z}}\left(mt\Sigma_{\mathcal{Z}}(\lambda^{\star})\right)$
$\displaystyle\leq\max_{z\in\mathcal{Z}}\;\|z\|^{2}\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|$
$\displaystyle\phantom{\leq}+f_{\mathcal{Z}}\left(\mathbf{A}^{\star}_{t}\right).$

Hence, one needs a bound on the maximum eigenvalue of
$\mathbf{A}_{t}^{-1}-(\mathbb{E}\mathbf{A}_{t})^{-1}$. Simple linear algebra
leads to:

$\mathbf{A}_{t}^{-1}-(\mathbb{E}\mathbf{A}_{t})^{-1}=\mathbf{A}_{t}^{-1}(\mathbb{E}\mathbf{A}_{t}-\mathbf{A}_{t})(\mathbb{E}\mathbf{A}_{t})^{-1}.$

Thus, in addition to bounding the maximum eigenvalue of $\mathbf{A}_{t}^{-1}$,
which is equal to the minimum eigenvalue of $\mathbf{A}_{t}$, we need a bound
on $\|\mathbf{A}_{t}-\mathbb{E}\mathbf{A}_{t}\|$. It may be derived from
concentration results on sum of random matrices derived in Tropp et al.,
(2015). We now state the result controlling the relative error obtained with
our randomized sampling allocation. The proof can be found in the Appendix C.

######  Theorem 5.2.

Let $\lambda^{\star}$ be a solution of the optimization problem
(G-relaxed-$\mathcal{Z}$). Let $0\leq\delta\leq 1$ and let $t_{0}$ be such
that

$t_{0}=2Ld^{2}\log(2d^{2}/\delta)/\nu_{\mathrm{min}},$

where $L=\max_{z\in\mathcal{Z}}\|z\|^{2}$ and $\nu_{\mathrm{min}}$ is the
smallest eigenvalue of the covariance
matrix$\frac{1}{K^{2}}\sum_{z\in\mathcal{Z}}zz^{\top}$. Then, at each round
$t\geq t_{0}$ with probability at least $1-\delta$, the randomized
G-allocation strategy for graphical bilinear bandit in Algorithm 2 produces a
matrix $\mathbf{A}_{t}$ such that:

$\displaystyle
f_{\mathcal{Z}}(\mathbf{A}_{t})\leq(1+\alpha)f_{\mathcal{Z}}(\mathbf{A}^{\star}_{t})$

where

$\displaystyle\alpha=\frac{Ld^{2}}{m\nu_{\mathrm{min}}^{2}}\sqrt{\frac{2v}{t}\log\left(\frac{2d^{2}}{\delta}\right)}+o\left(\frac{1}{\sqrt{t}}\right),$

and
$v\triangleq\mathbb{E}\big{[}(\mathbf{A}_{1}-\mathbb{E}\mathbf{A}_{1})^{2}\big{]}$.

We have just shown that the approximation value
$\max_{z\in\mathcal{Z}}z^{\top}\mathbf{A}_{t}^{-1}z$ converges to the optimal
value with a rate of $O\left(\sqrt{v}/(m\sqrt{t})\right)$. In Section 6, we
show that the best case graph implies a $v=O\left(m\right)$ matching the
convergence rate $O\left(1/\sqrt{mt}\right)$ of a linear bandit algorithm
using randomized sampling to pull $mt$ edge-arms without (graphical)
constraints. Moreover, we will see that the worst case graph implies that
$v=O\left(m^{2}\right)$.

Since we filled the gap between our constraint objective and the problem of
best arm identification in linear bandits, thanks to Theorem 5.1 and 5.2, we
are able to extend known results for best arm identification in linear bandits
on the sample complexity and its associated lower bound.

######  Corollary 5.3 (Soare et al., (2014), Theorem 1).

If the G-allocation is implemented with the random strategy of Algorithm 2,
resulting in an $\alpha$-approximation, then with probability at least
$1-\delta$, the best arm obtained with $\hat{\theta}_{t}$ is $z_{\star}$ and

$\displaystyle
t\leq\frac{128\sigma^{2}d^{2}(1+\alpha)\log{\left(\frac{6m^{2}t^{2}K^{4}}{\delta\pi}\right)}}{m\Delta_{\mathrm{min}}^{2}}\enspace,$

where
$\Delta_{\mathrm{min}}=\min_{z\in\mathcal{Z}\setminus\\{z_{\star}\\}}(z_{\star}-z)^{\top}\theta_{\star}$.

Moreover, let $\tau$ be the number of rounds sufficient for any algorithm to
determine the best arm with probability at least $1-\delta$. A lower bound on
the expectation of $\tau$ can be obtained from the one derived for the problem
of best arm identification in linear bandits (see e.g., Theorem 1 in Fiez et
al., (2019)):

$\mathbb{E}[\tau]\geq\min_{\lambda\in\mathcal{S}_{\mathcal{Z}}}\max_{z\in\mathcal{Z}\setminus\\{z_{\star}\\}}\log\left(\frac{1}{2.4\delta}\right)\frac{2\sigma^{2}\|z_{\star}-z\|^{2}_{\Sigma_{\mathcal{Z}}(\lambda)^{-1}}}{m\left(\left(z_{\star}-z\right)^{\top}\theta_{\star}\right)^{2}}\enspace.$

As observed in Soare et al., (2014) this lower bound can be upper bounded, in
the worst case, by $4\sigma^{2}d^{2}/(m\Delta_{\mathrm{min}}^{2})$ which
matches our bound up to log terms and the relative error $\alpha$.

##  6 Influence of the Graph Structure on $v$

The convergence bound in Theorem 5.2 depends on
$v=\mathbb{E}\big{[}(\mathbf{A}_{1}-\mathbb{E}\mathbf{A}_{1})^{2}\big{]}$. In
this section, we characterize the impact of the graph structure on this
quantity and, by extension, on the convergence rate. First of all, recall that

$\displaystyle\mathbf{A}_{1}=\sum_{(i,j)\in
E}\operatorname*{vec}{\left(X^{(i)}_{1}X^{(j)\top}_{1}\right)}\operatorname*{vec}{\left(X^{(i)}_{1}X^{(j)\top}_{1}\right)}^{\top}\enspace.$

Let denote
$\mathbf{A}_{1}^{(i,j)}=\operatorname*{vec}{(X^{(i)}_{1}X^{(j)\top}_{1})}\operatorname*{vec}{(X^{(i)}_{1}X^{(j)\top}_{1})}^{\top}$
such that $\mathbf{A}_{1}=\sum_{(i,j)\in E}\mathbf{A}_{1}^{(i,j)}$ and let
define for any random matrices $\mathbf{A}$ and $\mathbf{B}$ the operators
$\mathrm{Var}(\mathbf{A})\triangleq{\mathbb{E}\big{[}(\mathbf{A}-\mathbb{E}[\mathbf{A}])^{2}\big{]}}$
and
$\mathrm{Cov}(\mathbf{A},\mathbf{B})\triangleq{\mathbb{E}\big{[}(\mathbf{A}-\mathbb{E}[\mathbf{A}])(\mathbf{B}-\mathbb{E}[\mathbf{B}])\big{]}}$.
We can derive the variance of $\mathbf{A}_{1}$ as follows:

$\displaystyle\mathrm{Var}\left(\mathbf{A}_{1}\right)=$
$\displaystyle\sum_{(i,j)\in
E}\mathrm{Var}\left(\mathbf{A}^{(i,j)}_{1}\right)$
$\displaystyle+\sum_{(i,j)\in E}\sum_{\begin{subarray}{c}(k,l)\in E\\\
(k,l)\neq(i,j)\end{subarray}}\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(k,l)}_{1}).$

One can decompose the sum of the covariances into three groups: a first group
where $k\neq i,j$ and $l\neq i,j$ which means that the two edges do not share
any node and
$\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(k,l)}_{1})=\mathbf{0}$, and
two other groups where the edges share at least one node. For all edges
$(i,j)\in E$ we consider either the edges $(i,k)\in E$ where $k\neq j$,
yielding $\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(i,k)}_{1})$ or the
edges $(j,k)\in E$, yielding
$\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(j,k)}_{1})$.

(a)

(b)

Figure 1: Number of rounds $t$ needed to verify the stopping condition (3)
with respect to left: the number of edges $m$ where the dimension of the edge-
arm space $\mathcal{Z}$ is fixed and equal to $25$ and right: the dimension of
the edge-arm space $\mathcal{Z}$ where the number of edges is fixed and equal
to $156$. For both experiments we run 100 times and plot the average number of
rounds needed to verify the stopping condition.

Hence, one has

$\displaystyle\mathrm{Var}\left(\mathbf{A}_{1}\right)=$
$\displaystyle\sum_{(i,j)\in
E}\mathrm{Var}\left(\mathbf{A}^{(i,j)}_{1}\right)$
$\displaystyle+\sum_{i=1}^{n}\sum_{j\in\mathcal{N}(i)}\sum_{\begin{subarray}{c}k\in\mathcal{N}(i)\\\
k\neq
j\end{subarray}}\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(i,k)}_{1})$
$\displaystyle+\sum_{i=1}^{n}\sum_{j\in\mathcal{N}(i)}\sum_{k\in\mathcal{N}(j)}\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(j,k)}_{1})\enspace.$

Let $P\geq 0$ be such that for all $(i,j)\in E$,
$\mathrm{Var}\left(\mathbf{A}^{(i,j)}_{1}\right)\preceq P\times\mathbf{I}$ and
$M,N\geq 0$ such that for all $(i,j)\in E$:

$\displaystyle\forall
k\in\mathcal{N}(i),\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(i,k)}_{1})\preceq
M\times\mathbf{I}$ $\displaystyle\forall
k\in\mathcal{N}(j),\mathrm{Cov}(\mathbf{A}^{(i,j)}_{1},\mathbf{A}^{(j,k)}_{1})\preceq
N\times\mathbf{I}$

We want to compare the quantity $\|\mathrm{Var}(\mathbf{A}_{1})\|$ for
different types of graphs: star, complete, circle and a matching graph. To
have a fair comparison, we want graphs that reveal the same number of rewards
at each round of the learning procedure. Hence, we denote respectively
$n_{\mathrm{S}}$, $n_{\mathrm{Co}}$, $n_{\mathrm{Ci}}$ and $n_{\mathrm{M}}$
the number of nodes in a star, complete, circle and matching graph of $m$
edges and get:

Star graph:

$\displaystyle\|\mathrm{Var}(\mathbf{A}_{1})\|\leq
mP+n_{\mathrm{S}}^{2}(M+N).$

Complete graph:

$\displaystyle\|\mathrm{Var}(\mathbf{A}_{1})\|\leq
mP+n_{\mathrm{Co}}^{3}(M+N).$

Circle graph:

$\displaystyle\|\mathrm{Var}(\mathbf{A}_{1})\|\leq mP+n_{\mathrm{Ci}}(2M+4N).$

Matching graph:

$\displaystyle\|\mathrm{Var}(\mathbf{A}_{1})\|\leq mP+n_{\mathrm{M}}N.$

We refer the reader to Appendix D for more details on the given upper bounds.
Since the star (respectively, complete, circle and matching) graph of $m$
edges has a number of nodes $n_{\mathrm{S}}=m/2+1$ (respectively
$n_{\mathrm{Co}}=\left(1+\sqrt{4m+1}\right)/2$, $n_{\mathrm{Ci}}=m/2$ and
$n_{\mathrm{M}}=m$), we obtain the bounds stated in Table 1.

Graph | Upper bound on $\|\mathrm{Var}{\left(\mathbf{A}_{1}\right)}\|$ | $\alpha$
---|---|---
Star | $mP+(M+N)O\left(m^{2}\right)$ | $O\left(1/\sqrt{t}\right)$
Complete | $mP+(M+N)O\left(m\sqrt{m}\right)$ | $O\left(1/\left(m^{\frac{1}{4}}\sqrt{t}\right)\right)$
Circle | $mP+(M+N)O\left(m\right)$ | $O\left(1/\sqrt{mt}\right)$
Matching | $mP+mN$ | $O\left(1/\sqrt{mt}\right)$
Table 1: Upper bound on the variance and convergence rate of Algorithm 2 for
the star, complete, circle and matching graph with respect to the number of
edges $m$ and the number of rounds $t$.

These four examples evidence the strong dependency of the variance on the
structure of the graph. The more independent the edges are (i.e., with no
common nodes), the smaller the quantity $\|\mathrm{Var}(\mathbf{A}_{1})\|$ is.
For a fixed number of edges $m$, the best case is the matching graph where no
edge share the same node and the worst case is the star graph where all the
edges share a central node.

##  7 Experiments

In this section, we consider the modified version of a standard experiment
introduced by Soare et al., (2014) and used in most papers on best arm
identification in linear bandits (Xu et al.,, 2018; Tao et al.,, 2018; Fiez et
al.,, 2019; Zaki et al.,, 2019) to evaluate the sample complexity of our
algorithm on different graphs. We consider $d+1$ node-arms in
$\mathcal{X}\subset\mathbb{R}^{d}$ where $d\geq 2$. This node-arm set is made
of the $d$ vectors $(\mathbf{e}_{1},\dots,\mathbf{e}_{d})$ forming the
canonical basis of $\mathbb{R}^{d}$ and one additional arm
$x_{d+1}=(\cos(\omega),\sin(\omega),0,\dots,0)^{\top}$ with
$\omega\in]0,\pi/2]$. Note that by construction, the edge-arm set
$\mathcal{Z}$ contains the canonical basis
$(\mathbf{e}^{\prime}_{1},\dots,\mathbf{e}^{\prime}_{d^{2}})$ of
$\mathbb{R}^{d^{2}}$. The parameter matrix $\mathbf{M}_{\star}$ has its first
coordinate equal to 2 and the others equal to 0 which makes
$\theta_{\star}=\operatorname*{vec}{\left(\mathbf{M}_{\star}\right)}=\left(2,0,\dots,0\right)^{\top}\in\mathbb{R}^{d^{2}}$.
The best edge-arm is thus $z_{\star}=z^{(1,1)}=\mathbf{e}^{\prime}_{1}$. One
can note that when $\omega$ tends to 0, it is harder to differentiate this arm
from
$z^{(d+1,d+1)}=\operatorname*{vec}{\left(x_{(d+1)}x_{(d+1)}^{\top}\right)}$than
from the other arms. We set $\eta^{(i,j)}_{t}\sim\mathcal{N}(0,1)$, for all
edges $(i,j)$ and round $t$.

We consider the two cases where $\omega=0.1$ which makes the edge-arms
$z^{(1,1)}$ and $z^{(d+1,d+1)}$ difficult to differentiate, and $\omega=\pi/2$
which makes the edge-arm $z^{(1,1)}$ easily identifiable as the optimal edge-
arm. For each of these two cases, we evaluate the influence of the graph
structure, the number of edges $m$ and the edge-arm space dimension $d^{2}$ on
the sampling complexity. Results are shown in Figure 1.

When $\omega=0.1$, the type of the graph does not impact the number of rounds
needed to verify the stopping condition. This is mainly due to the fact that
the magnitude of its associated variance is negligible with respect to the
number of rounds. Hence, even if we vary the number of edges or the dimension,
we get the same performance for any type of graph including the matching
graph. This implies that our algorithm performs as well as a linear bandit
that draws $m$ edge-arms in parallel at each round. When $\omega=\pi/2$, the
number of rounds needed to verify the stopping condition is smaller and the
magnitude of the variance is no longer negligible. Indeed, when the number of
edges or the dimension increases, we notice that the star graph takes more
times to satisfy the stopping condition. Moreover, note that the sample
complexities obtained for the circle and the matching graph are similar. This
observation is in line with the dependency on the variance shown in Table 1.

##  8 Conclusion

We introduced a new graphical bilinear bandit setting and studied the best arm
identification problem with a fixed confidence. This problem being NP-Hard
even with the knowledge of the true parameter matrix $\mathbf{M}^{\star}$, we
first proposed an algorithm that provides a 1/2-approximation. Then, we
provided a second algorithm, based on G-allocation strategy, that uses
randomized sampling over the nodes to return a good estimate
$\hat{\mathbf{M}}$ that can be used instead of $\mathbf{M}_{\star}$. Finally,
we highlighted the impact of the graph structure on the convergence rate of
our algorithm and validated our theoretical results with experiments.
Promising extensions of the model include considering unknown parameters
$\mathbf{M}_{\star}^{(i,j)}$, different for each edge $(i,j)$ of the graph,
and investigating XY-allocation strategies.

## Acknowledgements

We thank anonymous reviewers, whose comments helped us improve the paper
significantly.

## References

  * Abbasi-Yadkori et al., (2011) Abbasi-Yadkori, Y., Pál, D., and Szepesvári, C. (2011).  Improved algorithms for linear stochastic bandits.  In Advances in Neural Information Processing Systems, pages 2312–2320.
  * Amin et al., (2011) Amin, K., Kearns, M., and Syed, U. (2011).  Graphical models for bandit problems.  In Proceedings of the Twenty-Seventh Conference on Uncertainty in Artificial Intelligence, page 1–10.
  * Audibert and Bubeck, (2010) Audibert, J.-Y. and Bubeck, S. (2010).  Best arm identification in multi-armed bandits.  In Proceedings of the 23th Annual Conference on Learning Theory, pages 41–53.
  * Audibert et al., (2011) Audibert, J.-Y., Bubeck, S., and Lugosi, G. (2011).  Minimax policies for combinatorial prediction games.  In Proceedings of the 24th Annual Conference on Learning Theory, pages 107–132.
  * Bargiacchi et al., (2018) Bargiacchi, E., Verstraeten, T., Roijers, D., Nowé, A., and Hasselt, H. (2018).  Learning to coordinate with coordination graphs in repeated single-stage multi-agent decision problems.  In International conference on machine learning, pages 482–490.
  * Bubeck et al., (2009) Bubeck, S., Munos, R., and Stoltz, G. (2009).  Pure exploration in multi-armed bandits problems.  In International conference on Algorithmic learning theory, pages 23–37. Springer.
  * Cao and Krishnamurthy, (2019) Cao, T. and Krishnamurthy, A. (2019).  Disagreement-based combinatorial pure exploration: Sample complexity bounds and an efficient algorithm.  In Proceedings of the Thirty-Second Conference on Learning Theory, volume 99, pages 558–588.
  * Cesa-Bianchi et al., (2013) Cesa-Bianchi, N., Gentile, C., and Zappella, G. (2013).  A gang of bandits.  In Advances in Neural Information Processing Systems, pages 737–745.
  * Cesa-Bianchi and Lugosi, (2012) Cesa-Bianchi, N. and Lugosi, G. (2012).  Combinatorial bandits.  Journal of Computer and System Sciences, 78(5):1404 – 1422.
  * Chen et al., (2014) Chen, S., Lin, T., King, I., Lyu, M. R., and Chen, W. (2014).  Combinatorial pure exploration of multi-armed bandits.  In Advances in Neural Information Processing Systems, volume 27, pages 379–387.
  * Chen et al., (2013) Chen, W., Wang, Y., and Yuan, Y. (2013).  Combinatorial multi-armed bandit: General framework and applications.  In International Conference on Machine Learning, pages 151–159.
  * Damla Ahipasaoglu et al., (2008) Damla Ahipasaoglu, S., Sun, P., and Todd, M. J. (2008).  Linear convergence of a modified frank–wolfe algorithm for computing minimum-volume enclosing ellipsoids.  Optimisation Methods and Software, 23(1):5–19.
  * Degenne et al., (2020) Degenne, R., Ménard, P., Shang, X., and Valko, M. (2020).  Gamification of pure exploration for linear bandits.  arXiv preprint arXiv:2007.00953.
  * Du et al., (2020) Du, Y., Kuroki, Y., and Chen, W. (2020).  Combinatorial pure exploration with full-bandit or partial linear feedback.  arXiv e-prints, pages arXiv–2006.
  * Erdos, (1975) Erdos, P. (1975).  Problems and results on finite and infinite graphs.  In Recent advances in graph theory (Proc. Second Czechoslovak Sympos., Prague, 1974), pages 183–192.
  * Fiez et al., (2019) Fiez, T., Jain, L., Jamieson, K. G., and Ratliff, L. (2019).  Sequential experimental design for transductive linear bandits.  In Advances in Neural Information Processing Systems, pages 10667–10677.
  * Frank et al., (1956) Frank, M., Wolfe, P., et al. (1956).  An algorithm for quadratic programming.  Naval research logistics quarterly, 3(1-2):95–110.
  * Guestrin et al., (2002) Guestrin, C., Lagoudakis, M. G., and Parr, R. (2002).  Coordinated reinforcement learning.  In Proceedings of the Nineteenth International Conference on Machine Learning, page 227–234.
  * Hassin and Khuller, (2001) Hassin, R. and Khuller, S. (2001).  z-approximations.  Journal of Algorithms, 41(2):429–442.
  * Jedra and Proutiere, (2020) Jedra, Y. and Proutiere, A. (2020).  Optimal best-arm identification in linear bandits.  arXiv preprint arXiv:2006.16073.
  * Jourdan et al., (2021) Jourdan, M., Mutý, M., Kirschner, J., and Krause, A. (2021).  Efficient pure exploration for combinatorial bandits with semi-bandit feedback.  In Proceedings of the 31st International Conference on Algorithmic Learning Theory.
  * Jun et al., (2019) Jun, K.-S., Willett, R., Wright, S., and Nowak, R. (2019).  Bilinear bandits with low-rank structure.  In International Conference on Machine Learning, pages 3163–3172.
  * Kazerouni and Wein, (2019) Kazerouni, A. and Wein, L. M. (2019).  Best arm identification in generalized linear bandits.  arXiv preprint arXiv:1905.08224.
  * Kiefer and Wolfowitz, (1960) Kiefer, J. and Wolfowitz, J. (1960).  The equivalence of two extremum problems.  Canadian Journal of Mathematics, 12:363–366.
  * Mannor and Shamir, (2011) Mannor, S. and Shamir, O. (2011).  From bandits to experts: On the value of side-observations.  In Advances in Neural Information Processing Systems, pages 684–692.
  * Perrault et al., (2020) Perrault, P., Boursier, E., Valko, M., and Perchet, V. (2020).  Statistical efficiency of thompson sampling for combinatorial semi-bandits.  In Advances in Neural Information Processing Systems.
  * Petersen and Pedersen, (2012) Petersen, K. B. and Pedersen, M. S. (2012).  The matrix cookbook, nov 2012.  URL http://www2. imm. dtu. dk/pubdb/p. php, 3274:14.
  * Pukelsheim, (2006) Pukelsheim, F. (2006).  Optimal Design of Experiments.  Society for Industrial and Applied Mathematics.
  * Rizk et al., (2019) Rizk, G., Colin, I., Thomas, A., and Draief, M. (2019).  Refined bounds for randomized experimental design.  NeurIPS Workshop on Machine Learning with Guarantees.
  * Robbins, (1952) Robbins, H. (1952).  Some aspects of the sequential design of experiments.  Bulletin of the American Mathematical Society, 58(5):527–535.
  * Sagnol, (2010) Sagnol, G. (2010).  Optimal design of experiments with application to the inference of traffic matrices in large networks: second order cone programming and submodularity.  PhD thesis, École Nationale Supérieure des Mines de Paris.
  * Siomina et al., (2006) Siomina, I., Varbrand, P., and Yuan, D. (2006).  Automated optimization of service coverage and base station antenna configuration in UMTS networks.  IEEE Wireless Communications, 13(6):16–25.
  * Soare et al., (2014) Soare, M., Lazaric, A., and Munos, R. (2014).  Best-arm identification in linear bandits.  In Advances in Neural Information Processing Systems, pages 828–836.
  * Tao et al., (2018) Tao, C., Blanco, S., and Zhou, Y. (2018).  Best arm identification in linear bandits with linear dimension dependency.  In Proceedings of the 35th International Conference on Machine Learning, volume 80, pages 4877–4886.
  * Tropp et al., (2015) Tropp, J. A. et al. (2015).  An introduction to matrix concentration inequalities.  Foundations and Trends® in Machine Learning, 8(1-2):1–230.
  * Valko, (2020) Valko, M. (2020).  Bandits on graphs and structures.
  * Valko et al., (2014) Valko, M., Munos, R., Kveton, B., and Kocák, T. (2014).  Spectral bandits for smooth graph functions.  In Xing, E. P. and Jebara, T., editors, International conference on machine learning, volume 32 of Proceedings of Machine Learning Research, pages 46–54.
  * Welch, (1982) Welch, W. (1982).  Algorithmic complexity: Three np-hard problems in computational statistics.  Journal of Statistical Computation and Simulation - J STAT COMPUT SIM, 15:17–25.
  * Xu et al., (2018) Xu, L., Honda, J., and Sugiyama, M. (2018).  A fully adaptive algorithm for pure exploration in linear bandits.  In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84, pages 843–851.
  * Zaki et al., (2019) Zaki, M., Mohan, A., and Gopalan, A. (2019).  Towards optimal and efficient best arm identification in linear bandits.  arXiv preprint arXiv:1911.01695.
  * Zaki et al., (2020) Zaki, M., Mohan, A., and Gopalan, A. (2020).  Explicit best arm identification in linear bandits using no-regret learners.  arXiv preprint arXiv:2006.07562.
  * Çivril and Magdon-Ismail, (2009) Çivril, A. and Magdon-Ismail, M. (2009).  On selecting a maximum volume sub-matrix of a matrix and related problems.  Theoretical Computer Science, 410(47):4801 – 4811.

##  Appendix A An NP-Hard Problem

###  A.1 Proof of Theorem 4.1

######  Theorem A.1.

Consider a given matrix $\mathbf{M}_{\star}\in\mathbb{R}^{d\times d}$ and a
finite arm set $\mathcal{X}\subset\mathbb{R}^{d}$. Unless P=NP, there is no
polynomial time algorithm guaranteed to find the optimal solution of

$\displaystyle\max_{\left(x^{(1)},\ldots,x^{(n)}\right)\in\mathcal{X}^{n}}\sum_{(i,j)\in
E}x^{(i)\top}\mathbf{M}_{\star}\ x^{(j)}\enspace.$

###### Proof.

We prove the statement by reduction to the Max-Cut problem. Let
$\mathcal{G}=(V,E)$ be a graph with $V=\\{1,\ldots,n\\}$. Let
$\mathcal{X}=\left\\{e_{0},e_{1}\right\\}$, where
$e_{0}=\left(1,0\right)^{\top}$ and $e_{1}=\left(0,1\right)^{\top}$. Let
$\mathbf{M}_{\star}=\left[\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right]$.
For any joint arm assignment $\left(x^{(1)}\ldots
x^{(n)}\right)\in\mathcal{X}^{n}$, let $F\subseteq E$ be defined as
$F=\left\\{i:x^{(i)}=e_{1}\right\\}$. Note that

$\sum_{(i,j)\in E}x^{(i)\top}\mathbf{M}_{\star}x^{(j)}=\sum_{(i,j)\in
E}\mathbf{1}\left[x^{(i)}\neq x^{(j)}\right]=2\times\sum_{(i,j)\in
E}\mathbf{1}\left[i\in F,j\notin F\right],$

where $\mathbf{1}[\cdot]$ is the indicator function. The assignement
$\left(x^{(1)},\ldots,x^{(n)}\right)$ induces a cut $(F,V\backslash F)$, and
the value of the assignment is _precisely_ twice the value of the cut. Thus,
if there was a polynomial time algorithm solving our problem, this algorithm
would also solve the Max-Cut problem. ∎

###  A.2 Proof of Theorem 4.2

######  Theorem A.2.

Let us consider the graph $\mathcal{G}=(V,E)$, a finite arm set
$\mathcal{X}\subset\mathbb{R}^{d}$ and the matrix $\mathbf{M}_{\star}$ given
as input to Algorithm 1. Then, the expected global reward $r=\sum_{(i,j)\in
E}x^{(i)\top}\mathbf{M}_{\star}x^{(j)}$ associated to the returned allocation
$\mathbf{x}=\left(x^{(1)},\dots,x^{(n)}\right)\in\mathcal{X}^{n}$ verifies:

$\displaystyle\frac{r-r_{\mathrm{min}}}{r_{\star}-r_{\mathrm{min}}}\geq\frac{1}{2}\enspace.$

where $r_{\star}$ and $r_{\mathrm{min}}$ are respectively the highest and
lowest global reward one can obtain with the appropriate joint arm. Finally,
the complexity of the algorithm is in $\mathcal{O}(K^{2}+n)$.

###### Proof.

Given the matrix $\mathbf{M}_{\star}$, the algorithm obtains the two node-arms
$(x_{\star},x_{\star}^{\prime})\in\mathcal{X}$ solution of

$\displaystyle\max_{(x,x^{\prime})\in\mathcal{X}}x^{\top}\mathbf{M}_{\star}x^{\prime}\enspace.$

Note that it is equivalent to obtain $z_{\star}$ solution of

$\displaystyle\max_{z\in\mathcal{Z}}\;z^{\top}\operatorname*{vec}{\left(\mathbf{M}_{\star}\right)}\enspace.$

Let us analyze a round of Algorithm 1 where we assign the arm of a node in
$V$. For sake of simplicity, we assume that node $i$ is assigned at round $i$.
At round $i$, we count the number $n_{1}^{(i)}$ of neighbors of $i$ that have
already been assigned the arm $x_{\star}$ and we count the number
$n_{2}^{(i)}$ of neighbors of $i$ that have already been assigned the arm
$x_{\star}^{\prime}$. Then, node $i$ is assigned the arm least represented
among its neighbors, that is arm $x_{\star}$ if $n_{2}^{(i)}\geq n_{1}^{(i)}$
and $x^{\prime}_{\star}$ otherwise. Eventually, the optimal edge-arm
$z_{\star}$ has been assigned $\max(n_{1}^{(i)},n_{2}^{(i)})$ times among node
$i$’s neighborhood. Hence, for each node $i$, if we denote $r_{i}$ the sum of
all the rewards obtained with the edge-arms constructed only during the round
$i$, we have

$\displaystyle r_{i}$
$\displaystyle=\max\left(n_{1}^{(i)},n_{2}^{(i)}\right)z_{\star}^{\top}\theta_{\star}+\min\left(n_{1}^{(i)},n_{2}^{(i)}\right)z^{\top}\theta_{\star}$
$\displaystyle\geq\frac{n_{1}^{(i)}+n_{2}^{(i)}}{2}(z_{\star}^{\top}\theta_{\star}+z^{\top}\theta_{\star})\enspace.$

One can notice that the arm $z$ can only be equal to
$\operatorname*{vec}{\left(x_{\star}x_{\star}^{\top}\right)}$ or
$\operatorname*{vec}{\left(x_{\star}^{\prime\top}x_{\star}^{\prime\top}\right)}$.
Let assume that
$\operatorname*{vec}{\left(x_{\star}x_{\star}^{\top}\right)}^{\top}\theta_{\star}\leq\operatorname*{vec}{\left(x_{\star}^{\prime}x_{\star}^{\prime\top}\right)}^{\top}\theta_{\star}$
without loss of generality and let consider the worst case where $z$ is always
equal to $\operatorname*{vec}{\left(x_{\star}x_{\star}^{\top}\right)}$. Since
z is constructed with the same node-arm $x_{\star}$, the allocation that
constructs at each edge the edge-arm $z$ exists (which is allocating
$x_{\star}$ to all the nodes), thus $m\times z^{\top}\theta_{\star}\geq
r_{\mathrm{min}}$.

Moreover one can also notice that $m\times z_{\star}^{\top}\theta_{\star}\geq
r_{\star}$. We thus have,

$\displaystyle r_{i}$
$\displaystyle\geq\frac{n_{1}^{(i)}+n_{2}^{(i)}}{2m}(r_{\star}+r_{\mathrm{min}})\enspace.$

Now let us sum all the rewards obtained with the constructed edge-arms at each
round of the algorithm, that is the global reward $r$ of the graph allocation
returned by the proposed algorithm:

$\displaystyle r$ $\displaystyle=\sum_{i=1}^{n}r_{i}$
$\displaystyle\geq\sum_{i=1}^{n}\frac{n_{1}^{(i)}+n_{2}^{(i)}}{2m}(r_{\star}+r_{\mathrm{min}})$
$\displaystyle=\frac{1}{2}(r_{\star}+r_{\mathrm{min}})$
$\displaystyle=\frac{1}{2}(r_{\star}-r_{\mathrm{min}})+r_{\mathrm{min}}\enspace.$

Moreover, the algorithm does $K^{2}$ estimation to find the best couple
$(x_{\star},x_{\star}^{\prime})\in\mathcal{X}^{2}$, and each of the $n$ rounds
of the algorithm is in $O(1)$. Hence the complexity is equal to $O(K^{2}+n)$.
∎

##  Appendix B Deriving the stopping condition

In this section, we remind key results to derive the stopping condition. We
refer the reader to Soare et al., (2014) and references therein for additional
details. Let $\mathcal{Z}\subset\mathbb{R}^{d^{2}}$ be the set of edge-arms
and let $K^{2}=|\mathcal{Z}|$. For $m,t>0$, we consider a sequence of edge-
arms $\mathbf{z}_{t}=(z_{1},\dots,z_{mt})\in\mathcal{Z}^{mt}$ and the
corresponding noisy rewards $(r_{1},\dots,r_{mt})$. We assume that the noise
terms in the rewards are i.i.d., following a $\sigma$-sub-Gaussian
distribution. Let
$\hat{\theta}_{t}=\mathbf{A}_{t}^{-1}b_{t}\in\mathbb{R}^{d^{2}}$ be the
solution of the ordinary least squares problem with
$\mathbf{A}_{t}=\sum_{i=1}^{mt}z_{i}z_{i}^{\top}\in\mathbb{R}^{d^{2}\times
d^{2}}$ and $b_{t}=\sum_{i=1}^{k}z_{i}r_{i}\in\mathbb{R}^{d^{2}}$. We first
recall the following property.

######  Proposition B.1 (Proposition 1 in Soare et al., (2014)).

Let $c=2\sigma\sqrt{2}$. For every fixed sequence $\mathbf{z}_{t}$, with
probability $1-\delta$, for all $t>0$ and for all $z\in\mathcal{Z}$, we have

$\displaystyle\left|z^{\top}\theta_{\star}-z^{\top}\hat{\theta}_{t}\right|\leq
c\|z\|_{\mathbf{A}_{t}^{-1}}\sqrt{\log\left(\frac{6m^{2}t^{2}K^{2}}{\delta\pi}\right)}\enspace.$

Our goal is to find the arm $z_{\star}$ that has the optimal expected reward
$z_{\star}^{\top}\theta_{\star}$. In other words, we want to find an arm
$z\in\mathcal{Z}$, such that for all $z^{\prime}\in\mathcal{Z}$,
$(z-z^{\prime})^{\top}\theta_{\star}\geq 0$. However, one does not have access
to $\theta_{\star}$, so we have to use its empirical estimate.

Let us consider a confidence set $\hat{S}(\mathbf{z}_{t})$ centered at
$\hat{\theta}_{t}\in\hat{S}(\mathbf{z}_{t})$ and such that
$\mathbb{P}\left(\theta_{\star}\notin\hat{S}(\mathbf{z}_{t})\right)\leq\delta$,
for some $\delta>0$. Since $\theta_{\star}$ belongs to
$\hat{S}(\mathbf{z}_{t})$ with probability at least $1-\delta$, one can stop
pulling arms when an arm has been found, such that the above condition is
verified for any $\theta\in\hat{S}(\mathbf{z}_{t})$. More formally, the best
arm identification task will be considered successful when an arm
$z\in\mathcal{Z}$ will verify the following condition for any
$z^{\prime}\in\mathcal{Z}$ and any $\theta\in\hat{S}(\mathbf{z}_{t})$:

$(z-z^{\prime})^{\top}(\hat{\theta}_{t}-\theta)\leq\hat{\Delta}_{t}(z,z^{\prime})\enspace,$

where
$\hat{\Delta}_{t}\left(z,z^{\prime}\right)=\left(z-z^{\prime}\right)^{\top}\hat{\theta}_{t}$
is the empirical gap between $z$ and $z^{\prime}$.

Using the upper bound in Proposition B.1, one way to ensure that
$\mathbb{P}\left(\theta_{\star}\in\hat{S}(\mathbf{z}_{t})\right)\geq 1-\delta$
is to define the confidence set $\hat{S}(\mathbf{z}_{t})$ as follows

$\displaystyle\hat{S}\left(\mathbf{z}_{t}\right)=\left\\{\theta\in\mathbb{R}^{d},\
\forall z\in\mathcal{Z},\ \forall
z^{\prime}\in\mathcal{Z},\left(z-z^{\prime}\right)^{\top}\left(\hat{\theta}_{t}-\theta\right)\leq
c\|z-z^{\prime}\|_{\left(\mathbf{A}_{t}\right)^{-1}}\sqrt{\log\left(\frac{6m^{2}t^{2}K^{4}}{\delta\pi}\right)}\right\\}\enspace.$

Then, the stopping condition can be reformulated as follows:

$\displaystyle\exists z\in\mathcal{Z},\ \forall z^{\prime}\in\mathcal{Z},\
c\|z-z^{\prime}\|_{\mathbf{A}_{t}^{-1}}\sqrt{\log\left(\frac{6m^{2}t^{2}K^{4}}{\delta\pi}\right)}\leq\hat{\Delta}_{t}\left(z,z^{\prime}\right)\enspace.$
(4)

##  Appendix C Estimation of the unknown parameter

###  C.1 Proof of Theorem 5.1

To prove Theorem 5.1, we first state some useful propositions and lemmas. For
any finite set $X\subset\mathbb{R}^{d}$, we define the function
$h_{X}:\mathcal{S}_{X}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ as follows: for
any $\lambda\in\mathcal{S}_{X}$,

$h_{X}(\lambda)=\begin{cases}\max_{x^{\prime}\in
X}x^{\prime\top}\Sigma_{X}(\lambda)^{-1}x^{\prime}&\text{if
}\Sigma_{X}(\lambda)\text{ is invertible}\\\
+\infty&\text{otherwise}\enspace.\end{cases}$

######  Lemma C.1.

Let $\mathcal{X}\subset\mathbb{R}^{d}$ be a finite set spanning
$\mathbb{R}^{d}$ and let
$\mathcal{Z}=\\{\operatorname*{vec}{(xx^{\prime\top})}\text{,
}(x,x^{\prime})\in\mathcal{X}^{2}\\}$. If
$\mu^{\star}\in\mathcal{S}_{\mathcal{X}}$ is a minimizer of $h_{\mathcal{X}}$,
then $\mu^{\star}$ is a solution of

$\displaystyle\min_{\mu\in\mathcal{S}_{\mathcal{X}}}\max_{z\in\mathcal{Z}}z^{\top}\left(\sum_{x\in\mathcal{X}}\sum_{x^{\prime}\in\mathcal{X}}\mu_{x}\mu_{x^{\prime}}\operatorname*{vec}{\left(xx^{\prime\top}\right)}\operatorname*{vec}{\left(xx^{\prime\top}\right)}^{\top}\right)^{-1}z\enspace.$

###### Proof.

First, let us notice that, for any $\mathcal{X}\subset\mathbb{R}^{d}$, one has
$h_{\mathcal{X}}\geq 0$. Thus, $\mu^{\star}$ is also a minimizer of
$h_{\mathcal{X}}^{2}$. In addition, $\mathcal{X}$ is spanning $\mathbb{R}^{d}$
so $h_{\mathcal{X}}(\mu^{\star})<+\infty$. Developing
$h_{\mathcal{X}}(\mu^{\star})^{2}$ yields:

$\displaystyle h_{\mathcal{X}}(\mu^{\star})\times
h_{\mathcal{X}}(\mu^{\star})$
$\displaystyle=\left(\max_{x\in\mathcal{X}}\;x^{\top}\Sigma_{X}(\mu^{\star})^{-1}x\right)\times\left(\max_{x\in\mathcal{X}}\;x^{\top}\Sigma_{X}(\mu^{\star})^{-1}x\right)$
$\displaystyle=\max_{x\in\mathcal{X}}\max_{x^{\prime}\in\mathcal{X}}\;x^{\top}\Sigma_{X}(\mu^{\star})^{-1}xx^{\prime\top}\Sigma_{X}(\mu^{\star})^{-1}x^{\prime}$
$\displaystyle=\max_{x\in\mathcal{X}}\max_{x^{\prime}\in\mathcal{X}}\;\operatorname*{vec}{\left(xx^{\prime\top}\right)}^{\top}\operatorname*{vec}{\left(\Sigma_{X}(\mu^{\star})^{-1}xx^{\prime\top}\Sigma_{X}(\mu^{\star})^{-1}\right)}$
$\displaystyle=\max_{x\in\mathcal{X}}\max_{x^{\prime}\in\mathcal{X}}\;\operatorname*{vec}{\left(xx^{\prime\top}\right)}^{\top}\left(\Sigma_{X}(\mu^{\star})^{-1}\otimes\Sigma_{X}(\mu^{\star})^{-1}\right)\operatorname*{vec}{\left(xx^{\prime\top}\right)}$
$\displaystyle=\max_{z\in\mathcal{Z}}\;z^{\top}\left(\Sigma_{X}(\mu^{\star})^{-1}\otimes\Sigma_{X}(\mu^{\star})^{-1}\right)z\enspace,$

where $\otimes$ denotes the Kronecker product. We can now focus on the central
term:

$\displaystyle\Sigma_{X}(\mu^{\star})^{-1}\otimes\Sigma_{X}(\mu^{\star})^{-1}$
$\displaystyle=\left(\sum_{x\in\mathcal{X}}\mu^{\star}_{x}xx^{\top}\right)^{-1}\otimes\left(\sum_{x\in\mathcal{X}}\mu^{\star}_{x}xx^{\top}\right)^{-1}$
$\displaystyle=\left(\sum_{x\in\mathcal{X}}\mu^{\star}_{x}xx^{\top}\otimes\sum_{x\in\mathcal{X}}\mu^{\star}_{x}xx^{\top}\right)^{-1}$
$\displaystyle=\left(\sum_{x\in\mathcal{X}}\sum_{x^{\prime}\in\mathcal{X}}\mu^{\star}_{x}\mu^{\star}_{x^{\prime}}\left(xx^{\top}\otimes
x^{\prime}x^{\prime\top}\right)\right)^{-1}$
$\displaystyle=\left(\sum_{x\in\mathcal{X}}\sum_{x^{\prime}\in\mathcal{X}}\mu^{\star}_{x}\mu^{\star}_{x^{\prime}}\operatorname*{vec}{\left(xx^{\prime\top}\right)}\operatorname*{vec}{\left(xx^{\prime\top}\right)}^{\top}\right)^{-1}\enspace,$

and the result holds. ∎

######  Theorem C.2.

Let $\mu^{\star}\in\mathcal{S}_{\mathcal{X}}$ be a minimizer of
$h_{\mathcal{X}}$. Let $\lambda^{\star}\in\mathcal{S}_{\mathcal{Z}}$ be the
distribution defined from $\mu^{\star}$ such that, for all
$z=\operatorname*{vec}{(xx^{\prime\top})}$,
$\lambda^{\star}_{z}=\mu^{\star}_{x}\mu^{\star}_{x^{\prime}}$. Then
$\lambda^{\star}$ is a minimizer of $h_{\mathcal{Z}}$.

###### Proof.

From Kiefer and Wolfowitz, (1960), we know that
$\min_{\lambda\in\mathcal{S}_{\mathcal{Z}}}h_{\mathcal{Z}}(\lambda)=d^{2}$ and
$\min_{\mu\in\mathcal{S}_{\mathcal{X}}}h_{\mathcal{X}}(\mu)=d$. Then, using
Proposition C.1, one has

$\displaystyle d^{2}$ $\displaystyle=h_{\mathcal{X}}(\mu^{\star})\times
h_{\mathcal{X}}(\mu^{\star})$
$\displaystyle=\max_{z\in\mathcal{Z}}\;z^{\top}\left(\sum_{x\in\mathcal{X}}\sum_{x^{\prime}\in\mathcal{X}}\mu^{\star}_{x}\mu^{\star}_{x^{\prime}}\operatorname*{vec}{\left(xx^{\prime\top}\right)}\operatorname*{vec}{\left(xx^{\prime\top}\right)}^{\top}\right)^{-1}z\enspace.$

This result implies that $h_{\mathcal{Z}}(\lambda^{\star})=d^{2}$. Since
$\min_{\lambda\in\mathcal{S}_{\mathcal{Z}}}h_{Z}(\lambda)=d^{2}$,
$\lambda^{\star}$ is a minimizer of $h_{\mathcal{Z}}$.

∎

###  C.2 Proof of Theorem 5.2

To prove our confidence bound, we need the two following proposition. The
first one is from (Tropp et al.,, 2015).

######  Proposition C.3 (Tropp et al., (2015), Chapter 5 and 6).

Let $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{t}$ be i.i.d. positive semi-definite
random matrices in $\mathbb{R}^{d^{2}\times d^{2}}$, such that there exists
$L>0$ verifying $\mathbf{0}\preceq\mathbf{Z}_{1}\preceq mL\mathbf{I}$. Let
$\mathbf{A}_{t}$ be defined as
$\mathbf{A}_{t}\triangleq\sum_{s=1}^{t}\mathbf{Z}_{s}$. Then, for any
$0<\varepsilon<1$, one can lowerbound $\lambda_{\mathrm{min}}(\mathbf{A}_{t})$
as follows:

$\mathbb{P}(\lambda_{\mathrm{min}}(\mathbf{A}_{t})\leq(1-\varepsilon)\lambda_{\mathrm{min}}(\mathbb{E}\mathbf{A}_{t}))\leq
d^{2}e^{-\displaystyle\frac{t\varepsilon^{2}\lambda_{\mathrm{min}}(\mathbb{E}\mathbf{Z}_{1})}{2mL}}.$

If in addition, there exists some $v>0$, such that
$\|\mathbb{E}\big{[}(\mathbf{Z}_{1}-\mathbb{E}\mathbf{Z}_{1})^{2}\big{]}\|\leq
v$, then for any $u>0$, one has

$\mathbb{P}\left(\left\|\mathbf{S}_{t}\right\|\geq u\right)\leq
2d^{2}e^{\displaystyle-\frac{u^{2}}{2mLu/3+2tv}},$

From the second inequality, (Rizk et al.,, 2019) derived a slightly different
inequality that we use here :

######  Proposition C.4 (Rizk et al., (2019), Appendix A.3).

Let $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{t}$ be $t$ i.i.d. random symmetric
matrices in $\mathbb{R}^{d^{2}\times d^{2}}$ such that there exists $L>0$ such
that $\|\mathbf{Z}_{1}\|\leq mL$, almost surely. Let
$\mathbf{A}_{t}\triangleq\sum_{i=1}^{t}\mathbf{Z}_{i}$. Then, for any $u>0$,
one has:

$\mathbb{P}\left(\left\|\mathbf{A}_{t}-\mathbb{E}\mathbf{A}_{t}\right\|\geq\sqrt{2tvu}+\frac{mLu}{3}\right)\leq
d^{2}e^{-u}\enspace.$

where
$v\triangleq\left\|\mathbb{E}\big{[}(\mathbf{Z}_{1}-\mathbb{E}\mathbf{Z}_{1})^{2}\big{]}\right\|$.

Finally, to prove our main theorem, we need the following lemma.

######  Lemma C.5.

One has
$\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|\leq\frac{d^{2}}{\nu_{\mathrm{min}}}$,
where $\nu_{\mathrm{min}}$ is the smallest eigenvalue of the covariance matrix
$\frac{1}{K^{2}}\sum_{z\in\mathcal{Z}}z^{\top}z$.

###### Proof.

Define $\mathcal{B}=\left\\{z\in\mathbb{R}^{d^{2}}:\|z\|=1\right\\}$. First,
for any semi-definite matrix $\mathbf{A}\in\mathbb{R}^{d^{2}\times d^{2}}$, we
have $\left\|\mathbf{A}\right\|=\max_{z\in\mathcal{B}}\;z^{\top}\mathbf{A}z$.
Because $\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}$ is positive definite and
symmetric, and by Rayleigh-Ritz theorem,

$\displaystyle\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|$
$\displaystyle=\max_{z\in\mathcal{B}}\;\frac{z^{\top}\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}z}{z^{\top}z}=\max_{z\in\mathcal{B}}\;z^{\top}\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}z\enspace.$

Let $\mathbf{Z}\in\mathbb{R}^{K^{2}\times d^{2}}$ be the matrix whose rows are
vectors of $\mathcal{Z}$ in an arbitrary order. Notice that $\mathcal{Z}$
spans $\mathbb{R}^{d^{2}}$, since $\mathcal{X}$ spans $\mathbb{R}^{d}$. Now
for any $z\in\mathcal{B}$, define $\beta^{(z)}\in\mathbb{R}^{K^{2}}$ as a
vector such that $z=\mathbf{Z}^{\top}\beta^{(z)}$ . Then,

$\displaystyle\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|$
$\displaystyle=\max_{z\in\mathcal{B}}\;\beta^{(z)^{\top}}\mathbf{Z}\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\mathbf{Z}^{\top}\beta^{(z)}$
$\displaystyle=\max_{z\in\mathcal{B}}\;\sum_{i=1}^{d^{2}}\sum_{j=1}^{d^{2}}\beta_{i}^{(z)}\beta_{j}^{(z)}z_{i}^{\top}\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}z_{j}$
$\displaystyle\leq\max_{z\in\mathcal{B}}\;\left\|\beta^{(z)}\right\|_{1}^{2}\times\max_{i,j}\;z_{i}^{\top}\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}z_{j}\enspace.$

Define
$\tilde{z}_{i}=\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-\frac{1}{2}}z_{i}$.
Clearly,
$\max_{i,j}\;z_{i}^{\top}\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}z_{j}=\max_{i,j}\;\tilde{z}_{i}^{\top}\tilde{z}_{j}=\max_{i}\;\tilde{z}_{i}^{2}$.
So we have

$\displaystyle\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|$
$\displaystyle\leq\max_{z\in\mathcal{B}}\;\left\|\beta^{(z)}\right\|_{1}^{2}\times\max_{z^{\prime}\in\mathcal{Z}}\;z^{\prime\top}\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}z^{\prime}$
$\displaystyle\leq\max_{z\in\mathcal{B}}\;\left\|\beta^{(z)}\right\|_{1}^{2}d^{2}\enspace.$

The last inequality comes from Kiefer and Wolfowitz equivalence theorem
(Kiefer and Wolfowitz,, 1960). Now observe that $\beta^{(z)}$ can be obtained
by least square regression :
$\beta^{(z)}=\left(\mathbf{Z}\mathbf{Z}^{\top}\right)^{-1}\mathbf{Z}z=\left(\mathbf{Z}^{\top}\right)^{\dagger}z$
where $\left(\cdot\right)^{\dagger}$ is the Moore-Penrose pseudo-inverse. Note
that $\mathbf{Z}\mathbf{Z}^{\top}$ is a Gram matrix. It is known that for a
matrix having singular values $\left\\{\sigma_{i}\right\\}_{i}$, its pseudo-
inverse has singular values $\begin{cases}\frac{1}{\sigma_{i}}&\text{if
}\sigma_{i}\neq 0\\\ 0&\text{otherwise}\end{cases}$ for all $i$. So for
$z\in\mathcal{B}$, we have:

$\displaystyle\left\|\beta^{(z)}\right\|_{1}^{2}$ $\displaystyle\leq
K^{2}\left\|\beta^{(z)}\right\|_{2}^{2}\leq
K^{2}\left\|\left(\mathbf{Z}^{\top}\right)^{\dagger}\right\|^{2}\leq\frac{K^{2}}{\sigma_{\mathrm{min}}\left(\mathbf{Z}\right)^{2}}\enspace,$

where $\sigma_{\mathrm{min}}\left(\cdot\right)$ refers to the smallest
singular value. Let $\nu_{\mathrm{min}}\left(\cdot\right)$ refer to the
smallest eigenvalue. Noting that

$\sigma_{\mathrm{min}}\left(\mathbf{Z}\right)^{2}=\nu_{\mathrm{min}}\left(\mathbf{Z}^{\top}\mathbf{Z}\right)=K^{2}\nu_{\mathrm{min}}\left(\frac{1}{K^{2}}\sum_{z\in\mathcal{Z}}zz^{\top}\right)\enspace,$

yields the desired result.

∎

We are now ready to state the bound on the random sampling error, relatively
to the objective value $\Sigma_{\mathcal{Z}}(\lambda^{\star})$ of the convex
relaxation solution.

######  Theorem C.6.

Let $\lambda^{\star}\in\mathcal{S}_{\mathcal{Z}}$ be a minimizer of
$h_{\mathcal{Z}}$. Let $0\leq\delta\leq 1$ and let $t_{0}>0$ be such that

$t_{0}=2Ld^{2}\log(2d^{2}/\delta)/\nu_{\mathrm{min}}\enspace,$

where $L=\max_{z\in\mathcal{Z}}\|z\|^{2}$ and $\nu_{\mathrm{min}}$ is the
smallest eigenvalue of the covariance matrix
$\frac{1}{K^{2}}\sum_{z\in\mathcal{Z}}z^{\top}z$. Then, at each round $t\geq
t_{0}$, with probability at least $1-\delta$, the randomized G-allocation
strategy for graphical bilinear bandit in Algorithm 2 produces a matrix
$\mathbf{A}_{t}$ such that:

$\displaystyle
h_{\mathcal{Z}}(\mathbf{A}_{t})\leq(1+\alpha)h_{\mathcal{Z}}(mt\times\Sigma_{\mathcal{Z}}(\lambda^{\star}))\enspace,$

where

$\displaystyle\alpha=\frac{Ld^{2}}{m\nu_{\mathrm{min}}^{2}}\sqrt{\frac{2v}{t}\log\left(\frac{2d^{2}}{\delta}\right)}+o\left(\frac{1}{\sqrt{t}}\right),$

and
$v\triangleq\|\mathbb{E}\big{[}(\mathbf{A}_{1}-\mathbb{E}\mathbf{A}_{1})^{2}\big{]}\|$.

###### Proof.

Let $(X^{(1)}_{s})_{s=1,\dots,t},\ldots,(X^{(n)}_{s})_{s=1,\ldots,t}$ be $nt$
i.i.d. random vectors in $\mathbb{R}^{d}$ such that for all $x\in\mathcal{X}$,
$\mathbb{P}\left(X^{(1)}_{1}=x\right)=\mu^{\star}_{x}$. For $(i,j)\in E$ and
$1\leq s\leq t$, we define the random matrix $\mathbf{Z}^{(i,j)}_{s}$ by

$\displaystyle\mathbf{Z}^{(i,j)}_{s}=\operatorname*{vec}{\left(X^{i}_{s}X^{j\top}_{s}\right)}\operatorname*{vec}{\left(X^{i}_{s}X^{j\top}_{s}\right)}^{\top}\enspace.$

Finally, let us define for all $1\leq s\leq t$, the edge-wise sum
$\mathbf{Z}_{s}\in\mathbb{R}^{d^{2}\times d^{2}}$, that is

$\displaystyle\mathbf{Z}_{s}=\sum_{(i,j)\in E}\mathbf{Z}^{(i,j)}_{s}\enspace.$

One can easily notice that $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{t}$ are i.i.d.
random matrices. We define the overall sum
$\mathbf{A}_{t}=\sum_{s=1}^{t}\mathbf{Z}_{s}$ and our goal is to measure how
close $f_{\mathcal{Z}}(\mathbf{A}_{t})$ is to
$f_{\mathcal{Z}}(mt\times\Sigma_{\mathcal{Z}}(\lambda^{\star}))$, where $mt$
corresponds to the total number of sampled arms $z\in\mathcal{Z}$ during the
$t$ rounds of the learning procedure. By definition of $\mathbf{A}_{t}$, one
has

$\displaystyle\max_{z\in\mathcal{Z}}\;z^{\top}\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}z$
$\displaystyle=\max_{z\in\mathcal{Z}}\;z^{\top}\left(\sum_{s=1}^{t}\sum_{(i,j)\in
E}\mathbb{E}\left[\mathbf{Z}_{s}^{(i,j)}\right]\right)^{-1}z$
$\displaystyle=\max_{z\in\mathcal{Z}}\;z^{\top}\left(\sum_{s=1}^{t}\sum_{(i,j)\in
E}\sum_{x,x^{\prime}\in\mathcal{X}}\mu^{\star}_{x}\mu^{\star}_{x^{\prime}}\operatorname*{vec}{(xx^{\prime\top})}\operatorname*{vec}{(xx^{\prime\top})}^{\top}\right)^{-1}z$
$\displaystyle=\max_{z\in\mathcal{Z}}\;z^{\top}\left(\sum_{s=1}^{t}\sum_{(i,j)\in
E}\sum_{z^{\prime}\in\mathcal{Z}}\lambda^{\star}_{z^{\prime}}z^{\prime}z^{\prime\top}\right)^{-1}z$
$\displaystyle=f_{\mathcal{Z}}(mt\Sigma_{\mathcal{Z}}(\lambda^{\star}))\enspace.$

This allows us to bound the relative error as follows:

$\displaystyle\alpha$
$\displaystyle=\frac{f_{\mathcal{Z}}(\mathbf{A}_{t})}{f_{\mathcal{Z}}(mt\times\Sigma_{\mathcal{Z}}(\lambda^{\star}))}-1$
$\displaystyle=\frac{\max_{z\in\mathcal{Z}}\;z^{\top}\left(\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}+\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right)z}{f_{\mathcal{Z}}(mt\times\Sigma_{\mathcal{Z}}(\lambda^{\star}))}-1$
$\displaystyle\leq\frac{\max_{z\in\mathcal{Z}}\;z^{\top}\left(\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right)z}{f_{\mathcal{Z}}(mt\times\Sigma_{\mathcal{Z}}(\lambda^{\star}))}\enspace.$

Using the fact that
$f_{\mathcal{Z}}(mt\Sigma_{\mathcal{Z}}(\lambda^{\star}))=d^{2}/mt$ (Kiefer
and Wolfowitz,, 1960), we obtain

$\displaystyle\alpha$
$\displaystyle\leq\frac{mt}{d^{2}}\times\max_{z\in\mathcal{Z}}\;z^{\top}\left(\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right)z$
$\displaystyle\leq\frac{mt}{d^{2}}\times\max_{z\in\mathcal{Z}}\;\|z\|^{2}\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|$
$\displaystyle\leq\frac{mtL}{d^{2}}\times\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|\enspace.$

Therefore, controlling the quantity
$\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|$ will
allow us to provide an upper bound on the relative error. Notice that

$\displaystyle\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|$
$\displaystyle=\|\mathbf{A}_{t}^{-1}\left(\mathbb{E}\mathbf{A}_{t}-\mathbf{A}_{t}\right)\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|$
$\displaystyle\leq\|\mathbf{A}_{t}^{-1}\|\
\|\mathbb{E}\mathbf{A}_{t}-\mathbf{A}_{t}\|\
\|\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|\enspace.$

Using Proposition C.3, we know that for any
$d^{2}e^{-\frac{t\lambda_{\mathrm{min}}(\mathbb{E}\mathbf{Z}_{1})}{mL}}<\delta_{h}<1$,
the following holds:

$\displaystyle\|\mathbf{A}_{t}^{-1}\|\leq\frac{\|\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|}{1-\sqrt{\frac{2mL}{t}\|\left(\mathbb{E}\mathbf{Z}_{1}\right)^{-1}\|\log(d^{2}/\delta_{h})}}\enspace,$

with probability at least $1-\delta_{h}$. Similarly, using Proposition C.4,
for any $0<\delta_{b}<1$, we have

$\|\mathbf{A}_{t}-\mathbb{E}\mathbf{A}_{t}\|\leq\frac{mL}{3}\log\frac{d^{2}}{\delta_{b}}+\sqrt{2tv^{2}\log\frac{d^{2}}{\delta_{b}}}\enspace,$

with probability at least $1-\delta_{b}$. Combining these two results with a
union bound leads to the following bound, with probability
$1-(\delta_{b}+\delta_{h})$:

$\displaystyle\left\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right\|\leq\left\|\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right\|^{2}\frac{(mL/3)\log(d^{2}/\delta_{b})+\sqrt{2tv\log(d^{2}/\delta_{b})}}{1-\sqrt{(2mL/t)\left\|\left(\mathbb{E}\mathbf{Z}_{1}\right)^{-1}\right\|\log(d^{2}/\delta_{h})}}\enspace.$

In order to obtain a unified bound depending on one confidence parameter
$1-\delta$, one could optimize over $\delta_{b}$ and $\delta_{h}$, subject to
$\delta_{b}+\delta_{h}=\delta$. This leads to a messy result and a negligible
improvement. One can use simple values $\delta_{b}=\delta_{h}=\delta/2$, so
the overall bound becomes, with probability $1-\delta$:

$\displaystyle\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\|\leq\frac{1}{tm^{2}}\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|^{2}\sqrt{\frac{2v}{t}\log\left(\frac{2d^{2}}{\delta}\right)}\left(\frac{1+\sqrt{\frac{m^{2}L^{2}\log(2d^{2}/\delta)}{18vt}}}{1-\sqrt{\frac{2L\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\|\log(2d^{2}/\delta)}{t}}}\right)\enspace.$

This can finally be formulated as follows:

$\displaystyle\left\|\mathbf{A}_{t}^{-1}-\left(\mathbb{E}\mathbf{A}_{t}\right)^{-1}\right\|\leq\frac{1}{tm^{2}}\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|^{2}\sqrt{\frac{2v}{t}\log\left(\frac{2d^{2}}{\delta}\right)}+o\left(\frac{1}{t\sqrt{t}}\right)\enspace.$

Using the obtained bound on
$\|\mathbf{A}_{t}^{-1}-\mathbb{E}(\mathbf{A}_{t})^{-1}\|$ yields

$\displaystyle\frac{f_{\mathcal{Z}}(\mathbf{A}_{t})}{f_{\mathcal{Z}}(mt\times\Sigma_{\mathcal{Z}}(\lambda^{\star}))}-1$
$\displaystyle\leq\frac{mtL}{d^{2}}\times\left(\frac{1}{tm^{2}}\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|^{2}\sqrt{\frac{2v}{t}\log\left(\frac{2d^{2}}{\delta}\right)}+o\left(\frac{1}{t\sqrt{t}}\right)\right)$
$\displaystyle\leq\frac{L}{md^{2}}\left\|\Sigma_{\mathcal{Z}}(\lambda^{\star})^{-1}\right\|^{2}\sqrt{\frac{2v}{t}\log\left(\frac{2d^{2}}{\delta}\right)}+o\left(\frac{1}{\sqrt{t}}\right)\enspace,$

By noticing that
$f_{\mathcal{Z}}(mt\times\Sigma_{\mathcal{Z}}(\lambda^{\star}))\leq
f_{\mathcal{Z}}(\mathbf{A}^{\star}_{t})$ and by using Lemma C.5, the result
holds. ∎

##  Appendix D Variance analysis

#### Star graph.

The covariance matrix of the star graph can be bounded as follows:

$\displaystyle\mathrm{Var}(\mathbf{A}_{1})\preceq m\times
P\cdot\mathbf{I}+(n_{\mathrm{S}}-1)(n_{\mathrm{S}}-2)M\cdot\mathbf{I}+n_{\mathrm{S}}(n_{\mathrm{S}}-1)N\cdot\mathbf{I}\enspace.$

Since the star graph of $m$ edges has a number of nodes
$n_{\mathrm{S}}=m/2+1$, we have

$\displaystyle\|\mathrm{Var}(\mathbf{A}_{1})\|\leq m\times P+(M+N)\times
O\left(m^{2}\right)\enspace.$

#### Complete graph.

As for the star graph,

$\displaystyle\mathrm{Var}(\mathbf{A}_{1})\preceq m\times
P\cdot\mathbf{I}+n_{\mathrm{Co}}(n_{\mathrm{Co}}-1)(n_{\mathrm{Co}}-2)M\cdot\mathbf{I}+n_{\mathrm{Co}}(n_{\mathrm{Co}}-1)(n_{\mathrm{Co}}-1)N\cdot\mathbf{I}\enspace.$

Since the complete graph of $m$ edges has a number of nodes
$n_{\mathrm{Co}}=\left(1+\sqrt{4m+1}\right)/2$, we have

$\displaystyle\|\mathrm{Var}(\mathbf{A}_{1})\|\leq m\times P+(M+N)\times
O\left(m\sqrt{m}\right)\enspace.$

#### Circle graph.

Again,

$\displaystyle\mathrm{Var}(\mathbf{A}_{1})\preceq m\times
P\cdot\mathbf{I}+2n_{\mathrm{Ci}}M\cdot\mathbf{I}+4n_{\mathrm{Ci}}N\cdot\mathbf{I}\enspace.$

Since the circle graph of $m$ edges has a number of nodes
$n_{\mathrm{Ci}}=m/2$, we have

$\displaystyle\|\mathrm{Var}(\mathbf{Z}_{1})\|\leq m\times P+(M+N)\times
O\left(m\right)\enspace.$

#### Matching graph.

Finally,

$\displaystyle\mathrm{Var}(\mathbf{A}_{1})\preceq m\times
P\cdot\mathbf{I}+n_{\mathrm{M}}N\cdot\mathbf{I}\enspace.$

Since the matching graph of $m$ edges has a number of nodes
$n_{\mathrm{M}}=m$, we have

$\displaystyle\|\mathrm{Var}(\mathbf{A}_{1})\|\leq m\times P+m\times
N\enspace.$

##  Appendix E Generalization

In this section, we provide some insights into the generalization to broader
reward settings.

###  E.1 When $\mathbf{M}_{\star}$ is not symmetric

Consider the same graphical bilinear bandit setting as the one explained in
the paper with the only difference that $\mathbf{M}_{\star}$ is not symmetric.
We recall here that in the graph $\mathcal{G}=(V,E)$ associated to the
graphical bilinear bandit setting, $(i,j)\in E$ if and only if $(j,i)\in E$.
Hence, for a given allocation $(x^{(1)},\dots,x^{(n)})\in\mathcal{X}^{n}$, one
can write the associated expected global reward as follows :

$\displaystyle\sum_{(i,j)\in E}x^{(i)\top}\mathbf{M}_{\star}x^{(j)}$
$\displaystyle=\sum_{i=1}^{n}\sum_{\begin{subarray}{c}j\in\mathcal{N}(i)\\\
j>i\end{subarray}}x^{(i)\top}\mathbf{M}_{\star}x^{(j)}+x^{(j)\top}\mathbf{M}_{\star}x^{(i)}$
$\displaystyle=\sum_{i=1}^{n}\sum_{\begin{subarray}{c}j\in\mathcal{N}(i)\\\
j>i\end{subarray}}x^{(i)\top}\mathbf{M}_{\star}x^{(j)}+\left(x^{(j)\top}\mathbf{M}_{\star}x^{(i)}\right)^{\top}$
$\displaystyle=\sum_{i=1}^{n}\sum_{\begin{subarray}{c}j\in\mathcal{N}(i)\\\
j>i\end{subarray}}x^{(i)\top}\mathbf{M}_{\star}x^{(j)}+x^{(i)\top}\mathbf{M}_{\star}^{\top}x^{(j)}$
$\displaystyle=\sum_{i=1}^{n}\sum_{\begin{subarray}{c}j\in\mathcal{N}(i)\\\
j>i\end{subarray}}x^{(i)\top}\left(\mathbf{M}_{\star}x^{(j)}+\mathbf{M}_{\star}^{\top}x^{(j)}\right)$
$\displaystyle=\sum_{i=1}^{n}\sum_{\begin{subarray}{c}j\in\mathcal{N}(i)\\\
j>i\end{subarray}}x^{(i)\top}\left(\mathbf{M}_{\star}+\mathbf{M}_{\star}^{\top}\right)x^{(j)}\enspace.$

Let us denote
$\bar{\mathbf{M}}_{\star}=\mathbf{M}_{\star}+\mathbf{M}_{\star}^{\top}$. One
can notice that $\bar{\mathbf{M}}_{\star}$ is symmetric. Solving the graphical
bilinear bandit with the matrix $\bar{\mathbf{M}}_{\star}$ is exactly what we
propose throughout the main paper.

###  E.2 When the reward captures more information than the interactions
between agents

Consider the real world problems introduced in the paper, but with the
difference that instead of a reward only related to the interaction between
two neighboring agents/nodes, there is an additional term that informs about
the absolute quality of the arm chosen by the agent itself. More formally we
consider the following reward $r^{(i,j)}_{t}$ for the node $i$:

$\displaystyle
r^{(i,j)}_{t}=x^{(i)\top}_{t}\mathbf{M}_{\star}x^{(j)}_{t}+x^{(i)\top}_{t}\beta_{\star}+\eta^{(i,j)}_{t}\enspace.$

where $\beta_{\star}\in\mathbb{R}^{d}$ is a second unknown parameter that
allows to capture the quality of the arm chosen by the node $i$ independently
of its neighbors.

In order to add a constant term in the reward, let us construct the set
$\tilde{\mathcal{X}}\subset\mathbb{R}^{d+1}$ such that each arm
$x\in\mathcal{X}$ is associated to a new arm $\tilde{x}\in\tilde{\mathcal{X}}$
defined as $\tilde{x}^{\top}=(x^{\top},1)$. Moreover, let us define the matrix
$\tilde{\mathbf{M}}^{\star}\in\mathbb{R}^{(d+1)\times(d+1)}$ as follows:

$\displaystyle\tilde{\mathbf{M}}_{\star}=\begin{pmatrix}\left[\begin{matrix}&&\\\
&\mathbf{M}_{\star}&\\\ &&\end{matrix}\right]\left[\begin{matrix}\\\
\beta_{\star}\\\ \\\ \end{matrix}\right]\\\\[17.07182pt]
\left[\begin{matrix}0\ \ \ &\cdots&\ \ \
0\end{matrix}\right]\end{pmatrix}\enspace.$

One can easily verify that for any edge $(i,j)\in E$ and any time step $t$,
the reward $r^{(i,j)}_{t}$ can now be written as follows:

$\displaystyle r^{(i,j)}_{t}$
$\displaystyle=\tilde{x}^{(i)\top}_{t}\tilde{\mathbf{M}}_{\star}\tilde{x}^{(j)}_{t}+\eta^{(i,j)}_{t}\enspace,$

which leads to the same graphical bilinear bandit setting explained in Section
3, this time in dimension $d+1$ instead of $d$. Hence, all the previous
results hold for this more general graphical bilinear bandit problem, provided
any dependence in $d$ is modified to $d+1$.

##  Appendix F Computing $\mu^{\star}$

In Algorithm 2, we need to find the solution $\mu_{\star}$ of
$\min_{\mu\in\mathcal{S}_{\mathcal{X}}}h_{\mathcal{X}}(\mu)$. In fact we need
$\mu_{\star}$ to sample from it. We show that for any set $X$, the function
$h_{X}$ is convex and we use the Frank-Wolfe algorithm (Frank et al.,, 1956)
to compute $\mu_{\star}$ and $\lambda_{\star}$. The convergence of the
algorithm has been proven in Damla Ahipasaoglu et al., (2008). Note that one
can only compute $\mu_{\star}$ or $\lambda_{\star}$ to obtain the other one
thanks to C.2.

######  Proposition F.1.

Let $d>0$, for any set $X\subset\mathbb{R}^{d}$, $h_{X}$ is convex.

###### Proof.

Let $(\lambda,\lambda^{\prime})\in\mathcal{S}_{X}^{2}$ be two distributions in
$\mathcal{S}_{X}$. If either $\Sigma_{X}(\lambda)$ or
$\Sigma_{X}(\lambda^{\prime})$ are not invertible, then for any $t\in[0,1]$
one has

$h_{X}(t\lambda+(1-t)\lambda^{\prime})\leq
th_{X}(\lambda)+(1-t)h_{X}(\lambda^{\prime})=+\infty\enspace.$

Otherwise, for $t\in[0,1]$, we define the positive definite matrix
$\mathbf{Z}(t)\in\mathbb{R}^{d\times d}$ as follows:

$\mathbf{Z}(t)=t\Sigma_{X}(\lambda)+(1-t)\Sigma_{X}(\lambda^{\prime})\enspace.$

Simple linear algebra (Petersen and Pedersen,, 2012) yields

$\frac{\partial\mathbf{Z}(t)^{-1}}{\partial
t}=\mathbf{Z}(t)^{-1}\frac{\partial\mathbf{Z}(t)}{\partial
t}\mathbf{Z}(t)^{-1}\enspace.$

Using this result and the fact that $\partial^{2}\mathbf{Z}(t)/\partial
t^{2}=0$, we obtain

$\frac{\partial^{2}\mathbf{Z}(t)^{-1}}{\partial
t^{2}}=2\mathbf{Z}(t)^{-1}\frac{\partial\mathbf{Z}(t)}{\partial
t}\mathbf{Z}(t)^{-1}\frac{\partial\mathbf{Z}(t)}{\partial
t}\mathbf{Z}(t)^{-1}\enspace.$

Therefore, for any $x\in X$,

$\displaystyle\frac{\partial^{2}x^{\top}\mathbf{Z}(t)^{-1}x}{\partial t^{2}}$
$\displaystyle=2x^{\top}\mathbf{Z}(t)^{-1}\frac{\partial\mathbf{Z}(t)}{\partial
t}\mathbf{Z}(t)^{-1}\frac{\partial\mathbf{Z}(t)}{\partial
t}\mathbf{Z}(t)^{-1}x$
$\displaystyle=2\left(\frac{\partial\mathbf{Z}(t)}{\partial
t}\mathbf{Z}(t)^{-1}x\right)^{\top}\mathbf{Z}(t)^{-1}\left(\frac{\partial\mathbf{Z}(t)}{\partial
t}\mathbf{Z}(t)^{-1}x\right)$ $\displaystyle\geq 0\enspace,$

which shows convexity for any fixed $x\in X$. The final results yields from
the fact that $h_{X}$ is a maximum over convex functions.

∎

##  Appendix G Additional experiment and information

We define the set of arms $\mathcal{X}\subset\mathbb{R}^{5}$ that is made of
$|\mathcal{X}|=100$ node-arms randomly sampled from a multivariate
5-dimensional Gaussian distribution $\mathcal{N}(0,I)$ and then normalized so
that $\|x\|=1$ for all $x\in\mathcal{X}$. In all the figures the results are
averaged over 100 random repetitions of the experiments.

We propose to validate our insight and compute the evolution of
$\|\mathrm{Var}\left(\mathbf{A}_{1}\right)\|$ for the three types of graphs
(star, complete and circle) and different number of edges. The results are
shown in Figure 2. One can notice that we retrieve the $O(m^{2})$ dependence
of the variance for the star graph, the $O(m\sqrt{m})$ for the complete graph
and the linear dependence $O(m)$ for the circle graph.

Figure 2: Evolution of the variance according to the number of edges and the
type of graph (star, complete, circle), the variance being averaged over 100
repetitions.

#### Machine used for all the experiments.

Intel(R) Xeon(R) CPU E5-2667 v4 @ 3.20GHz - 24 CPUs used.