To see the last equivalence, assume $z \in S^*(\theta)\setminus S^*(\lambda)$, then there exists $x \in \mathcal{Z}$ such that $E_{z-\lambda}vec(\lambda) \prec E_{x-\lambda}vec(\lambda)$, i.e. $E_{z-x}vec(\lambda) \prec 0_d$. If $x \in S^*(\theta)$ then (75) follows.

Further assume $x \notin S^*(\theta)$, then, either there exists $x' \in S^*(\theta)$ such that $E_{x-\lambda}vec(\lambda) \prec E_{x'-\lambda}vec(\lambda)$ then by transitivity $E_{z-\lambda}vec(\lambda) \prec E_{x'-\lambda}vec(\lambda)$ and $z, x' \in S^*(\theta)$ or for all $x' \in S^*(\theta)$, $E_{x-\lambda}vec(\lambda) \not\prec E_{x'-\lambda}vec(\lambda)$. Put together we have $(z \in S^*(\theta)\setminus S^*(\lambda))$ implies (75). Similarly, for $z \in S^*(\lambda)\setminus S^*(\theta)$, it is direct by Pareto-optimality (in the instance $\lambda$) that for all $x \in S^*(\theta)$, $E_{z-\lambda}vec(\lambda) \not\prec E_{x-\lambda}vec(\lambda)$.

Reversely, if (75) holds, then in the case a), there exists $z \in S^*(\theta)\setminus S^*(\lambda)$. If b) holds, then one cannot have $S^*(\theta) = S^*(\lambda)$. Indeed if we had $S^*(\theta) = S^*(\lambda)$, then it would follow that there exists $z \in \mathcal{Z}\setminus S^*(\theta) = \mathcal{Z}\setminus S^*(\lambda)$ such that for all $z \in S^*(\lambda)$, $E_{z-\lambda}vec(\lambda) \not\prec E_{z-\lambda}vec(\lambda)$, that in the instance $\lambda$, there would be a sub-optimal arm that is not dominated by any other Pareto-optimal arm, which is not possible. This concludes the reverse inclusion. Therefore,

$$ \mathrm{Alt}(\theta) := \underbrace{\{\lambda \mid z, x \in S^*(\theta) : E_{z-x} \mathrm{vec}(\lambda) \prec 0_d\}}_{A_\theta} \cup \underbrace{\{\lambda \mid \exists z \in \mathcal{Z} \setminus S^*(\theta) : \forall x \in S^*(\theta), E_{z-x} \mathrm{vec}(\lambda) \not\prec 0_d\}}_{B_\theta}, \quad (76) $$

then observe that

$$ A_{\theta} := \bigcup_{z \in S^*(\theta)} \bigcup_{x \in S^*(\theta) \setminus \{z\}} \{\lambda : E_{z-x} \mathrm{vec}(\lambda) \prec 0_d\}, \quad (77) $$

and for $z, x$ fixed, $\{\lambda \mid E_{z-x}vec(\lambda) \prec 0_d\}$ is a convex set. On the other side, for $z$ fixed, if $(E_{z-x}vec(\lambda) \not\prec 0_d, \forall x \in S^*(\theta))$ holds then there exists $\bar{d}^z \in \{1, \dots, d\}^{S^*(\theta)}$ such that for all $x \in S^*(\theta)$, $\langle E_{z,x}(d^z(x)), \mathrm{vec}(\lambda) \rangle \ge 0$. Putting these displays together yields

$$ B_{\theta} = \bigcup_{z \in \mathcal{Z} \setminus S^*(\theta)} \bigcup_{\bar{d}^z \in [d]^{S^*(\theta)}} \{\lambda : \forall x \in S^*(\theta) : \langle E_{z,x}(d^z(x)), \mathrm{vec}(\lambda) \rangle \ge 0\}, \quad (78) $$

then remark that for $z, \bar{d}^z$ fixed, the set $\{\lambda : \forall x \in S^*(\theta) : \langle E_{z,x}(d^z(x)), \mathrm{vec}(\lambda) \rangle \ge 0\}$ is convex. Finally, combining (78) with (77) yields the claimed result. $\square$

The following lemma allows to upper bound the probability that a multi-variate normal vector belongs to a convex set. This is proven in Lu and Li (2009) where the authors attributed it to Kuelbs et al. (1994).

**Lemma 27.** Let $\Sigma$ be a variance-covariance matrix, $\mu \in \mathbb{R}^n$, let $X \sim N(\mu, \Sigma)$ and $C \subset \mathbb{R}^n$, a convex set. Then, it holds that

$$ \mathbb{P}(X \in C) \le \frac{1}{2} \exp \left( - \inf_{x \in C} \frac{1}{2} \|x - \mu\|_{\Sigma^{-1}}^2 \right). $$

The following lemma can be found in algebra books, we prove it for the sake of self-containedness. For two $n \times n$ matrices $A, B$ we write $A \ge B$ iff $A - B$ is positive semi-definite i.e for all $x \in \mathbb{R}^n$, $x^\top(A - B)x \ge 0$.

