Geometry Meets Incentives: Sample-Efficient Incentivized Exploration with Linear Contexts

Thank you for your helpful comments and feedback. Below we address your specific questions:

First, we will happily provide more discussion on the related work in the final version of the paper and apologize for the initial confusion.

Can you provide more motivation and any potential applications of the problem?

Many modern online platforms recommend products to customers, for example an online retailer recommends certain products for a consumer to buy. The customers can then choose to accept the recommendation and buy the product or ignore the recommendation and buy a different product. The goal of the platform is to learn the optimal product to recommend based on customer feedback. However, if customers ignore the platform’s recommendation based on a prior belief that another option is better, then the platform may be unable to learn efficiently. This motivates our studying of Bayesian Incentive Compatible algorithms, i.e. algorithms that make recommendations in a way that incentivizes the customers to always follow the platform’s recommendation. Other similar instances of this problem are in mobile health applications that make recommendations to patients that patients can choose to follow or ignore. As mentioned at the start of [Sellke 2023], the BIC property also implies that "exploration helps everyone"; this may be especially relevant for reliability in health applications as a codification of the Hippocratic Oath. We will add further discussion of this point in the paper.

Do you have any results of corresponding lower bounds? Lower bounds can make the paper sound and more complete.

There is of course a trivial $\lambda d$ lower bound on the number of samples needed to do $\lambda$-exploration, so we also know that $poly(d)$ is the right rate. Closing the gap between our upper bound and the linear in $d$ lower bound is a very interesting question for future work!

Furthermore, the examples discussed in Appendix B actually imply lower bounds on the number of samples necessary for exploration in terms of the parameters of Assumption 2. More formally, these examples imply the following three lower bounds. For proof sketches of these three lower bounds, please see the rebuttal to Reviewer V6tS (they are omitted here due to space constraints).

Lower bound $c_v:$ There exist instances (i.e. prior distributions) that require $\Omega(1/c_v)$ samples for $1$-spectral exploration.

Lower bound $c_d:$ There exist instances that require $\Omega(c_d)$ samples for $1$-spectral exploration.

Lower bound $\epsilon_d:$ There exist instances that require $\Omega(\log(1/\epsilon_d))$ samples for $1$-spectral exploration.

The three lower bounds described above give a more complete picture of how are algorithm is tight in a polynomial sense for these problem parameters. Furthermore, these results combined with Propositions 4-5 imply a stronger lower bound than just linear in $d$ for the distributions discussed in these propositions. We are very happy to include these more detailed discussions about problem-parameter lower bounds in the final version of the paper.

The results heavily rely on the agents being myopic. However, in practice, some agents may consider future revenue and misreport now. In other words, they may sacrifice current revenue to mislead the planner. Could you discuss the importance of this assumption and any potential difficulty extending the existing results?

When there is a steady arrival of new agents interacting with the platform, we generally expect rational agents to act myopically (without trying to mislead the planner). This is because even if a single agent does interact with the platform multiple times, that agent is still one of many agents that the platform observes, and therefore that agent will have a very small affect on the planner. Intuitively, this suggests it should be hard for that single agent to act strategically in a way that is very advantageous for their own future interactions.

We do also agree that the constant arrival of new agents is a strong assumption in works on Bayesian Incentive Compatible algorithms. Indeed there is a substantial literature on more realistic models of incentivized exploration, incorporating e.g. agent heterogeneity ([IMSW 19]) and partial data disclosure ([IMSW 20]) and other references from the introduction. Actually the fact that we consider linear bandits is another step in this direction, as it makes our work more realistic than the finite-armed bandit models originally used to study incentive compatibility. In any case, for theoretical understanding it is also valuable to understand simple models, as they can inform a wide range of more realistic models which each have their own quirks. For example, covariate diversity among agents is generally understood to aid exploration (see e.g. the paper Mostly Exploration-Free Algorithms for Contextual Bandits). We thus expect the lack of agent-specific contexts in our setting to only make the problem harder, i.e. our algorithm could be potentially finetuned and/or streamlined to take further advantage of such information. We are happy to add more discussion of these points, including potential caveats and further citations to related works that study various problem settings.

