Moment Constrained Optimal Transport for Control Applications

## Abstract

Abstract Optimal transport is now a standard tool for solving many problems in statistics and machine learning. The optimal “transport of probability measures” is also a recurring theme in stochastic control and distributed control, where in the latter application the probability measure corresponds to an empirical distribution associated with a large collection of distributed agents, subject to local and global control. The goal of this paper is to make precise these connections, which inspires new relaxations of optimal transport for application in new and traditional domains. The proposed relaxation replaces a target measure with a “moment class”: a set of probability measures defined by generalized moment constraints. This is motivated by applications to control, outlier detection, and to address computational complexity. The main conclusions are (i) A characterization of the solution is obtained, similar to Kantorovich duality, in which one of the dual functions in the classical theory is replaced by a linear combination of the features defining the generalized moments. Hence the dimension of the optimization problem coincides with the number of constraints, even with an uncountable state space; (ii) By introducing regularization in the form of relative entropy, the solution can be interpreted as replacing a maximum with a soft-max in the dual; (iii) In applications such as control for which it is not known a-priori if the moment class is non-empty, a relaxation is proposed whose solution admits a similar characterization; (iv) The gradient of the dual function can be expressed in terms of the expectation of the features under a tilted probability measure, which motivates Monte-Carlo techniques for computation.

## 1 Introduction

#### Background and goals

Optimal Transport (OT) theory first emerged in the 18th century, and more recently has become a significant tool in the machine learning toolbox Villani (); Peyré et al. (). The goal is simply described: given two random variables $X$ and $Y$ on a common state space $\mathcal{X}$, with given distributions $\mu_{1}$ and $\mu_{2}$, find a joint probability measure $\gamma^{*}$ for the pair $(X,Y)$ that preserves the marginals, and solves some optimality criterion that attempts to couple the two random variables.

The Monge-Kantorovich formulation is expressed as follows. Let
$\Upgamma(\mu_{1},\mu_{2})$ set of probability measures on ${\cal B}(\mathcal{X}\times\mathcal{X})$ with specified marginals: $\gamma\in\Upgamma(\mu_{1},\mu_{2})$ if and only if $\gamma_{1}=\mu_{1}$ and $\gamma_{2}=\mu_{2}$.
Given a continuous function $c\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}_{+}$, the optimal transport problem is formulated as the minimum

$$ $\gamma^{*}\in\mathop{\rm arg\,min}_{\gamma\in\Upgamma(\mu_{1},\mu_{2})}\langle\gamma,c\rangle\,,$ (1) $$

