Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
102,072 characters · 18 sections · 67 citation commands
Semidiscrete optimal transport with unknown costs
{.5\baselineskip}
In the semidiscrete optimal transport problem PeCu19, we are given a continuous random vector $X$ with known density $f$, and a discrete random variable $Y$ with known probability mass function $p$. Our task is to choose a joint distribution for $\left(X,Y\right)$ to minimize an expected cost. Formally, we write
subject to
where $\mathcal{X}\subseteq \mathds{R}^d$ is the support of $X$ and $\left\{1,...,K\right\}$ is the support of $Y$. The marginal distributions are fixed by the problem inputs $p_k > 0$ and $f$. The objective function can also be written as $\mathbb{E}\left[c\left(X,Y\right)\right]$ under the assumption that $P\left(X\in dx,Y=k\right) = h\left(x,k\right)dx$. In other words, $h$ is the mixed joint likelihood of $\left(X,Y\right)$.
The name “optimal transport” has the same origin as the well-known “transportation problem” in linear programming FoFu56, which can be interpreted as a special case of ((ref)) where both $X$ and $Y$ have finite support. In that case, the joint distribution of $\left(X,Y\right)$ is described by a finite set of probabilities $h_{jk}$, whose cost coefficients are given by $c_{jk}$ for all possible values $\left(j,k\right)$ of $\left(X,Y\right)$. The optimization problem then reduces to a linear program, whose solution can be interpreted as a minimum-cost plan for matching supply from one set of locations (or facilities) indexed by $j$ with demand at a different set of locations indexed by $k$. Such plans have numerous applications in logistics. In semidiscrete optimal transport, $Y$ remains discrete but $X$ becomes continuous, which makes the problem more difficult and generally not solvable using LP methods. The semidiscrete version also has applications in logistics, specifically in geographical partitioning problems CaCaDe16,HaSc20 where demand can arise anywhere on the map, and the goal is to design service zones that assign spatial regions to supply facilities.
The optimal transport literature universally assumes that the cost function $c$ is known: in fact, most papers only work with some particular form for $c$, most commonly a Euclidean distance or $p$-norm Sa15. Our paper is the first to focus on a situation where $c$ is unknown, but can be estimated in a sequential manner from data. Our motivation for this setting is a class of applications where $k\in\left\{1,...,K\right\}$ represents a “choice” or “alternative” available to a decision-maker. One observes a sample of $X$ from the density $f$, then selects an alternative $k$ and incurs a cost $c\left(X,k\right)$. The goal is to assign $X$ values to alternatives in a way that minimizes expected cost while obeying the constraints ((ref)) on how often each alternative can be chosen.
Consider the following application to online advertising. An advertising platform has $K$ clients. A guaranteed contract Bharadwaj10 ensures that client $k$ receives a certain number of impressions through the platform, which can be expressed as a proportion $p_k$, agreed upon ahead of time, of the platform's total user base. Users have heterogeneous attributes which influence how they respond to advertisements. The platform observes the attributes $X$ of each randomly arriving user and selects a client, whose advertisement is then shown to the user in real time. Only one client can be selected per user, and the platform has the freedom to make individual assignments, provided that the long-term proportion of users assigned to client $k$ is equal to $p_k$. An assignment of a user with attributes $x$ to the $k$th client can be evaluated using a performance metric such as the clickthrough rate, which can be represented using, e.g., the linear model $c\left(x,k\right) = \beta^\top_k x$ LiChLaSc10 or the logistic model $c\left(x,k\right) = \frac{1}{1+e^{-\beta^\top_k x}}$ ChCa10. We can maximize clickthrough rates by minimizing the negative of these functions in ((ref)). The coefficients $\beta_k \in \mathbb{R}^d$ are unknown (note that each client has a different set of coefficients), reflecting the fact that the platform cannot predict clickthrough rates with perfect accuracy. However, they can be inferred over time by observing how users behave. The goal is to learn $\beta_k$ efficiently while optimizing the expected performance.
More generally, our paper considers parametric cost functions of the form $c\left(x,k\right) = \chi\left(\beta^\top_k x\right)$, where $\chi$ can be any link function associated with a generalized linear statistical model McNe89. We do not know $\beta_k$, but when an observed $x$ is assigned to the $k$th alternative, we collect an observation of the form $c\left(x,k\right)+\varepsilon$, where $\varepsilon$ is an additive noise. We then use an appropriate statistical procedure (e.g., linear or logistic regression) to estimate $\beta_k$ from the observations for that $k$ value. Only one $k$ can be observed at a time, creating a tradeoff: one more observation of the $k$th cost function means one less for the others.
This is essentially the same tradeoff that arises in classical optimal learning problems, such as ranking and selection ChChLePu15 and multi-armed bandits GiGlWe11, where one similarly observes noisy samples of the values of various “alternatives” with the goal of efficiently identifying the best. Our problem belongs to the class of contextual learning, where the value of an alternative additionally depends on a random vector $X$; see, e.g., HaRyDe16, AlChZo21 and ShHoZh21 from the ranking and selection literature, or AbPaSz11 and BaBaKh21 from the bandit community. However, the presence of the target constraint ((ref)) in our problem is a significant departure from existing models. In optimal learning, there is a clear notion of “correct selection,” when the alternative with the smallest true value also has the smallest estimated value, and of “regret” when it does not. The proportional assignments $p_k$ in those problems vary between algorithms GlJu04, and many papers aim to learn proportions that maximize the probability of correct selection Ru20,ChRy22. However, under ((ref)), we sometimes have to make choices that do not achieve the lowest cost, even if the costs are known. Thus, the standard definitions of correct selection or regret do not apply in this setting.
In this paper, we bridge optimal transport theory and optimal learning to formalize a notion of correct selection, and we create a sampling algorithm that makes correct selections at an optimal rate. As a byproduct, we also obtain new results for the more standard setting of semidiscrete optimal transport with known costs. Although the known-cost setting is not the main focus of this research, we use it as a starting point to develop the analytical technique that we eventually use to handle unknown costs. We first explain the essential elements of this technique in Section (ref) below, and then discuss the main results in Section (ref).
There is a vast literature on optimal transport with known $c$. Vi21 gives an introduction to the deep theory of this field. Computational approaches include methods based on discretization, such as the Sinkhorn Cu13 and Greenkhorn LiHoJo22 algorithms or the proximal splitting procedure of PaPeOu14, as well as methods based on PDEs BeFrOb14 and difference of convex programming BoPe24.
The specific setting of semidiscrete optimal transport has distinct structure that greatly aids computation. Problem ((ref))-((ref)) admits an optimal policy $\pi^*$, a decision rule that assigns an observed $X$ to some $\pi^*\left(X\right)\in\left\{1,...,K\right\}$, in the presence of probabilistic targets. This characterization gives us a notion of correct selection that we will later use in the unknown-cost setting: given $X$, we want to choose the same $k$ as the optimal policy. The crucial distinguishing feature of the semidiscrete problem is that $\pi^*$ is completely characterized by a finite-dimensional vector of parameters, namely, the shadow prices of ((ref)) that optimally solve the Kantorovich dual of ((ref))-((ref)). As long as we know $c$ and are able to observe samples from $f$, we can learn these parameters using stochastic approximation (SA) techniques. SA is very simple to implement, performs well empirically GeCuPeBa16, and is guaranteed to converge to the optimal parameters due to the concavity of the Kantorovich dual problem. An important detail is that SA does not require discretization, unlike many general-purpose optimal transport methods.
We build on this approach and use SA under both known and unknown costs. Our main focus is on theoretical performance guarantees. In the literature on SA, the strongest possible guarantee BaMo11 is that the sequence of iterates converges in $L^2$ at the canonical rate $\mathcal{O}\left(\frac{1}{n}\right)$, where $n$ is the number of samples drawn from $f$. In general, however, this rate is only achievable with Polyak averaging PoJu92 included into the SA procedure and with the assumption of strong convexity (concavity) on the SA objective function. Even in later work (e.g., BaMo13), which has yielded sharper rates, strong convexity is still needed to bound the $L^2$ error of the SA iterates. The main technical issue in our case is that the objective of the Kantorovich dual is concave, but not strongly concave. The optimal transport literature is aware of this difficulty, and has developed a workaround (see, e.g., GeCuPeBa16 or TaShKu22) that replaces the objective by a smoothed approximation. This approach recovers the canonical rate at the expense of one of the most attractive features of the algorithm, namely, the ability to solve the original problem exactly.
In marked contrast, we do not modify the problem. Instead, we adopt a perspective from the recent bandit literature (e.g., GoZe13, BaBa20, and BaBaKh21), which has found that convergence rates improve when the data-generating distribution exhibits a sufficient degree of random variation (“richness”). With comparable assumptions on the distribution of $X$, the objective is still not strongly concave, but it has a kind of “local” strong concavity. Our paper is the first to identify this property and show that it enables recovery of the canonical rate without smoothing or averaging. This part of our analysis holds generally, without assuming any particular structure on $c$, as long as it is known. However, our analysis of unknown costs also draws upon this result.
When $c$ is unknown, we still use SA to learn the dual variables, but the situation becomes much more complicated because the gradient estimates are no longer unbiased (they now have to be computed using estimated cost coefficients $\beta_k$). The problem couples optimal learning together with SA in such a way that both aspects create difficulties for each other: statistical error in the estimated costs creates bias in the gradient estimator, but the gradient estimator drives the decision rule and thus determines which alternatives are sampled.
We solve this problem using a simple, computationally efficient online algorithm with a “semi-myopic” structure. Most of the time, the algorithm approximates the optimal policy using plug-in estimators of both the costs and the dual variables, but occasionally it is forced to sample each individual cost function. These periods of forced exploration become less frequent as time goes on, to reduce their impact on performance. Semi-myopic structure has a long history in optimal learning, dating back to the work of BrRu12 on dynamic pricing, and followed by extensions to ranking and selection GaKa16 and contextual bandits GoZe13. Virtually all of these papers require assumptions on the richness of the data-generating distribution in order to obtain optimal rates. Semi-myopic algorithms also vary in how often they need to conduct forced exploration. The gap between exploration periods grows linearly in BrRu12; quadratically in KeZe14, GaKa16, and BaKe21; and exponentially in GoZe13 and BaBa20. The gap for our procedure, in the more complex setting of semidiscrete optimal transport, grows sub-exponentially.
Our paper is the first to couple this type of exploration with stochastic approximation. As in SA, we update the iterates using (biased) gradient estimates; as in optimal learning, information is collected sequentially in a semi-myopic fashion. Building on our results for the known-cost setting, we prove that SA still achieves the canonical rate under this setup. As a corollary, the expected number of incorrect decisions made by the semi-myopic policy (the most relevant analog of “regret” in our setting) grows at a rate $\mathcal{O}\left(\sqrt{n}\right)$, the best possible. This performance is achieved in an online manner, without requiring the sample size to be known in advance. Numerical experiments show that, in some situations, cost uncertainty produces only minor performance loss relative to an “ideal” SA procedure with unbiased gradients.
In summary, our work bridges and contributes to three distinct streams of literature. We contribute to optimal transport by being the first to formulate and study the semidiscrete problem with unknown costs. This problem couples optimal learning and stochastic approximation in a novel and challenging way. We contribute to optimal learning by studying a novel contextual problem with probabilistic targets, and by proving that semi-myopic exploration can learn the optimal policy efficiently. We contribute to SA by showing that optimal convergence rates can be achievable without strong concavity, as long as the data-generating distribution is sufficiently rich.
Section (ref) presents preliminaries on the semidiscrete optimal transport problem, most importantly the structure of the optimal policy. We then state the algorithm and key theoretical results: Section (ref) presents stochastic approximation (Section (ref)) and semi-myopic learning (Section (ref)) together with our main assumptions and rate results. Then, Section (ref) presents the main steps of our technical analysis for SA (local strong concavity under known costs), while Section (ref) does the same for semi-myopic learning of unknown costs. Section (ref) presents numerical examples and insights, and Section (ref) concludes.
We begin by reviewing the structure of the semidiscrete optimal transport problem. The literature has characterized the structure of the optimal solution to ((ref))-((ref)) when $c$ is known; these results can be found in, e.g., CaCaDe16 or GeCuPeBa16, but we summarize them here to provide background, sketching the key steps in the derivation. In brief, the problem ((ref))-((ref)) is solved by a policy, or decision rule for assigning any observed $X$ to a $k$ value. In other words, $Y$ should be set equal to a certain deterministic function of $X$. It is necessary to understand this function before we can formally define a performance metric for the problem that we will eventually tackle.
For notational convenience, we write $c_k\left(x\right) = c\left(x,k\right)$, to visually distinguish between costs $c_k$ corresponding to different “alternatives,” from which we will eventually make a selection. We do not impose any parametric structure on $c_k$. By Kantorovich duality (Theorem 1.3 in Vi21), the problem ((ref))-((ref)) has the same optimal value as the functional optimization problem
subject to the constraint
In ((ref)), the marginal density of $X$ is $f$, and the marginal pmf of $Y$ is $p$, as in the primal problem. To handle ((ref)), we take
Letting $g = \left(\psi\left(1\right),...,\psi\left(K\right)\right)^\top$, we can rewrite ((ref)) as
The minimum of a finite number of linear functions is concave, and expectations preserve concavity, so the objective in ((ref)) is concave in $g$. An important distinguishing characteristic of semidiscrete optimal transport is that the dual is both concave and finite-dimensional, and thus it can be solved without the use of discretization or nonconvex programming (see, e.g., BoPe24).
In fact, there are no constraints on $g$, so we can simply take the gradient and set it equal to zero. Under some mild conditions, the gradient can be interchanged with the expectation, an instance of “infinitesimal perturbation analysis,” a well-known technique in the simulation community Fu06,Ki06. We show this formally for completeness, as GeCuPeBa16 did not do so. The proof is deferred to the Appendix.
Under the conditions of Proposition (ref), we combine ((ref)) with ((ref)) and find that any optimal $g^*$ solves the system of equations
In words, $g^*$ is a vector of bonuses and penalties (since its elements can be either positive or negative) which are subtracted from the costs $c_j\left(X\right)$ to ensure that each element is the smallest with precisely the target probability.
These same bonuses and penalties can be used to construct an optimal solution for the primal problem ((ref))-((ref)). By weak duality, we have
where $h$ is the joint likelihood of $\left(X,Y\right)$ satisfying ((ref))-((ref)), as before. Now define $Y^* = \arg\min_j c_j\left(X\right) - g^*_j$. By ((ref)), we have $P\left(Y^* = k\right) = p_k$, so the joint likelihood of $\left(X,Y^*\right)$ satisfies ((ref))-((ref)). We then write
where the last line follows by ((ref)). Then, letting
it follows that the joint distribution of $\left(X,\pi^*\left(X\right)\right)$ optimally solves ((ref))-((ref)).
The policy $\pi^*$ is a natural benchmark against which other policies may be compared. Suppose that $\left\{X^n\right\}^{\infty}_{n=1}$ is a sequence of i.i.d. samples from $f$, and $\left\{\pi^n\right\}^{\infty}_{n=1}$ is a sequence of random variables taking values in $\left\{1,...,K\right\}$. We say that a correct selection is made at time $n$ if $\pi^n = \pi^*\left(X^n\right)$. We may use the expected total number $\mathbb{E}\left(\sum^n_{m=1} 1_{\left\{\pi^m\neq\pi^*\left(X^m\right)\right\}}\right)$ of incorrect selections to evaluate the performance of the sequence up to time $n$. The probability of incorrect selection (PICS), written as $P\left(\pi^n\neq \pi^*\left(X^n\right)\right)$, plays an important role in the analysis of this performance metric.
Many papers on optimal learning present their results in terms of the regret, which in our context would be the difference in cost achieved by $\pi^n$ vs. $\pi^*\left(X^n\right)$. However, because of ((ref)), regret relative to $\pi^*$ can be negative, unlike in traditional bandit problems. We therefore focus on the expected number of incorrect selections, which has a very clear, straightforward interpretation.
Our overarching goal is to learn $\pi^*$ iteratively based on a sequence $\left\{X^n\right\}^{\infty}_{n=1}$ of i.i.d. samples from the data-generating density $f$. Ultimately, we wish to do this in a setting where the cost functions $c_k$ are unknown and have to be estimated. Our algorithmic approach combines two elements: stochastic approximation (for learning the optimal dual variables $g^*$) and semi-myopic exploration (for learning the costs).
Integrating these elements is quite challenging. Therefore, we begin by examining stochastic approximation separately, in a setting where the costs are known. Section (ref) states the assumptions and main results needed to establish convergence rates for SA in semidiscrete optimal transport. These results, particularly the local strong concavity property, will continue to play a vital role when we switch to the unknown-cost setting. Section (ref) presents our algorithm and main results for that case.
It was observed by GeCuPeBa16 that, when the cost functions $c_k$ are known, ((ref)) can be approached as a stochastic root-finding problem PaKi11. We can straightforwardly apply the stochastic approximation (SA) algorithm
where $\left\{\alpha_{n}\right\}^{\infty}_{n=0}$ is a stepsize sequence satisfying the usual conditions $\sum_n \alpha_n = \infty$, $\sum_n \alpha^2_n < \infty$. The gradient is bounded, so the convergence $g^n\rightarrow g^*$ easily follows from classical SA theory KuYi03. The only input to ((ref)) is the stepsize $\alpha_n$. The algorithm requires knowledge of the targets $p_k$ and the cost functions $c_j$, but does not require us to know the density $f$ from which each $X^{n+1}$ is sampled. This is an attractive feature of SA, as data-generating distributions are extremely difficult to estimate. If the goal is to maximize expected costs, as in the motivating example from Section (ref), we can use ((ref)) with a trivial sign change.
Although our main focus is on unknown costs, the convergence rate of ((ref)) is an important prerequisite for that analysis. As mentioned in Section (ref), the main technical issue is the fact that the objective function $g \mapsto \mathbb{E}\left(F\left(g,X\right)\right)$ is concave, but not strongly concave. We will explain in Section (ref) how this issue is resolved; for the moment, let us state the assumptions and results. It is important to note that these assumptions do not require any specific parametric structure on $c_k$. Such structure will only be imposed in Section (ref) when we consider unknown costs.
From ((ref)), we see that adding a constant to every component of $g^*$ does not change the optimal decision. Therefore, we only need to learn the differences $g^*_k-g^*_1$ for $k > 1$, so we may assume without loss of generality that $g^*_1 = 0$. Our main assumption concerns the “richness” of the data-generating distribution. It is based on the recent optimal learning literature, which has found, in various settings, that desirable convergence rates require the observed values of $X$ to exhibit a sufficient degree of random variation. This is made more precise as follows.
Assumption (ref)(i) is a so-called “margin condition,” requiring sufficient random variation near the decision boundary (where multiple components of $Z$ are close to zero). Conditions of this type in the existing literature include Assumption A2 of GoZe13, Assumption 2 of BaBa20, and Assumption 2 of BaBaKh21, all of which also assume boundedness as in Assumption (ref)(ii). The main difference between these papers and our work is that they impose margin conditions for pairs of choices (in our setting, these would be the pairwise differences $Z_j - Z_k$), while we make the slightly stronger assumption that the margin condition holds on the entire multivariate density of $Z$. The distinctions between assumptions are due to the fact that, in the above-cited papers, the penalties $g^*_j$ are not present at all, and the ultimate goal is only to learn the cost functions, which is done in a “marginal” way, by observing one $k$ at a time. In our setting, even when the costs are known, we still have the problem of learning $g^*_j$, which is done “jointly” by updating our estimates of all the penalties in each iteration of ((ref)). It therefore becomes necessary to consider the decision boundary in $K$ dimensions.
When Assumption (ref) is imposed, one then obtains certain useful technical properties for the pairwise differences. These properties will be useful for our analysis later. We state them below; the proof is deferred to the Appendix. Note that the second property considers the pairwise differences jointly by minimizing over them, which is one reason why Assumption (ref) holds on the multivariate density of $Z$.
Under known costs, our algorithm for learning $g^*$ is given by ((ref)) using a simple stepsize $\alpha_{n} = \frac{\alpha}{n+1}$ with $\alpha>0$ being a pre-specified constant. Our main result is a rate bound on the convergence in $L^2$ of the sequence $\left\{\delta^n\right\}^{\infty}_{n=0}$, with $\delta^n = g^n - g^*$.
Thus, under Assumption (ref), classical SA achieves the canonical $\mathcal{O}\left(\frac{1}{n}\right)$ rate. It is interesting to note that this does not happen in general: one usually resorts to Polyak averaging PoJu92 to recover the optimal convergence rate, which is indeed done in BaMo11. However, in the specific setting of semidiscrete optimal transport with known costs, neither averaging nor strong concavity are necessary as long as the data-generating process is rich in the sense of Assumption (ref). This result is completely new to our paper, and is valuable for preserving one of the most attractive characteristics of the SA approach to semidiscrete optimal transport, namely its ability to learn (in the limit) the exact optimal policy without resorting to approximation or discretization of the problem.
The result of Theorem (ref) depends on two parameters, the threshold $\bar{\alpha}$ used by the stepsize, and the constant $D$ in the bound. We briefly clarify their dependence on the constants $\kappa_1,...,\kappa_4$ obtained from Assumption (ref). The threshold $\bar{\alpha}$ is the main source of this dependence, because $D$ is essentially an increasing function of $\alpha > \bar{\alpha}$. Our proofs indicate that $\bar{\alpha}$ is decreasing in $\kappa_1$ and $\kappa_2$ (that is, $\bar{\alpha}$ is larger when these constants are smaller) and increasing in $\kappa_3$ and $\kappa_4$. It is not possible to precisely quantify each relationship because, for example, $\bar{\alpha}$ depends on the support of $Z$ (which is controlled by $\kappa_4$) indirectly, through the support of the pairwise differences $Z_{j,k}$. However, roughly speaking, expanding the support of $X$ generally increases the support of $Z_{j,k}$ and thus the threshold.
Now, let us turn to the case where the cost functions are unknown. In each iteration, we still observe $X^{n+1}$, but we do not know $c_k\left(X^{n+1}\right)$. We only have access to approximate cost functions $\hat{c}^n_k$, which now have to be used to update $g^n$. We can, however, select one particular $k$ and collect a noisy observation of the form $W^{n+1}_k = c_k\left(X^{n+1}\right) + \varepsilon^{n+1}_k$, where $\left\{\varepsilon^n_k\right\}$ are i.i.d. noise terms with mean zero, assumed to be independent of $\left\{X^n\right\}$. Additionally, we can only observe one function at a time: if we choose $k$ at iteration $n$, we cannot see $W^{n+1}_j$ for $j\neq k$.
The techniques developed in our paper apply to cost functions of the form $c_k\left(X\right) = \chi\left(\beta^\top_k X\right)$, where $\chi$ is a link function associated with some generalized linear model McNe89 used in statistics. The simplest link function is the identity, in which case the cost functions are linear. Other possible $\chi$ include the logistic link $\chi\left(t\right) = \frac{1}{1+e^{-t}}$ and the Poisson link $\chi\left(t\right) = e^t$. In all of these cases, $\chi$ is known; the choice of link function depends on the application at hand, with, e.g., the logistic link modeling bounded costs, the Poisson link modeling integer-valued costs, and so on. However, in all cases, the parameters $\beta_k$ are unknown, and the approximate cost functions are given by $\hat{c}^n_k\left(X\right) = \chi\left(\left(\hat{\beta}^n_k\right)^\top X\right)$, where $\hat{\beta}^n_k$ is a maximum likelihood estimator of $\beta_k$ computed using all past observations (up to the $n$th time stage) of the $k$th cost function.
The cost function to be observed, formally denoted as $\pi^{n+1}$, is selected in an online manner: it depends on the available estimators $\hat{\beta}^n_k$ and $g^n$, but also on the new data $X^{n+1}$. Thus, $\pi^{n+1}$ affects the approximate costs in the next iteration, and through these approximations, the next gradient update as well. This algorithmic framework, laid out in Algorithm (ref), can be viewed as a coupling of the SA procedure ((ref)) with optimal learning (e.g., bandit problems). As in optimal learning, we can only sample one “alternative” (cost function) at a time, but the policy used to do this is linked to the bias of the SA update.
As discussed in Section (ref), our analysis uses a particular construction of $\pi^{n+1}$ known as a semi-myopic policy. This approach approximates ((ref)) using plug-in estimates of $g^*$ and $\beta_k$. Most decisions are made based on this approximation, but we occasionally force ourselves to sample different $k$ to prevent the estimates from stalling. Forced exploration occurs less frequently as time goes on.
Formally, we first define $\hat{\pi}^{n+1} = \arg\min_j \hat{c}^n_j\left(X^{n+1}\right) - g^n_j$ to be the plug-in estimate of the optimal decision, based on estimates $\hat{\beta}^n$ and $g^n$. We also define, for each $k = 1,...,K$, the set
of time periods where we will force ourselves to explore $k$. Any $a>0$ works for the theory; in practice, this value is tunable (for instance, $a \geq 5$ produces fairly infrequent forced exploration). Then, the $\left(n+1\right)$st decision is made according to the policy
The estimate $g^n_j$ is updated using ((ref)), and the regression coefficients are computed using a suitable generalized linear model. Thus, in order to implement the algorithm, we require only an exploration parameter $a > 0$ and a stepsize sequence; as in Section (ref), we will use $\alpha_n = \frac{\alpha}{n+1}$, so only a single parameter $\alpha$ needs to be specified.
We now proceed to the assumptions and main results. Here and in Section (ref), we present our main results for the special case of linear costs $c_k\left(X\right) = \beta^\top_k X$, for which we use the ridge-like estimator
where $\rho^n = 1 + \left(\log n\right)^3$. A rigorous extension to more general $\chi$ is provided in the Appendix. The results and analysis for that case are very similar, with the estimator $\hat{\beta}^{n+1}_k$ being a modification of maximum likelihood, just as ((ref)) is a modification of ordinary least squares. The reason to focus on linear costs here is for ease of presentation: the analysis becomes somewhat more streamlined, avoiding technical nuisances and allowing the reader to understand the main concepts. We emphasize, however, that the results hold for general $\chi$ and the extension does not add much difficulty to what we present here.
As before, we let $\delta^n = g^n - g^*$. We also define $\Delta^n_k = \hat{\beta}^n_k - \beta_k$. Recall the distinction between $\hat{\pi}^n$, which is the plug-in estimate of the optimal decision, and $\pi^n$, which is the semi-myopic policy combining that estimate with forced exploration as defined in ((ref)).
We require two additional assumptions on the data-generating distribution. The first is a standard assumption on the residual noise terms, used in virtually all proofs in contemporary learning theory. The second is very similar to Assumption A3 in GoZe13, and is weaker than Assumption 3 of BaBaKh21.
The bulk of the analysis focuses on bounding the convergence rate of the estimation error for both $g^n$ and $\hat{\beta}^n_k$. These results are stated in Theorem (ref). Then, Theorem (ref) obtains a rate for the accuracy of the plug-in estimate $\hat{\pi}^{n}$, in the sense of the PICS metric defined in Section (ref).
As in the case of Theorem (ref), we may briefly discuss the dependence of the threshold $\bar{\alpha}$ on the constants $\kappa_1,...,\kappa_6$. As before, the threshold is decreasing in $\kappa_1,\kappa_2$ and increasing in $\kappa_3$. The dependence on $\kappa_4$ is more complex because the support of $X$ indirectly influences various constants in the proof, but it is possible to say that the threshold becomes increasing in $\kappa_4$ for large enough values of this constant. The new constants $\kappa_5,\kappa_6$ that are specific to the unknown-cost setting do not influence $\bar{\alpha}$, only the rate constants $C,C'$. Unfortunately, the dependence of $C,C'$ on $\kappa_1,...,\kappa_6$ is non-monotonic and very difficult to characterize.
The square-root rate in ((ref)) is due to the fact that some $X$ values can make it arbitrarily difficult to distinguish between choices. Faster rates (e.g., exponential rates, as in GaDuCh19) can only be obtained if one assumes some minimal amount of separation between objective values, for example by making the support of $X$ finite. If this is not the case, it will be impossible to improve on ((ref)).
While we focus on PICS rather than regret in this work, our result is comparable to an optimal regret bound. For example, GoZe13 obtains a single-period regret rate of $\mathcal{O}\left(\frac{1}{n}\right)$ (a known universal lower bound) from a bound that multiplies together two terms, one involving the estimation error $\|\Delta^n_j\|_2$ of the regression coefficients, and one that looks similar to PICS. Individually, both terms converge at a rate of $\mathcal{O}\left(\frac{1}{\sqrt{n}}\right)$, as they do in our paper. One can view that problem as a special case of ours with known $g^*\equiv 0$, which implies that it is not possible to guarantee faster rates than these.
A straightforward consequence of Theorem (ref) is that the expected number of incorrect selections made by $\hat{\pi}^{n+1}$ (essentially adding up ((ref)) over $n$) grows according to $\mathcal{O}\left(\sqrt{n}\right)$. This rate does not change if we count incorrect selections made by the semi-myopic policy (incorporating the forced exploration periods) rather than simply the plug-in estimate, because forced exploration occurs only in $\mathcal{O}\left(\left(\log n\right)^9\right)$ time periods.
Earlier, in Section (ref), we separated the presentation of the main results according to the core components of our algorithm: first, stochastic approximation under known costs (Theorem (ref)), and then semi-myopic learning under unknown costs (Theorems (ref)-(ref)). We make the same distinction when presenting the proofs. In this section, we confine ourselves to the setting of Section (ref) and prove Theorem (ref). The unknown-cost case will be treated later in Section (ref).
As was discussed previously, the main technical challenge in establishing Theorem (ref) is the lack of strong concavity in the objective function. A crucial first step is the derivation (in Section (ref)) of a kind of “local strong concavity” property that is weaker than strong concavity, but still sufficient to prove (in Section (ref)) the main rate result. We will later use this property in our analysis of unknown costs, but it holds more generally and does not require $c_k$ to have any parametric structure. For that reason, the results in this section are of some stand-alone interest, because we are the first to show that semidiscrete optimal transport can be solved at an optimal rate without the need for discretization or smoothing, as long as the data-generating distribution is sufficiently rich.
We establish two technical results concerning the smoothness of the objective $\mathbb{E}\left(F\left(g,X\right)\right)$. The proof techniques will be used again later on, in Section (ref), where these properties will be generalized to the setting of unknown cost functions. We will assume Assumption (ref) throughout, denoting by $\kappa'_1,...,\kappa'_4$ the constants obtained from Lemma (ref).
\noindentProof: To avoid notational clutter, we write $g^*$ as simply $g$ in this proof, since no other $g$ will be considered. By definition of $F$, we can write
where
Clearly, $\mathbb{E}\left(T_1\right) = \delta^\top p$. We will work toward a lower bound on $\mathbb{E}\left(T_2\right)$, which will require some algebra on ((ref)). Because the indices minimizing $F\left(g+\delta,X\right)$ and $F\left(g,X\right)$ are unique by Assumption (ref), we can write
The mapping $z\mapsto \left[\left(c_k\left(X\right)-g_k-\delta_k\right)-z\right]\cdot 1_{\left\{c_k\left(X\right)-g_k = F\left(g,X\right),c_k\left(X\right) - g_k-\delta_k > z\right\}}$ is decreasing in $z$. Because $F\left(g+\delta,X\right)\leq c_j\left(X\right)-g_j-\delta_j$ for any $j$ by definition, we then have
Consequently,
We now seek to bound the expectation of each term in ((ref)). For convenience, let $Z_{j,k} = \left(c_j\left(X\right)-g_j\right)-\left(c_k\left(X\right)-g_k\right)$ for any $j,k$. Similarly, let $\bar{\delta}_{j,k} = \delta_j-\delta_k$. We write
Taking expectations, we obtain
where ((ref)) follows by the tower property, and ((ref)) follows because $\bar{\delta}_{j,k}-z\geq 0$ for $z \in \left[0,\bar{\delta}_{j,k}\right]$.
Using Lemma (ref)(ii), and taking $0 < \kappa<\kappa'_1$, we write
and, since $P\left(Z_{j,k}\in dz\right) \geq \kappa'_2$ on $0 \leq z \leq \kappa'_1$, ((ref)) yields
Then, returning to ((ref)), we consider two cases. First, suppose that $0< \bar{\delta}_{j,k} < \kappa'_1$. Then, ((ref)) can be applied directly, yielding
In the second case, we consider $\bar{\delta}_{j,k} \geq \kappa'_1$. The mapping
is increasing on the interval $\left[\kappa'_1,\infty\right)$. Therefore, repeating the arguments used to derive ((ref)), and applying ((ref)) again, we obtain
Combining ((ref)) and ((ref)), we have
for any $\bar{\delta}_{j,k}> 0$. Since the left-hand side of ((ref)) is always positive, we have
for any value (positive or negative) of $\bar{\delta}_{j,k}$.
We may now derive the desired result. Combining ((ref)) and ((ref)), we have
Because, in ((ref)), all of the terms in the double sum are positive, we may take another lower bound in which we only keep those terms in which either $j =1$ or $k=1$. However, for any $j,k$, either $\bar{\delta}_{j,k}\geq 0$ or $\bar{\delta}_{k,j}\geq 0$. Therefore, we have
whence
as required. $\Box$
The second property represents, again, a kind of local strong concavity. In ordinary concavity, the right-hand side of inequality ((ref)) below is simply zero. If we had strong concavity, the right-hand side would be $\|\delta\|^2_2$. Thus, we have that property only when $\|\delta\|_2$ is small, i.e., $g$ is in a neighborhood of $g^*$. For notational convenience, we let $p^g$ be the vector whose components are given by $p^g_k = P\left(c_k\left(X\right)-g_k = F\left(g,X\right)\right)$, noting that $p^{g^*} = p$.
\noindentProof: To avoid notational clutter, we write $g^*$ as simply $g$ in this proof, since no other $g$ will be considered. By the concavity of the function $g\mapsto F\left(g,X\right)$, we have
We then apply Lemma (ref). It remains only to show that
To see this, consider two cases. First, suppose that there is some $j$ for which $\left|\delta_j\right| \geq \kappa'_1$. Then,
Second, suppose that $\left|\delta_k\right| < \kappa'_1$ for all $k$. Then,
This completes the proof. $\Box$
We rewrite ((ref)) in vector form as $g^{n+1} = g^n - \alpha_{n}\left(\zeta^{n+1}-p\right)$, where $\zeta^{n+1}_k = 1_{\left\{k=\arg\min_j c_j\left(X^{n+1}\right)-g^n_j\right\}}$. We may also write a similar recursion
for the error sequence. Before proceeding to the convergence rate, we first prove two intermediate results: an almost sure (but non-vanishing) bound on $\|\delta_n\|_2$, and a concentration inequality showing that the tail probabilities of $\|\delta^n\|_2$ vanish quickly outside a certain bounded range. The proofs are deferred to the Appendix (Section (ref)).
Finally, we prove Theorem (ref). Let $C_0 = 2$ and define $p^n = p^{g^*+\delta^n}$. Using ((ref)), we may write
where ((ref)) uses Lemma (ref) with an appropriate $C_1 > 0$. Letting $\mathcal{F}^n$ be the sigma-algebra generated by $X^1,...,X^n$, we have
On the event $\left\{\|\delta^n\|_2 \leq C_0\right\}$, we have
Therefore,
where $C_2 > 0$ is a constant obtained from Lemma (ref). Taking unconditional expectations and applying Lemma (ref), we obtain
Therefore, for some suitable $C_3,C_4 > 0$, we have
Letting $a_n = n\mathbb{E}\left(\|\delta^{n+1}\|^2_2\right)$, we have
A technical lemma in the Appendix (Lemma (ref)) proves that a non-negative sequence $\left\{a_n\right\}$ satisfying ((ref)) is uniformly bounded if the coefficient of $a_n$ on the right-hand side is positive. Thus, we obtain the desired result as long as $\alpha > \frac{1}{C_3}$.
Finally, we may prove our main result (Theorems (ref)-(ref)). As in Section (ref), we will focus on the case of linear costs $c_k\left(X\right) = \beta^\top_k X$ where ((ref)) is used for estimation. The Appendix presents a fully rigorous extension to the case of more general link functions $\chi$. The analysis in that extension will be very similar, with some additional technical complications in several steps.
The order of the arguments is similar to Section (ref). Section (ref) proves analogous properties to those in Section (ref) for the case of unknown costs. Section (ref) outlines the proof of the main rate results, with some intermediate steps moved to the Appendix. The core of the proof is the derivation of new concentration inequalities for both the SA error and the estimation error. The rates of these inequalities allow us to place a sub-exponential gap between periods of forced exploration.
We proceed largely the same way as in Section (ref). This time, however, we generalize the notation to allow for different choices of the regression coefficients. Let $F^b\left(g,x\right) = \min_j b^\top_j x - g_j$, and let $p^{b,g}$ be the vector whose components are given by $p^{b,g}_k = P\left(b^\top_k X - g_k = F^b\left(g,X\right)\right)$. Note that $F^{\beta} = F$ and $p^{\beta,g^*} = p$.
The first two results are extensions of Lemmas (ref)-(ref). Their main purpose is to prove local strong concavity in the presence of estimation error for both $g^*$ and $\beta$. The overall structure of the proofs closely follows that of their counterparts in Section (ref), but some additional work is required to handle the error of the regression coefficients. The proofs are deferred to the Appendix.
The final property is completely new to the case of unknown costs, and bounds the estimation error of the vector $p$ of target probabilities in a setting where we know $g^*$, but not $\beta$.
\noindentProof: To avoid notational clutter, we write $g^*$ as simply $g$ in this proof, since no other $g$ will be considered. Note, however, that we do distinguish between $b$ and $\beta$.
As in the proof of Lemma (ref), we define $Z^{\beta}_{j,k} = \beta^\top_j X-g_j-\left(\beta^\top_k X-g_k\right)$. For additional convenience, denote $d_k = \Delta^\top_k X$. Then, $b^\top_k X = \beta^\top_k X + d_k$. Observe that, for any $k$,
By Assumption (ref)(ii), $\left|d_k\right| \leq \kappa_4\Delta^{\max}$. We focus on the first term of ((ref)) and derive
Using very similar arguments, the same bound can be shown for the second term of ((ref)). Thus, we have
Consider the case $\Delta^{\max}\leq \frac{\kappa'_1}{2\kappa_4}$. By Lemma (ref)(i), ((ref)) yields
At the same time, since we always have $\left|p^{b,g^*}_k - p_k\right| \leq 2$, we may continue ((ref)) as
We observe the elementary inequality that, for any $a_1,a_2,a_3 > 0$ and $\rho\geq 0$,
Applying this to ((ref)) yields the desired result. $\Box$
Finally, we provide an outline of the steps made in proving Theorems (ref)-(ref). The notation from Section (ref) is carried over. The notation $\|\cdot\|_{\text{sp}}$ will be used to denote the spectral norm (largest singular value).
The first result relates the probability of incorrect selection (for an arbitrary policy) to the estimation error of both $g^*$ and $\beta$. It is used in several places. Theorem (ref) will be a direct consequence of it once Theorem (ref) is established.
\noindentProof: We observe that
so it is sufficient to consider a single $k$. Let $d_k = \left(b_k-\beta_k\right)^\top X$. By repeating some arguments in the proof of Lemma (ref) (starting from eq. ((ref)) onwards), we obtain
The arguments in the proof of Lemma (ref) can then be repeated to yield the desired bound on ((ref)). $\Box$
Next, we state three intermediate technical results related to the information matrix of the observations collected using the proposed semi-myopic policy. The proofs of these results are deferred to the Appendix.
We now begin to study the properties of the semi-myopic algorithm. We first establish concentration inequalities on the two types of estimation error (that is, the error of estimating $\beta$ and $g^*$, respectively). The proofs are deferred to the Appendix, but we note that this is essentially the same approach as in Section (ref); in particular, the statement and proof of Lemma (ref) are similar in structure to Lemma (ref). The difference is that we now require separate concentration inequalities for two interrelated objects. The result of Lemma (ref) is used inside the proof of Lemma (ref).
Finally, we complete the proof. Theorem (ref) is the combination of Propositions (ref) and (ref) below. The proof of Proposition (ref) is highly technical and deferred to the Appendix. Proposition (ref) is an extension of Theorem (ref) incorporating estimation error from regression (and using the result of Proposition (ref)). As mentioned earlier, Theorem (ref) follows straightforwardly from these results together with Lemma (ref).
\noindentProof: Let $\Delta^{\max,n} = \max_k \|\Delta^n_k\|_2$. As in the proof of Theorem (ref), we first write
where ((ref)) is obtained by applying Lemmas (ref)-(ref) with appropriate $D_1,D_2,D_3 > 0$. Letting $\mathcal{F}^n$ be the sigma-algebra generated by $\pi^1,...,\pi^{n}$, $W^1_{\pi^1},...,W^n_{\pi^{n}}$, and $X^1,...,X^n$, we write
Let $D_0 = \frac{\kappa_6}{4C_1}$, where $C_1$ is the constant obtained from Lemma (ref). On the event $\left\{\|\delta^n\|_2 \leq D_0\right\}$, we obtain
by repeating arguments from the proof of Theorem (ref). Then, in general, we have
where $D_4>0$ in ((ref)) is a constant obtained from Lemma (ref). We now take the unconditional expectation, obtaining
where ((ref)) uses Lemma (ref), and ((ref)) uses Proposition (ref). We thus have
for suitable $D_5,D_6,D_7>0$. Letting $a_n = n\mathbb{E}\left(\|\delta^n\|^2_2\right)$, we have
As in the proof of Theorem (ref), we apply Lemma (ref) and obtain the desired result as long as $\alpha > \frac{1}{D_5}$. $\Box$
We illustrate our framework on two numerical examples with linear costs. Section (ref) considers synthetic test instances with $X$ generated uniformly on the unit sphere. Section (ref) applies our framework to learn facility locations in a geographical partitioning problem (a common application area of semidiscrete optimal transport). Both types of experiments seek to explore the performance of the semi-myopic policy in settings that lie outside the assumptions used in our theoretical analysis.
We consider a set of $1000$ synthetic test problems. In all of these problems, the data $X$ are generated uniformly on the unit sphere in $d=10$ dimensions; since $X$ does not have a density in $\mathbb{R}^d$, Assumption (ref) does not hold. Each problem uses $K = 5$ different vectors $\beta_k$, also generated uniformly on the unit sphere, and a different $p$ generated uniformly on the $d$-dimensional simplex.
For ease of implementation, we modified the semi-myopic policy to conduct forced exploration probabilistically: that is, in the $n$th iteration, we select $k$ uniformly at random with probability $\frac{\left(\frac{1}{a}\log n\right)^9}{n}$, where we used $a = 4.5$. The rest of the time, we select a cost function according to the myopic policy $\hat{\pi}^n$. We set $\rho^n\equiv 0.001$ for simplicity, and in order to enable use of recursive least squares. At each $n$, we evaluate the policy using the correct selection indicator $1_{\left\{\hat{\pi}^n=\pi^*\left(X^n\right)\right\}}$. By averaging this quantity across $10$ runs, each containing $1000$ iterations and starting from $g^0 \equiv 0$, $\hat{\beta}^0_k \equiv 0$, we obtain a trajectory of the probability $P\left(\hat{\pi}^n = \pi^*\left(X^n\right)\right)$ of correct selection. We then average this trajectory over $1000$ test problems.
The optimal policy $\pi^*$ for each test problem is precomputed via brute force by running ((ref)) with the true cost coefficients for a very large number of iterations. As a benchmark, we also implement a policy that knows the true costs (and, therefore, does not need to conduct any exploration or collect any samples), but does not know $g^*$, and runs ((ref)) in an online manner. We also calculate the average probability of correct selection for this policy, which allows us to evaluate any loss incurred by having to learn the coefficients.
Figure (ref) shows average performance for the semi-myopic policy with stepsize parameter $\alpha = 5$ and two levels of Gaussian noise: $Var\left(\varepsilon\right) = 0.02^2$ and $Var\left(\varepsilon\right)=0.2^2$. Since the benchmark policy knows the costs, it is unaffected by the noise variance. We see that, when the variance is low, there is virtually no performance loss caused by unknown costs. When the variance is high, the average PCS of the semi-myopic policy lags behind that of the known-cost policy by $0.039$ (and has narrowed this gap relative to the earlier iterations). Performance also does not appear to be very sensitive to the frequency of forced exploration: we also tested a much higher value $a = 6$, for which the probability of forced exploration is only $0.0036$ by the end of the time horizon, but obtained very similar trajectories. This also indicates that, in an online implementation, we could allow forced exploration to happen quite infrequently while maintaining good performance.
Figure (ref) uses the higher variance level and varies the stepsize parameter $\alpha$ over several orders of magnitude. We observe some sensitivity with respect to this parameter, but, interestingly, this sensitivity is partially smoothed out by the presence of noise: under known costs, the difference between $\alpha = 0.5$ and $\alpha = 50$ is far more pronounced than under unknown costs. By the end of the time horizon, there is virtually no difference between the two best-performing choices, and the third-best is relatively close as well. Only the largest choice $\alpha = 500$ lags behind the others, perhaps reflecting the fact that the multiplicative constant in Theorem (ref) is increasing in $\alpha$.
Figure (ref) repeats the analysis of Figure (ref) for the same experimental setup, with the sole difference that the noise follows a Laplace (rather than Gaussian) distribution with mean zero and $Var\left(\varepsilon\right)=0.2^2$. Thus, this experiment considers a situation where the assumption of sub-Gaussianity is violated (instead, the noise is sub-exponential). However, despite the heavier tails, the results are virtually identical to the normal case, indicating that our approach is fairly robust to violations of this assumption. We also compared noise levels similarly to Figure (ref), but the results were again nearly identical to the Gaussian case and so we omit them here.
Finally, Figure (ref) uses the same setup for $X$ and $\beta_k$, but now we assume that each observation $W^{n+1}_k$ is generated from a Bernoulli distribution with success probability $\chi\left(\beta^\top_k X^{n+1}\right) = \frac{1}{1+e^{-\beta^\top_k X^{n+1}}}$. In other words, the true cost structure is an instance of logistic regression with coefficients $\beta_k$. However, we do not use logistic regression to estimate the coefficients, due to the high computational cost of estimating a GLM. Instead, we apply recursive least squares directly to the binary observations, treating the error $W^{n+1}_k - \chi\left(\left(\hat{\beta}^n_k\right)^\top X^{n+1}\right)$ as the residual. Thus, the semi-myopic policy works with a partially misspecified statistical model, though ChRy20 showed that recursive least squares is still consistent in this situation (i.e., we will still learn the true $\beta_k$ values given enough time). Since estimation is more difficult with misspecification, we extend the time horizon to $10^4$ iterations for this experiment.
Because the benchmark policy knows the costs, misspecification is not a concern for it, and its performance only depends on $\alpha$. The semi-myopic policy has much more difficulty in the misspecified setting than in the previous two, but makes consistent improvement throughout the majority of the time horizon. Because cost estimation is now the most difficult aspect of the problem, the performance of the semi-myopic policy is now almost entirely unaffected by $\alpha$.
It is well-known in the stochastic approximation literature that the performance of SA can be sensitive to the stepsize rule. It may be possible to speed up convergence by using other choices of $\alpha_n$, many of which are surveyed in GePo06. Other, more sophisticated procedures can be found in BrCiZe11 or ScZhLe13, but every one of these methods involves at least one tunable parameter, so there is no way to avoid the tuning issue entirely. However, it is noteworthy that the presence of noise in the observed costs has the effect of reducing sensitivity to $\alpha$ in all three of the settings considered here.
For each $k = 1,...,K$, let $x_k \in \left[0,1\right]^2$ be a fixed location (“facility”) in the unit square. A user appears at some random location $X \in \left[0,1\right]^2$ and is assigned to a facility. The cost incurred by assigning the user to facility $k$ is $c_k\left(X\right) = \|X-x_k\|_2$, the Euclidean distance between the facility and the user. The goal is to choose an assignment rule that minimizes the expected distance subject to the probabilistic targets (which, in this context, require each facility to serve a certain proportion of the population). It is well-known in computational geometry that the optimal policy can be visualized as an additively weighted Voronoi diagram CaCaDe16 whose weights are the values $g^*_k$. The diagram partitions the unit square into regions, each containing a single facility, to which all users appearing in that region are assigned.
We now suppose that the locations $x_k$ are unknown, but a user appearing at location $X$ can obtain a noisy observation $W_k = \|X-x_k\|_2 + \varepsilon$ of the distance to any one facility of his or her choice. This observation is not linear in $x_k$, but we can still apply a linear regression model using the transformation
The last two terms in ((ref)) have zero mean and can be treated as “residual error.” Because $X$ is observed, the quantity $X^\top X$ is known. Thus, we can learn $x_k$ using a linear regression model in which $\left[1,-2X^\top\right]^\top$ is the vector of features, and $W^2_k - X^\top X$ is the response variable. The coefficients to be learned are $\beta_k = \left[\mathbb{E}\left(\varepsilon\right)^2 + x^\top_k x_k,x^\top_k\right]^\top$. Clearly, $x_k$ is identifiable in this model.
We considered two instances of this model where $\varepsilon \sim \mathcal{N}\left(0,0.02^2\right)$ and $X$ is uniformly distributed on $\left[0,1\right]^2$. Note that, even with the normality assumption, the residual error in ((ref)) is not sub-Gaussian. The instances used different facility locations $x_k$ and targets $p$. The settings of the semi-myopic policy were the same as in Section (ref).
In this discussion, we aim to take advantage of the unique visual interpretation offered by the geographical partitioning setting. Thus, we do not present average trajectories of PCS over a large number of runs, as we did in Section (ref). Rather, we will visualize the final partition learned by the semi-myopic policy after a large number $10^6$ of iterations, and compare this to the optimal partition induced by $\pi^*$. Unlike in Section (ref), we do not need to separately present a policy that knows $x_k$ but attempts to learn $g^*$ in an online manner, because the number of iterations is so large that such a policy will have converged to the true optimal partition already.
Figure (ref) presents the true and estimated partitions for both instances, with the facility locations $x_k$ shown as black dots. The difference between locations turns out to exert a very significant impact on performance. In Instance 1, the facility close to the right edge of the unit square is also assigned a small $p_k$ value (specifically, $0.0258$). Thus, only a small portion of observations will lead us to sample this facility, causing us to rely more on forced exploration. Even after $1$ million iterations, the PCS achieved by the semi-myopic policy (i.e., the proportion of users assigned to the correct facility) is only $0.8944$. In Instance 2, this is less of a concern, and the semi-myopic policy achieves a PCS of $0.9542$. The estimated partition is also visibly closer to the true one in Instance 2.
In general, the estimation error of $\beta_k$ is much more of an issue in the partitioning problem than it was in Section (ref). Even after $1$ million iterations, the semi-myopic policy noticeably lags behind SA with known costs (which was used to find the true partition). We found that our policy was effective in learning the regression coefficients, with $\hat{\beta}^n_k$ being within three decimal places of $\beta_k$ by the end of the time horizon. However, by the time this happened, the stepsize in the gradient update had become so small that $g^n$ had fallen behind (though it is still making progress, as PCS continues to increase if we keep running the policy). One reason may be that the regression noise in this problem is not sub-Gaussian, but has heavier tails than it did in Section (ref), so the gradient update is more heavily biased. Even so, this example shows that the semi-myopic policy is able to learn the true partition even when some of the underlying model assumptions are violated.
We have presented a novel variant of the semidiscrete optimal transport problem in which the costs are unknown, but can be estimated through sequential observations. This process is subject to the usual tradeoff of optimal learning, where sampling one cost function means one less sample available for the others. The problem thus brings together elements of optimal transport, optimal learning, and (through the characterization of the optimal policy) stochastic approximation. We have developed a simple and provably efficient algorithm based on the well-established notion of semi-myopic exploration; however, our analysis couples this structure with stochastic approximation in a novel way. In the process, we are able to resolve several known issues with applying SA to semidiscrete optimal transport, having to do with the non-smoothness of the stochastic gradient and the lack of strong concavity in the objective function.
There are many opportunities for future work. First, it may be interesting to investigate the theoretical performance of other types of learning algorithms, such as Thompson sampling AgGo13, in this particular context. Second, it is worth taking a closer look at practical issues such as stepsize selection, which is known to impact the performance of SA procedures. Third, there may be potential for further connections between optimal transport and optimal learning beyond the semidiscrete version studied here. We hope that our work will spark new interest in bridging these fields.