The authors don't provide any experiments. Even a simulation comparing their algorithms and existing algorithms will improve the paper significantly.

First, we want to emphasize that the main contribution and algorithm of this paper are theoretical and show that it is possible to explore every dimension in $poly(d)$ time. Furthermore, there are no existing algorithms for this problem. That being said, the algorithm itself is relatively simple to implement. Due to the interest from reviewers, we implemented our main algorithm (Algorithm 1) along with the two helper algorithms (Algorithms 2 and 3) in Python (which only takes ~100 lines of code). We then tested the algorithm on synthetic data with varying dimensions when the prior distribution is a $d$-dimensional Gaussian that is independent across dimensions and has mean $1$ in the first dimension and mean $0$ in all other dimensions. Note that our algorithm also can be applied to prior distributions with dependencies between coordinates, however we ran experiments for the independent case as this simplifies the code significantly (and also matches the independent priors in [Sellke-Slivkins 2023]). For this simple setting, we run our algorithm with $\epsilon_d = 0.1$, $c_d = 1$, $K= 1$, and $c_v = 1$.

Of course, the bound in Theorem 3 is a worst-case bound on the sample complexity, and therefore we would expect that on most instances of the problem the number of steps is significantly less than that upper bound. Furthermore, we note that the constants in our algorithm are certainly not optimal for this specific instance of the problem, yet these (relatively large) constants are what allows our theoretical results to hold for any prior distribution. While we are unable to load a graph in this rebuttal, we ran our algorithm for values of $d$ ranging from $d = 1$ to $d = 24$, and sampled $\ell^*$ from the prior distribution 100 times for each value of $d$ and then calculated the number of samples necessary to achieve $\lambda$-spectral exploration. The results show a clear quadratic dependence of the sample complexity on the dimension, and quartic scaling for the running time.

We do want to emphasize that for a specific instance of our problem (i.e. a specific prior), a user could tune the worst-case constants used in the algorithm to get faster $\lambda$-spectral exploration. However, the above table just uses the default constants of the algorithm with the goal of showing that

• the algorithm is implementable in relatively few lines of python code
• the algorithm is efficient to run for reasonably large values of $d$
• in practice, the number of steps grows polynomially in $d$ at a much better rate than the worst-case bound in our theoretical results (in this experiment we see a roughly quadratic relationship between the dimension and the number of steps). We are happy to add this result to the final version of our paper with more details about the implementation and with graphs showing the trends.

Thank you for your helpful comments and feedback. Below we address your specific questions:

The paper is only theoretical. It would have been very useful to get some experimental results of the proposed algorithm [...]

First, we want to emphasize that the main contribution and algorithm of this paper are theoretical and show that it is possible to explore every dimension in $poly(d)$ time. That being said, the algorithm itself is relatively simple to implement. Due to the interest from reviewers, we implemented our main algorithm (Algorithm 1) along with the two helper algorithms (Algorithms 2 and 3) in Python (which only takes ~100 lines of code). We then tested the algorithm on synthetic data with varying dimensions when the prior distribution is a $d$-dimensional Gaussian that is independent across dimensions and has mean $1$ in the first dimension and mean $0$ in all other dimensions. Note that our algorithm also can be applied to prior distributions with dependencies between coordinates, however we ran experiments for the independent case as this simplifies the code significantly (and also matches the independent priors in [Sellke-Slivkins 2023]). For this simple setting, we run our algorithm with $\epsilon_d = 0.1$, $c_d = 1$, $K= 1$, and $c_v = 1$.

Of course, the bound in Theorem 3 is a worst-case bound on the sample complexity, and therefore we would expect that on most instances of the problem the number of steps is significantly less than that upper bound. Furthermore, we note that the constants in our algorithm are certainly not optimal for this specific instance of the problem, yet these (relatively large) constants are what allows our theoretical results to hold for any prior distribution. While we are unable to load a graph in this rebuttal, we ran our algorithm for values of $d$ ranging from $d = 1$ to $d = 24$, and sampled $\ell^*$ from the prior distribution 100 times for each value of $d$ and then calculated the number of samples necessary to achieve $\lambda$-spectral exploration. The results show a clear quadratic dependence of the sample complexity on the dimension, and quartic scaling for the running time.