in which the inner product denotes expectation:
$\langle\gamma,c\rangle\mathbin{:=}{\sf E}[c(X,Y)]$ when $(X,Y)\sim\gamma$.
In this paper we take $\mathcal{X}$ a closed subset of $\mathbb{R}^{N}$ with $N\geq 1$.
It is always assumed that $c(x,x)=0$ for each $x$;
a typical choice is $c(x,y)={\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|x-y\|^{2}$.

Letting $\gamma^{*}$ denote an optimizer, we can write

$$ $\gamma^{*}(dx,dy)=\mu_{1}(dx)T^{*}(x,dy)$ (2) $$

where $T^{*}$ is a transition kernel. The optimization problem is a linear program, which is why we can expect degeneracy: there is a mapping $\uptau\colon\mathcal{X}\to\mathcal{X}$ such that $Y=\uptau(X)$ with probability one under $\gamma^{*}$.

It will be useful to recall the convex dual of ([1](#S1.E1)) since it is similar to those constructed to obtain solutions to the optimization problems posed here.
Letting $\psi_{i}\colon\mathcal{X}\to\mathbb{R}$, $i=1,2$, denote continuous functions representing Lagrange multipliers for the two marginal constraints, the dual becomes the maximum over all such functions:

$$ $\displaystyle\langle\gamma^{*},c\rangle$ $\displaystyle=\max_{\psi_{1},\psi_{2}}\bigl{\{}\langle\mu_{1},\psi_{1}\rangle+\langle\mu_{2},\psi_{2}\rangle:\psi_{1}(x)+\psi_{2}(y)\leq c(x,y)\ \textit{for all $x,y$}\bigr{\}}$ (3) $$

Letting $\{\psi_{i}^{*}\}$ denote optimizers, the solution to the primal $\gamma^{*}$ is supported on the set $\{(x,y):\psi_{1}^{*}(x)+\psi_{2}^{*}(y)=c(x,y)\}$, which is an instance of complementary slackness. Conditions for the existence of optimizers and justification of ([3](#S1.E3)) may be found in standard texts.

There are many other academic communities interested in transforming probability measures cheaply. Examples include the fully probabilistic control design of Kárný () and related linearly-solvable Markov decision framework introduced in Todorov ().
The area of Mean Field Games begins with a multi-objective control problem, but the final solution technique amounts to transporting a probability measure on a high dimensional space in such a way as to minimize some objective function. Similar to mean field games is ensemble control, with applications ranging from power systems to medicine Hochberg et al. (); Chertkov and Chernyak (); relaxations of this technique are found in Cammardella et al. (); Bušić and Meyn (). More examples may be found in the recent survey Garrabe and Russo ().

The goal of the research summarized here is to build bridges between these fields, and in particular:
1.  Obtain new techniques for efficiently approximating the solution to the OT problem, and 2.  Investigate how computational techniques from OT theory might improve computation of optimal control solutions.

#### Feature projected optimal transport

The relaxations of OT considered in this paper are related to concepts introduced by A. Markov in the 19th century, but left unpublished until the survey of Kreĭn ().
Given a family of continuous, real-valued functions $\{f_{i}:1\leq i\leq M\}$ on $\mathcal{X}$ and scalars $\{r_{i}:1\leq i\leq M\}$, consider the set of probability measures satisfying the generalized moment constraints:

$$ ${\cal P}_{f,r}=\{\mu\in{\cal B}(\mathcal{X}):\langle\mu,f_{i}\rangle=r_{i}\ :\ 1\leq i\leq M\}$ (4) $$

The respective canonical distributions are defined as
the extreme points of this subset of the simplex of probability measures on ${\cal B}(\mathcal{X})$. Subject to assumptions on $\mathcal{X}$ (such as compactness), any canonical distribution $\mu^{*}$ has at most $M+1$ points of support—see Kemperman ().

This theory lends itself naturally to function approximation; e.g., Micchelli and Rivlin ().
Applications to hypothesis testing can be found in Pandit and Meyn (); Unnikrishnan et al. () and to channel coding in
Lapidoth (); Abbe et al. () (the connection between canonical distributions and mismatched decoding was pointed out in Abbe et al. ()). Inspired by this history, in this paper we explore three relaxations of the OT problem ([1](#S1.E1)). The functions $\{f_{i}\}$ will be called features. In control applications the random vector $f(Y)$ might represent observed costs or rewards. In signal processing, features might be analogous to a basis as used in wavelet transforms, for which $\langle\mu,f_{i}\rangle$ represents the $i$ coefficient.

#### Problem 1:

Feature projected optimal transport (OT-FP):
With given $\mu_{1}$, features $\{f_{i}\}$, and scalars $\{r_{i}\}$,

$$ $\displaystyle d(\mu_{1},{\cal P}_{f,r})=\min\bigl{\{}\langle\gamma,c\rangle:\gamma\in\Upgamma(\mu_{1},\mu)\,,\ \mu\in{\cal P}_{f,r}\bigr{\}}$ (5) $$

This is a relaxation of ([1](#S1.E1))
if we are given $\mu_{2}$, and set $r_{i}=\langle\mu_{2},f_{i}\rangle$ for each $i$.

We also consider the regularized optimization problem with convex regularizer

$$ $\mathcal{C}(\gamma)=D(\gamma\|\mu_{1}\otimes\mu_{2})$ (6) $$

This is the functional used to define the Sinkhorn distance in Cuturi (), which is a major part of the algorithms introduced to efficiently approximate solutions to ([1](#S1.E1)).

In this paper only the first marginal is constrained: $\gamma\in\Upgamma(\mu_{1},\mu)$, with $\mu\in{\cal B}(\mathcal{X})$ regarded as a variable. The probability measure $\mu_{2}\in{\cal P}$ may be chosen based on intuition regarding the form of $\gamma_{2}^{*}$, or chosen for ease of computation. In applications to control proposed in Section [3](#S3) we take $\mu_{2}=\mu_{1}$, since we seek a solution for which $\gamma^{*}_{2}$ is not far from $\mu_{1}$.

#### Problem 2:

Feature projected optimal transport with regularization (OT-FPR):
With $\mu_{1}$, features $\{f_{i}\}$, scalars $\{r_{i}\}$, and $\varepsilon>0$ given,

$$ $\displaystyle d_{\varepsilon}(\mu_{1},{\cal P}_{f,r})=\min_{\mu,\gamma}\bigl{\{}\langle\gamma,c\rangle+\varepsilon\mathcal{C}(\gamma):\gamma\in\Upgamma(\mu_{1},\mu)\,,\ \mu\in{\cal P}_{f,r}\bigr{\}}$ (7) $$

In some applications, we may not know a-priori if the moment class ${\cal P}_{f,r}$ defined in ([4](#S1.E4)) is non-empty, or we may find that matching moments exactly is too costly. Let ${\cal R}\colon\mathbb{R}^{M}\to\mathbb{R}_{+}$ be a convex function that vanishes only at the origin, and assume its convex dual ${\cal R}^{*}\colon\mathbb{R}^{M}\to\mathbb{R}_{+}$ is everywhere finite. Recall this is defined as

$$ ${\cal R}^{*}(x)\mathbin{:=}\sup\{\delta^{\intercal}x-{\cal R}(\delta):\delta\in\mathbb{R}^{M}\}\,,\qquad x\in\mathbb{R}^{M}\,.$ (8) $$

For example, the convex dual of the quadratic ${\cal R}(\delta)=\tfrac{\kappa}{2}\|\delta\|^{2}$, with $\kappa>0$, is also quadratic:
${\cal R}^{*}(x)=\tfrac{1}{2\kappa}\|x\|^{2}$.

#### Problem 3:

Feature projected optimal transport with regularization and penalty (OT-FPRP): with $\mu_{1}$, features $\{f_{i}\}$, scalars $\{r_{i}\}$, $\varepsilon>0$, and ${\cal R}$ given,

$$ $\displaystyle d_{\varepsilon,{\cal R}}(\mu_{1},{\cal P}_{f,r})=\min_{\mu,\gamma,\delta}\bigl{\{}\langle\gamma,c\rangle+\varepsilon\mathcal{C}(\gamma)+{\cal R}(\delta)\ :\ \gamma\in\Upgamma(\mu_{1},\mu)\,,\ \delta_{i}=\langle\mu,\tilde{f}_{i}\rangle\,,\ 1\leq i\leq M\bigr{\}}$ (9) $$

with $\tilde{f}_{i}=f_{i}-r_{i}$.

#### Contributions

The dual of the optimization problem ([5](#S1.E5)) is similar to what is found in the classical OT theory, and leads to the following representation under mild conditions:

$$ $\displaystyle d(\mu_{1},{\cal P}_{f,r})=\max_{\psi,\lambda}\bigl{\{}\langle\mu_{1},\psi\rangle:\psi(x)+\lambda^{\intercal}\tilde{f}(y)\leq c(x,y)\ \textit{for all $x,y\in\mathcal{X}$}\bigr{\}}$ (10) $$

This optimization problem appears very similar to the classical dual ([3](#S1.E3)), with two important differences: in
([10](#S1.E10)) the objective function involves $\mu_{1}$ only (recall that $\gamma_{2}$ is not constrained in OT-FP), and the function $\psi_{2}$ appearing in ([3](#S1.E3)) is replaced with $\lambda^{\intercal}\tilde{f}$.

Regularization leads to an approximation that is much more easily computed: for both OT-FPR and OT-FPRP we find that a convex dual can be constructed of dimension $M$. The optimal transport kernel is determined by a parameter $\lambda^{*}\in\mathbb{R}^{M}$ solving the respective convex optimization problems:
for a concave function $\varphi^{*}\colon\mathbb{R}^{M}\to\mathbb{R}$,

$$ OT-FPR: $\displaystyle\displaystyle\lambda^{*}=\mathop{\rm arg\,max}\varphi^{*}(\lambda)$ (11a) OT-FPRP: $\displaystyle\displaystyle\lambda^{*}=\mathop{\rm arg\,max}\{\varphi^{*}(\lambda)-{\cal R}^{*}(-\lambda)\}$ (11b) $$

It is shown in Proposition [2.3](#S2.Thmtheorem3) that the solution to OT-FP can be approximated arbitrarily closely using OT-FPR by choosing $\varepsilon>0$ sufficiently small.
This is used to establish that there is no duality gap for OT-FP, which provides the proof of ([10](#S1.E10)).

The structure of $\varphi^{*}$ lends itself to Monte-Carlo methods. Each $\lambda\in\mathbb{R}^{M}$ determines a probability measure $\gamma^{\lambda}$.
The gradient of $\varphi^{*}$ at $\lambda$ coincides with the negative of the mean of $\tilde{f}$ under this probability measure, and its Hessian is a negative scaling of a conditional autocovariance of $\tilde{f}(Y)$ under $\gamma^{\lambda}$.

These results suggest many extensions that can be addressed using similar analysis.
If there is motivation to relax the constraint $\gamma_{1}=\mu_{1}$, we can take a pair of moment classes $\{{\cal P}_{f^{i},r^{i}}:i=1,2\}$ and consider the resulting convex programs for any of Problems 1–3. OT-FPRP becomes

$$ $\displaystyle\!\!\!\!d_{\varepsilon,{\cal R}}({\cal P}_{f^{1},r^{1}},{\cal P}_{f^{2},r^{2}})\mathbin{:=}\min_{\mu_{1},\mu_{2},\gamma}\Bigl{\{}\langle\gamma,c\rangle+\varepsilon\mathcal{C}(\gamma)+{\cal R}(\delta)\ :\ \gamma\in\Upgamma(\mu_{1},\mu_{2})\,,\ \mu_{i}\in{\cal P}_{f^{i},r^{i}}\Bigr{\}}$ (12) $$

Duality theory for this problem follows the same path as in analysis of OT-FPRP; theory is simplified substantially because there are only a finite number of constraints.

When $\varepsilon=0$ and ${\cal R}$ is eliminated we omit the subscript in $d_{\varepsilon,{\cal R}}({\cal P}_{f^{1},r^{1}},{\cal P}_{f^{2},r^{2}})$.
In this special case, the optimization problem ([12](#S1.E12))
fits squarely in the theory of Markov’s canonical distributions: a linear program over the space of probability measures, with a finite number of constraints. If the state space is compact and all functions are continuous, the theory in this paper extends easily to obtain

$$ $\displaystyle d({\cal P}_{f^{1},r^{1}},{\cal P}_{f^{2},r^{2}})$ $\displaystyle=\max_{\lambda^{1},\lambda^{2}}\Bigl{(}\min_{x,y}\bigl{\{}c(x,y)-{\lambda^{1}}^{\intercal}\tilde{f}^{1}(x)-{\lambda^{2}}^{\intercal}\tilde{f}^{2}(y)\bigr{\}}\Bigr{)}\,,\quad\tilde{f}^{i}=f^{i}-r^{i}$ (13) $$

#### Literature review

Standard references for optimal transport include Villani (); Ollivier et al. (), and much more history may be found in the recent survey Peyré and Cuturi ().

It was already made clear that the framework for the optimization problems considered in this paper draws from the 2013 paper of Cuturi (), which was greatly expanded soon after publication. More recent analysis of Cuturi’s approach may be found in
Peyré and Cuturi (); Eckstein and Nutz (); Nutz and Wiesel ().

Relaxations of marginals in optimal transport was proposed in Li and Lin (). Their OT-RMC is similar to OT-FPRP in the form ([12](#S1.E12)) in which the set $\mathcal{X}=\{x^{i}:1\leq i\leq L\}$ is finite, $f_{i}(x)=\mathbb{I}\{x=x^{i}\}$ for each $1\leq i\leq L$, and with penalty function
${\cal R}$ a weighted $\ell_{1}$ norm. Theory is illustrated with application to image processing, which seems an ideal area of application for the methods introduced here.

Relaxations of marginals is proposed in Balaji et al. (); Le et al. () to improve numerical performance. The motivation and analysis is very different from what is reported here; it will be worthwhile to see if the approaches can be combined to obtain more efficient algorithms to solve OT-FPRP.

The applications to stochastic control in Section [3](#S3) apply OT-FPR in the form

$$ $\min\Bigl{\{}\sum_{i=1}^{M}{\sf E}[(X_{i}-Y_{i})^{2}]+\varepsilon\mathcal{C}(\gamma)\ :\ {\sf E}[{\cal U}(Y_{i})]=r_{i}\ \textit{for each $i$}\Bigr{\}}$ $$

This is designed for distributed control applications in which ${\sf E}[{\cal U}(Y_{i})]=r_{i}$ is interpreted as perfect tracking of the ensemble of agents, and the quadratic cost ${\sf E}[\|X-Y\|^{2}]$ represents a penalty for deviation from nominal behavior.

Chertkov and Chernyak () consider a similar objective, with $\varepsilon\to\infty$, and a different regularizer:

$$ $\min_{\mu}D(\mu\|\mu_{1})\,,\ \ \text{s.t.}\ {\sf E}_{\mu}[{\cal U}(Y_{i})]=r_{i}\ \textit{for each $i$}$ $$

This approach is similar to (and inspired by) Todorov ().
The regularizer is smaller:

$$ $D(\mu\|\mu_{1})=D(\mu_{1}T\|\mu_{1})\leq\int_{x\in\mathcal{X}}\mu_{1}(dx)D(T(x,\,\cdot\,)\|\mu_{1})=\mathcal{C}(\gamma)\,,\quad\gamma=\mu_{1}\otimes T\in\Upgamma(\mu_{1},\mu)$ $$

The algorithms in Cammardella et al. (); Bušić and Meyn () address similar control objectives, but the optimality criterion is more closely related to OT-FPRP using a quadratic penalty.

## 2 Optimal Transport And Relaxations

Solutions to each version of OT-FP will involve a family of transition kernels $\{T^{\lambda}:\lambda\in\mathbb{R}^{M}\}$. For each $\lambda$ we denote $\gamma^{\lambda}=\mu_{1}\odot T^{\lambda}$, and let $\mu^{\lambda}=\gamma^{\lambda}_{2}$ denote the second marginal:

$$ $\mu^{\lambda}(A)\mathbin{:=}\int\mu_{1}(dx)T^{\lambda}(x,A)\,,\qquad A\in{\cal B}(\mathcal{X})$ $$

For measurable $g\colon\mathcal{X}\to\mathbb{R}$ we adopt the operator-theoretic notation,

$$ $T^{\lambda}g\,(x)\mathbin{:=}\int T^{\lambda}(x,dy)g(y)\,,\qquad x\in\mathcal{X}$ $$

It is sometimes useful to use probabilistic notation: for measurable $g\colon\mathcal{X}\to\mathbb{R}$ and $h\colon{\sf X}\times{\sf X}\to\mathbb{R}$,

$$ ${\sf E}^{\lambda}[g(Y)\mid X=x]\mathbin{:=}T^{\lambda}g\,(x)\,,\ \ x\in\mathcal{X}\,,\qquad{\sf E}^{\lambda}[h(X,Y)]\mathbin{:=}\langle\gamma^{\lambda},h\rangle$ (14) $$

Proofs of all of the technical results of the paper are contained in the supplementary material.
Assumptions are required for existence of optimizers and desirable properties of the dual for each optimization problem.

Assumptions:
(A1)
$c\colon\mathcal{X}\times\mathcal{X}\to\mathbb{R}_{+}$ and $f\colon\mathcal{X}\to\mathbb{R}^{M}$ are continuous, and there is an open neighborhood $N(r)\subset\mathbb{R}^{M}$ containing $r$ such that
${\cal P}_{f,r^{\prime}}$ is non-empty for all $r^{\prime}\in N(r)$.

(A2)
$\mu_{1}$ and $\mu_{2}$ have compact support, and Problem 2 is feasible under perturbations: for any $r^{\prime}\in N(r)$,
there is $\gamma$ and $\mu$ satisfying $\mu\in{\cal P}_{f,r^{\prime}}$ and
$\gamma\in\Upgamma(\mu_{1},\mu)$.

(A3)
$\Sigma^{0}\mathbin{:=}\text{\rm Cov}\,(Y)$ is positive definite when $Y\sim\mu_{2}$.

### 2.1 Dual for OT-FP

Characterization of a solution to Problem 1 is based on a Lagrangian relaxation.
Introduce two classes of Lagrange multipliers for ([5](#S1.E5)):
$\psi$ is for the first marginal constraint, a real-valued measurable function on $\mathcal{X}$, and $\lambda\in\mathbb{R}^{M}$ for the moment constraints. The dual functional is defined as the infimum,

$$ $\displaystyle\varphi^{*}(\psi,\lambda)$ $\displaystyle\mathbin{:=}\inf_{\gamma}\ \ \langle\gamma,c\rangle-\langle\gamma_{1}-\mu_{1},\psi\rangle-\langle\gamma_{2},\lambda^{\intercal}\tilde{f}\rangle$ (15) $\displaystyle=\langle\mu_{1},\psi\rangle+\inf_{x,y}\{c(x,y)-\psi(x)-\lambda^{\intercal}\tilde{f}(y)\}$ $$

The convex dual of ([5](#S1.E5)) is defined to be the supremum of $\varphi^{*}(\psi,\lambda)$ over all $\psi$ and $\lambda$. The dual optimization problem admits a familiar representation.
Compactness is assumed in Proposition [2.1](#S2.Thmtheorem1) (ii), as in prior work such as Kemperman ()
concerning canonical distributions.

###### Proposition 2.1

Suppose that (A1) and (A2) hold. Then,

- (i)
With $\varphi^{*}$ defined in ([15](#S2.E15)), the dual convex program admits the representation
$\displaystyle d^{*}\mathbin{:=}\sup_{\psi,\lambda}\varphi^{*}(\psi,\lambda)=\sup_{\psi,\lambda}\bigl{\{}\langle\mu_{1},\psi\rangle:\psi(x)+\lambda^{\intercal}\tilde{f}(y)\leq c(x,y)\ \ \textit{for all $x,y$}\bigr{\}}$
(16)
On replacing $\psi$ with $\psi^{\lambda}(x)\mathbin{:=}\inf_{y}\{c(x,y)-\lambda^{\intercal}\tilde{f}(y)\}$ we obtain the equivalent max-min problem
$d^{*}=\sup_{\lambda}\int\inf_{y}[c(x,y)-\lambda^{\intercal}\tilde{f}(y)]\mu_{1}(dx)$
(17)
- (ii)
Suppose in addition the set $\mathcal{X}$ is compact. Then the supremum in ([16](#S2.E16))
is achieved, and there is no duality gap: for a vector $\lambda^{*}\in\mathbb{R}^{M}$,
$d(\mu_{1},{\cal P}_{f,r})=d^{*}=\int\min_{y}\{c(x,y)-{\lambda^{*}}^{\intercal}\tilde{f}(y)\}\mu_{1}(dx)$

We present here the proof of part (i). The proof of (ii) is contained in the supplementary material, based on approximation with solutions to OT-FPR. A summary of the approach is contained in Proposition [2.3](#S2.Thmtheorem3).

###### Proposition 2.1

#### Proof of Proposition 2.1 (i)

The dual function is invariant under a constant shift in $\psi$, so we may assume the infimum in ([15](#S2.E15))
is precisely zero by adding a constant to $\psi$.
This gives

$$ $\displaystyle\max_{\psi,\lambda}\varphi^{*}(\psi,\lambda)$ $\displaystyle=\max_{\psi,\lambda}\bigl{\{}\langle\mu_{1},\psi\rangle:\inf_{x,y}[c(x,y)-\psi(x)-\lambda^{\intercal}\tilde{f}(y)]=0\bigr{\}}$ $$

The value of the maximum is unchanged with the equality constraint replaced by inequality: $\inf_{x,y}\{c(x,y)-\psi(x)-\lambda^{\intercal}\tilde{f}(y)\}\geq 0$, which then yields the representation ([16](#S2.E16)).
$\sqcap$$\sqcup$

Once we solve ([16](#S2.E16)),
we obtain $\gamma^{*}$ through complementary slackness:

$$ $0=\sum_{x,y}\gamma^{*}(x,y)\{\psi^{*}(x)+{\lambda^{*}}^{\intercal}f(y)-c(x,y)\}$ $$

which means that $\gamma^{*}$ is supported on the set $\{(x,y):{\lambda^{*}}^{\intercal}f(y)+\psi^{*}(x)=c(x,y)\}$.

### 2.2 Regularization

Recall that the functional $\mathcal{C}(\gamma)\mathbin{:=}D(\gamma\|\mu_{1}\otimes\mu_{2})$ is used to define the Sinkhorn distance in Cuturi (),
and coincides with mutual information when the marginals of $\gamma$ agree with the given probability measures $\mu_{1}$ and $\mu_{2}$. In the present paper, the marginal $\mu_{2}$ is a design parameter. The general examples in Section [3](#S3)
illustrate how it might be selected. In particular, it is explained why the special case $\mu_{2}=\mu_{1}$ is useful in applications to stochastic control.

#### OT-FPR geometry and duality

A close cousin to OT-FPR uses $\mathcal{C}$ as a constraint rather than penalty, which is the approach taken in Cuturi (). Consider for fixed $\delta>0$,

$$ $\displaystyle d_{\delta}^{c}(\mu_{1},{\cal P}_{f,r})=\min$ $\displaystyle\ \ \langle\gamma,c\rangle\,,\quad\text{s.t.}\ \ \gamma\in\Upgamma(\mu_{1},\mu)\,,\ \mu\in{\cal P}_{f,r}\,,\ \mathcal{C}(\gamma)\leq\delta$ (18) $$

The parameter $\varepsilon>0$ in ([7](#S1.E7)) may be regarded as a Lagrange multiplier corresponding to the constraint $\mathcal{C}(\gamma)\leq\delta$. Under general conditions there is $\delta(\varepsilon)$ such that the optimizers of
([18](#S2.E18)) and ([7](#S1.E7)) coincide.

In considering the dual of ([7](#S1.E7)) we choose a relaxation of the moment constraints only: letting $\lambda\in\mathbb{R}^{M}$ denote the Lagrange multiplier as before,

$$ $\varphi^{*}(\lambda)\mathbin{:=}\inf_{\gamma}\{\langle\gamma,c\rangle+\varepsilon\mathcal{C}(\gamma)-\langle\gamma_{2},\lambda^{\intercal}\tilde{f}\rangle:\gamma_{1}=\mu_{1}\}$ (19) $$

The convex dual of OT-FPR is by definition the supremum of the concave function $\varphi^{*}$ over all $\lambda\in\mathbb{R}^{M}$. The optimizer, when it exists, is denoted $\gamma^{\lambda}$.

Figure: Figure 1: Dual geometry for OT-FPR
Refer to caption: /html/2208.01958/assets/x1.png

Introducing the notation

$$ $\ell_{0}^{\lambda}(x,y)=\lambda^{\intercal}\tilde{f}(y)-c(x,y)\,,\quad x,y\in\mathcal{X}$ (20) $$

the dual function may be expressed

$$ $\varphi^{*}(\lambda)=-\max_{\gamma}\{\langle\gamma,\ell_{0}^{\lambda}\rangle-\varepsilon\mathcal{C}(\gamma):\gamma_{1}=\mu_{1}\}$ $$

The dual of ([18](#S2.E18)) with $d=d(\varepsilon)$ yields better geometric insight. If the maximum above exists, then the maximizer $\gamma^{\lambda}$ solves

$$ $\gamma^{\lambda}\in\mathop{\rm arg\,max}\{\langle\gamma,\ell_{0}^{\lambda}\rangle:\mathcal{C}(\gamma)\leq\delta\,,\ \gamma_{1}=\mu_{1}\}$ $$

The structure of the solution is illustrated in Fig. [1](#S2.F1). The convex region shown is the set of all $\gamma$ for which $\gamma_{1}=\mu_{1}$ and $\mathcal{C}(\gamma)\leq\delta$ (recall that $\mathcal{C}(\gamma)=D(\gamma\|\mu_{1}\otimes\mu_{2})$). The optimizer $\gamma^{\lambda}$ lies on the intersection of this region and the hyperplane shown in the figure, indicated with a dashed line: $\{\gamma:\langle\gamma,\ell_{0}^{\lambda}\rangle=\langle\gamma^{\lambda},\ell_{0}^{\lambda}\rangle\}$. This value of $\lambda$ does not optimize $\varphi^{*}$ because the hyperplane is not the boundary of the
half-space shown in the figure.

For computation it is convenient to make a change of variables:
since $\gamma_{1}=\mu_{1}$ is constrained, the infimum is over all transition kernels:

$$ $\varphi^{*}(\lambda)\mathbin{:=}\inf_{T}\{-\langle\mu_{1}\odot T,\ell_{0}^{\lambda}\rangle+\varepsilon\mathcal{C}(\mu_{1}\odot T)\}\,,\quad\lambda\in\mathbb{R}^{M}$ (21) $$

For each $\lambda\in\mathbb{R}^{M}$ and $\varepsilon>0$, denote

$$ $\displaystyle B_{\lambda,\varepsilon}(x)$ $\displaystyle=\varepsilon\log\int_{y\in\mathcal{X}}\exp\bigl{(}\varepsilon^{-1}\ell_{0}^{\lambda}(x,y)\bigr{)}\mu_{2}(dy)\quad$ $\displaystyle x\in\mathcal{X}$ (22) $$

###### Proposition 2.2

Suppose that (A1)–(A3) hold.

- (i)
The function ([21](#S2.E21))
may be expressed
$\varphi^{*}(\lambda)=-\langle\mu_{1},B_{\lambda,\varepsilon}\rangle$.
- (ii)
The transition kernel maximizing ([21](#S2.E21)) is
$T^{\lambda}(x,dy)=\mu_{2}(dy)\exp(L^{\lambda}(x,y))\,,\qquad L^{\lambda}(x,y)=\varepsilon^{-1}\{\ell_{0}^{\lambda}(x,y)-B_{\lambda,\varepsilon}(x)\}$
(23a)
- (iii)
There is no duality gap: a unique $\lambda^{*}\in\mathbb{R}^{M}$ exists satisfying
$\varphi^{*}(\lambda^{*})=d_{\varepsilon}(\mu_{1},{\cal P}_{f,r})$
(23b)

The similarity between Proposition [2.2](#S2.Thmtheorem2) and Proposition [2.1](#S2.Thmtheorem1) is found through examination of
([17](#S2.E17)), and the recognition that $-B_{\lambda,\varepsilon}(x)$ is a ($\mu_{2}$-weighted) soft minimum of $-\ell_{0}^{\lambda}(x,y)=c(x,y)-\lambda^{\intercal}\tilde{f}(y)$ over $y\in\mathcal{X}$. Subject to this interpretation, the convex dual of OT-FPR can be expressed in a form entirely analogous to ([17](#S2.E17)):

$$ $\max_{\lambda}\varphi^{*}(\lambda)=\max_{\lambda}\int\mathop{\rm softmin}_{y}\{c(x,y)-\lambda^{\intercal}\tilde{f}(y)\}\mu_{1}(dx)$ (24) $$

###### Proposition 2.2

#### OT-FP approximation

Consider the following procedure to obtain a solution to OT-FP (without regularization), but with $\mathcal{X}$ compact. Choose $\mathcal{C}$ so that the support of $\mu_{1}$ and $\mu_{2}$ are each equal to all of $\mathcal{X}$. Let $\{\gamma^{\varepsilon},\lambda^{\varepsilon}:\varepsilon>0\}$ denote primal-dual solutions to OT-FPR, where $\varepsilon>0$ is the scaling in ([7](#S1.E7)). Hence for each $\varepsilon>0$,

$$ $d_{\varepsilon}(\mu_{1},{\cal P}_{f,r})=\langle\gamma^{\varepsilon},c\rangle+\varepsilon\mathcal{C}(\gamma^{\varepsilon})=-\langle\mu_{1},B_{\lambda^{\varepsilon},\varepsilon}\rangle$ $$

###### Proposition 2.3

Suppose that the assumptions of Proposition [2.1](#S2.Thmtheorem1) (ii) hold, so in particular $\mathcal{X}$ is compact. Then, any weak subsequential limit of $\{\gamma^{\varepsilon},\lambda^{\varepsilon}:\varepsilon>0\}$ as $\varepsilon\downarrow 0$ defines a pair $(\gamma^{0},\lambda^{0})$ for which $\gamma^{0}$ solves OT-FP and $\lambda^{0}$ achieves the supremum in ([17](#S2.E17)).

###### Proposition 2.3

#### OT-FPR calculus

We turn next to representation of the derivatives of the dual function.
The quantity $\varepsilon^{-1}B_{\lambda,\varepsilon}(x)$ is a log moment generating function for each $x$; for this reason, it is not difficult to obtain suggestive expressions for the first and second derivatives with respect to $\lambda$. This of course is what is needed to apply gradient ascent to maximize the function $\varphi^{*}$.

It is convenient to make the change of variables $\zeta=\varepsilon^{-1}\lambda$, and consider

$$ ${\cal J}(\zeta)\mathbin{:=}-\varepsilon^{-1}\varphi^{*}(\varepsilon\zeta)=\int_{x\in\mathcal{X}}\mu_{1}(dx)\log\int_{y\in\mathcal{X}}\mu_{2}(dy)\exp\bigl{(}\{\zeta^{\intercal}\tilde{f}(y)-\varepsilon^{-1}c(x,y)\}\bigr{)}$ $$

###### Proposition 2.4

The function ${\cal J}$ is convex and continuously differentiable.
The first and second derivatives of ${\cal J}$ admit the following representations:

$$ $\displaystyle\nabla{\cal J}(\zeta)=m^{\lambda}\,,\qquad\nabla^{2}{\cal J}(\zeta)=\Sigma^{\lambda}$ (25a) $$

in which $m^{\lambda}_{i}=\langle\mu^{\lambda},\tilde{f}_{i}\rangle={\sf E}^{\lambda}[\tilde{f}_{i}(Y)]$ for each $i$,
and the Hessian ([25a](#S2.E25.1)) coincides with the conditional covariance:

$$ $\Sigma^{\lambda}={\sf E}^{\lambda}[f(Y)f(Y)^{\intercal}]-{\sf E}^{\lambda}\bigl{[}{\sf E}^{\lambda}[f(Y)\mid X]{\sf E}^{\lambda}[f(Y)\mid X]^{\intercal}\bigr{]}$ (25b) $$

It follows that ${\cal J}$ is strictly convex:

###### Lemma 2.5

Suppose that (A1)–(A3) hold. Then, the covariance
$\Sigma^{\lambda}$ defined in ([25b](#S2.E25.2)) is full rank for any $\lambda\in\mathbb{R}^{M}$.

###### Proposition 2.4

###### Lemma 2.5

### 2.3 Soft Projection

The construction of a dual to ([9](#S1.E9)) proceeds as in the analysis of ([7](#S1.E7)), in which we relax only the moment constraints:

$$ $\displaystyle\varphi^{*}(\lambda)\mathbin{:=}\min_{\mu,\gamma,\delta}\Bigl{\{}\langle\gamma,c\rangle+\varepsilon\mathcal{C}(\gamma)+{\cal R}(\delta)+\sum_{i=1}^{M}\lambda_{i}\{\delta_{i}-\langle\mu,\tilde{f}_{i}\rangle\}:\gamma\in\Upgamma(\mu_{1},\mu)\,,\ \delta\in\mathbb{R}^{M}\Bigr{\}}$ $$

###### Proposition 2.6

The dual may be expressed
$\varphi^{*}(\lambda)=-\langle\mu_{1},B_{\lambda,\varepsilon}\rangle-{\cal R}^{*}(-\lambda)$, $\lambda\in\mathbb{R}^{M}$,
in which the function $B_{\lambda,\varepsilon}$ is defined in ([22](#S2.E22)).

###### Proposition 2.6

### 2.4 Monte-Carlo approximation

Provided we can efficiently draw samples from $\gamma^{\lambda}$ for each $\lambda$, we can construct a stochastic approximation algorithm to approximate $\lambda^{*}$ using stochastic gradient descent. Suppose that $\{X_{n}\}$ is i.i.d. $[\mu_{1}]$, and given an estimate $\lambda^{n}\in\mathbb{R}^{M}$ and the observation $x=X_{n+1}$, we draw $Y_{n+1}\sim T^{\lambda_{n}}(x,\,\cdot\,)$ independently of $\{(X_{k},Y_{k}):k\leq n\}$.

The algorithms described below are inspired by Proposition [2.4](#S2.Thmtheorem4), and hence we adopt the scaling $\zeta^{n}=\varepsilon^{-1}\lambda^{n}$.
We consider the OT-FPRP problem whose dual function is given in Proposition [2.6](#S2.Thmtheorem6), and replace maximization of $\varphi^{*}$ with minimization of

$$ ${\cal J}(\zeta)=\varepsilon^{-1}\{\langle\mu_{1},B_{\varepsilon\zeta,\varepsilon}\rangle+{\cal R}^{*}(-\varepsilon\zeta)\}$ $$

Assume that $\nabla{\cal R}^{*}$ is Lipschitz continuous.

Given an initial condition $\zeta^{0}\in\mathbb{R}^{M}$, a non-negative step-size sequence $\{\alpha_{n}\}$, and positive definite matrices $\{G^{n}\}$, the stochastic gradient descent algorithm is defined by the recursion

$$ $\zeta^{n+1}=\zeta^{n}-\alpha_{n+1}G^{n+1}\{\widetilde{m}^{n+1}+\nabla{\cal R}^{*}(-\varepsilon\zeta^{n})\}$ (26) $$

in which ${\sf E}[\widetilde{m}^{n+1}_{i}\mid\lambda_{n}]=\langle\gamma^{\lambda^{n}}_{2},\tilde{f}_{i}\rangle$ for each $i$. For example, $\widetilde{m}^{n+1}=\tilde{f}(Y_{n+1})$.

If $M$ is not large we might opt for Zap stochastic approximation, in which $G^{n}$ approximates the inverse of $R(\lambda^{n})\mathbin{:=}\nabla^{2}{\cal J}\,(\lambda^{n})$. This can be achieved using $G^{n}=[\overline{\Sigma}^{n}]^{-1}$ for each $n$, where the estimates evolve according to

$$ $\overline{\Sigma}^{n+1}=\overline{\Sigma}^{n}+\beta_{n+1}\{\widetilde{\Sigma}^{n+1}-\overline{\Sigma}^{n}\}$ (27) $$

initialized with $\overline{\Sigma}^{0}>0$ (positive definite), and
where $\widetilde{\Sigma}^{n+1}$ is random with conditional mean ${\sf E}[\widetilde{\Sigma}^{n+1}\mid\lambda_{n}]=\Sigma^{\lambda_{n}}$ (see ([25b](#S2.E25.2))). The stepsize sequence is chosen with $\beta_{n}\gg\alpha_{n}$ for large $n$ (see Devraj et al. () for motivation).

Two choices for the construction of $\{\widetilde{\Sigma}^{n+1}\}$ are summarized in the following:

1. Conditional computation plus sampling:
In each of the general examples described in Section [3](#S3) it is not computationally expensive to compute the conditional means $m_{i}\mathbin{:=}T^{\lambda}f_{i}$ and second moments $m_{ij}^{2}\mathbin{:=}T^{\lambda}[f_{i}f_{j}]$ for each $i,j$. In this case we draw $\{X_{n}\}$ i.i.d., and take for each $n$,

$$ $\widetilde{m}^{n+1}_{i}=m_{i}(X_{n+1})\,,\quad\widetilde{\Sigma}^{n+1}_{ij}=m_{ij}^{2}(X_{n+1})-\widetilde{m}^{n+1}_{i}[\widetilde{m}^{n+1}_{j}]^{\intercal}\qquad 1\leq i,j\leq M$ $$

2. Split Sampling: This approach requires that we obtain $\{X_{n}\}$ i.i.d., and also be able to draw a sample from the transition kernel.
With $K\geq 2$ an integer:
draw $x=X_{n+1}$ from $\mu_{1}$, and then draw $K$ independent samples $\{Y_{n+1}^{k}:1\leq k\leq K\}$ from $T^{\lambda}(x,\,\cdot\,)$, independently of $\{(X_{k},Y_{k}):k\leq n\}$.
We then have for $k\neq j$,

$$ ${\sf E}^{\lambda}\bigl{[}{\sf E}^{\lambda}[f(Y)\mid X]{\sf E}^{\lambda}[f(Y)\mid X]^{\intercal}\bigr{]}={\sf E}[f(Y_{n+1}^{k})f(Y_{n+1}^{j})^{\intercal}]={\sf E}[f(Y_{n+1}^{j})f(Y_{n+1}^{k})^{\intercal}]$ $$

This justifies the choice $\displaystyle\widetilde{m}^{n+1}=\frac{1}{K}\sum_{k=1}^{K}\tilde{f}(Y_{n+1}^{k})$ and

$$ $\displaystyle\widetilde{\Sigma}^{n+1}$ $\displaystyle=\frac{1}{K}\sum_{k=1}^{K}f(Y_{n+1}^{k})f(Y_{n+1}^{k})^{\intercal}-\frac{1}{K^{2}-K}\sum_{k=1}^{K}\sum_{j\neq k}f(Y_{n+1}^{k})f(Y_{n+1}^{j})^{\intercal}$ $$

## 3 Examples

#### Quadratic moment constraints with Gaussian regularizer

Consider the following special case in which the function $f$ is designed to specify all first and second moments for $Y$. To solve Problem 2 we adopt the following notational conventions for the Lagrange multiplier:
${\sf E}[Y_{i}]=m^{1}_{i}\longleftrightarrow\lambda^{1}_{i}$ and ${\sf E}[Y_{i}Y_{j}]=m^{2}_{ij}\longleftrightarrow\lambda^{2}_{i,j}$.
Of course we have $m^{2}_{ij}=m^{2}_{ji}$ for each $i,j$. The total number of constraints is thus $M=n+n(n+1)/2$.
For purposes of calculation it is useful to introduce the symmetric matrices $M_{Y}^{2}$ and $\Lambda^{2}$ with respective entries $\{m^{2}_{ij}\}$ and $\{\lambda^{2}_{ij}\}$; similarly notation is used for $m_{Y}$ and $\lambda^{1}$, the $n$-dimensional vectors with entries $\{m^{1}_{i}\}$ and $\{\lambda^{1}_{i}\}$.

Recalling the notation ([22](#S2.E22)) gives $\ell_{0}^{\lambda}(x,y)=\lambda^{\intercal}\tilde{f}(y)-c(x,y)$ with

$$ $\lambda^{\intercal}\tilde{f}(y)=y^{\intercal}\Lambda^{2}y-\langle\Lambda^{2},M_{Y}^{2}\rangle+y^{\intercal}\lambda^{1}-m_{Y}^{\intercal}\lambda^{1}$ (28) $$

An explicit solution to problem 2 is obtained when $c$ is quadratic and $\mu_{2}$ is Gaussian:

###### Proposition 3.1

Consider the OT-FPR optimization problem ([7](#S1.E7))
in the following special case: $c(x,y)={\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|x-y\|^{2}$, with $\mu_{2}=N(0,1)$ in the choice of regularizer ([6](#S1.E6)). Suppose moreover that the target covariance $\Sigma_{Y}\mathbin{:=}M_{Y}^{2}-m_{Y}m_{Y}^{\intercal}$ is positive definite.

$$ The following conclusions then hold: $$

- (i)
For each $\lambda$ for which $\Lambda^{2}<{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}(1+\varepsilon)I$, the transition kernel $T^{\lambda}$ is Gaussian:
conditioned on $X=x$, the distribution of $Y$ is Gaussian $N(m_{T^{\lambda}}^{x},\Sigma_{T^{\lambda}})$ with
$m_{T^{\lambda}}^{x}=\varepsilon^{-1}\Sigma_{T^{\lambda}}[x+\lambda^{1}]\qquad\Sigma_{T^{\lambda}}=\bigl{[}I+\varepsilon^{-1}[I-2\Lambda^{2}]\bigr{]}^{-1}$
(29a)
The representation ([29a](#S3.E29.1)) holds for any choice of $\mu_{1}$.
In the following it is assumed that $\mu_{1}=\mu_{2}$ so that $X\sim N(0,I)$:
- (ii)
When $\varphi^{*}(\lambda)$ is finite, the bivariate distribution $\gamma^{\lambda}$ is a two dimensional Gaussian random vector satisfying ${\sf E}^{\lambda}[Y]=\varepsilon^{-1}\Sigma_{T^{\lambda}}\lambda^{1}$, and with $(2n)\times(2n)$ covariance matrix
$\Sigma^{\lambda}=\begin{bmatrix}I+\varepsilon^{-2}\Sigma_{T^{\lambda}}&\varepsilon^{-1}I\\
\varepsilon^{-1}I&\Sigma_{T^{\lambda}}^{-1}\end{bmatrix}^{-1}=\begin{bmatrix}I&-\varepsilon^{-1}\Sigma_{T^{\lambda}}\\
-\varepsilon^{-1}\Sigma_{T^{\lambda}}&\Sigma_{T^{\lambda}}+\varepsilon^{-2}\Sigma_{T^{\lambda}}^{2}\end{bmatrix}$
(29b)
- (iii)
The optimizer $\lambda^{*}$ determines $T^{\lambda^{*}}$ as follows: Write
$\Sigma_{Y}=U^{-1}DU$ where $U$ is an orthogonal transformation and $D$ is diagonal. Then,
$\Sigma_{T^{\lambda^{*}}}=U^{-1}ZU$
(29c)
where $Z$ is the diagonal matrix with entries
$Z_{i,i}=\frac{1}{2}\bigl{\{}-\varepsilon^{2}+\sqrt{\varepsilon^{4}+4\varepsilon^{2}D_{i,i}}\bigr{\}}$.
This combined with ([29a](#S3.E29.1)) give $\lambda^{*}$, and
$\gamma^{*}=\mu_{1}\odot T^{\lambda^{*}}$.

###### Proposition 3.1

#### Applications to stochastic control

We consider here a tracking problem, based on the following conventions:

- 1.
The random variable $X$ is expressed as a vector of length $M+1$, with $X_{k}\in{\sf X}$ the state of the system at time $k$, so that $\mathcal{X}={\sf X}^{M+1}$. We write $X=(X_{0},\dots,X_{M})$, with $X_{0}$ interpreted as the initial condition, and
$M$ is considered a time horizon. $X\sim\mu_{1}$ represents nominal behavior, and $Y\sim\gamma_{2}^{*}$ behavior under optimal control.
- 2.
For a fixed function ${\mathcal{U}}\colon{\sf X}\to{\sf X}$ and any $y=(y_{0},\dots,y_{M})\in\mathcal{X}$, we take
$f_{i}(y)={\mathcal{U}}(y_{i})$ for $1\leq i\leq M$, and $r\in\mathbb{R}^{M}$ is a reference signal to be tracked.
In OT-FPR we demand perfect tracking, ${\sf E}[{\mathcal{U}}(Y_{k})]=r_{k}$ for $1\leq k\leq M$.
- 3.
The regularizer $\mathcal{C}$ is used with $\mu_{2}=\mu_{1}$. The objective function $\langle\gamma,c\rangle+\varepsilon\mathcal{C}(\gamma)$ is
regarded as a control cost: a large penalty is incurred if the distribution of $Y$ is far from $X$.
Moreover, the initialization $Y_{0}=X_{0}$ is enforced, whose distribution on ${\sf X}$ is denoted $\nu_{0}$.

Suppose that $X$ is Markovian, so that $\mu_{1}$ takes the form

$$ $\mu_{1}(dx)=\nu_{0}(dx_{0})\prod_{i=1}^{M}P_{i}(x_{i-1},dx_{i})$ $$

where $P_{i}$ is a transition kernel for each $i$.
Note that $\mu_{1}$ may have finite support, in which case each $P_{i}$ is a transition matrix.

The following is the main conclusion of Proposition [A.2](#A1.Thmtheorem2) in the supplementary material.

###### Proposition 3.2

The conditional distribution defined in ([23a](#S2.E23.1)) is Markovian: for a collection of transition kernels $\{{\check{P}}_{i}^{\lambda}\}$ parameterized by $x$,

$$ $T^{\lambda}(x,dy)=\nu_{0}(dy_{0})\prod_{i=1}^{M}{\check{P}}_{i}^{\lambda}(y_{i-1},dy_{i};x)$ (30) $$

The proposition motivates a mixed Monte-Carlo approach. Consider the stochastic gradient descent algorithm ([26](#S2.E26)) in which $\widetilde{m}^{n+1}$ is constructed in two phases. First $x=X_{n+1}$ is drawn from $\mu_{1}$, and then for each $k$ define a sequence of conditional marginals inductively: $\nu_{0}(dy_{0};x)=\nu_{0}(dy_{0})$, and for $k\geq 1$,

$$ $\displaystyle\nu_{k}(dy_{k};x)$ $\displaystyle=\int_{\sf X}\nu_{k-1}(dy;x){\check{P}}_{k}^{\lambda}(y,dy_{k};x)\,,\qquad\widetilde{m}^{n+1}_{k}=\int_{\sf X}{\mathcal{U}}(y)\nu_{k}(dy;x)-r_{k}$ $$

This is practical when ${\sf X}$ is discrete so that each integral reduces to a sum.

###### Proposition 3.2

## 4 Conclusions

Relaxation of marginal constraints is a natural approach to address complexity in complex applications of OT, and prior work motivates relaxations for applications to decentralized control and outlier detection. It is expected that the general framework of this paper will find application in other areas.

Remaining open questions concern computation. In particular: (i) It is likely that relaxation of $\mu_{1}$ will reduce computational complexity, but may lead to undesirable side effects in some applications; (ii) We have established that the solution to OT-PF (without regularization) can be obtained via solution to OT-PFR, with $\varepsilon>0$ treated as a variable. An ODE approach will be considered in future research:
using the change of variable $t=\varepsilon^{-1}$, let $\gamma^{t}$ denote the solution to OT-PFR using $\varepsilon=1/t$, and seek a vector field $F\colon\mathbb{R}^{M+1}\to\mathbb{R}^{M}$ for which

$$ ${\mathchoice{\genfrac{}{}{}{1}{d}{dt}}{\genfrac{}{}{}{1}{d}{dt}}{\genfrac{}{}{}{3}{d}{dt}}{\genfrac{}{}{}{3}{d}{dt}}}\lambda^{t}=F(\lambda^{t},t)\,,\qquad\lambda^{0}=\mu_{1}\otimes\mu_{2}$ $$

This approach was adopted in Bušić and Meyn () for computation in stochastic control.

## Appendix A Supplementary Material

Much of the analysis that follows is based on convex duality between relative entropy and log moment generating functions. For any probability measure $\mu$ on ${\cal B}(\mathcal{X})$ and function $g\colon\mathcal{X}\to\mathbb{R}$, the log moment generating function is denoted,

$$ $\Lambda_{\mu}(g)=\log\langle\mu,e^{g}\rangle$ $$

With $\mu$ fixed, this is viewed as an extended-valued, convex functional on the space of Borel measurable functions. Lemma [A.1](#A1.Thmtheorem1) is a standard a standard tool in information theory
Dembo and Zeitouni (), and a reason that relative entropy is popular for use as a regularizer in optimization.

###### Lemma A.1

Relative entropy and the log moment generating function are related via convex duality:

- (i)
For any probability measure $\pi$ we have
$D(\pi\|\mu)=\sup_{g}\{\langle\pi,g\rangle-\Lambda_{\mu}(g)\}$
(31a)
If $D(\pi\|\mu)<\infty$ then the supremum is achieved, with optimizer equal to the log likelihood ratio, $g^{*}=\log(d\pi/d\mu)$.
- (ii)
For Borel measurable $g\colon\mathcal{X}\to\mathbb{R}$,
$\Lambda_{\mu}(g)=\sup_{\pi}\{\langle\pi,g\rangle-D(\pi\|\mu)\}$
(31b)
If $\Lambda_{\mu}(g)<\infty$ then the supremum is achieved, where the optimizer $\pi^{*}$ has log likelihood ratio $\log(d\pi^{*}/d\mu)=g-\Lambda_{\mu}(g)$.
$\sqcap$$\sqcup$

###### Lemma A.1

#### Proof of Proposition 2.2

For each $\lambda$ we have by definition,

$$ $\displaystyle\varphi^{*}(\lambda)$ $\displaystyle=\min_{T}\int_{x\in\mathcal{X}}\mu_{1}(dx)\Bigl{\{}\varepsilon D(T(x,\,\cdot\,)\|\mu_{2})-\int_{y\in\mathcal{X}}T(x,dy)\ell_{0}^{\lambda}(x,y)\Bigr{\}}$ (32) $\displaystyle=-\varepsilon\max_{T}\int_{x\in\mathcal{X}}\mu_{1}(dx)\Bigl{\{}\varepsilon^{-1}\int_{y\in\mathcal{X}}T(x,dy)\ell_{0}^{\lambda}(x,y)-D(T(x,\,\cdot\,)\|\mu_{2})\Bigr{\}}$ $$

For each $x$ we have an optimization problem of the form ([31b](#A1.E31.2)). Applying
Lemma [A.1](#A1.Thmtheorem1) (ii) gives the representation ([23a](#S2.E23.1))
and by substitution (or applying ([31b](#A1.E31.2))) we obtain

$$ $\displaystyle\Bigl{\{}\varepsilon^{-1}\int_{y\in\mathcal{X}}T^{\lambda}(x,dy)\ell_{0}^{\lambda}(x,y)-D(T^{\lambda}(x,\,\cdot\,)\|\mu_{2})\Bigr{\}}=\varepsilon^{-1}B_{\lambda,\varepsilon}(x)$ (33) $$

Integrating with respect to $\mu_{1}$ and applying ([32](#A1.E32)) completes the proof.
$\sqcap$$\sqcup$

#### Proof of Lemma 2.5

Suppose that $v\in\mathbb{R}^{M}$ is in the null space: $\Sigma^{\lambda}v=0$. From the definition ([25b](#S2.E25.2)) it follows that

$$ $0=v^{\intercal}\Sigma^{\lambda}v={\sf E}^{\lambda}\bigl{[}\bigl{\{}v^{\intercal}\bigl{(}f(Y)-{\sf E}^{\lambda}[f(Y)\mid X]\bigl{)}\bigl{\}}^{2}\bigl{]}$ $$

Equivalently, there is a function $g\colon\mathcal{X}\to\mathbb{R}$ such that

$$ $v^{\intercal}f(Y)=g(X)\qquad a.s.\ [\gamma^{\lambda}]$ $$

The probability measures $\gamma^{\lambda}$ and $\gamma^{0}\mathbin{:=}\mu_{1}\otimes\mu_{2}$ are mutually absolutely continuous, so the same equation holds under a.s. $[\gamma^{0}]$. Independence gives

$$ $v^{\intercal}f(Y)={\sf E}^{0}[v^{\intercal}f(Y)\mid Y]={\sf E}^{0}[g(X)\mid Y]=\langle\mu_{1},g\rangle\qquad a.s.\ [\gamma^{0}]$ $$

That is, the variance of $v^{\intercal}f(Y)$ is equal to zero.
Under (A3) this is possible only if $v=0$.
$\sqcap$$\sqcup$

#### Proof of Proposition 2.4

Recall the notation $\mu^{\lambda}=\mu_{1}T^{\lambda}$, which is the second marginal of $\gamma^{\lambda}$, and the probabilistic notation ([14](#S2.E14)). Also, by definition we have ${\cal J}(\zeta)=\varepsilon^{-1}\langle\mu_{1},B_{\varepsilon\zeta,\varepsilon}\rangle$.

We have for each $i$,

$$ $\varepsilon^{-1}\frac{\partial}{\partial\zeta_{i}}B_{\varepsilon\zeta,\varepsilon}(x)=\frac{\int_{y\in\mathcal{X}}\mu_{2}(y)\exp\bigl{(}\{\zeta^{\intercal}\tilde{f}(y)-\varepsilon^{-1}c(x,y)\}\bigr{)}\tilde{f}_{i}(y)}{\int_{y\in\mathcal{X}}\mu_{2}(y)\exp\bigl{(}\{\zeta^{\intercal}\tilde{f}(y)-\varepsilon^{-1}c(x,y)\}\bigr{)}}=T^{\lambda}\tilde{f}_{i}\,(x)$ $$

Integrating each side over $\mu_{1}$ gives ([25a](#S2.E25.1)) (recall that $\mu^{\lambda}=\mu_{1}T^{\lambda}$).

To obtain the second derivative of ${\cal J}(\zeta)$ requires the first derivative of the log-likelihood:

$$ $L^{\varepsilon\zeta}_{j}(x,y)\mathbin{:=}\frac{\partial}{\partial\zeta_{j}}L^{\varepsilon\zeta}(x,y)=\frac{\partial}{\partial\zeta_{j}}\bigl{[}\zeta^{\intercal}\tilde{f}(y)-\varepsilon^{-1}B_{\varepsilon\zeta,\varepsilon}(x)\bigr{]}=\tilde{f}_{j}(y)-T^{\lambda}\tilde{f}_{j}\,(x)$ $$

From this we obtain,

$$ $\displaystyle\frac{\partial^{2}}{\partial\zeta_{i}\partial\zeta_{j}}B_{\varepsilon\zeta,\varepsilon}(x)$ $\displaystyle=\frac{\partial}{\partial\zeta_{j}}T^{\varepsilon\zeta}\tilde{f}_{i}\,(x)$ $\displaystyle=\int T^{\varepsilon\zeta}(x,dy)\{L^{\varepsilon\zeta}_{j}(x,y)\tilde{f}_{i}(y)\}$ $\displaystyle=\int T^{\varepsilon\zeta}(x,dy)\tilde{f}_{j}(y)\tilde{f}_{i}(y)-T^{\lambda}\tilde{f}_{j}(x)\int T^{\varepsilon\zeta}(x,dy)\tilde{f}_{j}(y)$ $\displaystyle={\sf E}^{\lambda}[\tilde{f}_{j}(Y)\tilde{f}_{i}(Y)\mid X=x]-{\sf E}^{\lambda}[\tilde{f}_{i}(Y)\mid X=x]{\sf E}^{\lambda}[\tilde{f}_{j}(Y)\mid X=x]$ $$

Integrating each side over $\mu_{1}$ gives
([25b](#S2.E25.2)).
$\sqcap$$\sqcup$

#### Proofs of Proposition 2.1 and Proposition 2.3

Part (i) of
Proposition [2.1](#S2.Thmtheorem1) was established following the proposition. The proofs of part (ii) and Proposition [2.3](#S2.Thmtheorem3) are established simultaneously.

Let $(\gamma^{\varepsilon},\lambda^{\varepsilon})$ denote the solution to OT-FPR, with $\varepsilon>0$ regarded as a variable.

We let
$(\gamma^{0},\lambda^{0})$ denote any weak sub-sequential limit: for a sequence $\{\varepsilon_{i}\downarrow 0\}$,

$$ $\gamma^{\varepsilon_{i}}\buildrel\rm w\over{\longrightarrow}\gamma^{0}\,,\qquad\lambda^{\varepsilon_{i}}\to\lambda^{0}\,,\qquad i\to\infty.$ $$

Optimality of $\gamma^{0}$ is established in the following steps:

- (i)
Subject to (A1) and (A2) we know that $\gamma^{0}\in\Upgamma_{f,r}(\mu_{1})$.
- (ii)
For any $\gamma\in\Upgamma_{f,r}(\mu_{1})$ for which $\mathcal{C}(\gamma)<\infty$ and any $\varepsilon>0$ we have
$\displaystyle\langle\gamma^{0},c\rangle=\lim_{i\to\infty}\langle\gamma^{\varepsilon_{i}},c\rangle$
$\displaystyle\leq\lim_{i\to\infty}\{\langle\gamma^{\varepsilon_{i}},c\rangle+\varepsilon_{i}\mathcal{C}(\gamma^{\varepsilon_{i}})\}$
$\displaystyle\leq\lim_{i\to\infty}\{\langle\gamma,c\rangle+\varepsilon_{i}\mathcal{C}(\gamma)\}=\langle\gamma,c\rangle$
- (iii)
Under the support assumption we can approximate in the weak topology any $\gamma\in\Upgamma_{f,r}(\mu_{1})$ by a probability measure $\gamma^{\delta}$ satisfying $\mathcal{C}(\gamma^{\delta})<\infty$ and
$\langle\gamma^{0},c\rangle\leq\langle\gamma^{\delta},c\rangle\leq\langle\gamma,c\rangle-\delta$
Since $\delta>0$ is arbitrary this establishes optimality.

We next show $\lambda^{0}$ provides a solution to ([17](#S2.E17)). Proposition [2.2](#S2.Thmtheorem2) gives for any $\lambda$,

$$ $\langle\gamma^{0},c\rangle\geq-\lim_{i\to\infty}\langle\mu_{1},B_{\lambda,\varepsilon_{i}}\rangle=\int\inf_{y}\{c(x,y)-\lambda^{\intercal}\tilde{f}(y)\}\mu_{1}(dx)$ $$

The lower bound is achieved using $\lambda^{0}$ by allowing $\lambda$ to depend on $i$:

$$ $\langle\gamma^{0},c\rangle\leq\lim_{i\to\infty}\{\langle\gamma^{\varepsilon_{i}},c\rangle+\varepsilon_{i}\mathcal{C}(\gamma^{\varepsilon_{i}})\}=-\lim_{i\to\infty}\langle\mu_{1},B_{\lambda^{\varepsilon_{i}},\varepsilon_{i}}\rangle=\int\inf_{y}\{c(x,y)-{\lambda^{0}}^{\intercal}\tilde{f}(y)\}\mu_{1}(dx)$ $$

$\sqcap$$\sqcup$

#### Proof of Proposition 2.6

Minimizing over $T$ in the representation $\gamma=\mu_{1}\odot T$ we obtain

$$ $\displaystyle\varphi^{*}(\lambda)=$ $\displaystyle\min_{\delta}\bigl{\{}{\cal R}(\delta)+\lambda^{\intercal}\delta\bigr{\}}$ (34) $\displaystyle+\min_{T}\int_{x\in\mathcal{X}}\mu_{1}(dx)\Bigl{\{}\varepsilon D(T(x,\,\cdot\,)\|\mu_{2})-\int_{y\in\mathcal{X}}T(x,dy)\ell_{0}^{\lambda}(x,y)\Bigr{\}}$ $$

The solution to the minimization over $T$ may be found in the proof of Proposition [2.2](#S2.Thmtheorem2).

The minimum over $\delta$ can be transformed to a maximization of the form ([8](#S1.E8)):

$$ $\min_{\delta}\bigl{\{}{\cal R}(\delta)+\lambda^{\intercal}\delta\bigr{\}}=-\max_{\delta}\bigl{\{}(-\lambda)^{\intercal}\delta-{\cal R}(\delta)\bigr{\}}=-{\cal R}^{*}(-\lambda)$ $$

$\sqcap$$\sqcup$

#### Proof of Proposition 3.1

From ([28](#S3.E28)) and using $c(x,y)={\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|x-y\|^{2}$
we obtain an expression for the likelihood $L^{\lambda}$ appearing in ([23a](#S2.E23.1)):

$$ $L^{\lambda}(x,y)=\varepsilon^{-1}\bigl{\{}y^{\intercal}\Lambda^{2}y+y^{\intercal}\lambda^{1}-\kappa^{\lambda}-B_{\lambda,\varepsilon}(x)\}-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}(\|x\|^{2}-2x^{\intercal}y+\|y\|^{2})\bigr{\}}$ $$

with $\kappa^{\lambda}=\langle\Lambda^{2},M_{Y}^{2}\rangle+m_{Y}^{\intercal}\lambda^{1}$.
The expression for $T^{\lambda}$ in ([23a](#S2.E23.1)) using $\mu_{2}=N(0,I)$ then implies that for any $x$, $T^{\lambda}(x,dy)$ admits the Gaussian density

$$ $\uptau^{\lambda}(y\mid x)=\frac{1}{n^{\lambda}(x)}\exp\bigl{(}-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|y\|^{2}\bigr{)}\exp\bigl{(}\varepsilon^{-1}\{-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}y^{\intercal}[I-2\Lambda^{2}]y+y^{\intercal}[x+\lambda^{1}]\}\bigr{)}$ $$

where $n^{\lambda}(x)=(2\pi)^{n/2}\exp\bigl{(}\varepsilon^{-1}\{\kappa^{\lambda}+B_{\lambda,\varepsilon}(x)+{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|x\|^{2}\}\bigr{)}$ may be regarded as a normalizing constant. This proves (i).

To establish (ii) we begin with the representation for the density $\uptau^{\lambda}$ of $\gamma^{\lambda}$:

$$ $\uptau^{\lambda}(x,y)=\frac{1}{(2\pi)^{n/2}}\exp(-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|x\|^{2})\uptau^{\lambda}(y\mid x)$ $$

To compute the right hand side requires that we compute the normalizing constant $n^{\lambda}$. We can avoid direct computation by exploiting the fact that $T^{\lambda}$ is constructed to preserve the marginal constraint $X\sim\mu_{1}=N(0,1)$. It can be shown that $\log(n^{\lambda})$ is quadratic in $x$, and hence for a symmetric matrix $A$ and vector $z$,

$$ $\uptau^{\lambda}(x,y)\propto\exp\bigl{(}-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}x^{\intercal}Ax+x^{\intercal}z-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|y\|^{2}+\varepsilon^{-1}\{-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}y^{\intercal}[I-2\Lambda^{2}]y+y^{\intercal}[x+\lambda^{1}]\}\bigr{)}$ $$

This implies part of the desired representation for the covariance:

$$ $\Sigma^{\lambda}=\begin{bmatrix}\Sigma^{\lambda}_{1,1}&\Sigma^{\lambda}_{1,2}\\ \Sigma^{\lambda}_{2,1}&\Sigma^{\lambda}_{2,2}\end{bmatrix}=\begin{bmatrix}A&-\varepsilon^{-1}I\\ -\varepsilon^{-1}I&\Sigma_{T^{\lambda}}^{-1}\end{bmatrix}^{-1}$ $$

To obtain $A$ we impose the constraint $\Sigma^{\lambda}_{1,1}=I$ and apply the block inverse formula Bernstein ():

$$ $\begin{bmatrix}{A}&{B}\\ {C}&{D}\end{bmatrix}^{-1}={\begin{bmatrix}\left({A}-{BD}^{-1}{C}\right)^{-1}&-\left({A}-{BD}^{-1}{C}\right)^{-1}{BD}^{-1}\\ -{D}^{-1}{C}\left({A}-{BD}^{-1}{C}\right)^{-1}&\quad D^{-1}+{D}^{-1}{C}\left({A}-{BD}^{-1}{C}\right)^{-1}{BD}^{-1}\end{bmatrix}}$ $$

giving $I=\Sigma^{\lambda}_{1,1}=[A-\varepsilon^{2}\Sigma_{T^{\lambda}}]^{-1}$, so that
$A=I+\varepsilon^{-2}\Sigma_{T^{\lambda}}$. Given $\left({A}-{BD}^{-1}{C}\right)^{-1}=I$, the matrix inverse formula simplifies to give

$$ $\displaystyle\Sigma^{\lambda}_{1,2}$ $\displaystyle=-{BD}^{-1}=\varepsilon^{-1}\Sigma_{T^{\lambda}}$ $\displaystyle\Sigma^{\lambda}_{2,2}$ $\displaystyle=D^{-1}+{D}^{-1}{C}{BD}^{-1}=\Sigma_{T^{\lambda}}+\varepsilon^{-2}\Sigma_{T^{\lambda}}^{2}$ $$

This completes the proof of ([29b](#S3.E29.2))
and (ii).

The proof of (iii) is obtained by matching second moments:
For $\lambda=\lambda^{*}$ we must achieve the target covariance $\Sigma_{Y}$, giving the Riccati equation:

$$ $-\Sigma_{T^{\lambda^{*}}}-\varepsilon^{-2}\Sigma_{T^{\lambda^{*}}}^{2}+\Sigma_{Y}=0$ $$

Setting $Z=U\Sigma_{T^{\lambda^{*}}}U^{-1}$ gives
$Z+\varepsilon^{-2}Z^{2}-D=0$,
whose unique positive definition solution is diagonal, with entries given in (iii).
$\sqcap$$\sqcup$

#### Computation for non-Gaussian 𝝁 𝟏 subscript 𝝁 1 \mu_{1}

In this case it is necessary to compute the normalizing constant in the definition of $T^{\lambda}$:

$$ $\displaystyle n^{\lambda}(x)$ $\displaystyle=n^{\lambda}(x)\int\uptau^{\lambda}(y\mid x)\,dy$ (35) $\displaystyle=\int\exp\bigl{(}-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}y^{\intercal}\Sigma_{T^{\lambda}}^{-1}y+\varepsilon^{-1}y^{\intercal}[x+\lambda^{1}]\bigr{)}\,dy$ $\displaystyle=\sqrt{{(2\pi)^{d}\,{\mathop{\rm det}}(\Sigma_{T^{\lambda}})}}\exp\bigl{(}{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\varepsilon^{-2}[x+\lambda^{1}]^{\intercal}\Sigma_{T^{\lambda}}[x+\lambda^{1}]\bigr{)}$ $$

Monte-Carlo methods can be used to estimate $\lambda^{*}$. Denote for each $x$,

$$ $q^{\lambda}(x)=\int T^{\lambda}(x,dy)f(y)\,,\quad m^{\lambda}(x)=\int T^{\lambda}(x,dy)f(y)f(y)^{\intercal}$ $$

Each have polynomial entries: $q^{\lambda}_{i}$ is a quadratic function of $x$ and
$m^{\lambda}_{i,j}(x)$ is a fourth order polynomial in $x$ for each $i,j$. In applying any of the algorithms described in Section [2.4](#S2.SS4) you might take

$$ $\displaystyle\widetilde{m}^{n+1}$ $\displaystyle=q^{\lambda_{n}}(X_{n+1})\,,\qquad\widetilde{\Sigma}^{n+1}=m^{\lambda_{n}}(X_{n+1})-\widetilde{m}^{n+1}[\widetilde{m}^{n+1}]^{\intercal}$ $$

These functions will have finite means provided ${\sf E}[\|X\|^{4}]$ is finite under $\mu_{1}$.

#### Proof of Proposition 3.2

The proof reduces to justifying ([30](#S3.E30)), which
is one component of Proposition [A.2](#A1.Thmtheorem2)
that follows.

Write $L_{i}^{\lambda}(x_{i},y_{i})=\varepsilon^{-1}\{\lambda_{i}({\mathcal{U}}(y_{i})-r_{i})-{\mathchoice{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{1}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}{\genfrac{}{}{}{3}{1}{2}}}\|x_{i}-y_{i}\|^{2}\}$, and for each $i$ consider the positive kernel,

$$ ${\widehat{P}}_{i}^{\lambda}(y_{i-1},dy_{i})=P_{i}(y_{i-1},dy_{i})\exp\bigl{(}L_{i}^{\lambda}(x_{i},y_{i})\bigr{)}$ $$

###### Proposition A.2

The conditional distribution defined in ([23a](#S2.E23.1)) can be expressed

$$ $T^{\lambda}(x,dy)=\nu_{0}(dy_{0})\exp\bigl{(}-\varepsilon^{-1}B_{\lambda,\varepsilon}(x)\bigr{)}\prod_{i=1}^{M}{\widehat{P}}_{i}^{\lambda}(y_{i-1},dy_{i})$ (36) $$

Consequently, conditioned on $X=x$, the process $Y$ is of the form ([30](#S3.E30)),
in which each kernel in the product takes the form,

$$ ${\check{P}}_{i}^{\lambda}(y_{i-1},dy_{i};x)=\frac{1}{g_{i-1}(y_{i-1};x)}{\widehat{P}}_{i}^{\lambda}(y_{i-1},dy_{i})g_{i}(y_{i};x)$ $$

The functions $\{g_{i}:0\leq i\leq M\}$ are defined inductively: $g_{M}(y_{M};x)\equiv 1$, and for $1\leq i\leq M$,

$$ $g_{i-1}(y;x)\mathbin{:=}\int{\widehat{P}}_{i}^{\lambda}(y,dy_{i})g_{i}(y_{i};x)\,,\quad y\in{\sf X}$ $$

This results in $g_{0}(y_{0},x)=\exp\bigl{(}\varepsilon^{-1}B_{\lambda,\varepsilon}(x)\bigr{)}$.

###### Proposition A.2

#### Proof

The representation ([36](#A1.E36)) follows from the definition ([23a](#S2.E23.1)) and the structure imposed on $\tilde{f}$ and $\mu_{1}$. It is then immediate that ([36](#A1.E36)) can be transformed to ([30](#S3.E30)): by construction,

$$ $\prod_{i=1}^{M}{\check{P}}_{i}^{\lambda}(y_{i-1},dy_{i};x)=\frac{1}{g_{0}(y_{0};x)}\prod_{i=1}^{M}{\widehat{P}}_{i}^{\lambda}(y_{i-1},dy_{i})$ $$

Since $y_{0}=x_{0}$ by construction, it also follows that

$$ $\exp\bigl{(}\varepsilon^{-1}B_{\lambda,\varepsilon}(x)\bigr{)}=g_{0}(x_{0};x)$ $$

$\sqcap$$\sqcup$

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: - Some comparison against existing EV charging algorithms (there are a lot) would be useful.
soundness: 2
strengths: - The problem being considered is an interesting one. 
- The idea of controlling an distribution to look like another (more optimal) distribution is certainly useful.
summary: This paper considers using optimal transport (OT) in mean field control, with constraints on the moments of the distributions. An algorithm (Sinkhorn algorithm) is proposed and an example of EV charging is considered.
weaknesses: It's important to note that I'm not that familiar with the field of optimal transport. But I felt the following are the weaknesses:
- The abstract, intro conclusion seem to promise a lot more than what the math actually delivers? The algorithm relies on Gibbs kernels, which feels pretty standard. How broadly applicable is this? 
- The EV problem presented is somewhat strange. The paper seems to say that the controllable variables are the EV arrival time and state-of-charge? But these are typically the main source of randomness in EV problems. The controllable knobs are typically the charging profiles. 
- How many EVs does there need to be for a mean field approximation to be valid? At any single charging station, there won't be that many EVs.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 1
questions: Please see my comments about the weakness.
soundness: 2
strengths: The EV Charging problem has received much attention recently. This work aims to optimize the consumption while satisfying the grid constraints.
summary: This paper studies the application of optimal transport to mean-field control. The main contribution is a representation of the mean-field control, Moment Constrained Optimal Transport for Control (MCOT-C), that aims to coordinate the agents and enforce constraints. The authors propose a variant of the Sinkhorn algorithm and apply the proposed algorithm to an EV charging application.
weaknesses: The theoretical part of this paper is hard for me to follow. The introduction starts with the mathematical problem settings without sufficient discussion about the background, significant challenges, and the motivation for using optimal transport in mean-field control. Besides, the theoretical problem setting, assumptions, and propositions are difficult to interpret. I suggest the authors add more discussions about the connections between the general framework and a specific example (e.g., the EV charging problem) that is easier to understand. Discussing the intuition/significance after stating each proposition or lemma is also helpful.

For the experiment part, I encourage the authors to compare the problem setting and the performance with related works on EV charging (e.g., [1]). Such comparisons can help the readers understand the advantages/limitations of the proposed approach.

[1] B. Alinia, M. H. Hajiesmaili, Z. J. Lee, N. Crespi, and E. Mallada, "Online EV Scheduling Algorithms for Adaptive Charging Networks with Global Peak Constraints," in IEEE Transactions on Sustainable Computing, vol. 7, no. 3, pp. 537-548, 1 July-Sept. 2022.

---

confidence: 2
contribution: 2
initial_rating: 3
presentation: 1
questions: - What is $l$ in equation (7)?
- Could you provide what the core observations and messages are that you want to share through the experiment sections?
- Could the authors reorganize the contributions they want to claim in an itemized format? It would be helpful if the reorganized contributions included a clear explanation of the new approaches, the challenges faced, and the advantages compared to previous approaches. Additionally, if the authors feel there are any points that the reviewer may have overlooked, highlighting those would be welcome.
soundness: 2
strengths: The approach of leveraging computational techniques from optimal control theory for control problems and the observations obtained from experiments applying the approach can be interesting. They provide the background theoretical derivation of such an approach. The approach focuses on a finite set of moments, so it could be more tractable in practice.
summary: This paper introduces moment constrained optimal transport for control (MCOT-C), which leverages computational techniques from optimal control theory for control problems. They provide an algorithm obtained by modifying the Sinkhorn algorithm by replacing the update on the second marginal with gradient descent on the dual. Then, they provide how their proposed approaches apply to mean field control applications, further providing an online version of MCOT-C.
weaknesses: **Writing:**
The reviewer thinks the writing of this paper needs to be improved. The reviewer was confused by the abstract and couldn't understand what contributions were made in this paper at first. The authors barely use phrases like 'this paper' or 'we,' so the actions taken in the paper were not distinguished clearly. It seems this issue occurs throughout the paper as well. The reviewer feels that the authors didn't clearly articulate the prior approaches, what they did new, and what the advantages are. This makes it challenging to understand the contributions they are claiming.

**Contribution:**
To the best of the reviewer's understanding, the contribution of this paper is that they are establishing the theoretical background to leverage computational techniques from optimal control theory to the mean field control problem, as presented in sections 2.1 and 2.2. In Proposition 2, they provide the calculation of derivatives for their problem and introduce a gradient descent-based algorithm in section 2.3. Then they directly jump to the experiments, and sections 3 and 4 consist of explanations of their experiments. However, the reviewer believes they should have provided more discussion about the method before proceeding to the experiments, for example, the motivation behind the specific design of the algorithm, or theoretical analysis, or some justifications for why they expect it would work well.

Additionally, the reviewer couldn't grasp what the authors were trying to claim with the experiments. The reviewer believes the authors should have provided guidance on interpreting their experimental results and offered clear conclusions or messages derived from the experiments. However, most of the content in the experiment section seems to focus on the details of the experiments.

**Minor issues:**
- Typo in line 072; a space between "S.It" is required.
- Typo in line 077; "litterature" should be corrected.
- Typo in line 419; "(i)" should be changed to "(ii)."

---

confidence: 2
contribution: 2
initial_rating: 3
presentation: 1
questions: Why is the complexity in section 3.1 $N_t^3 \times N_b$, when the state consists of only two variables with time, the arrival time and the plugging time?
soundness: 2
strengths: The authors propose an interesting and relevant control application.
They bridge the fields of optimal transport and mean field control by modifying the Sinkhorn algorithm, where the update on the second marginal is replaced by gradient descent.
summary: The authors propose a solution for an application concerning the charging of a fleet of electric vehicles, where a central planner decides on the plugging time of the vehicles. Their proposed method applies a modified version of the Sinkhorn algorithm, typically used in optimal transport problems, to the considered control problem.
weaknesses: The main weakness of the paper is its structure.
Instead of introducing the considered problem in the introduction, the authors start with their definitions of optimal transport and mean field control. The considered problem is formulated later. This makes it hard for the reader to follow.
When presenting a paper which mainly focuses on one application, it would be better to first clearly describe the application and then introduce the math and methods needed for the solution.

Another major issue is the lack of consistency in the notation. Some examples:
In Eq. 7 there is the variable $l$ which has not been introduced before.
In Eq. 16 there is $h$ on the left hand side and $f$ on the right hand side. Are these equivalent?
In section 3.1 it say the gradient is calculated on $\mathcal X \times \mathcal W$. Should this be $\mathcal S \times \mathcal W$?

Additionally there are many typos, missing punctuation and missing comma placement, further reducing the readability of the paper.

Finally the experiments do not compare to a baseline, apart from a naive decision rule, making it hard to evaluate the efficacy of the proposed method.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: I am not sure I have specific questions; perhaps one is "is the control specification in your EV charging case study truly realistic?". My issues are not necessarily fixable by a small change in the paper.
soundness: 3
strengths: - the paper seeks to connect two important areas of applied mathematics, and does it in an elegant fashion that is amenable to an analytical solution
- the paper is largely well-written and I suspect that it would be considered fairly easily readable by experts in the field
- the rigorous approach of the paper is refreshing and exposes the mathematical meat of the problem
summary: This is a half-theoretical, half-application paper which initially seeks to propose and solve the moment-constrained optimal transport for control problem, a variation of the constrained optimal transport problem for mean field control. After solving its first objective, the paper seeks to illustrate, and partially adapt, the developed results on an example of charging a large fleet of electric vehicles.
weaknesses: In short, I like this paper a lot, but I am just not convinced that ICLR is the right venue for it. I urge the authors to consider once again whether the ICLR audience is what this paper is shooting for (without any offense at all to either the paper or the ICLR audience). For instance, the L in ICLR stands for "learning". Yet, any connection to learning -- if there is any -- remains unexplained. My issues largely stem from this perceived mismatch:

- relegating the proofs, including the central results (Proposition 1/2), to the appendix probably does not make the paper more readable to the non-experts (who will find notation burdensome to start with), and only signals (possibly incorrectly) that the authors do not see these results as their primary contribution. In that case, I am not sure what is the primary contribution

- there seems to be some dissonance between the claimed contributions (which speak about control of an ensemble of agents), the actual technical contributions (which can indeed serve this goal, but the connection is not really described in detail), and the application (which is farfetched at best: why should the power consumption exactly track a reference trajectory?)

- the "second" technical part of the paper, coming within the application section, applies largely domain-agnostic mathematical theory to a problem that is tailored to the application (which is, as I already mentioned, farfetched). My suggestion would really be to split this paper into two: the first one dealing with the general problem of MCOT-C (and its online version in as much generality as possible), with perhaps only a small academic example if there is no room for more, and the second one applied, with a realistic case study and a detailed description of the algorithm implementation and possible minor modifications

- partly because of decoupling of the proof (and a weird ordering of the proofs in the Appendix), it is not clear to me how challenging the main result actually is: the duality seems rather straightforward. If this is not the case, that should be emphasized. If it is, perhaps it would be good to try to develop online results on MCOT-C in more generality

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 1
questions: Please see the "Weaknesses" section for my questions.
soundness: 1
strengths: -	The paper tackles the important, contemporary problem of charging EVs and uses a dataset to do so.
-	The combination of mean field control and optimal transport to optimize EV charging appears to be interesting.
-	The paper contains theoretical results to complement empirical findings.
summary: This paper combines the fields of mean field control and optimal transport to solve multi agent control problems. When doing so, the authors constrain some of the marginal distributions to specific moment classes and adopt the Sinkhorn algorithm to their setup. The theoretical and algorithmic considerations are complemented by an extensive example of EV charging in the Netherlands, based on a real dataset.
weaknesses: -	Introduction: In my opinion, the introduction should be less technical in the sense that there shouldn’t be several extensive mathematical expressions. In this way, especially non-expert readers can get a first impression of the paper’s contributions without being confronted with mathematical details.
-	Line 74: EVs are mentioned for the first time here (excluding the abstract). In my opinion the authors should lead with a motivating example and then move on to the technical tools like MFC. In in its current form, the introduction reads like a random collection of mathematical concepts. The authors should put more emphasis on their goals and high-level ideas.
-	Line 84-90: The contributions are formulated far too vague, for example, ‘‘Coordination of an ensemble of agents to achieve a desired goal’’ basically describes any cooperative multi-agent problem. Similarly, a discussion and comparison to the existing literature is missing. What sets the contributions of this paper apart from the existing literature?
-	Assumptions (A1) to (A3): The assumptions are neither explained nor is there a discussion of how realistic or restrictive they are.
-	Proposition 1, Proposition 2, Lemma 1: Like the assumptions, the theoretical results are just stated but not discussed or explained. This presentation style makes it very hard to follow the train of thought in this paper.
-	Section 2.2: It would be helpful for first time readers to explain why the dual problem can be useful for solving these types of problems. Just stating that it is “needed for the algorithm” (line 140) does not provide any intuitive insight.
-	Section 2.3: In this section it is hard for me to understand the algorithmic contributions of the paper. If the contribution, compared to the existing Sinkhorn algorithm, is just the update of $\zeta^k$, the authors should explain when and why this update makes an important difference.
-	Section 3: I am not very familiar with the EV charging literature, but I wonder if there are not any existing papers that focus on similar use cases. Since there is not a single reference in Section 3, it seems like this model is completely new and has no connections to existing work. Is this really the case?
-	Section 4 (like the previous concern): Aren’t there any existing methods to compare against? What exactly are the advantages of the proposed approach?


Minor Comments:

-	Line 35: Is it supposed to be “… common state space $\mathcal{X}$ …’’? How is this state space defined? Are there any restrictions on $\mathcal{X}$ or is it completely arbitrary? (Line 101 seems to contain the precise definition)
-	Line 36: the extensive mathematical definitions should appear later in the paper, but not in first paragraph of the introduction. That aside, I think that $\mu_1$ and $\mu_2$ are not properly defined here.
-	Line 58: What values can the variables $S_k$ and $W_k$ take?
-	Line 73: space missing at “… $S$.It“
-	Lines 79-82: Although the sentence “Inspired by … of optimal control solutions” is somewhat vague, it is nevertheless informative about the goals of this paper. I think it needs to appear earlier in the introduction.
-	Line 103: Are the marginals defined correctly? Shouldn’t the two $dx$ for the first marginal just be $x$? The same question applies to the second marginal and $dy$.
-	Line 104: there is one word “problem” too many.
-	Line 104: How are the probability kernels $T^\lambda$ in this family defined?
-	Line 118: I think the notation is wrong here. It should be “… \leq 0 for all 1 \leq m \leq M}” instead of “… \leq 0: 1 \leq m \leq M}”
-	Lines 119-120: While I do understand that an equality can be equivalently written as two inequalities, I am unsure how this applies to the previously defined moment class. Does that mean that for equalities, we would define a moment class with more than $M$ inequalities?
-	Line 383: Why is a quadratic penalization chosen? If it is standard in the literature, corresponding references should be added.
-	Line 397: Why is the infinite norm a good candidate?
-	Line 419: It should be “(ii)” instead of “(i)”, right?

Meta-review: While the reviewers did not find any disqualifying technical errors for the submission, and are appreciative of its general contributions, its connections with machine learning remain tenuous, and the authors have not sufficiently endeavored to make the connection to problems involving distribution matching or constraints that arise in ML problems.

Decision: Reject