Entropic Regularization

Due the convex structure of the Kantorovich problem, it is natural to turn to methods that, upon discretizing, exploit the fact that we aim to solve a linear programming problem. Upon replacing our measures $\mu ,\nu$ with $\mu =\sum _{i=1}^{N}a_{i}\delta _{x_{i}}$ and $\nu =\sum _{j=1}^{M}b_{j}d_{y_{j}}$ for some appropriately determined $a_{i},b_{j}$ and fixed $N,M.$ The resulting minimization problem becomes:

$\min \left(\sum _{i,j}c_{ij}\alpha _{ij}:\alpha _{ij}\geq 0,\sum _{i}\alpha _{ij}=b_{j},\sum _{j}\alpha _{ij}=a_{i}\right)$

which we may approach using classical methods that exploit the structure of the resulting linear system.

The Entropic Regularization

In studying a slight variant to this problem, for fixed $\varepsilon >0$ we may now instead focus on tackling

$\min \left(\sum _{i,j}(c_{ij}\alpha _{ij}+\varepsilon \alpha _{ij}\log(\alpha _{ij})):\alpha _{ij}\geq 0,\sum _{i}\alpha _{ij}=b_{j},\sum _{j}\alpha _{ij}=a_{i}\right)$

which converges as $\varepsilon \rightarrow 0.$ Rewriting the objective function at hand by using

$c_{ij}\alpha _{ij}+\varepsilon \alpha _{ij}\log(\alpha _{ij})=\varepsilon \alpha _{ij}\log({\frac {\alpha _{ij}}{\eta _{ij}}}),\ \ \ \sum _{i,j}(c_{ij}\alpha _{ij}+\varepsilon \alpha _{ij}\log(\alpha _{ij}))=\varepsilon {\textbf {KL}}(\alpha |\eta ),$

where $\eta _{ij}=e^{-c_{ij}/\varepsilon }$ and \textbf{KL} denotes the ${\textit {Kullback-Leibler}}$ divergence (based on the relative entropy).

Thus, from here, this minimization problem reads as a projection from this divergence, of the point $\eta \in \mathbb {R} ^{N\times N}$ on the set of constraints

$C=C^{x}\cap C^{y}=\{\alpha \in \mathbb {R} ^{N\times N}:\sum _{j}\alpha _{ij}=a_{i}\}\bigcap \{\alpha \in \mathbb {R} ^{N\times N}:\sum _{i}\alpha _{ij}=b_{j}\}$

In particular, this Kullback-Leibler divergence has the property that the projection onto the intersection of two linear subspaces such as $C^{x},C^{y}$ can be achieved by alternate projections. This is done by defining $\alpha ^{2k+1}=Proj_{C^{x}}(\alpha ^{2k})$ and $\alpha ^{2k+2}=Proj_{C^{y}}(\alpha ^{2k+1})$ with $\alpha ^{0}=\eta$ and looking at the limit as we take $k\rightarrow \infty$ .

As such, the problem now becomes of projecting onto one of the linear subspaces, which can be done. If we are given a $\beta$ , we wish to solve

$\min \left(\sum _{i,j}\alpha _{ij}\log \left({\frac {\alpha _{ij}}{\beta _{ij}}}\right):\sum _{j}\alpha _{ij}=a_{i}\right)$

which we can solve an approximated problem of with the following algorithm.

Algorithm

Start from $\alpha ^{0}=\eta$
define $\alpha _{ij}^{2k+1}=p_{i}\alpha _{ij}^{2k}$ with $p_{i}=a_{i}/\sum _{j}\alpha _{ij}^{2k}$
then define $\alpha _{ij}^{2k+2}=q_{j}\alpha _{ij}^{2k+1}$ with $q_{j}=b_{j}/\sum _{i}\alpha _{ij}^{2k+1}$

The limit we observe as $k\rightarrow \infty$ provides us with the solution to the regularized problem. This idea is called the iterative proportional fitting procedure ( ${\textbf {IPFP}}$ ). In utilizing this entropic regularization, this algorithm lends itself to the following advantages:

we may compute the projections explicitly at every step the moment we need only project to $C^{x}$ or $C^{y}$
it is quite memory-efficient since the cumulating product of the coefficients $p_{i},q_{j}$ appear at each step, thus allowing us to only store $2N$ as opposed to $N^{2}$
it can be easily parallelized given the separability of the structure at any given step.