We do want to emphasize that for a specific instance of our problem (i.e. a specific prior), a user could tune the worst-case constants used in the algorithm to get faster $\lambda$-spectral exploration. However, the above table just uses the default constants of the algorithm with the goal of showing that

• the algorithm is implementable in relatively few lines of python code
• the algorithm is efficient to run for reasonably large values of $d$
• in practice, the number of steps grows polynomially in $d$ at a much better rate than the worst-case bound in our theoretical results (in this experiment we see a roughly quadratic relationship between the dimension and the number of steps). We are happy to add this result to the final version of our paper with more details about the implementation and with graphs showing the trends.
• Please refer to the rebuttal to Reviewer Rjs1 for a code snippet (not included here due to char limit).

Lack of pedagogy: the setting is not clearly / formally defined. For readers that do not know the incentivized exploration problem, it is very hard to capture what are the objectives of this kind of problem. [...] Beyond the lack of formal presentation of the problem, the paper also uses many notations that are never defined [...]

We sincerely apologize for this. Paralleling [Sellke 2023], we will add clearer presentation, motivation and other discussion for the problem formulation, as well as definitions of these notations, to the final version of the paper.

While I find the problem theoretically interesting, I am not fully convinced by its real interest for pratical applications, as users can interact multiple times, have their own interests, their own history, the knowledge and understanding of the planner's algorithm, etc)

We agree that the constant arrival of new agents is a strong assumption in works on Bayesian Incentive Compatible algorithms. Firstly, there is a substantial literature on more realistic models of incentivized exploration, incorporating e.g. agent heterogeneity ([IMSW 19]) and partial data disclosure ([IMSW 20]) and other references from the introduction. Actually the fact that we consider linear bandits makes our work more realistic than the original finite-armed bandit models used to study incentive compatibility. In any case, for theoretical understanding it is also valuable to understand simple models, as they can inform a wide range of more realistic models which each have their own quirks. For example, covariate diversity among agents is generally understood to aid exploration (see e.g. the paper Mostly Exploration-Free Algorithms for Contextual Bandits). Thus we expect the simple nature of our setting to only make the problem harder, meaning that our algorithms could be potentially finetuned and/or streamlined if the model were refined to be more realistic. We are happy to add more discussion of these points, including potential caveats and further citations to related works that study various problem setting.

Regarding applications, many modern online platforms recommend products to customers, for example an online retailer recommends certain products for a consumer to buy. The customers can then choose to accept the recommendation and buy the product or ignore the recommendation and buy a different product. The goal of the platform is to learn the optimal product to recommend based on customer feedback. However, if customers ignore the platform’s recommendation based on a prior belief that another option is better, then the platform may be unable to learn efficiently. This motivates our studying of Bayesian Incentive Compatible algorithms, i.e. algorithms that make recommendations in a way that incentivizes the customers to always follow the platform’s recommendation. Other similar instances of this problem are in mobile health applications that make recommendations to patients that patients can choose to follow or ignore. As mentioned at the start of [Sellke 2023], the BIC property also implies that "exploration helps everyone"; this may be especially relevant for reliability in health applications as a codification of the Hippocratic Oath. We will add further discussion of this point in the paper.

The central assumption of the algorithm is that actions live in a d-dimentionnal unit ball. Isn't it limitative for some applications ? And more importantly, in that setting what would justify a prior that would not be uniform over every vector from that ball ? This can be very limitative for the paper if no reasonnable pratical justification for this can be found since taking a fully non informative prior is a trivial setting, for which TS is BIC without any initial exploration phase. Or am I wrong on this point ?

First, we want to clarify: it is definitely not known that Thompson Sampling is BIC on its own for spherical action sets and uniform priors, even though all arms have the same mean reward under the prior. The only bandit setting where this kind of result has been proved is for discrete action sets with independent priors across the arms (in [Sellke-Slivkins 2023]). The proof there applies correlation inequalities to the increments of two Doob martingales, and strongly makes use of this independence. The work [HNSW22] follows up on [Sellke-Slivkins 2023] considers combinatorial semi-bandits (where actions are sets of multiple arms) but does not include a result of this type. [Sellke 2023] considers linear bandits and provides a counterexample showing the Thompson sampling may be BIC at time 1 but not time 2. It is plausible that a linear bandit extension of this result is true under some rotational symmetry conditions, but it would be a new, very interesting theorem.