**Lemma 28.** Let $A, B$ be two symmetric order $n$ non-singular matrices. If $A - B \le 0$ then $B^{-1} - A^{-1} \le 0$.

*Proof.* As $A$ is symmetric, there exists $P$ such that $A = P^2$ which is denoted by $A^{1/2} = P$. We have

$$ \begin{align*} A - B \le 0 &\iff \forall x \in \mathbb{R}^n, x^\top A x - x^\top B x \le 0 \\ &\iff \forall x \in \mathbb{R}^n, (A^{1/2} x)^\top (A^{1/2} x) - (A^{1/2} x)^\top A^{-1/2} B A^{-1/2} (A^{1/2} x) \le 0 \\ &\iff I_n - \underbrace{A^{-1/2} B A^{-1/2}}_{M} \le 0 \quad (\text{as } A^{1/2} \text{ is invertible}) \\ &\iff \forall x \in \mathbb{R}^n, (M^{1/2} x)^\top M^{-1} (M^{1/2} x) - (M^{1/2} x)^\top (M^{1/2} x) \le 0 \quad (\text{as } M \text{ is symmetric and invertible}) \\ &\iff M^{-1} - I_n \le 0, \end{align*} $$

then by plugging in the expression of $M$, we have

$$ \begin{align*} A - B \le 0 &\iff (A^{-1/2} B A^{-1/2})^{-1/2})^{-1} - I_n \le 0 \\ &\iff \forall x \in \mathbb{R}^n, x^\top A^{1/2} B^{-1} A^{1/2} x - (A^{1/2} x)^\top A^{-1} (A^{1/2} x) \le 0 \\ &\iff B^{-1} - A^{-1} \le 0, \end{align*} $$

where the last inequality follows again as $A^{1/2}$ is invertible which follows as $A$ is symmetric non-singular. $\square$
In the result below, we show that the Euclidean norm is $\alpha$-exp concave over a bounded domain.

**Lemma 29.** The function $\lambda \mapsto \exp(-\alpha \|\text{vec}(\theta' - \lambda)\|_{\Sigma^{-1} \otimes aa^T}^2)$ defined a bounded domain $D$ is concave for $\alpha \in (0, \frac{1}{2 \max_{\lambda \in D} \|\text{vec}(\theta' - \lambda)\|_{\Sigma^{-1} \otimes aa^T}^2}]$.

