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.
71,042 characters · 10 sections · 53 citation commands
Quantile optimization in semidiscrete optimal transport
{.5\baselineskip}
Optimal transport Vi21 is a rich and challenging branch of optimization focusing on the efficient transformation of one probability distribution into another. Given two random variables, $X$ and $Y$, with known marginal distributions, the goal is to design a joint distribution that minimizes some cost function depending on both $X$ and $Y$. In the special case where $X$ and $Y$ are both discrete, optimal transport reduces to the well-known “transportation problem” FoFu56, in which linear programming is used to match supply to demand. However, when one or both of the random variables are continuous, the problem becomes much more difficult. Applications of these more general settings arise in facility logistics CaCaDe16, computer vision BoDi23, text mining Yurochkin19, portfolio optimization Nguyen25, and algorithmic fairness Chzhen20. Optimal transport also interfaces with more theoretical areas of operations research, e.g., the literature on distributionally robust optimization, which often uses very similar representations of cost BlMuZh22,GaKl23. The purely theoretical literature on optimal transport is surveyed in multiple monographs such as Sa15, PeCu19, Vi21, and Ma23.
The objective optimized in this extensive literature is, invariably, an expected cost. In other words, if $\pi$ is the joint distribution of $\left(X,Y\right)$, the problem is written as $\inf_{\pi} \mathbb{E}\left[c\left(X,Y\right)\right]$ for some cost function $c$, subject to constraints that preserve the marginal distributions. In fact, the vast majority of the optimal transport literature is focused on one specific cost function, namely $c\left(x,y\right) = \|x-y\|^2_2$, in which case the objective becomes the well-known Wasserstein distance metric CuPe18. Some theoreticians have examined other choices, such as the $1$-norm XiLiZhZh24, the $\infty$-norm ChDeJu08, and more generally, cost functions satisfying the “twist condition” McRi16, a class that includes $p$-norms. Even these special cases are somewhat rarely seen, and in any case, they are also minimized in expectation.
The present paper is, to our knowledge, the first study of optimal transport where the objective function is not the expected cost, but a quantile of it, that is,
where $\alpha \in \left(0,1\right)$ and
The rest of the optimal transport framework remains unchanged: we still choose a joint distribution under pre-specified marginals. In general, quantile optimization is used when the underlying distributions are heavy-tailed or skewed, so that the means either do not exist or are heavily affected by extreme values. Such situations can arise in energy KiPo11, finance CaGaMoOl22, and even traditional OR application areas like inventory management GlLi96, where rare events such as supply chain disruptions may yield unusually high costs. One may then prefer to optimize the median or an upper quantile of the cost to produce more risk-averse decisions. With regard to optimal transport in particular, we may consider the problem of geographical partitioning CaDe13, where a territory is divided between facilities with given capacities; here, the decision-maker may prefer to design the service zones to, e.g., reduce the incidence of very high travel times.
Our study takes place in the context of semidiscrete optimal transport, meaning that $X$ has a general (continuous) distribution, but $Y$ is discrete. While this is not the most general setting, there are compelling reasons to examine it. First, semidiscrete OT is, in and of itself, highly relevant to applications, including the aforementioned geographical partitioning problem, as well as economic choice modeling BoPe24, artificial intelligence Chen19, data analytics Gao24, and computational biology SePrDa22. Second, semidiscrete OT is far more flexible regarding the choice of $c$ than the general case, and does not have to be restricted to squared distances or $p$-norms. It is worth noting that, in geographical partitioning, the most commonly used cost function is the unsquared Euclidean distance HaSc20, which already falls outside much of the OT literature. Third, as TaShKu22 points out, many numerical algorithms for general optimal transport work by discretizing at least one of the distributions anyway; for example, the celebrated Sinkhorn algorithm Cu13 is based on discretization.
In fact, outside of very narrow special cases, it is exceedingly difficult to compute exact solutions for general OT, and the relevant literature offers mainly theoretical characterizations GaSw98. Semidiscrete OT is the exception to the rule. GeCuPeBa16 was the first to show that, in the classical formulation where the expected cost is minimized, it is possible to learn the exact optimal solution using very simple, fast, simulation-based iterative methods (essentially a one-line stochastic approximation algorithm; see, e.g., KuYi03) without having to discretize $X$. This tractability hinges on the fact that, in the semidiscrete setting, although $\pi$ itself is infinite-dimensional, the optimal solution can be completely characterized using only finitely many parameters. Thus, by studying quantile optimization in semidiscrete OT, we may sacrifice some nominal theoretical generality (partially compensated by the ability to handle a much greater variety of cost functions), but we gain a very clean characterization of the optimal solution that (as we show) can be computed very efficiently.
We provide a comprehensive analysis of the problem, proceeding in several steps:
The results of this analysis leave us with two computational challenges: solving the root-finding problem, and computing the tie-breaking rule. Both problems can be approached using simulation, i.e., we generate data from the distribution of $X$ and use these samples to estimate various difficult expectations. We do not even need to know the distribution of $X$; we only require a means of sampling from it. Thus, we may leverage techniques used in the simulation community. As laid out by KiHe08, simulation-based methods generally fall into two categories: stochastic approximation (SA, also known as stochastic gradient descent), and sample average approximation (SAA). Both approaches have complementary strengths, and in fact, we have reason to use both in this work. Although stochastic approximation is often used to solve stochastic root-finding problems (even for some types of quantile optimization problems; see, e.g., HuPeZhZh22), in this particular case that problem is coupled with a difficult convex program whose value is not differentiable in the decision variable. Therefore, we instead use SAA to learn the optimal quantile, because it requires us to solve a relatively simple linear program and has strong convergence guarantees. On the other hand, when it comes to learning the tie-breaking rule, SAA is unwieldy because there are too many constraints, while SA turns out to be extremely easy to compute. We give a convergence proof showing that the true optimal transport plan is learned at the canonical $\mathcal{O}\left(N^{-\frac{1}{2}}\right)$ rate, where $N$ is the sample size.
Finally, we present numerical examples illustrating the computational tractability of our approach in the setting of geographical partitioning, and in the process obtain new insights into that application. For these problems, traditional (mean-minimizing) semidiscrete OT is solved by a classical geometric structure known as the additively weighted Voronoi diagram Au91. Quantile optimization gives rise to a completely different, partially randomized structure, which, to our knowledge, has never been observed before in the literature on computational geometry. We visualize and contrast both mean- and quantile-optimizing solutions and show their impact on the cost distribution.
To summarize, our paper makes the following contributions. 1) We are the first to optimize a quantile objective for any kind of optimal transport. 2) For the semidiscrete setting, we develop a novel analytical framework that produces a clean characterization of the optimal solution. Two key aspects of the solution, both novel relative to the traditional setting, are that the optimal quantile solves a root-finding problem, and that the optimal $\pi$ requires a probabilistic tie-breaking rule. 3) We show how both aspects may be addressed using simulation-based methods, which are both computationally tractable and provably efficient. In particular, the estimated choice probabilities characterizing the optimal transport plan converge to their true values at the canonical square-root rate. 4) Using these methods, we compute and illustrate a novel geometric structure for the partitioning problem. In this way, our work provides cross-disciplinary insights that, we hope, will be of interest to operations researchers and computational geometers, as well as to the community working on optimal transport.
Formally, our problem has the objective ((ref)) subject to the constraints
where $\mathcal{X}\subseteq \mathds{R}^d$ is the support of $X$, and $\left\{1,...,K\right\}$ is the support of $Y$. We assume that $c$ satisfies $P\left(c\left(X,k\right)<\infty\right) = 1$ for any $k$. The decision variable in ((ref))-((ref)) is $\pi$, which can be viewed as the mixed joint likelihood of $\left(X,Y\right)$. In other words, $P\left(X\in dx,Y=k\right) = \pi\left(x,k\right)dx$. The density $f$ and the probability vector $p$ (with $p_k > 0$ for all $k$) are fixed problem inputs representing the marginal density of $X$ and the marginal pmf of $Y$, respectively. Constraints ((ref))-((ref)) ensure that $\pi$ respects the pre-specified marginals.
Our technical analysis consists of three parts. First, in Section (ref), we formulate a simpler feasibility determination problem and characterize its solutions. Section (ref) then connects this problem to quantile optimization; essentially, the optimal quantile is related to a particular instance of the feasibility determination problem with appropriately specified inputs. Finally, Section (ref) characterizes “tie-breaking rules,” which are needed to design an implementable optimal transport plan that achieves the desired quantile.
Before proceeding to quantile optimization, let us first consider a simpler feasibility problem. For fixed $t$ and $\alpha$, we wish to determine if there exists a joint distribution $\pi$ satisfying ((ref))-((ref)) as well as
The answer is not obvious because $\alpha$ and $t$, as well as $f$ and $p$, are arbitrary. If a feasible distribution does not always exist, we wish to characterize the conditions under which feasibility occurs.
We may rewrite ((ref)) as
a linear constraint in $\pi$. Thus, ((ref)) together with ((ref))-((ref)) is an infinite-dimensional system of linear equations. With the addition of the nonnegativity constraints ((ref)), we may apply the following infinite-dimensional analog of Farkas' lemma. The proof is given in the Appendix.
Thus, determining the feasibility of the original system is equivalent to solving the infinite-dimensional linear program
subject to ((ref)), and verifying that the optimal value is positive. We may simplify this problem by rewriting ((ref)) as
and setting
for every $x\in\mathcal{X}$. Substituting ((ref)) into ((ref)) yields the nonlinear, but finite-dimensional, unconstrained, and convex optimization problem
As the next result shows, the feasibility of the system ((ref))-((ref)) is determined entirely by whether ((ref)) has an optimal solution.
\noindentProof: If ((ref)) is unbounded, ((ref))-((ref)) are feasible by Theorem (ref). Suppose now that ((ref)) is not unbounded. Then, because this problem is unconstrained and thus cannot be infeasible, it must have an optimal solution $\left(\zeta^*,\psi^*\right)$.
Suppose that the optimal value of ((ref)) is strictly negative, i.e.,
The inequality is preserved if we multiply both sides by arbitrarily large $a > 0$. Thus, we can find $a$ for which $\left(a\zeta^*,a\psi^*\right)$ has a better objective value than $\left(\zeta^*,\psi^*\right)$, contradicting the fact that $\left(\zeta^*,\psi^*\right)$ is optimal. Therefore, the optimal value of ((ref)) cannot be strictly negative, whence ((ref))-((ref)) are infeasible by Theorem (ref). Finally, because $\zeta = 0$ and $\psi = 0$ are feasible for ((ref)), the optimal value of ((ref)) must be zero. $\Box$
Optimal solutions of ((ref)), if they exist (they need not be unique), also provide important structural insight into the original feasibility problem, i.e., the solution of ((ref)) and ((ref))-((ref)). Suppose that some nontrivial $\left(\zeta^*,\psi^*\right)$ is optimal for ((ref)), and consider a random variable $Y$ satisfying
for all $\omega$. There is some indeterminacy in ((ref)) because the argmin may not be unique. However, if there exists a (possibly random) tie-breaking rule that causes $Y$ to have the correct marginal distribution, then such a $Y$ will be feasible for the original problem. We formally prove this result as follows.
\noindentProof: Constraints ((ref))-((ref)) are satisfied by the properties of $Y$. It remains to show ((ref)).
First, by Proposition (ref), the optimal value of ((ref)) must be zero. Therefore,
At the same time, we can derive
where ((ref)) follows from ((ref)), and ((ref)) holds because $P\left(Y=k\right) = p_k$. Setting ((ref)) equal to ((ref)) yields ((ref)), as required. $\Box$
It is worth noting that the lack of a unique argmin in ((ref)) constitutes a profound difference between our setting and classical optimal transport. Much of the OT literature is focused on situations where $Y$ can be written as a deterministic function of $X$, known as a “Monge map.” This is the well-known Kantorovich-Monge equivalence property Pr07,Be21. In our setting, however, the mass $\pi\left(x,\cdot\right)$ must be divided between different $\pi\left(x,k\right)$, and the tie-breaking rule provides the correct splitting.
We have not yet said anything about how one may compute a tie-breaking rule satisfying the conditions of Proposition (ref). In fact, thus far we have not yet established that such a rule even exists. We will return to these issues later, in Section (ref). At the moment, however, we can see how the existence of an optimal solution ((ref)) guarantees the feasibility of the original system and allows (potentially) a feasible joint distribution to be constructed.
We now use the insights from Section (ref) to develop an approach for quantile optimization. Again, by Theorem (ref), the original system described by ((ref)) and ((ref))-((ref)) is feasible if
for all $\zeta,\psi$. Note that, if $\zeta = 0$, then ((ref)) reduces to $\sum_k p_k\psi_k - \min_k \psi_k \geq 0$, which always holds for all $\psi$. If $\zeta \neq 0$, we may rewrite ((ref)) as
Since we may rescale $\psi$ without changing the inequality, it is sufficient to verify ((ref)) for $\zeta \in\left\{-1,1\right\}$. In other words, for all $\psi$, we must have
Therefore,
must hold for all $\psi$. Equivalently, we may write
Using a change of variables $\psi' = -\psi$ for the upper bound, we have
By the max-min inequality, the lower bound in ((ref)) is always less than or equal to the upper bound. Both the maximization problem in the lower bound and the minimization problem in the upper bound are unconstrained, and therefore feasible. Consequently, if we are given a value of $t$, ((ref)) describes a nonempty interval of $\alpha$ values for which it is possible to make that $t$ the $\alpha$-quantile of the cost.
If our goal is quantile optimization, we should make $\alpha$ as large as possible for fixed $t$. In other words, given a fixed threshold for the cost, we would prefer to have as much of the cost distribution below this threshold as possible. We can connect the upper bound in ((ref)) back to ((ref)) through the following result.
\noindentProof: Problem ((ref)) is not unbounded because $\alpha^*$ is, by construction, an element of the interval ((ref)). Plugging $\left(\zeta^*,\psi^*\right)$ into the objective of ((ref)) yields
where the last line is due to the fact that $\tilde{\psi}$ is a minimizer in ((ref)). The desired conclusion follows by Proposition (ref). $\Box$
We also have additional structure on the optimization problem in ((ref)), namely, that the right-hand side can always be minimized by binary-valued variables. The proof of this property is deferred to the Appendix so as not to distract from the presentation.
We now have a correspondence between $t$ and $\alpha^*\left(t\right)$, the largest feasible quantile for that cost threshold. Proposition (ref) then provides us with a corresponding solution of ((ref)), and Proposition (ref) links that solution to a choice of $Y$ (modulo a suitable tie-breaking rule, which still remains to be discussed) that solves ((ref)) and ((ref))-((ref)) for $t$ and $\alpha^*\left(t\right)$.
The quantile optimization problem that we originally set out to solve is then readily seen to be an inversion of this correspondence: given $\alpha$, we search for $t$ such that $\alpha^*\left(t\right) = \alpha$. It is obvious that $\alpha^*\left(t\right)$ is increasing in $t$, and with some mild regularity conditions, it becomes continuous and strictly increasing, so a suitable $t$ always exists and is unique. The following results examine some simple conditions under which these properties hold; the proofs are deferred to the Appendix.
In the following, suppose that we already have the desired value of $t$ (the corresponding $\alpha^*$ is simply the target quantile $\alpha$). Thus, we have the optimal value of the quantile optimization problem ((ref)). For practical purposes, however, we wish to also know the optimal solution, i.e., the joint distribution $\pi$ of $\left(X,Y\right)$ that achieves the infimum.
We abuse notation slightly by letting $\psi^*_k$ be a minimizer of the right-hand side of ((ref)); from Proposition (ref), we know that any such minimizer can be converted into an optimal solution of ((ref)), so it is no longer necessary to distinguish between them in the notation. By combining Propositions (ref) and (ref), we know that a suitable $\pi$ can be obtained from any $Y$ that satisfies
with ties in the argmax broken in a (possibly randomized) manner that satisfies $P\left(Y = k\right) = p_k$. We will now prove that a randomized tie-breaking rule satisfying these conditions always exists, and give a characterization that will later enable its tractable computation.
First, let us observe that, by definition,
which may be rewritten as
This means precisely that $\psi^*$ is optimal if and only if $p$ is a subgradient of the function
at $\psi^*$. In order to eventually establish the existence of an optimal tie-breaking rule, it is necessary to describe all of the subgradients of $h$.
For notational convenience, let $H_k\left(\psi,X\right) = 1_{\left\{c\left(X,k\right)\leq t\right\}} + \psi_k$ and $H\left(\psi,X\right) = \max_k H_k\left(\psi,X\right)$. Then, $h\left(\psi\right) = \mathbb{E}\left(H\left(\psi,X\right)\right)$. Let us first focus on the subgradients of $H\left(\psi,x\right)$, holding $x$ fixed, and then pass to the expectation. We will also use the notation
so that we do not have to keep writing out the argmax.
For any $k$ and $x$, $H_k\left(\psi,x\right)$ is differentiable with $\nabla_\psi H_k\left(\psi,x\right) = e_k$, where $e_k$ is a vector of zeroes with the $k$th component equal to $1$. By Lemma 3.1.10 of Ne04, the subdifferential $\partial H\left(\psi,x\right)$ of $H\left(\psi,x\right)$ at $\psi$ is given by
In words, $\partial H\left(\psi,x\right)$ is the probability simplex on only those coordinates $k$ that are elements of the argmax. To put it differently, $v \in \partial H\left(\psi,x\right)$ represents a discrete probability distribution on the elements of $S\left(\psi,x\right)$. Note that the subdifferential only depends on $x$ through $S\left(\psi,x\right)$.
By the main theorem of RoWe82, the subdifferential of $h$ is given by $\partial h\left(\psi\right) = \mathbb{E}\left(\partial H\left(\psi,X\right)\right)$. Here, we are referring to the Aumann expectation Au65 of a random set $\partial H\left(\psi,X\right)$ whose realization is determined by $S\left(\psi,X\right)$; in other words, the expectation may be taken over the discrete distribution of $S\left(\psi,X\right)$, rather than the continuous distribution of $X$. Formally, $w \in \partial h\left(\psi\right)$ if and only if we can write
where
and $v\left(S\right) \in \mathbb{R}^K_+$ satisfies $\sum_k v_k\left(S\right) = 1$, $v_k = 0$ for $k {\not \in} S$. We may restrict the summation set in ((ref)) to only those $S$ for which $Q\left(\psi,S\right) > 0$.
Now recall from before that $p \in \partial h\left(\psi^*\right)$. Therefore, by ((ref)), for every $S$ with $Q\left(\psi^*,S\right) > 0$, we may find $v_S$ such that $p = \sum_S Q\left(\psi,S\right) v_S$. Each vector $v_S$ is a probability distribution, and this distribution is precisely our randomized tie-breaking rule: we first observe $X$ and therefore $S\left(\psi^*,X\right)$, and then let $P\left(Y = k\mid S\left(\psi^*,X\right)\right) = v_{S\left(\psi^*,X\right),k}$. Thus, ((ref)) is satisfied because $v_{S\left(\psi^*,X\right),k} > 0$ only for $k \in S\left(\psi^*,X\right)$, and
where the last equality follows by ((ref)).
We have now established that a suitable tie-breaking rule always exists. However, in order to compute it, we must find a probability distribution $v_S$ for each $S$ such that $p = \sum_S Q\left(\psi^*,S\right) v_S$. This is difficult, in part because the number of distributions to be computed is combinatorially large, and in part because the probabilities $Q\left(\psi^*,S\right)$ are computationally challenging; however, in Section (ref), we will present a computationally tractable scheme that does not require knowledge of $Q$ or enumeration of the possible $S$.
Our problem poses two algorithmic challenges. The first is to find the root $\alpha^*\left(t\right) = \alpha$ for given $\alpha$, as described in Section (ref), and the second is to compute the tie-breaking rule using the characterization in Section (ref). Both involve the computation of difficult integrals: in the first case, we must deal with a difficult expectation in ((ref)), and in the second case, we must deal with a combinatorially large number of difficult probabilities in ((ref)).
As discussed in Section (ref), we use a simulation-based approach where independent data $\left\{X^n\right\}^N_{n=1}$ are generated from the distribution of $X$. We will use sample average approximation (SAA) to solve the root-finding problem and learn the optimal quantile, and we will use stochastic approximation (SA) to learn the tie-breaking rule. These algorithms are presented and studied in Sections (ref) and (ref), respectively.
As discussed in Section (ref), our task is to solve $\alpha^*\left(t\right) = \alpha$ given some prespecified $\alpha$. Because the computation of $\alpha^*\left(t\right)$ involves a difficult expectation, we can view this as a stochastic root-finding problem PaKi11. Such problems are often solved using the stochastic approximation (SA) method KuYi03. SA has previously been applied GeCuPeBa16,ZhRy26 to solve semidiscrete optimal transport problems with the standard objective of minimizing expected costs. However, in the setting of quantile optimization, we have two “layers” of root-finding rather than just one, because $\alpha^*\left(t\right)$ for fixed $t$ is itself the optimal value of a stochastic optimization problem. While variants of SA do exist for such settings Bo97, their theoretical guarantees Do22 are weaker, or require much stronger assumptions, than what is available for standard SA. Moreover, our case has the additional complication that, in ((ref)), the expression inside the expected value is not pathwise differentiable due to the possibility of ties between different indices $k$.
For all of these reasons, we do not attempt to learn the optimal quantile through SA. Instead, we use the sample average approximation (SAA) method KiPaHe14,ShDeRu21, which replaces the expectation in ((ref)) by its empirical equivalent given a sample $\left\{X^n\right\}^N_{n=1}$. The optimization problem over $\psi$ can then be reformulated as a linear program by introducing additional decision variables $z_n$, $n=1,...,N$. For fixed $t$, we solve
subject to
This LP has $K+N$ decision variables and $K\cdot N$ constraints, and can be tractable for a fairly large sample size due to the availability of large-scale LP solvers.
However, we must solve a problem like this for every $t$ value that we might consider. Since $\psi$ is bounded by Proposition (ref) and $\hat{\alpha}^N$ only involves an expectation over a finite distribution, the optimal value is bounded. Then, since $\hat{\alpha}^N\left(t\right)$ is increasing in $t$, we can solve $\hat{\alpha}^N\left(t\right) = \alpha$ through the bisection method. Because we are using a sample average instead of the true expectation, $\hat{\alpha}^N\left(t\right)$ may have jumps and a root may not exist; in this case, the procedure will converge to the value
Note that every iteration of the bisection method works with the same sample, i.e., we may approximate $\hat{t}^N$ to arbitrary precision without increasing $N$. The consistency of $\hat{\alpha}^N$ can be obtained through standard techniques ShDeRu21. However, we can go further and characterize the convergence rates of $\hat{\alpha}^N$ and $\hat{t}^N$ under the same assumptions as in Propositions (ref)-(ref). The proofs of the following results can be found in the Appendix.
A consequence of Theorem (ref) is that, with high probability, the SAA solution will give us $\psi$ that is exact for the estimated root $\hat{t}^N$, in the sense of minimizing the right-hand side of ((ref)) with $t = \hat{t}^N$. In other words, all of the SAA estimation error will come from $\hat{t}^N$, not the $\psi$ value obtained from the LP.
Recall from ((ref)) that the joint distribution $\pi$ is designed by computing $Y$ according to a criterion that depends on both $t$ and $\psi$. If $t=\hat{t}^N$, but $\psi\in\Psi^*\left(\hat{t}^N\right)$ with a proper tie-breaking rule, we will obtain a joint distribution that is feasible and has $\hat{t}^N$ as its $\alpha^*\left(\hat{t}^N\right)$-quantile. The only source of error will then be $\hat{t}^N$. Theorem (ref) tells us that we find $\psi$ that is “correct” for the estimated quantile $\hat{t}^N$ with high probability.
If we slightly strengthen the conditions of Theorem (ref) (essentially making $\alpha^*\left(t\right)$ Lipschitz as per Proposition (ref)), we can further show that $\hat{\Psi}^N\left(\hat{t}^N\right) \subseteq \Psi^*\left(t^*\right)$ with high probability. We will use this result in Section (ref) to assess how close our estimated $\pi$ is to the optimal transport plan.
For the moment, let us suppose, as in Section (ref), that the desired $t^*$, and its corresponding $\psi^*$, have already been found. In practice, $t^*$ would be replaced by the SAA estimate $\hat{t}^N$. However, from Theorem (ref), we would still have the correct $\psi^*$ value for $\hat{t}^N$ with high probability. Thus, ((ref)) combined with a suitable tie-breaking rule will produce a cost distribution whose $\alpha^*\left(\hat{t}^N\right)$-quantile is $\hat{t}^N$. The developments below will produce a rule that achieves $\alpha^*\left(t\right)$ for any $t$ as long as $\psi^*$ minimizes ((ref)) for that $t$.
From Section (ref), we know that a valid tie-breaking rule is represented by a collection $\left\{v_S\right\}$ of probability distributions. Each $S\subseteq\left\{1,...,K\right\}$ satisfies $Q\left(\psi^*,S\right) > 0$ and each $v_S \in \mathbb{R}^K_+$ satisfies $\sum_k v_{S,k} = 1$ and $v_{S,k} = 0$ for $k {\not \in} S$. Thus, any solution of the system
is a valid tie-breaking rule. It is crucial that any solution will yield an optimal transport plan through Proposition (ref). Thus, we are free to select any convenient objective function and optimize it under the constraints ((ref))-((ref)). The main difficulty is that the number of constraints is combinatorially large, making it prohibitive to apply SAA as in Section (ref).
We propose to use the expected entropy objective
This objective is strictly convex and thus ((ref)) will have a unique optimal solution. Under this choice, we obtain an unconstrained dual problem with only $K$ decision variables, as shown in the following result.
\noindentProof: By combining Example 3.21, Sec. 3.3.2, and eq. (5.11) in BoVa04, we obtain the dual function
By eq. (5.16) in BoVa04, the dual problem maximizes ((ref)) subject to $\lambda\geq 0$ and no constraints on $w,\theta$. It is optimal to set $\lambda \equiv 0$ and $w_S = Q\left(\psi^*,S\right)\log\left(\sum_{k\in S} e^{-\theta_k -1}\right)$. Plugging these values into ((ref)) leads to ((ref)). $\Box$
Moreover, an optimal solution of the dual can be immediately converted into a valid $v_S$ for any $S$. For implementation purposes, one can do this “on demand,” storing only the dual variables and creating a probability distribution only after $X$ and $S\left(\psi^*,X\right)$ have been observed.
\noindentProof: We first verify feasibility. Clearly $\sum_k v^*_{S,k} = 1$ for all $S$. The first-order optimality conditions of ((ref)) are given by
which is exactly ((ref)).
Lastly, we plug our choice of $v^*$ into the primal objective and find
where the last line follows by ((ref)). So, strong duality holds for the pair $\left(v^*,\theta^*\right)$, which completes the proof. $\Box$
Finally, we explain how to optimize the dual. Problem ((ref)) can be written as $\max_\theta p^\top \theta - \mathbb{E}\left(\log \sum_{k\in S\left(X\right)} e^{\theta_k}\right)$. The expression inside the expected value is pathwise differentiable (unlike in Section (ref)), giving rise to a simple stochastic approximation algorithm
where $\gamma_m$ is a stepsize and $X^{m+1}$ is an independent sample from the distribution of $X$. We use $m$ instead of $n$ in the superscript because, in a standard implementation of SA, a new independent sample (separate from the ones that were used to estimate $\hat{t}^N$) would be required in every iteration to ensure the validity of the procedure. However, this algorithm does not require us to know $Q$, to enumerate $S$, or to store $v_S$ for every $S$. The initial iterate $\theta^0$ can be arbitrary. A common variation known as Polyak averaging PoJu92, which often leads to improved convergence rates, adds an auxiliary iterate
which simply averages the iterates obtained through ((ref)). We then use $\bar{\theta}^m$ as our estimates of the dual variables. Whether or not Polyak averaging is used, the computational and storage cost of stochastic approximation is trivial.
There is an extensive theory on the convergence rates of SA; see, e.g., the seminal work of BaMo11,BaMo13. There are, however, some issues with applying these results to the sequence $\left\{\theta^m\right\}$ obtained using ((ref)). First, the theory often makes strong assumptions, such as strong convexity, on the objective function being optimized. When these assumptions are not met, one typically obtains a convergence rate for the value of the objective function ((ref)), not the decision variable. Second, and perhaps more importantly, the theory would require us to have a second sample from the distribution of $X$, independent of the sample $\left\{X^n\right\}^N_{n=1}$ used to compute $\hat{t}^N$. This would increase the burden of data collection.
For these reasons, we derive a performance guarantee in the following framework. First, since we will be working with SAA estimates of $t^*$ and $\psi^*$, we will abuse notation slightly to make the dependence of the choice on these parameters explicit, by writing $S\left(t,\psi,x\right) = \arg\max_k 1_{\left\{c\left(X,k\right)\leq t\right\}} + \psi_k$. Now suppose that, in ((ref)), we generate $X^{m+1}$ not by sampling from the density $f$, but by bootstrapping from the original sample $\left\{X^n\right\}^N_{n=1}$. The set $S\left(\psi^*,X^{m+1}\right)$ in ((ref)) will be replaced by $S\left(\hat{t}^N,\hat{\psi}^N,X^{m+1}\right)$. Then, ((ref)) will optimize, not ((ref)), but the sample analog
There are two advantages to doing this. First, we may work with one sample of size $N$ for all of our computations, without having to divide the total data-collection effort into two samples. Second, the computational cost of running SA becomes separated from $N$, because ((ref)) is very fast and can easily be run for much more than $N$ iterations. Thus, we may assume that the optimal solution of ((ref)) can be cheaply computed to arbitrary precision, and the only remaining task is to characterize the convergence of the optimal solution $\hat{\theta}^N$ of ((ref)), and the corresponding optimal transport plan, to the true optimal solutions.
First, we build on Theorem (ref) to characterize the probability that the estimated $\hat{t}^N$ leads to an error in the choice set (before tie-breaking). Note the distinction: Theorem (ref) studies the set of parameters characterizing the optimal assignment rule, whereas Proposition (ref) below concerns a specific assignment of a new $X$ not in the sample. In other words, we have moved from assessing the accuracy of the estimated parameters to evaluating their quality in making decisions for other $X$ that had not previously been observed. The probability that the choice set estimated for that $X$ under $\hat{t}^N$ and $\hat{\psi}^N$ does not match the true choice set (computed using $t^*$ and $\psi^*$) is shown to be of the canonical order $\mathcal{O}\left(N^{-\frac{1}{2}}\right)$.
Next, we study the accuracy of the solution of ((ref)). Since adding the same constant to each component $\theta$ will not change the value of the objective function, we may normalize $\theta_K = 0$ and think of both ((ref)) and ((ref)) as problems in $K-1$ dimensions. It is then fairly straightforward to show that ((ref)) has a unique solution $\theta^*$. The gradient of ((ref)) is given by
The next result establishes that, when ((ref)) is sufficiently small, we have an accurate approximation of $\theta^*$.
Note that $\hat{\theta}$ in ((ref)) need not be deterministic. In fact, if $\hat{\theta} = \hat{\theta}^N$, i.e., we have the exact optimal solution of ((ref)), the condition holds trivially since the left-hand side of ((ref)) is equal to zero. We allow some small amount of error, the idea being that SA can always be run long enough to make ((ref)) hold.
Our final result bounds the error of the choice probabilities computed from $\hat{t}^N$, $\hat{\psi}^N$, and $\hat{\theta}^N$ (or a close approximation), again for a new data point $X$. This completes our analysis, as the choice probabilities are precisely what characterizes the optimal transport plan for an arbitrary $X$. We see that the canonical $\mathcal{O}\left(N^{-\frac{1}{2}}\right)$ rate is obtained.
Geographical partitioning problems are a major application area of semidiscrete optimal transport CaCaDe16,CaPeRy24. In such problems, we consider a geographical region $\mathcal{X}\subseteq\mathbb{R}^2$ in which there are $K$ facilities with fixed locations $x_k \in\mathcal{X}$. The random variable $X$ represents the location of a randomly selected customer, with $f$ being the population density. The random variable $Y$ takes values in $\left\{1,...,K\right\}$ and represents a choice of facility. The $k$th facility handles a proportion $p_k$ of the total population.
The classical formulation of this problem chooses $\pi$ to minimize the expected cost $\mathbb{E}_\pi\left(c\left(X,Y\right)\right)$, where typically $c\left(x,k\right) = \|x-x_k\|_2$. Other reasonably common cost variants include squared Euclidean distance Oh07 and $1$-norms MaHoSoLe18. It is well-known that the optimal assignment rule, given a random location $X$, sets $Y\left(X\right) = \arg\min_j c\left(X,j\right) - g^*_j$ for a certain specific choice of $g^*\in\mathbb{R}^K$, where the argmin is a.s. unique under mild regularity conditions. Thus, the $k$th facility serves all customers in the set $A_k = \left\{x:Y\left(x\right) = k\right\}$, and the sets $A_1,...,A_K$ comprise a partition of $\mathcal{X}$. This partitioning structure is known in the literature as an “additively weighted Voronoi diagram”. The weights $g^*$ can be easily computed (for virtually any cost function) using stochastic approximation GeCuPeBa16.
In the following, we consider an illustrative example of the partitioning problem with $\mathcal{X} = \left[0,1\right]^2$ and $f$ being the uniform density on $\mathcal{X}$. There are six facilities, and the inputs $p$ and $x_k$ are pre-generated randomly. In Figure (ref), we take $c\left(x,k\right) = \|x-x_k\|_2$ and use our approach to solve $\inf_\pi \mathbb{Q}^{\alpha}_\pi\left(c\left(X,Y\right)\right)$ for various choices of $\alpha$. For reference, we also include the additively weighted Voronoi diagram (i.e., the partition that minimizes expected cost) in Figure (ref). The Voronoi diagram shows six disjoint regions whose areas correspond to the inputs $p_k$. Starting with the upper-leftmost facility and moving clockwise, these areas are $0.0522$ (green), $0.2140$ (orange), $0.1864$ (cyan), $0.1091$ (yellow), $0.2296$ (purple), and $0.2087$ (blue).
Figure (ref) shows the optimal division of the region when the objective is to minimize the median travel distance (that is, $\alpha = 0.5$). At a high level, the geometric structure of the optimal assignment is the same for all values of $\alpha$, but let us focus on the median to make the discussion more concrete. Using a sample size $N = 10^4$ for SAA, and the CVX package cvx with the Mosek solver, we can compute an estimate $t^* \approx 0.2099$. For each facility $k$, we draw a circle $E_k = \left\{x:\|x-x_k\|_2 \leq t^*\right\}$ of radius $t^*$ centered at $x_k$. The optimal $\psi^*_k$ in this example are equal to $1$ for all $k$ except the upper-leftmost (green) facility, which has $\psi^*_k = 0$. The main reason for this difference is that $p_k$ is very small for the green facility, requiring us to select it less frequently than the others.
For every $x$, let $B\left(x\right) = \left\{k:x\in E_k\right\}$. We assign $x$ to a facility as follows:
Note, however, that we still need to run the algorithms in Section (ref) even if we know this general structure, because that is the only way to obtain the correct $t^*$, $\psi^*$, and $\theta^*$, or a close approximation thereof.
We see the same structure in Figures (ref)-(ref), where $\alpha$ is set to $0.25$ and $0.95$, respectively. Although one typically prefers to optimize upper quantiles in cost minimization, we include $\alpha = 0.25$ purely for visual insight. The main difference between these two cases and $\alpha = 0.5$ is that the radius $t^*$ of each circle (i.e., the optimal quantile) is smaller in the first case (relative to $\alpha = 0.5$) and larger in the second case. For $\alpha = 0.25$, we have $\psi^*_k = 1$ for all $k$ and so the green circle uses only one color. For $\alpha = 0.95$, we again have $\psi^*_k = 1$ for all facilities except the green one, which has $\psi^*_k = 0$. A consequence is that the locations close to the orange facility are assigned randomly, because they exist in the overlap between the orange and cyan circles, but the locations close to the green facility are deterministically assigned to the orange one, because they exist only in the overlap between orange and green.
Figure (ref) shows the empirical densities of $c\left(X,Y\right)$ when $Y$ is chosen according to the solution of quantile optimization. We compare $\alpha = 0.5$ (median-minimizing), $\alpha = 0.95$, and also include the mean-minimizing distribution for reference. We omit $\alpha = 0.25$ from the comparison for brevity. The mean, median, and $95$th percentile of $c\left(X,Y\right)$ under each distribution are summarized in Table (ref). In Figure (ref), we see that the median-minimizing solution prefers to reduce larger costs without much regard for how far out they are in the tails; thus, there are quite a few extreme values, but their incidence is low. On the other hand, in Figure (ref), the incidence of very high costs is very sharply reduced, at the expense of the median.
It is also interesting to see how the visualization of the optimal assignment changes depending on the cost function. Figure (ref) presents a second example with different $p$ and $x_k$ as well as a different $c\left(x,k\right) = \|x-x_k\|_1$. The $1$-norm changes the shape of the regions in the Voronoi diagram (Figure (ref)). However, when the objective is to optimize the median of the $1$-norm, the high-level structure is the same as in the previous example, with the only difference being that the circles are now diamonds. Once more, the lower-leftmost facility is the only one with $\psi^*_k = 0$, because its low $p_k$ value requires us to select it less frequently. We do not show results for other quantiles here, since they generally repeat the insights from the previous example, but it is certainly possible to optimize them if desired.
To our knowledge, the geometric structure in Figures (ref)-(ref) has never before been observed in the literature on either optimal transport or computational geometry. These examples show that, in the setting of quantile optimization, not only can the optimal assignments be characterized theoretically, they are also computationally tractable using our methods.
To our knowledge, this is the first paper to consider quantile optimization in any form of optimal transport. The semidiscrete setting is a natural first choice to develop these ideas because it enables cleaner theoretical characterization and efficient computation. While some of the structure in our paper may be generalizable, it is highly unlikely that optimal plans would be computable under general continuous distributions. Even under the classical objective of minimizing expected cost, it is virtually impossible to obtain computationally tractable plans without resorting to discretization. Therefore, from a practical point of view, the semidiscrete setting may in fact be the most relevant to study.
Even in the semidiscrete case, quantile optimization has shown itself to be more complex than mean minimization because the optimal plan does not have the Monge property. In other words, given $X$, the optimal choice of $Y$ is not deterministic, but may require randomization via a tie-breaking rule. Such rules are quite under-studied in the optimal transport literature, where one usually works with assumptions under which Kantorovich-Monge equivalence holds. We believe that our approach to characterizing and computing the tie-breaking rule is one of the most novel aspects of our work. This aspect may be useful outside the quantile optimization setting, in any situation where unusual or non-standard cost functions preclude the existence of an optimal deterministic map. We hope that our work will lead to increased interest in such scenarios.