Regarding the prior, note that in our setting the prior encapsulates the knowledge that a fresh agent enters the system with. Thus in modeling, the prior should take into account all publicly available knowledge. For example, the publication of a scientific paper or news article related to some of the actions should certainly influence the prior. While it might be challenging to describe exactly the correct prior to use in a particular real-life setting, we believe this illustrates the need for flexibility in handling rather general priors. Note that using a completely wrong prior means that the BIC guarantees we pursue are rendered invalid. This is in contrast with ordinary bandit problems, where the only downside to using a simpler prior is potentially slower learning in the early stages.

Thank you for your helpful comments and feedback. Below we address your specific questions:

Can authors discuss the tightness of dependence on d?

There is of course a trivial $\lambda d$ lower bound on the number of samples needed to do $\lambda$-exploration, so we also know that $poly(d)$ is the right rate. Closing the gap between our upper bound and the linear in $d$ lower bound is a very interesting question for future work!

Furthermore, the examples discussed in Appendix B actually imply lower bounds on the number of samples necessary for exploration in terms of the parameters of Assumption 2. More formally, these examples imply the following three lower bounds, and we include a brief proof sketch of each below.

Lower bound $c_v:$ There exist instances (i.e. prior distributions) that require $\Omega(1/c_v)$ samples for $1$-spectral exploration.

Proof sketch Consider the following example with $d=2$ and where the coordinates of $\ell^\ast$ are independent (and assume that $c_d < 1, \epsilon_d < 1, c_v < 1$):

\\ell^\\ast_1 = \\begin{cases} -c_d & \\text{w.p. $\epsilon_d$} \newline 1 & \\text{w.p. $1-\\epsilon_d$} \\end{cases},\\quad\\quad \\ell^\ast_2 = \\begin{cases} -c_v & \\text{w.p. $1/2$} \newline c_v & \\text{w.p. $1/2$} \\end{cases}

Note that if $E[\ell_1^* \mid \psi] = 1-O(\epsilon_d) \ge 0.5$ for sufficiently small $\epsilon_d$, then no optimal action will ever put weight more than $O(c_v)$ on the second coordinate because the magnitude of $\ell_2^*$ is bounded by $c_v$. This implies that we must need $1/c_v$ steps in order to guarantee $1$-spectral exploration.

$\blacksquare$

Lower bound $c_d:$ There exist instances that require $\Omega(c_d)$ samples for $1$-spectral exploration.

Proof Sketch

Consider the following example

$\ell^\ast_1 = \begin{cases} -c_d & \text{w.p. $\epsilon_d$}\newline 2c_d & \text{w.p. $2\epsilon_d$} \newline 1 & \text{w.p. $1-3\epsilon_d$} \end{cases},\quad\quad \ell^\ast_2 = \begin{cases} -c_v & \text{w.p. $1/2$}\newline c_v & \text{w.p. $1/2$} \end{cases}$

Once again, the first action must be $e_1$. Furthermore, in this case we need $O(poly(1/c_d))$ actions of $e_1$ in order to decrease the conditional expectation of $\ell^\ast_1$ to be $0$, as we need this many samples to be able to effectively distinguish between $\ell_1^* = -c_d$ and $\ell_1^* = 2c_d$. This implies that we require $O(poly(1/c_d))$ samples to explore the second dimension in this example. $\blacksquare$

Lower bound $\epsilon_d:$ There exist instances that require $\Omega(\log(1/\epsilon_d))$ samples for $1$-spectral exploration.

Proof sketch:

This proof is more subtle and requires inductive control of the information gain on the 2nd coordinate. We provide a preparatory lemma and corollary.

Lemma: let $L=±c$ be uniformly random for some $|c|≤1$. Suppose we receive noisy observations $r_i = s_i L + Z_i$ for a sequence $r_1,...$ that is adapted to the filtration $\mathcal F_t$ generated by $(s_1,r_1,s_2,r_2,...,s_t,r_t)$. (I.e. the signal strengths $s_i$ may depend on the past.) Let $T_t=\sum_{i=1}^t s_i^2$. Then the expected information gain on $L$ is at most $O(E[T_t])$.

