Semidiscrete Optimal Transport

Semidiscrete optimal transport refers to situations in optimal transport where two input measures are considered, and one measure is a discrete measure and the other one is absolutely continuous with respect to Lebesgue measure. Hence, because only one of the two measures is discrete, we arrive at the appropriate name "semidiscrete."

Voronoi Cells. These partitions of the plane are Voronoi cells, which are used in the semidiscrete dual optimal transport problem. The cell structures are determined by the cost function used, here being the Euclidean metric. The bold black lines represent the cell structure.

Formulation of the semidiscrete dual problem

In particular, we will examine semidiscrete optimal transport in the case of the dual problem. The general dual problem for continuous measures can be stated as

$$\max _{(\psi ,\varphi )\in R(c)}{\Big \{}\int _{X}\psi d\mu +\int _{Y}\varphi d\nu {\Big \}}$$

where $\mu ,\nu$ denote probability measures on domains $X,Y$ respectively, and $c(x,y)$ is a cost function defined over $X\times Y$. $R(c)$ denotes the set of possible dual potentials, and the condition $\varphi (x)+\psi (y)\leq c(x,y)$ is satisfied. It should also be noted that $\mu$ has a density such that $\mu =f(x)dx$. Now, we would like to extend this notion of the dual problem to the semidiscrete case. Such a case can be reformulated as

$$\max _{\varphi \in \mathbb {R} ^{m}}{\Big \{}{\mathcal {E}}(\varphi )=\int _{X}\varphi ^{c}d\mu +\sum _{j}\varphi _{j}b_{j}{\Big \}}.$$

Aside from using a discrete measure in place of what was originally a continuous one, there are a few other notable distinctions within this reformulation. The first is that $\varphi ^{c}$ denotes the c-transform of $\varphi$. The c-transform can be defined as $\varphi ^{c}(x):=\min _{j}\{c(x,y_{j})-\varphi _{j}\}$. $\varphi _{j}$ is used to denote $\varphi (y_{j})$. Furthermore, we note that our original measure $\nu$ is a sum of Dirac masses evaluated at locations $y_{j}$ with weights $b_{j}$, i.e., $\nu =\sum _{j=1}^{N}b_{j}\delta _{y_{j}}$.

Voronoi cells to find weights

Now, we will establish the notion of Voronoi cells. The Voronoi cells refer to a special subset of $X$, and the reason we are interested in such a subset is because we can use the Voronoi cells to find the regions that are sent to each $y_{j}$. In particular, if we denote the set of Voronoi cells as $V_{\varphi }(j)$, we can find our values of $\varphi _{j}$ using the fact $b_{j}=\int _{V_{\varphi }(j)}f(x)dx$. Recall that $f(x)$ refers to a density of the measure $\mu$, i.e., $\mu =f(x)dx$. We define the Voronoi cells with

$$V_{\varphi }(j)={\Big \{}x\in X:\forall j'\neq j,{\frac {1}{2}}|x-y_{j}|^{2}-\varphi _{j}\leq {\frac {1}{2}}|x-y_{j'}|^{2}-\varphi _{j'}{\Big \}}.$$

We use the specific cost function $c(x,y)={\frac {1}{2}}|x-y|^{2}$ here. This is a special case and we may generalize to other cost functions if we desire. When we have this special case, the decomposition of our space $X$ is known as a "power diagram." Using our power diagram as a domain of integration, we can successfully find the weights $b_{j}$.

Finding the weights via the gradient

Finding the weights via the above method is equivalent to maximizing ${\mathcal {E}}(\varphi )$, and we may do this by taking the partial derivatives of this function with respect to $\varphi _{j}$. This is the same as taking the gradient of ${\mathcal {E}}(\varphi )$. In partial derivative form, we have

$${\frac {\partial {\mathcal {E}}}{\partial \varphi _{j}}}=-\int _{V_{\varphi }(j)}f(x)dx+b_{j}$$

,

and in gradient form, we have

$$\nabla {\mathcal {E}}(\varphi )_{j}=-\int _{V_{\varphi }(j)}f(x)dx+b_{j}.$$

Since $\nabla {\mathcal {E}}(\varphi )_{j}=0$ when it attains a maximum, we have the relation between the weights and the measure density that we established in the previous section. Note that the maximum is taken and not the minimum because our function ${\mathcal {E}}(\varphi )$ is a concave function. The discrete summation contained within this function is linear, but an infimum of a linear function is evaluated for the integration part, making the overall function concave.

Algorithm discussion

We may search for a maximum of value of ${\mathcal {E}}(\varphi )$ by means of certain algorithms, the first being gradient ascent. Whether or not such an algorithm is capable of being implemented effectively is contingent on the ability to find our power diagram in a practical way. Certain computational geometry algorithms allow the cells to be found efficiently. A second suitable algorithm is Newton's method. Using Newton's method to find the zeros of ${\frac {\partial {\mathcal {E}}}{\partial \varphi _{j}}}$, one must compute the second derivative of ${\mathcal {E}}$, as well as justify certain constraints are met, such as bounded eigenvalues of the Hessian.