*Proof.* Let $\tilde{V} = \Sigma^{-1} \otimes aa^T$ and define $g(\lambda) := \exp(-\alpha \|\theta' - \lambda\|_{\Sigma^{-1} \otimes aa^T}^2)$ for $\lambda \in \text{vec}(D)$, with $\theta' := \text{vec}(\theta')$. Direct calculation of the Hessian of $g$ gives

$$ \nabla^2 g(\lambda) = (-2\alpha \tilde{V} + 4\alpha^2 \tilde{V}(\lambda - \theta')(\theta' - \theta')^T \tilde{V}) g(\lambda) $$

further, observe that,

$$ \begin{aligned} u^T (-2\alpha \tilde{V} + 4\alpha^2 \tilde{V}(\lambda - \theta')(\lambda - \theta')^T \tilde{V}) u &= -2\alpha \|u\|_{\tilde{V}}^2 + 4\alpha^2 (u^T \tilde{V}(\lambda - \theta'))^2 \\ &\leq -2\alpha \|u\|_{\tilde{V}}^2 + 4\alpha^2 \|u\|_{\tilde{V}}^2 \|\lambda - \theta'\|_{\tilde{V}}^2, \end{aligned} $$

which proves the claimed statement as the Hessian is negative semi-definite for $\alpha \in (0, 1/(2 \max_{\lambda \in D} \|\text{vec}(\theta' - \lambda)\|_{\Sigma^{-1} \otimes aa^T}^2))$. $\square$

The following lemma expresses the difference in estimation of the empirical mean between two consecutive rounds.

**Lemma 30.** Let $\tilde{V}_t := \Sigma^{-1} \otimes V_t$, then, it holds that

$$ \tilde{V}_t \text{vec}(\hat{\theta}_t - \hat{\theta}_{t-1}) = [\Sigma^{-1} \otimes a_t] (\varepsilon_t + [I_d \otimes a_t^T] \text{vec}(\theta - \hat{\theta}_{t-1})). $$

*Proof.* By Sherman-Morrison's formula, with $S_t = \sum_{s \in [t]} \text{vec}(a_s X_s^T)$

$$ \begin{aligned} \text{vec}(\hat{\theta}_t) &= [I_d \otimes (V_{t-1} + a_t a_t^T)^{-1}] S_t \\ &= \left[ I_d \otimes \left( V_{t-1}^{-1} - \frac{V_{t-1}^{-1} a_t a_t^T V_{t-1}^{-1}}{1 + a_t^T V_{t-1}^{-1} a_t} \right) \right] S_t \\ &= \text{vec}(\hat{\theta}_{t-1}) - \left[ I_d \otimes \left( \frac{V_{t-1}^{-1} a_t a_t^T V_{t-1}^{-1}}{1 + a_t^T V_{t-1}^{-1} a_t} \right) \right] S_t + [I_d \otimes V_{t-1}^{-1}] \text{vec}(a_t X_t^T) \\ &= \text{vec}(\hat{\theta}_{t-1}) - \frac{[I_d \otimes (V_{t-1}^{-1} a_t a_t^T)][I_d \otimes V_{t-1}^{-1}]}{1 + a_t^T V_{t-1}^{-1} a_t} S_t + [I_d \otimes V_{t-1}^{-1}] \text{vec}(a_t X_t^T) \\ &= \text{vec}(\hat{\theta}_{t-1}) - \frac{[I_d \otimes (V_{t-1}^{-1} a_t a_t^T)] \text{vec}(\hat{\theta}_{t-1}) + [I_d \otimes (V_{t-1}^{-1} a_t a_t^T V_{t-1}^{-1})] \text{vec}(a_t X_t^T)}{1 + a_t^T V_{t-1}^{-1} a_t} + [I_d \otimes V_{t-1}^{-1}] \text{vec}(a_t X_t^T) \end{aligned} $$

we recall that for two vector $a, X$,

$$ \text{vec}(aX^T) = [I_d \otimes a]X \quad \text{and} \quad [A \otimes B][C \otimes D] = (AC) \otimes (BD) $$

then

$$ \begin{aligned} [I_d \otimes (V_{t-1}^{-1} a_t a_t^T V_{t-1}^{-1})] \text{vec}(a_t X_t^T) &= [I_d \otimes (V_{t-1}^{-1} a_t a_t^T V_{t-1}^{-1})][I_d \otimes a_t] X_t \\ &= (a_t^T V_{t-1}^{-1} a_t) [I_d \otimes (V_{t-1}^{-1} a_t)] X_t. \end{aligned} $$

Proceeding similarly yields

$$ \begin{aligned} [I_d \otimes V_{t-1}^{-1}] \text{vec}(a_t X_t^T) &= [I_d \otimes V_{t-1}^{-1}][I_d \otimes a_t] X_t \\ &= [I_d \otimes V_{t-1}^{-1} a_t] X_t, \end{aligned} $$

therefore,

$$ \text{vec}(\hat{\theta}_t) = \text{vec}(\hat{\theta}_{t-1}) + \frac{[I_d \otimes (V_{t-1}^{-1} a_t)] X_t - [I_d \otimes (V_{t-1}^{-1} a_t a_t^T)] \text{vec}(\hat{\theta}_{t-1})}{1 + a_t^T V_{t-1}^{-1} a_t}. $$
By Kronecker product on the left, we have

$$
\begin{align*}
[\Sigma^{-1} \otimes V_t] \mathrm{vec}(\hat{\theta}_t - \hat{\theta}_{t-1})
  &= \frac{[\Sigma^{-1} \otimes V_t][I_d \otimes V_{t-1}^T a_t] X_t - [\Sigma^{-1} \otimes V_t][I_d \otimes (V_{t-1}^T a_t a_t^T)] \mathrm{vec}(\hat{\theta}_{t-1})}{1 + a_t^T V_{t-1}^T a_t} \\
  &= \frac{[\Sigma^{-1} \otimes (V_t V_{t-1}^T a_t)] X_t - [\Sigma^{-1} \otimes (V_t V_{t-1}^T a_t a_t^T)] \mathrm{vec}(\hat{\theta}_{t-1})}{1 + a_t^T V_{t-1}^T a_t}.
\end{align*}
$$

Then observe that

$$
\begin{align*}
[\Sigma^{-1} \otimes (V_t V_{t-1}^{-1} a_t)] X_t
  &= [\Sigma^{-1} \otimes ((I_h + a_t a_t^T V_{t-1}^{-1}) a_t)] X_t \\
  &= ([\Sigma^{-1} \otimes a_t] + (a_t^T V_{t-1}^{-1} a_t) [\Sigma^{-1} \otimes a_t]) X_t \\
  &= (1 + a_t^T V_{t-1}^{-1} a_t) [\Sigma^{-1} \otimes a_t] X_t,
\end{align*}
$$

and similar developments yields

$$
\begin{align*}
[\Sigma^{-1} \otimes (V_t V_{t-1}^{-1} a_t a_t^T)] \mathrm{vec}(\hat{\theta}_{t-1})
  &= [\Sigma^{-1} \otimes ((I_h + a_t a_t^T V_{t-1}^{-1}) a_t a_t^T)] \mathrm{vec}(\hat{\theta}_{t-1}) \\
  &= ([\Sigma^{-1} \otimes (a_t a_t^T)] + (a_t^T V_{t-1}^{-1} a_t) [\Sigma^{-1} \otimes (a_t a_t^T)]) \mathrm{vec}(\hat{\theta}_{t-1}) \\
  &= (1 + a_t^T V_{t-1}^{-1} a_t) [\Sigma^{-1} \otimes (a_t a_t^T)] \mathrm{vec}(\hat{\theta}_{t-1})
\end{align*}
$$

so that we finally have

$$
\begin{align*}
\tilde{V}_t \mathrm{vec}(\hat{\theta}_t - \hat{\theta}_{t-1})
  &= [\Sigma^{-1} \otimes V_t] \mathrm{vec}(\hat{\theta}_t - \hat{\theta}_{t-1}) \\
  &= [\Sigma^{-1} \otimes a_t] X_t - [\Sigma^{-1} \otimes (a_t a_t^T)] \mathrm{vec}(\hat{\theta}_{t-1}) \\
  &= [\Sigma^{-1} \otimes a_t] X_t - ([\Sigma^{-1} \otimes a_t] [I_d \otimes a_t^T]) \mathrm{vec}(\hat{\theta}_{t-1}) \\
  &= [\Sigma^{-1} \otimes a_t] (X_t - (I_d \otimes a_t^T) \mathrm{vec}(\hat{\theta}_{t-1})).
\end{align*}
$$

We recall that we have $X_t := (I_d \otimes a_t^T)\mathrm{vec}(\theta) + \epsilon_t$ from the linear model. Thus, all put together, we recover the claimed result. □

**Lemma 9.** Let $\omega \in \mathbb{R}_+^K$. For any $\mu' := (\mu'_1 \dots \mu'_K)^\top \in \mathbb{R}^{K \times d}$, the following statement holds in the unstructured setting

$$
\inf_{\lambda \in \operatorname{Alt}(S^*(\mu'))} \|\operatorname{vec}(\lambda - \mu')\|_{\Sigma^{-1} \otimes \operatorname{diag}(\omega)}^2 = \inf_{\lambda \in \operatorname{Alt}(S^*(\mu')) \cap \{\lambda \mid \max_i \|\mu'_i - \lambda_i\|_{\Sigma^{-1}} < \epsilon\}} \|\operatorname{vec}(\lambda - \mu')\|_{\Sigma^{-1} \otimes \operatorname{diag}(\omega)}^2, \quad (16)
$$

$$
\text{where } \epsilon := \max(2 \max_{i \notin S} \max_{j \in S} \|\mu'_i - \mu'_j\|_{\Sigma^{-1}}, \max_{i,j \in S^2} \|\mu'_i - \mu'_j\|_{\Sigma^{-1}}), \text{ and } S = S^*(\mu').
$$

This result shows that the best response always exists in the unstructured setting as an inf over a compact subset of the alternative. Although its closed form is unknown in PSI, we show that it belongs to a ball centered at $\boldsymbol{\mu}'$ and whose radius depends also on $\boldsymbol{\mu}'$.

*Proof.* The idea of the proof is to show that when there is a ball $\mathbb{B}$ such that if an alternative parameter of $\boldsymbol{\mu}'$ does not belong to $\mathbb{B}$, then there is a parameter in $\operatorname{Alt}(S^*(\boldsymbol{\mu}')) \cap \mathbb{B}$ for which the transportation cost will be smaller. Let $\boldsymbol{\mu}' := (\mu'_1 \dots \mu'_K)^\top$, where $\mu'_i$ denotes the vector mean of arm $i$ and we let $\boldsymbol{\lambda} := (\lambda_1 \dots \lambda_K)^\top$ where $\lambda_i \in \mathbb{R}^d$. To ease notation, we let $S = S^*(\boldsymbol{\mu}')$. Introducing

$$
\operatorname{Alt}^{-}(S) := \bigcup_{i,j \in S^2 : i \neq j} W_{i,j} \quad \text{and} \quad \operatorname{Alt}^{+}(S) := \bigcup_{i \in S^c} V_i, \qquad (79)
$$

where $\Theta := \mathbb{R}^{K \times d}$, we define

$$
W_{i,j} := \{ \boldsymbol{\lambda} := (\lambda_1 \dots \lambda_K)^\top \in \Theta \mid \lambda_i \le \lambda_j \} \text{ and } V_i := \{ \boldsymbol{\lambda} := (\lambda_1 \dots \lambda_K)^\top \in \Theta \mid \exists i \in S \text{, } \forall j \in S, \lambda_i \ne \lambda_j \}.
$$

By simply expanding the expression, we have in the unstructured case,

$$
D(\omega, \lambda; \mu') := \|\mathrm{vec}(\lambda - \mu')\|_{\Sigma^{-1} \otimes \mathrm{diag}(\omega)}^2 = \sum_{i=1}^{K} \omega_i \|\lambda_i - \mu'_i\|_{\Sigma^{-1}}^2.
$$
Let $(i, j) \in S^2$ be fixed and $\lambda \in W_{i,j}$. Let $\alpha_{i,j} = \|\mu'_i - \mu'_j\|_{\Sigma^{-1}}$. If $\|\lambda_i - \mu'_i\|_{\Sigma^{-1}} > \alpha_{i,j}$ then define the instance $\tilde{\lambda}$ as

$$ \tilde{\lambda}_k := \begin{cases} \mu'_k & \text{if } k \notin \{i, j\}, \\ \mu'_j & \text{else,} \end{cases} $$

which satisfies $\tilde{\lambda}_i \le \tilde{\lambda}_j$, so $\tilde{\lambda} \in W_{i,j}$, and

$$ \begin{aligned} D(\omega, \tilde{\lambda}; \mu') &= \omega_i \|\mu'_i - \mu'_j\|_{\Sigma^{-1}}^2 + \sum_{k \notin \{i, j\}} \omega_k \|\lambda_k - \mu'_k\|_{\Sigma^{-1}}^2, \\ &= \omega_i \|\mu'_i - \mu'_j\|_{\Sigma^{-1}}^2 < \omega_i \|\lambda_i - \mu'_i\|_{\Sigma^{-1}}^2 \\ &< D(\omega, \lambda, \mu'), \end{aligned} $$

and further observe that $\max_{k \in [K]} \|\tilde{\lambda}_k - \mu'_k\|_{\Sigma^{-1}} \le \alpha_{i,j}$. We proceed similarly if $\|\lambda_j - \mu'_j\|_{\Sigma^{-1}} > \alpha_{i,j}$, by defining

$$ \tilde{\lambda}_k := \begin{cases} \mu'_k & \text{if } k \notin \{i, j\}, \\ \mu'_i & \text{else,} \end{cases} $$

and the same conclusion follows. So we have proved that there exists $\tilde{\lambda} \in W_{i,j}$ with $\max_{k \in [K]} \|\tilde{\lambda}_k - \mu'_k\|_{\Sigma^{-1}} \le \alpha_{i,j}$ and for which the transportation cost is not larger than that of $\lambda$. We prove a similar property for $\text{Alt}^+(S)$.

Now, fix $i \notin S$ and take $\lambda \in V_i$. Let $b_i = \max_{k \in S} \|\mu'_i - \mu'_k\|_{\Sigma^{-1}}$. If $\|\lambda_i - \mu'_i\|_{\Sigma^{-1}} > b_i$ then it suffices to define $\tilde{\lambda}$ as

$$ \tilde{\lambda}_p := \begin{cases} \mu'_p & \text{if } p \ne i, \\ \mu'_i & \text{else,} \end{cases} $$

for $\tilde{i} \in S$, to ensure that $\tilde{\lambda} \in V_i$ and

$$ \begin{aligned} D(\omega, \tilde{\lambda}; \mu') &= \omega_i \|\mu'_i - \mu'_k\|_{\Sigma^{-1}}^2 + \sum_{k \ne i} \omega_k \|\lambda_k - \mu'_k\|_{\Sigma^{-1}}^2 \\ &\le \omega_i b_i^2 \\ &< \omega_i \|\lambda_i - \mu'_i\|_{\Sigma^{-1}}^2 \\ &< D(\omega, \lambda; \mu'), \end{aligned} $$

where the second line also uses the fact that $\tilde{i} \in S$. Assume $\|\lambda_i - \mu'_i\|_{\Sigma^{-1}} < b_i$ and that $H_i := \{k \in S : \|\lambda_k - \mu'_k\|_{\Sigma^{-1}} > 2b_i\}$ is non-empty. Let us define the instance $\tilde{\lambda}$:

$$ \tilde{\lambda}_k := \begin{cases} \lambda_i & \text{if } k \in H_i, \\ \lambda_k & \text{if } k \in (S \cup \{i\}) \setminus H_i, \\ \mu'_k & \text{else.} \end{cases} $$

Since $\lambda \in V_i$, $\tilde{\lambda}$ as defined above satisfies $\tilde{\lambda}_i \not\le \tilde{\lambda}_k$, $k \in S$, that is $\tilde{\lambda} \in V_i$ and

$$ \begin{aligned} D(\omega, \tilde{\lambda}; \mu') &= \sum_{k \in (S \cup \{i\})} \omega_k \|\tilde{\lambda}_k - \mu'_k\|_{\Sigma^{-1}}^2 \\ &= \sum_{k \in (S \cup \{i\}) \setminus H_i} \omega_k \|\lambda_k - \mu_k\|_{\Sigma^{-1}}^2 + \sum_{k \in H_i} \omega_k \|\lambda_i - \mu'_k\|_{\Sigma^{-1}}^2 \\ &< \sum_{k \in (S \cup \{i\}) \setminus H_i} \omega_k \|\lambda_k - \mu'_k\|_{\Sigma^{-1}}^2 + \sum_{k \in H_i} \omega_k 4b_i^2, \end{aligned} $$

which follows since $\|\lambda_i - \mu'_i\|_{\Sigma^{-1}} \le \|\lambda_i - \mu'_i\|_{\Sigma^{-1}} + \|\mu'_i - \mu'_k\|_{\Sigma^{-1}} \le 2b_i$. Recalling that for $k \in H_i$,

$$ 4b_i < \|\lambda_k - \mu'_k\|_{\Sigma^{-1}}^2, $$
it follows that

$$D(\omega, \tilde{\lambda}; \mu') < D(\omega, \lambda; \mu').$$

Put together, we have proved that for all $\lambda \in V_i$, there exists $\tilde{\lambda} \in V_i$ whose transportation cost is not larger than that of $\lambda$ and which additionally satisfies: $\max_{k \in [K]} \|\tilde{\lambda}_k - \mu'_k\|_{\Sigma^{-1}} \le 2b_i$.

To conclude, let us define $\epsilon := \max(2 \max_{i \notin S} b_i, \max_{i,j \in S^2} \alpha_{i,j})$. Using what precedes, we have proved that

$$\inf_{\lambda \in \text{Alt}(S) \cap \{\lambda: \max_k \|\lambda_k - \mu'_k\|_{\Sigma^{-1}} < \epsilon\}} D(\omega, \lambda; \mu') \le \inf_{\lambda \in \text{Alt}(S)} D(\omega, \lambda; \mu'),$$

which proves the claimed result as we have $\text{Alt}(S) \cap \{\lambda: \max_k \|\lambda_k - \mu'_k\|_{\Sigma^{-1}} < \epsilon\} \subset \text{Alt}(S)$. $\square$

# J IMPLEMENTATION DETAILS AND ADDITIONAL EXPERIMENTS

After presenting the implementation details and computational cost in Appendix J.1, we display supplementary experiments in Appendix J.3. The datasets are described in Appendix J.2.

## J.1 Implementation Details and Complexity

### J.1.1 Setup

The algorithm is implemented mainly in standard C++17. Our complete code is available on the repository at https://github.com/cyrille-kone/psips with instructions to compile and run. We run our algorithms mainly on a personal laptop with 8GB RAM, 256GB SSD, and an Apple M2 CPU. The code is compiled with the GCC13. In the experiments, we include the “Oracle” algorithm, which computes the optimal weights of the underlying bandit instance. At each round, the oracle pulls an arm according to this optimal weights vector. Another algorithm included in our benchmark is the APE algorithm of Kone et al. (2023), which is an LUCB-type algorithm for PSI. GAPSI denotes the tracking algorithm of Crepon et al. (2024) for PSI.

In the following paragraphs, we analyze the time and memory complexity of the different parts of our algorithms, and we discuss its total time complexity.

### J.1.2 Time and memory complexity

**Compute Pareto Set** When $d = 2$, it is known that the Pareto set can be computed in an average $O(|Z| \log(|Z|))$ time complexity and a worst-case $O(|Z| \log(|Z|) + |Z| |p|)$ where $p$ is the size of the Pareto set. For a general dimension $d > 3$, the algorithm of Kung et al. (1975) achieves a time complexity of $O(|Z| \log(|Z|)^{d-2})$. Their algorithm is based on divide and conquer and in the case $d = 2$, it consists in sorting the arms along one coordinate, then by traversing the sorted items in decreasing order (for sorted dimension) each element is added to the set $S$ if it is note dominated by an element of $S$. The resulting set is then the Pareto set, and combining the initial sorting and the construction of the Pareto set, the worst case complexity will be $O(\log(|Z|) + |Z|p + p(1-p)/2)$, $p$ is the size of the Pareto set. The space complexity is $O(|Z|)$. This algorithm is described below.

To check if two parameters have the same Pareto set, in the worst case, one can compute their Pareto sets using the algorithms described earlier. However, it is possible (and more efficient) to check this condition without computing the Pareto sets of the parameters by simply checking from the definition of Alt.

**Update $\hat{\theta}_{t+1}$** To compute $\hat{\theta}_{t+1}$ in the linear setting, we used Sherman-Morrison formula to avoid computing the inverse of $V_t$. Indeed as $V_t = V_{t-1} + a_t a_t^\top$ we have by Sherman-Morrison formula

$$V_t^{-1} = V_{t-1}^{-1} - \frac{V_{t-1}^{-1} a_t (V_{t-1}^{-1} a_t)^\top}{1 + \|a_t\|_{V_{t-1}^{-1}}^2},$$

which can be computed in time $O(h^2)$. The computation of $Z_{t+1}$ is done in time $O(hd)$, so that the total cost of computing the least-squares estimate is dominated by the matrix product $V_{t+1}^{-1} Z_{t+1}$, which can be done naively in $O(h^2d)$.
**Generate Random Samples** To pull from $\Pi_t$ we generate samples from $\mathcal{N}(0_{dh}, I_{dh})$ (whose average cost is $\mathcal{O}(dh)$) and we compute a Cholesky decomposition of the $\Sigma \otimes V_t^{-1}$ which we do by computing the Cholesky of $\Sigma$ and that of $V_t^{-1}$. Note that in the unstructured setting, the Cholesky of $V_t^{-1}$ is trivial to compute. In summary, in the worst case, the total time complexity of this step is $\mathcal{O}(\max(d^2h^2, h^3))$ in the linear setting and $\mathcal{O}(d^2h)$ in the unstructured setting.

**Miscellaneous** For the other parts of the algorithm, the cost of a step of AdaHedge (De Rooij et al., 2014) is $\mathcal{O}(K)$, same for the memory complexity. When $\Theta$ is defined as a ball, the cost of the verification $\lambda \in \Theta$ is simply the cost of computing the norm of $\Theta$. For Euclidean norms, it will be $\mathcal{O}(dh)$. The computation of $M(t, \delta)$ is $\mathcal{O}(1)$, as we only need to compute $(q(\delta, i)_{i \in [K]})$ once at initialization.

**Summary** Combining the displays above, the time complexity of our algorithm at each step is

$$ \mathcal{O}(|\mathcal{Z}| \log(|\mathcal{Z}|)^{\max(1, d-2)}) + \mathcal{O}(\max(d^2h^2, h^3)). $$

In the unstructured setting, it can be reduced to

$$ \mathcal{O}(\max(K^2d, K \log(K)^{\max(1, d-2)})). $$

## J.2 Datasets

**COV-BOOST** This dataset has been publicly released in Munro et al. (2021) based on a phase 2 booster trial for Covid19 vaccine. Kone et al. (2023) further processed the data and extracted the parameters of a bandit instance with 20 arms and $d=3$. We use the values tabulated by the authors to simulate our algorithms; see Tables 3 and 4.

<table><caption>**Table 3:** Table of the empirical arithmetic mean of the log-transformed immune response for three immunogenicity indicators. Each acronym corresponds to a vaccine. There are two groups of arms corresponding to the first 2 doses: one with prime BNT/BNT (BNT as first and second dose) and the second with prime ChAd/ChAd (ChAd as first and second dose). Each row in the table gives the values of the 3 immune responses for an arm (i.e., a combination of three doses).</caption><thead><tr><th rowspan="2">Dose 1/Dose 2</th><th rowspan="2">Dose 3 (booster)</th><th colspan="3">Immune response</th></tr><tr><th>Anti-spike IgG</th><th>NT<sub>50</sub></th><th>cellular response</th></tr></thead><tbody><tr><td rowspan="10">Prime BNT/BNT</td><td>ChAd</td><td>9.50</td><td>6.86</td><td>4.56</td></tr><tr><td>NVX</td><td>9.29</td><td>6.64</td><td>4.04</td></tr><tr><td>NVX Half</td><td>9.05</td><td>6.41</td><td>3.56</td></tr><tr><td>BNT</td><td>10.21</td><td>7.49</td><td>4.43</td></tr><tr><td>BNT Half</td><td>10.05</td><td>7.20</td><td>4.36</td></tr><tr><td>VLA</td><td>8.34</td><td>5.67</td><td>3.51</td></tr><tr><td>VLA Half</td><td>8.22</td><td>5.46</td><td>3.64</td></tr><tr><td>Ad26</td><td>9.75</td><td>7.27</td><td>4.71</td></tr><tr><td>m1273</td><td>10.43</td><td>7.61</td><td>4.72</td></tr><tr><td>CVn</td><td>8.94</td><td>6.19</td><td>3.84</td></tr><tr><td rowspan="10">Prime ChAd/ChAd</td><td>ChAd</td><td>7.81</td><td>5.26</td><td>3.97</td></tr><tr><td>NVX</td><td>8.85</td><td>6.59</td><td>4.73</td></tr><tr><td>NVX Half</td><td>8.44</td><td>6.15</td><td>4.59</td></tr><tr><td>BNT</td><td>9.93</td><td>7.39</td><td>4.75</td></tr><tr><td>BNT Half</td><td>8.71</td><td>7.20</td><td>4.91</td></tr><tr><td>VLA</td><td>7.51</td><td>5.31</td><td>3.96</td></tr><tr><td>VLA Half</td><td>7.27</td><td>4.99</td><td>4.02</td></tr><tr><td>Ad26</td><td>8.62</td><td>6.33</td><td>4.66</td></tr><tr><td>m1273</td><td>10.35</td><td>7.77</td><td>5.00</td></tr><tr><td>CVn</td><td>8.29</td><td>5.92</td><td>3.87</td></tr></tbody></table>

**Network on Chip** We use the dataset studied by Almer et al. (2011) and publicly released by Zuluaga et al. (2013) and available at http://www.spiral.net/software/pal.html. We apply further preprocessing by normalizing the
<table><caption>**Table 4:** Pooled variance of each group.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Immune response</th></tr><tr><th>Anti-spike IgG</th><th>NT<sub>50</sub></th><th>cellular response</th></tr></thead><tbody><tr><th>Pooled sample variance</th><td>0.70</td><td>0.83</td><td>1.54</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Means of each arm in the NoC instance with normalized features.</figcaption></figure>

features to extract a linear instance with parameter $\theta^*$ given below.

$$ \theta^* = \begin{pmatrix} -3.08665453 & -3.35487744 \\ -3.66027623 & 0.19333635 \\ -2.68963781 & -1.39779755 \\ -7.90670356 & -4.44360318 \end{pmatrix}. $$

In Figure 7 we plot the resulting PSI instance.

**Reproducibility** We provide below the parameters used in the experiments reported in the main paper. As explained in the experiment section, since their calibration relies on overly conservative union bounds, we used $M(t, \delta) = \frac{1}{\delta} \log(t/\delta)$ and $c(t, \delta) = 1 + \frac{\log(\log(t))}{\log(1/\delta)}$ for which we noticed a negligible empirical error for $\delta \le 0.1$ in our experiments. For AdaHedge, we used the implementation of De Rooij et al. (2014), translating the pseudo-code provided in their paper. In the supplementary material, we provide our C++ implementation together with the scripts to run each experiment. The complete scripts to reproduce our experiments will be available on an anonymous repository with the instructions to compile and run the algorithms.

## J.3 Supplementary Experiments

We provide additional experimental results. In particular, we conduct experiments to evaluate the performance of our algorithm in more complex scenarios. These experiments include instances in higher dimension $d$ or with a larger number of arms $K$. Moreover, we also report the probability of error in the covid19 experiment. Each experiment is repeated on 100 independent runs, and the quantities reported are averaged statistics. Due to its high runtime, we omit GAPSI of Crepon et al. (2024) from these experiments.
<figure><img src="image_9.png" /><figcaption>**Figure 8:** Empirical stopping time on 10-armed random Gaussian instances in dimension $d \in \{3, 4, 5, 6\}$ (left to right).</figcaption></figure>

**Higher Dimension $d$** In this experiment, we evaluate our algorithm in the unstructured setting on 100 random Gaussian instances with higher dimension $d \in \{3, 4, 5, 6\}$, $K = 10$ and $\Sigma = I_d/2$. Each arm is drawn from $U([-1, 1]^d)$, and to have reasonable runtime, we reject instances whose complexity $H(\theta)$ is larger than 500. The average Pareto set size was respectively 4.42, 6.78, 8.57, 9.21 in dimensions $d \in \{3, 4, 5, 6\}$. Using $\delta = 0.01$, we observe a negligible empirical error.

In Figure 8, we see that PSIPS has competitive empirical performance for higher dimension $d$. It consistently outperforms APE and performs on par with uniform sampling. The good performance of uniform sampling in this experiment can be attributed to the large size of the Pareto set, which contains most of the arms. We conjecture that $\omega^*(\theta)$ is close to the uniform allocation on those instances.

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Empirical stopping time on 2-dimensional random Bernoulli instances with $K \in \{10, 15, 25, 40\}$ arms (left to right).</figcaption></figure>

**Large Number of Arms $K$** We benchmark the algorithms on instances with a larger number of arms. We report the results on 100 random Bernoulli instances with $K \in \{10, 15, 25, 40\}$ arms in dimension $d = 2$. The means of the instances are drawn from $U([0.2, 0.9]^{K \times d})$, and we use $\Sigma = I_2/4$, and we reject instances whose complexity $H(\theta)$ is larger than 500. The average Pareto set size was respectively 1.73, 1.64, 1.37, 1.28 for the number of arms $K \in \{10, 15, 25, 40\}$. Using $\delta = 0.01$, we observe a negligible empirical error.

Figure 9 shows that PSIPS has competitive empirical performance for higher dimension $d$. It significantly outperforms APE and slightly outperforms uniform sampling. The worsening of uniform sampling's performance stems from the smaller size of the Pareto set, which almost contains only one arm.
<figure><img src="image_11.png" /><figcaption>**Figure 10:** Average probability of misidentification of the Pareto set in the covid19 experiment. The results are averaged over 1000 runs.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11:** Theoretical impact of the correlation coefficient on the sample complexity for a 5-armed instance with $d = 2$.</figcaption></figure>

**Probability of Error** Both the recommendation rule and the sampling rule of PSIPS are $\delta$-independent and defined as any time $t$. Therefore, by disabling the stopping rule, it is possible to run PSIPS for $T = 5000$ time steps and report the empirical error $\mathbb{1}(\hat{S}_t \neq S^*)$ at any time $t \in [K, 5000]$ and averaged over 1000 runs. Uniform sampling enjoys the same anytime property and has the same recommendation rule. Therefore, we use it as a benchmark.

Figure 10 reveals that the sampling rule from PSIPS yields an empirical error which is close to the one achieved by uniform sampling. Even though PSIPS has only theoretical guarantees in the fixed-confidence setting, it performs empirically well in the anytime setting.

**Toward Understanding The Theoretical Impact of Correlation** As a step toward understanding the theoretical effect of correlation on the sample complexity, we picked an instance defined by $\mu_1 := (0.73, 1.20)^T$, $\mu_2 := (0.45, -0.63)^T$, $\mu_3 := (0.63, 1.28)^T$, $\mu_4 := (0.94, 2.31)^T$, $\mu_5 := (2.08, 1.48)^T$ with a correlation matrix $\Sigma_\rho$: unit diagonal and an off-diagonal term $\rho \in (-1, 1)$ (correlation coefficient as in the setup of Section 4). Using an approximate convex solver and the exponentiated gradient algorithm, we approximate the theoretical complexity $T^*(\mu, \Sigma_\rho)$ (as outlined in Lemma 1), which we report in Figure 11 w.r.t. $\rho$ and $\rho = 0$ (uncorrelated setting) is used as reference. In this instance, we note that negative correlation decreases the sample complexity (w.r.t. the same instance with uncorrected objectives).