Proof: this is a special case of observing Brownian motion with drift $L$ up to the random stopping time $T_t$. (Since observing $r_i = s_i L + Z_i$ is equivalent to observing Brownian motion with drift $L$ for time $s_i^2$, up to rescaling.) Let $Q_±(T_t)$ be the laws of Brownian motion with drifts $±c$ up to time $T_t$. By symmetry the expected information gain can be computed assuming $L=c$, and is bounded by

$ &E[KL(Q_+(T_t),[Q_+(T_t)+Q_-(T_t)]/2)] \\ &\stackrel{\text{Convexity of KL}}{\leq} E[KL(Q_+(T_t),Q_+(T_t))] + E[KL(Q_+(T_t),Q_{-}(T_t))] \\ &\leq E[T_t]. $

Here we used $E[KL(Q_+(T_t)|Q_-(T_t))]=c^2 E[T_t]≤E[T_t]$ by Girsanov's theorem. $\blacksquare$

Corollary: in the setting of the previous lemma, let $C_t = E[L|\mathcal F_t] = E[L|(s_1,...,s_t,r_1,...,r_t)].$ Then $E[C_t^2] \leq O(E[T_t])$.

Proof: This follows from the previous lemma: the expected relative entropy between the prior and posterior distributions of $L$ is precisely the information gain. In turn, the relative entropy is a strictly convex and even function of $C_t$ which is $\Omega(C_t^2)$. $\blacksquare$

We now return to the first example prior from the $c_v$ lower bound, and apply the Corollary to the observations of the 2nd coordinate. To make the application direct, whenever an action $(a_{1,t},a_{2,t})$ is played, we replace the noisy reward observation with separate observations for both $(a_{1,t},0)$ and $(0,a_{2,t})$ (with half the noise level, which only affects constant factors in the argument). This gives strictly more information since the two separate observations can be added to recover the original observation. We let $C_t,\mathcal F_t$ from above correspond to observations in the second coordinate.

At each time $t$, by Jensen's inequality on the convex function $f(x)=(1/3 - x)_+$, we see that for any signal $\psi$:

\leq O(E[f(E[\ell_1 | \psi])]) \leq O(E[f(\ell_1)]) \leq O(\epsilon_d).$$ On the main high-probability event that $E[\ell_1 | \psi]≥1/2$, the 2nd coordinate of $Exploit(\psi)$ has absolute value at most $O(|C_t|)$. Via the Lemma and Corollary above, it follows that $$\mathbb E[T_{t+1}] - \mathbb E[T_t] ≤ O(\mathbb E[T_t] + \epsilon_d).$$ Namely the case $\{E[\ell_1 | \psi]≤1/2\}$ contributes $O(\epsilon_d)$ while the remaining case contributes $O(\mathbb E[C_t^2])\leq O(\mathbb E[T_t])$. Since $T_0=0$, this implies by induction that $$\mathbb E[T_t] ≤ e^{O(t)} \epsilon_d.$$ Finally note that we need $T_t ≥ 1$ for $1$-spectral exploration. Indeed, $1$-spectral exploration requires that the actions $a_1,\dots,a_t$ satisfy

\sum_{i=1}^t \langle a_t^{\top}, M a_t\rangle \geq Tr(M)

for all positive semi-definite matrices $M$, and taking $M=\begin{pmatrix} 0 & 0 \newline 0 & 1\end{pmatrix}$ recovers the claim. In all this gives the desired $\log(1/\epsilon_d)$ lower bound. $\blacksquare$ The three lower bounds described above give a more complete picture of how are algorithm is tight in a polynomial sense for these problem parameters. Furthermore, these results combined with Propositions 4-5 imply a stronger lower bound than just linear in $d$ for the distributions discussed in these propositions. We are very happy to include these more detailed discussions about problem-parameter lower bounds in the final version of the paper. >*What's the computational cost and memory cost of the proposed algorithms?* The memory of the algorithms as-written scales linearly with the number of steps (which is of course polynomial in $d$), as the only information that needs to be stored are the action and reward from each step. The computational cost is also also polynomial in $d$, with the exception potentially of computing the function $f$ used in Line 4 of Algorithm 5. In our paper, we only prove the existance of such a function $f$. However, for many high-dimensional distributions (for example correlated Gaussians) this can be found in polynomial (in $d$) time using a simple linear program. In fact, for independent Gaussians (as in our experimental example), this $f$ can be directly computed in constant time without even needing a linear program. We mention that this situation is similar to [Sellke-Slivkins 2023] in the $K$-armed setting, where the main algorithm is also shown to be efficient only in special cases. Their examples involve Beta prior distributions (see Appendix G therein) rather than Gaussians. >*If the reward mean is not as a linear function but as a general function, can we get similar results under the same conditions? If not, what's could be the main challenges?* This is a great question, and indeed [Sellke 2023] briefly considers generalized linear bandits where a non-linear link function is applied to the inner product. The link function there is required to be strictly increasing with uniformly upper and lower bounded derivative. We do not see any conceptual barriers to extending our results to a similar class of functions, since the main property we rely on is the sensitivity of the optimal action with respect to the posterior mean reward. At the same time it is not straightforward and would require re-thinking many intermediate steps in the algorithm and analysis. (Both of these comments parallel the discussion at the end of page 3 of potentially generalizing to smooth non-spherical action sets.)

The paper did not include any experiments.

First, we want to emphasize that the main contribution and algorithm of this paper are theoretical and show that it is possible to explore every dimension in $poly(d)$ time. That being said, the algorithm itself is relatively simple to implement. Due to the interest from reviewers, we implemented our main algorithm (Algorithm 1) along with the two helper algorithms (Algorithms 2 and 3) in Python (which only takes ~100 lines of code). We then tested the algorithm on synthetic data with varying dimensions when the prior distribution is a $d$-dimensional Gaussian that is independent across dimensions and has mean $1$ in the first dimension and mean $0$ in all other dimensions. Note that our algorithm also can be applied to prior distributions with dependencies between coordinates, however we ran experiments for the independent case as this simplifies the code significantly (and also matches the independent priors in [Sellke-Slivkins 2023]). For this simple setting, we are able to run our algorithm with the constants of $\epsilon_d = 0.1$, $c_d = 1$, $K= 1$, and $c_v = 1$.

Of course, the bound in Theorem 3 is a worst-case bound on the sample complexity, and therefore we would expect that on most instances of the problem the number of steps is significantly less than that upper bound. Furthermore, we note that the constants in our algorithm are certainly not optimal for this specific instance of the problem, yet these (relatively large) constants are what allows our theoretical results to hold for any prior distribution. While we are unable to load a graph in this rebuttal, we ran our algorithm for values of $d$ ranging from $d = 1$ to $d = 24$, and sampled $\ell^*$ from the prior distribution 100 times for each value of $d$ and then calculated the number of samples necessary to achieve $\lambda$-spectral exploration. The results show a clear quadratic dependence of the sample complexity on the dimension, and quartic scaling for the running time.

We do want to emphasize that for a specific instance of our problem (i.e. a specific prior), a user could tune the worst-case constants used in the algorithm to get faster $\lambda$-spectral exploration. However, the above table just uses the default constants of the algorithm with the goal of showing that

• the algorithm is implementable in relatively few lines of python code
• the algorithm is efficient to run for reasonably large values of $d$
• in practice, the number of steps grows polynomially in $d$ at a much better rate than the worst-case bound in our theoretical results (in this experiment we see a roughly quadratic relationship between the dimension and the number of steps). We are of course happy to add this result to the final version of our paper with more details about the implementation and with graphs showing the trends.
• One reason why the sample complexity is so much better here than in our worst-case bound is that a factor of $d^4$ in our worst-case bound comes from applying Lemma 8 in our paper. This lemma is very much a worst-case bound for how much exploration we gain in each step in terms of the eigenvalues, and comes from subtleties of high-dimensional geometry. In most practical settings we would not expect to need as many steps as implied by that lemma.

We don't include all of the code due to space constraints, but the main algorithm corresponds to the algorithm in the paper pretty clearly as follows

def bic_exploration(d, lambda_val, eps_d, c_d, cLA2, cLA4, K, cv, ell_star, ell0_mean):
    kappa = max(1 / (lambda_val * cLA4), (4 * d * (K * np.sqrt(np.pi) + 1) * (1 + 1 / lambda_val)) / (cv ** 2 / (8 * np.pi)))
    v_list = [np.eye(d)[0]]  # Start with e1
    q_list = [[0]*int(kappa)]

    # Initial exploration with e1
    for i in range(int(kappa)):
        q_list[0][i] = np.dot(v_list[0], ell_star)+ np.random.normal()  # reward proxy

    t = int(kappa)
    j = 0
    M = lambda V: sum(np.outer(v, v) for v in V)
    while np.min(np.linalg.eigvals(M(v_list))) \le lambda_val:
        eigvals, eigvecs = np.linalg.eigh(M(v_list))
        idx = np.argsort(-eigvals)
        lambdas = eigvals[idx]
        w = eigvecs[:, idx].T
        ell_lambda = np.sum(lambdas >= lambda_val)

        a, t = initial_exploration(w, lambdas[:ell_lambda], v_list, q_list, t, cLA4, lambda_val, eps_d, c_d, K, ell_star, j+1)
        while np.linalg.norm(a - w[:ell_lambda].T @ (w[:ell_lambda] @ a)) \le np.sqrt(lambda_val):
            a, t = exponential_growth(a, w, lambdas[:ell_lambda], v_list, q_list, t, ell0_mean, cLA2, lambda_val, ell_star, j+1)
        v_list.append(a)
        q_list.append([0]*int(kappa))
        for i in range(int(kappa)):
            q_list[-1][i] = np.dot(a, ell_star)+ np.random.normal()
        t += int(kappa)
        j += 1
    return t

The paper focuses solely on the initial sampling exploration, and needs to be combined with Sellke 2023's result to have an end-to-end BIC algorithm. Thus, the contribution may seem a bit incremental.

We respectfully disagree that the contribution is incremental. The key technical gap left open by Sellke 2023 was that its Thompson‑sampling stage becomes BIC only after collecting an exponentially large amount of initial data in the dimension d. Consequently, there was no sample‑efficient end‑to‑end algorithm in high dimensions prior to our work, and it was not clear whether such guarantees should even be possible.

The paper did not explain the constants in assumption 2. Having more elaborate justification of what each constant means and how they should be chosen in practice would help with understanding these assumptions.

Propositions 4 and 5 both give examples of common distributions and associated constants for Assumption 2 that can be used in practice. We also include in Appendix B a more detailed discussion of the necessity for each of these constants and intuition for why we need these assumptions. While we unfortunately had to move this to the appendix for the submission, we are also happy to elaborate on this more in the body for the final version.

Furthermore, do the authors have any insight on whether this initial sample complexity is tight or can it still be improved in the future?

There is of course a trivial $\lambda d$ lower bound on the number of samples needed to do $\lambda$-exploration, so we also know that $poly(d)$ is the right rate. Closing the gap between our upper bound and the linear in $d$ lower bound is a very interesting question for future work!

Furthermore, the examples discussed in Appendix B actually imply lower bounds on the number of samples necessary for exploration in terms of the parameters of Assumption 2. More formally, these examples imply the following three lower bounds. For proof sketches of these three lower bounds, please see the rebuttal to Reviewer V6tS (they are omitted here due to space constraints).

Lower bound $c_v:$ There exist instances (i.e. prior distributions) that require $\Omega(1/c_v)$ samples for $1$-spectral exploration.

Lower bound $c_d:$ There exist instances that require $\Omega(c_d)$ samples for $1$-spectral exploration.

Lower bound $\epsilon_d:$ There exist instances that require $\Omega(\log(1/\epsilon_d))$ samples for $1$-spectral exploration.

The three lower bounds described above give a more complete picture of how are algorithm is tight in a polynomial sense for these problem parameters. Furthermore, these results combined with Propositions 4-5 imply a stronger lower bound than just linear in $d$ for the distributions discussed in these propositions. We are very happy to include these more detailed discussions about problem-parameter lower bounds in the final version of the paper.