A Sinkhorn-type Algorithm for Constrained Optimal Transport

## Abstract

Abstract Entropic optimal transport (OT) and the Sinkhorn algorithm have made it practical for machine learning practitioners to perform the fundamental task of calculating transport distance between statistical distributions. In this work, we focus on a general class of OT problems under a combination of equality and inequality constraints. We derive the corresponding entropy regularization formulation and introduce a Sinkhorn-type algorithm for such constrained OT problems supported by theoretical guarantees. We first bound the approximation error when solving the problem through entropic regularization, which reduces exponentially with the increase of the regularization parameter. Furthermore, we prove a sublinear first-order convergence rate of the proposed Sinkhorn-type algorithm in the dual space by characterizing the optimization procedure with a Lyapunov function. To achieve fast and higher-order convergence under weak entropy regularization, we augment the Sinkhorn-type algorithm with dynamic regularization scheduling and second-order acceleration. Overall, this work systematically combines recent theoretical and numerical advances in entropic optimal transport with the constrained case, allowing practitioners to derive approximate transport plans in complex scenarios.

## 1 Introduction

Obtaining the optimal transport (OT) (Villani et al., ; Linial et al., ; Peyré et al., ) plan between statistical distributions is an important subroutine in machine learning (Sandler and Lindenbaum, ; Jitkrittum et al., ; Arjovsky et al., ; Salimans et al., ; Genevay et al., ; Chen et al., ; Fatras et al., ). In this work, we focus on an optimal transportation problem with a combination of inequality and equality constraints. A typical example is an OT problem with one inequality constraint:

$$ $\min_{\begin{subarray}{c}P:P\bm{1}=r,P^{\top}\bm{1}=c,P\geq 0\end{subarray}}C \cdot P,\,\text{such that $D\cdot P\geq t$},$ (1) $$

where $\cdot$ stand for entry-wise inner product, $C\in\mathbb{R}^{n\times n}$ is the cost matrix, $t\in\mathbb{R},D\in\mathbb{R}^{n\times n}$ encode the inequality constraint, and $c,r\in\mathbb{R}^{n}$ are respectively the source and target density. As an illustration, in Figure 1, we plot the optimal transport plan between 1D distributions when the main cost is the transport cost induced by the $l_{1}$ Manhattan distance, but an inequality constraint is placed on the transport cost induced by Euclidean distance. One can see that solving inequality-constrained OT problems in equation 1 allows one to obtain transport maps under more complex geometric structures.

For the unconstrained OT problem, the most important recent breakthrough is the introduction of entropic regularization and the resultant Sinkhorn algorithm (Yule, ; Sinkhorn, ; Cuturi, ). With simple matrix scaling steps, the Sinkhorn algorithm gives an approximate OT solution in near-linear time (Altschuler et al., ), which fuels the wide adoption in the machine learning community. This prompts the following natural question:

Is there an extension of the Sinkhorn algorithm to optimal transport problems with constraints?

This paper answers the above question in the affirmative. The key insight underlying the development in this work is that the entropic optimal transport problem is an instance of entropic linear programming (Fang, ), from which it is natural to extend entropy regularization to the constrained case. For an illustration, the example in equation 1 has the following relaxation:

where $s$ is the slack variable, $\eta>0$ is the entropy regularization parameter, and $H(P,s):=\sum_{ij}p_{ij}\log{p_{ij}}+s\log(s)$ is the entropy regularization term. The primal-dual form of equation 2 leads to a direct generalization of the Sinkhorn algorithm to the constrained case, which we develop in Section 2 and Section 3.

Figure: Figure 1: Illustration of 1D optimal transport under different inequality constraints. By incrementally decreasing the upper bound for the Euclidean distance transport cost, the transport plan evolves from minimizing the Euclidean distance transport cost to minimizing the Manhattan distance ($l_{1}$) transport cost. The transport plan is computed with the Sinkhorn-type procedure in Algorithm 1.
Refer to caption: x1.png

### 1.1 Contributions

In addition to the introduction of a novel Sinkhorn-type algorithm for constrained OT, this paper systematically generalizes existing theoretical results as well as practical acceleration techniques from the existing entropic OT literature. In particular,

- •
As a corollary to the analysis in Weed (), we show that entropic optimal transport in the constrained case is exponentially close to the optimal solution.
- •
By extending the result in Altschuler et al. (), we show that a Sinkhorn-type algorithm reaches approximate first-order stationarity in polynomial time.
- •
We further improve practical performance by introducing dynamic regularization scheduling in Chen et al. () to the Sinkhorn-type algorithm.
- •
We extend the approximate sparsity analysis in Tang et al. () to the constrained OT case and introduce an accelerated second-order algorithm through sparse Newton iterations.
- •
We introduce a variational formulation of the task, which naturally enables the use of first-order methods such as in Dvurechensky et al. ().

### 1.2 Related literature

#### Constrained optimal transport

The numerical treatment for constrained OT focuses on special cases such as capacity constraints (Korman and McCann, , ), multi-marginal transport (Gangbo and Swiech, ; Buttazzo et al., ; Benamou et al., ; Pass, ; Khoo et al., ), and martingale optimal transport (Tan and Touzi, ; Beiglböck et al., ; Galichon et al., ; Dolinsky and Soner, ; Guo and Obłój, ). In terms of numerical methods for constrained OT, Benamou et al. () defines an iterative Bregman projection approach for equality constraints and a Bregman–Dykstra iteration for general constraints. In Appendix E, we show that the iterative Bregman projection is equivalent to our proposed update step in the equality case. Our contribution in relation to Benamou et al. () is threefold. First, for equality constraints, we propose a generalized Sinkhorn-type algorithm that implements the iterative Bregman projection, whereas previous work only defines the projection step in special cases. Second, for inequality constraint, our approach is novel and enjoys theoretical guarantee as seen in Theorem 1 and Theorem 3. Finally, the Bregman–Dykstra iteration in (Benamou et al., ) is only implemented when a simple closed-form solution exists for each projection step, whereas our proposed method applies to general cases.

#### OT in machine learning

A substantial amount of research literature exists on the use of optimal transport in different areas of machine learning. Notable applications include statistical learning (Kolouri et al., ; Vayer et al., ; Genevay et al., ; Luise et al., ; Oneto et al., ; Huynh et al., ), domain adaptation (Fernando et al., ; Redko et al., ; Courty et al., ; Alvarez-Melis et al., ; Nguyen et al., , ; Xu et al., ; Turrisi et al., ), and using optimal transport distance in designing training targets (Genevay et al., ; Bousquet et al., ; Sanjabi et al., ; Deshpande et al., ; Lei et al., ; Patrini et al., ; Onken et al., ). The usage of entropic optimal transport in the constrained case allows practitioners to query transport plans which satisfy more complex structure than in the unconstrained case, which are beneficial across all the applications above.

#### Acceleration for Sinkhorn

There is a considerable body of work in speeding up the runtime of the Sinkhorn algorithm, and a significant portion of our numerical treatment is devoted to the extension of these algorithms to the constrained case. In addition to (Chen et al., ; Tang et al., ) to be discussed in Section 3, we mention a few noteworthy methods for acceleration.
Randomized or greedy row/column scaling (Genevay et al., ; Altschuler et al., ) can be directly extended to the constrained case by the inclusion of the dual variables for the linear constraints. Likewise, methods based on the variational form of entropic OT, such as Nesterov acceleration (Dvurechensky et al., ) and mirror descent (Lin et al., ; Kemertas et al., ), can be used in the constrained case by considering the Lyapunov function introduced in Section 2. While the majority of existing techniques can be directly extended to the constrained case, one notable exception is the class of methods based on approximation to the kernel $K=\exp(-C\eta)$ (Deriche, ; Solomon et al., ; Bonneel et al., ; Altschuler et al., ; Scetbon and Cuturi, ; Lakshmanan et al., ; Huguet et al., ; Li et al., ). The main obstruction from performing kernel approximation in the constrained case is that the kernel is changing throughout the optimization process due to update in the constraint associated dual variables. As the kernel compression step is costly and often done in the offline stage, it poses a significant challenge if one has to update the kernel approximation with each dual variable update dynamically.

## 2 Theoretical foundations of constrained entropic optimal transport

This section summarizes the entropic optimal transport in the constrained case.

#### Notations

We use $M\cdot M^{\prime}:=\sum_{ij}m_{ij}m^{\prime}_{ij}$ to denote the entry-wise inner product. For a matrix $M$, the notation $\log{M}$ stands for entry-wise logarithm, and similarly $\exp(M)$ denotes entry-wise exponential. We use the symbol $\lVert M\rVert_{1}$ to denote the entry-wise $l_{1}$ norm, i.e. $\lVert M\rVert_{1}:=\lVert\mathrm{vec}(M)\rVert=\sum_{ij}\lvert m_{ij}\rvert$. The $\lVert M\rVert_{\infty}$ and $\lVert M\rVert_{2}$ norms are defined likewise as the entry-wise $l_{\infty}$ and $l_{2}$ norms, respectively. The notation $\mathbf{1}_{n\times n}$ is the all-one $n\times n$ matrix, and the notation $\mathbf{1}$ denotes the all-one vector of appropriate size.

#### Background

For simplicity, we assume the target and source density $r,c\in\mathbb{R}^{n}$ satisfies $\sum_{i}r_{i}=\sum_{j}c_{j}=1$.
For any inequality constraint of the form $D\cdot P\geq t$, note one can convert the condition to $(D-t\mathbf{1}_{n\times n})\cdot P\geq 0$, and similarly in the equality case. Let $K$ and $L$ respectively denote the number of inequality and equality constraints. By the construction above, for each optimal transport problem under linear constraint, there exists $D_{1},\ldots D_{K},D_{K+1},\ldots,D_{K+L}\in\mathbb{R}^{n\times n}$ so that the linear constraints are encoded by $\mathcal{I},\mathcal{E}$ where

$$ $\displaystyle\mathcal{I}:=$ $\displaystyle\bigcap_{k=1,\ldots,K}\{P\mid D_{k}\cdot P\geq 0\},$ (3) $\displaystyle\mathcal{E}:=$ $\displaystyle\bigcap_{l=1,\ldots,L}\{P\mid D_{l+K}\cdot P=0\}.$ $$

We summarize the *general form of constrained optimal transport* by the following the following linear program (LP):

$$ $\min_{\begin{subarray}{c}P:P\bm{1}=r,P^{\top}\bm{1}=c,P\geq 0\end{subarray}}C \cdot P,\,\text{such that $P\in\mathcal{S}$},$ (4) $$

where $\mathcal{S}:=\mathcal{I}\cap\mathcal{E}$ for $\mathcal{I},\mathcal{E}$ defined in equation 3.

#### Entropic optimal transport under constraint

Under entropic linear programming (Fang, ), one can write down the formulation for entropic optimal transport under general linear constraints. For $k=1,\ldots,K$, we define $s_{k}$ to be the slack variable corresponding to $D_{k}\cdot P$. The
*constrained entropic optimal transport* follows from the following equation:

$$ $\displaystyle\min_{\begin{subarray}{c}P,s:P\bm{1}=r,\,P^{\top}\bm{1}=c,\,\\ P\geq 0,\,s_{k}\geq 0\,\forall\,k\end{subarray}}$ $\displaystyle C\cdot P+\frac{1}{\eta}H(P,s_{1},\ldots,s_{K}),$ (5) subject to $\displaystyle D_{k}\cdot P=s_{k}\text{ for $k=1,\ldots,K$.}$ $\displaystyle D_{l+K}\cdot P=0\text{ for $l=1,\ldots,L$,}$ $$

where the entropy term is defined by $H(P,s_{1},\ldots,s_{K})=\sum_{ij}p_{ij}\log(p_{ij})+\sum_{k=1}^{K}s_{k}\log s_
{k}$.

We motivate the optimization task in equation 5 by Theorem 1, which shows that the entropy-regularized optimal solution is exponentially close to the optimal solution:

###### Theorem 1 .

For simplicity, assume that $\sum_{i}r_{i}=\sum_{j}c_{j}=1$, the LP in equation 4 has a unique solution $P^{\star}$, and assume that $\lVert D_{k}\rVert_{\infty}\leq 1$ for $k=1,\ldots,K$. Denote $P^{\star}_{\eta}$ as the unique solution to equation 5. There exists a constant $\Delta$, depending only on the LP in equation 4, so that the following holds for $\eta\geq\frac{(K+1)(1+\ln(4n^{2}(K+1)))}{\Delta}$:

$$ $\lVert P^{\star}_{\eta}-P^{\star}\rVert_{1}\leq 8n^{\frac{2}{K+1}}(K+1)\exp \left(-\eta\frac{\Delta}{K+1}\right).$ $$

The definition of $\Delta$ and the proof are deferred to Appendix A. In this work, we only consider examples for which $K=O(1)$, in which case the bound in Theorem 1 does not significantly differ from the unconstrained case. It should be noted that the dependency on $n$ goes down with an increasing number of inequality constraints $K$, but one can check the bound is monotonically worse as $K$ increases. Interpreting the bound could mean that adding more inequality constraints limits the degree of freedom but may also potentially amplify the error caused by entropy regularization.

###### Theorem 1 .

#### Variational formulation under entropic regularization

By introducing the Lagrangian variable and using the minimax theorem (for a detailed derivation, see Appendix D), we formulate the associated primal-dual problem to equation 5 as:

$$ $\displaystyle\max_{x,y,a}\min_{P,s}\big{[}L(P,s,x,y,a)$ $\displaystyle:=$ $\displaystyle\frac{1}{\eta}P\cdot\log P+C\cdot P-x\cdot(P\bm{1}-r)-y\cdot(P^{ \top}\bm{1}-c)$ $\displaystyle+$ $\displaystyle\frac{1}{\eta}\sum_{k=1}^{K}s_{k}\log{s_{k}}+\sum_{k=1}^{K}a_{k}s _{k}-\sum_{m=1}^{K+L}a_{m}(D_{m}\cdot P)\big{]},$ $$

where $s=(s_{1},\ldots,s_{K})$ is the shorthand for slack variable, $a=(a_{1},\ldots,a_{K},a_{K+1},\ldots,a_{K+L})$ is the shorthand for the constraint dual variables (excluding the original row/column constraints).

By eliminating $P,s$ (see Appendix D), the function $f(x,y,a):=\min_{P,s}L(P,s,x,y,a)$ admits the following form

$$ $\displaystyle f(x,y,a)$ (6) $\displaystyle=$ $\displaystyle\frac{-1}{\eta}\sum_{ij}\exp{\Big{(}\eta(-C_{ij}+\sum_{m=1}^{L+K} a_{m}\left(D_{m}\right)_{ij}+x_{i}+y_{j})-1\Big{)}}$ $\displaystyle+$ $\displaystyle\sum_{i}x_{i}r_{i}+\sum_{j}y_{j}c_{j}-\frac{1}{\eta}\sum_{k=1}^{K }{\exp(-\eta a_{k}-1)}.$ $$

As a consequence of the minimax theorem, maximizing over $f$ is equivalent to solving the problem defined in equation 5. We emphasize that $f$ is *concave*, allowing one to use routine convex optimization techniques.

Define $P=\exp{\big{(}\eta(-C+\sum_{m}a_{m}D_{m}+x\bm{1}^{\top}+\bm{1}y^{\top})-1\big{
)}}$ as the intermediate matrix corresponding to dual variables $x,y,a$. We write down the first derivative of the Lyapunov function $f$ :

$$ $\displaystyle\nabla_{x}f=r-P\mathbf{1},\quad\nabla_{y}f=c-P^{\top}\mathbf{1},$ (7) $\displaystyle\partial_{a_{k}}f=\exp(-\eta a_{k}-1)-P\cdot D_{k},\ \ \forall k \in[K],$ $\displaystyle\partial_{a_{l+K}}f=-P\cdot D_{l+K},\ \ \forall l\in[L]\,.$ $$

One can use any first-order method on $f$ using the derivatives provided in equation 7. For example, one option with a desirable performance guarantee would be accelerated gradient ascent (Nesterov, ; Dvurechensky et al., ). Overall, the formulation in equation 6 shows that the introduction of constraint does not pose a substantial change when compared to the optimization task in the unconstrained case.

## 3 Main algorithm

This section proposes an efficient Sinkhorn-type algorithm in the constrained case. We assume $K+L=O(1)$ to ensure the efficiency of the proposed approach.

#### Main idea of the algorithm

The variational formulation under entropic regularization shows that one can effectively solve for the entropic formulation in equation 5 by solving for the optimization task:

$$ $\max_{x\in\mathbb{R}^{n},y\in\mathbb{R}^{n},a\in\mathbb{R}^{K+L}}f(x,y,a)$ $$

for $f$ given in equation 6. We define a Sinkhorn-type algorithm by introducing three iteration steps:

- 1.
($x$ update) $x\leftarrow\operatorname*{arg\,max}_{\tilde{x}}f(\tilde{x},y,a)$,
- 2.
($y$ update) $y\leftarrow\operatorname*{arg\,max}_{\tilde{y}}f(x,\tilde{y},a)$,
- 3.
($a$ update) $a,t\leftarrow\operatorname*{arg\,max}_{\tilde{a},\tilde{t}}f(x+\tilde{t}
\mathbf{1},y,\tilde{a})$

The proposed alternating update approach is summarized in Algorithm 1.

#### Implementation detail

Let $(x,y,a)$ be the current dual variables and define

$$ $P=\exp{\left(\eta(-C+\sum_{m}a_{m}D_{m}+x\bm{1}^{\top}+\bm{1}y^{\top})-1\right)}$ $$

as the current intermediate transport plan defined by these dual variables.
By equation 7, one has $\nabla_{x}f=0\iff P\mathbf{1}=r$ and $\nabla_{y}f=0\iff P^{\top}\mathbf{1}=c$. Thus, the $x,y$ update steps amount to row/column scaling of the matrix $P$, which is identical to the Sinkhorn algorithm.

The $a$ update step constitutes the main novelty of our algorithm in the constrained case. In particular, the inclusion of the $t$ variable improves numerical stability by enforcing a normalization condition: Suppose $P$ is the intermediate transport plan formed by the dual variable $(x+t\mathbf{1},y,a)$, then

$$ $\partial_{t}f(x+t\mathbf{1},y,a)=0\iff\sum_{ij}P_{ij}=\sum_{i}r_{i},$ (8) $$

which ensures that terms of the form $M\cdot P$ can be bounded by $\lVert M\rVert_{\infty}\left(\sum_{i}r_{i}\right)$. Thus, optimality in the $t$ variable ensures boundedness in the derivatives of $f$ such as those in equation 7.

We propose to use Newton’s method for the $a$ update step. Namely, by directly computing the gradient and the Hessian, $\nabla_{at}f,\nabla_{at}^{2}f$, one uses the search direction $(\Delta a,\Delta t)=-\left(\nabla^{2}_{at}f\right)^{-1}\nabla_{at}f$. The learning rate is obtained through the standard backtracking line search scheme (Boyd and Vandenberghe, ).

Figure: Algorithm 1 Sinkhorn-type algorithm under linear constraint

#### Complexity analysis of Algorithm 1

The row/column scaling step is identical to the Sinkhorn algorithm and thus costs $O(n^{2})$ per iteration. For the constraint dual update step, the computation cost is dominated by the calculation of $\nabla_{at}^{2}f$. For the formula of the Hessian, one can refer to equation 37 in Appendix F and see that the cost is dominated by the computation of $\sum_{ij}P_{ij}\left(D_{m}\right)_{ij}\left(D_{m^{\prime}}\right)_{ij}$ for $m,m^{\prime}=1,\ldots,K+L$, which leads to a computation complexity of $O(N_{a}(K+L)^{2}n^{2}))$, where $N_{a}$ is the number of iterations taken in the Newton step. The computation complexity in the constraint dual update step is $O(n^{2})$ for two reasons. First, this algorithm assumes $K+L=O(1)$. Second, Newton’s method enjoys super-exponential convergence practically (Boyd and Vandenberghe, ), and thus one can set $N_{a}=O(1)$, which is usually sufficient to reach machine precision.

### 3.1 Acceleration techniques

#### Entropy regularization scheduling

An important feature of the Sinkhorn algorithm is that the iteration complexity heavily depends on the entropy regularization term $\eta$, the tuning of which plays a significant part in practical performance. To that end, one can aid acceleration by using the doubling entropy regularization scheduling technique introduced in (Chen et al., ). For a desired entropy regularization value $\eta_{\mathrm{final}}$, we take an initial regularization strength $\eta_{\mathrm{init}}$ and take $N_{\eta}=\lceil\mathrm{log}_{2}(\eta_{\mathrm{final}}/\eta_{\mathrm{init}})\rceil$. Then, one defines successively doubling regularization levels $\eta_{0},\ldots,\eta_{N_{\eta}}$ so that $\eta_{\mathrm{init}}=\eta_{0}<\ldots<\eta_{N_{\eta}}=\eta_{\mathrm{final}}$ and $\eta_{i}=2\eta_{i-1}$ for $i=1,\ldots,N_{\eta}-1$. For each step $i$, one runs the subroutine in Algorithm 1 at $\eta_{i-1}$, and the obtained dual variable is used as the initialization when calling the subroutine for $\eta_{i}$.

#### Second-order acceleration through sparsity

We further accelerate the Sinkhorn-type algorithm with Sinkhorn-Newton-Sparse (SNS) (Tang et al., ), a second-order method originally developed for the unconstrained case. Instead of alternating maximization as in Algorithm 1, one can instead introduce the combined variable $z=(x,y,a)$, where a naive strategy would be to optimize directly $f$ through a full Newton update of the type $z=z-\alpha(\nabla_{z}^{2}f)^{-1}\nabla_{z}f$, where $\alpha$ can be obtained through backtracking line search. However, a full Newton step would require an impractical $O((n+K+L)^{3})$ scaling.

The SNS algorithm introduces a practical second-order method through a sparse approximation of the Hessian matrix, and we show how one can extend the algorithm to the constrained OT case. The key observation of SNS is that the Hessian submatrix corresponding to variable $x,y$ has a special structure:

$$ $\nabla_{xy}^{2}f(x,y,a)=\eta\begin{bmatrix}\mathrm{diag}(P\bm{1})&P\\ P^{\top}&\mathrm{diag}(P^{\top}\bm{1})\end{bmatrix},$ (9) $$

which admits a sparse approximation as long as $P$ admits a sparse approximation. Moreover, the full Hessian matrix to $f$ admits a sparse approximation: As $K+L=O(1)$, the blocks of $\nabla^{2}f$ corresponding to $\nabla_{xy}\nabla_{a}f$ and $\nabla_{a}^{2}f$ lead to at most $O(n)$ nonzero entries to keep track of.

The rationale for approximate sparsity is simple: Under a mild uniqueness assumption, the optimal solution $P^{\star}$ to the LP in equation 4 has at most $2n-1+K+L=O(n)$ nonzero entries due to the fundamental theorem of linear programming (Luenberger et al., ).
Moreover, the exponential closeness result in Theorem 1 implies that the entries of $P^{\star}_{\eta}$ decay at a rate of $\exp{(-\frac{\eta\Delta}{K+1})}$ except at $2n-1+K+L=O(n)$ entries, which proves that approximate sparsity holds in the constrained case as well.

We propose an extended SNS algorithm by sparsifying the $\nabla_{xy}^{2}f$ block in which one keeps only $O(n)$ nonzero entries. The sparse Newton step is used in combination with Algorithm 1 as a warm start to achieve rapid acceleration. The proposed scheme leads to the same $O(n^{2})$ per-iteration complexity as SNS in the unconstrained case. Implementation details can be found in Algorithm 2 in Appendix F.

## 4 Convergence analysis

In this section, we present a convergence bound for a modified version of the proposed Sinkhorn-type algorithm in terms of the stationarity condition. Similar to the proof in (Altschuler et al., ) for the unconstrained case, our proof strategy relies on characterizing convergence through the Lyapunov function $f$ introduced in equation 6. For the rest of this section, we assume $\sum_{i}r_{i}=\sum_{j}c_{j}=1$.

The following Theorem characterizes the improvement of the Lyapunov function from one step of the proposed algorithm (proof is in Appendix B):

###### Theorem 2 .

Let $(x,y,a)$ be the current dual variable which has undergone at least one update step, and let $P$ be the associated transport matrix, i.e.,

$$ $P=\exp{\left(\eta(-C+\sum_{m}a_{m}D_{m}+x\bm{1}^{\top}+\bm{1}y^{\top})-1\right )}.$ $$

For the $x$ update step satisfying $x^{\prime}=\operatorname*{arg\,max}_{\tilde{x}}f(\tilde{x},y,a)$, one has

$$ $f(x^{\prime},y,a)-f(x,y,a)=\mathrm{KL}\left(r||P\mathbf{1}\right),$ (10) $$

where $\mathrm{KL}\left(\cdot||\cdot\right)$ stands for the Kullback-Leibler divergence.

For the $y$ update step satisfying $y^{\prime}=\operatorname*{arg\,max}_{\tilde{y}}f(x,\tilde{y},a)$, one has

$$ $f(x,y^{\prime},a)-f(x,y,a)=\mathrm{KL}\left(c||P^{\top}\mathbf{1}\right).$ (11) $$

Consider the $a$ update step satisfying $(a^{\prime},t^{\prime})=\operatorname*{arg\,max}_{\tilde{a},\tilde{t}}f(x+
\tilde{t}\mathbf{1},y,\tilde{a})$. For $m=1,\ldots,K+L$, define $d_{m}$ as the optimality condition corresponding to $a_{m}$, i.e.

$$ $\displaystyle d_{k}:=\exp(-\eta a_{k}-1)-P\cdot D_{k},$ $\displaystyle k=1,\ldots,K$ (12) $\displaystyle d_{l+K}:=-P\cdot D_{l+K},$ $\displaystyle l=1,\ldots,L.$ $$

Then, for $c_{d}:=\max_{m\in[K+L]}\lVert D_{m}\rVert_{\infty}$, one has

$$ $\displaystyle f(x+t^{\prime}\mathbf{1},y,a^{\prime})-f(x,y,a)$ (13) $\displaystyle\geq$ $\displaystyle\sum_{k=1}^{K}\lvert d_{k}\rvert\min\left(\frac{1}{8\eta},\frac{ \lvert d_{k}\rvert}{8\eta c_{d}+4\eta(K+L)c_{d}^{2}}\right)$ $\displaystyle+$ $\displaystyle\sum_{l=1}^{L}\frac{d_{l+K}^{2}}{2\eta(K+L)c_{d}^{2}}.$ $$

We now prove a bound on the first-order stationarity condition.
As one is performing alternating optimization for more than two variables, the convergence proof requires a modification to Algorithm 1, in which one examines the stationary condition on each variable and chooses the update step greedily. The following Theorem characterizes the approximate stationarity of the greedy version of Algorithm 1 (proof is in Appendix C):

###### Theorem 3 .

Let $(x,y,a)$ be the current variable, let $P$ be the associated transport matrix, and let $d_{m}$ be as in equation 12. Define $Q_{x},Q_{y},Q_{a}$ by

$$ $\displaystyle Q_{x}=$ $\displaystyle\mathrm{KL}\left(r||P\mathbf{1}\right),Q_{y}=\mathrm{KL}\left(c|| P^{\top}\mathbf{1}\right),$ $\displaystyle Q_{a}=$ $\displaystyle\sum_{k=1}^{K}\lvert d_{k}\rvert\min\left(\frac{1}{8\eta},\frac{ \lvert d_{k}\rvert}{8\eta c_{d}+4\eta(K+L)c_{d}^{2}}\right)$ $\displaystyle+$ $\displaystyle\sum_{l=1}^{L}\frac{d_{l+K}^{2}}{2\eta(K+L)c_{d}^{2}},$ $$

which are respectively the right-hand-side of equation 10, equation 11, and equation 13.

Consider a greedy version of Algorithm 1, in which only one update step is performed at each iteration, and an $x$ (resp. $y$, $a$) update step is chosen if $Q_{x}$ (resp. $Q_{y}$, $Q_{a}$) is the largest among $(Q_{x},Q_{y},Q_{a})$. Let $c_{d}$ be as in Theorem 2, and define $c_{g}$ as the gap term in the Lyapunov function $f$ at initialization, i.e.

$$ $c_{g}=\max_{\tilde{x},\tilde{y},\tilde{a}}f(\tilde{x},\tilde{y},\tilde{a})-f(x _{\mathrm{init}},y_{\mathrm{init}},a_{\mathrm{init}}).$ $$

Let $(x^{u},y^{u},a^{u})$ be the dual variable after $u$ iterations and let $P_{u}$ be the associated intermediate transport plan. Then, the greedy algorithm outputs dual variables $(x,y,a)$ which satisfies the below approximate stationarity condition in $O(c_{g}\epsilon^{-2})$ iterations:

$$ $\lVert\nabla f(x,y,a)\rVert_{1}\leq\epsilon.$ (14) $$

Moreover, the matrix $P$ associated with the outputted dual variable satisfies the following condition:

$$ $\displaystyle\epsilon\geq$ $\displaystyle\lVert P\mathbf{1}-r\rVert_{1}+\lVert P^{\top}\mathbf{1}-c\rVert_ {1}$ (15) $\displaystyle+$ $\displaystyle\sum_{k=1}^{K}\lvert\min\left(P\cdot D_{k},0\right)\rvert+\sum_{l =1}^{L}\lvert P\cdot D_{l+K}\rvert,$ $$

which shows that $P$ approximately satisfies the linear constraints in equation 4.

###### Theorem 2 .

###### Theorem 3 .

## 5 Numerical experiments

Figure: Figure 2: Random assignment problem. Plot of the proposed Sinkhorn-type algorithm in terms of assignment cost and constraint violation. Specifically, constraint violation is defined by $\mathrm{Violation}(P)=|\min(0,D_{I}\cdot\mathrm{Round}(P,\mathcal{U}_{r,c}))|+ |D_{E}\cdot\mathrm{Round}(P,\mathcal{U}_{r,c})|$.
Refer to caption: x2.png

We conduct numerical experiments to showcase the performance of the proposed Sinkhorn-type algorithm and its acceleration techniques. Let $\mathcal{U}_{r,c}$ be the set of transport matrices from $r$ to $c$. We use the rounding algorithm in Altschuler et al. () to obtain projection into $\mathcal{U}_{r,c}$. The performance is evaluated through the cost and constraint violation of the transport matrix one obtains through rounding. Specifically, one first uses the
dual variable $(x,y,a)$ to form the intermediate matrix $P$. Then, one uses the rounding algorithm to get a transport matrix, denoted $\mathrm{Round}(P,\mathcal{U}_{r,c})$. The cost or score of the transport is evaluated as

$$ $\displaystyle\mathrm{Cost}(P)=C\cdot\mathrm{Round}(P,\mathcal{U}_{r,c}),$ $\displaystyle\mathrm{Score}(P)=(-C)\cdot\mathrm{Round}(P,\mathcal{U}_{r,c}).$ $$

For constraint violation, we use the following metric:

$$ $\displaystyle\mathrm{Violation}(P)=$ $\displaystyle\sum_{k=1}^{K}\lvert\min\left(\mathrm{Round}(P,\mathcal{U}_{r,c}) \cdot D_{k},0\right)\rvert$ $\displaystyle+$ $\displaystyle\sum_{l=1}^{L}\lvert\mathrm{Round}(P,\mathcal{U}_{r,c})\cdot D_{l +K}\rvert.$ $$

In Proposition 1 in Appendix C, we give an upper bound on $\mathrm{Violation}(P)$ when one runs the greedy version of Algorithm 1 in Theorem 3.

#### Random assignment problem under constraints

In the first numerical test, we consider the random assignment problem (Mézard and Parisi, ; Steele, ; Aldous, ) with additional inequality and equality constraints. In this setting, we set the problem size of $n=500$ and an entropy regularization of $\eta=1200$. The source and target vectors are $c=r=\frac{1}{n}\mathbf{1}$. We consider three $n\times n$ matrices $C,D_{I},D_{E}$, which respectively encode the cost, the inequality constraint, and the equality constraint. We generate the entries of $C,D_{I},D_{E}$ by i.i.d. random variables following the distribution $\mathrm{Unif}([0,1])$. We then set two threshold variables $t_{I},t_{E}$ and consider the following optimal transport task:

$$ $\displaystyle\min_{\begin{subarray}{c}P:P\bm{1}=r,\,P^{\top}\bm{1}=c\end{ subarray}}$ $\displaystyle C\cdot P,$ (16) subject to $\displaystyle D_{I}\cdot P\leq t_{I},$ $\displaystyle D_{E}\cdot P=t_{E}.$ $$

The conversion of equation 16 to the general form (see equation 4) can be performed by the procedure detailed in Section 2.
Subsequently, we perform the proposed modified Sinkhorn-type algorithm detailed in Algorithm 1, as well as the accelerated Sinkhorn-Newton-Sparse algorithm in Algorithm 2 (see Appendix F). We set the threshold parameters to be $t_{I}=t_{E}=\frac{1}{2}$, and it has been verified that the tested instances of equation 16 is feasible. We test the performance for Algorithm 1 and Algorithm 2. For Algorithm 2, we use $N_{1}=20$ Sinkhorn steps for initialization.

In Figure 3, we plot the performance of Algorithm 1 and Algorithm 2.
One can see that the proposed modified Sinkhorn-type algorithm quickly converges to an approximately optimal solution that approximately satisfies the additional constraints. Algorithm 2 achieves similar performance under a vastly smaller number of iterations. We remark that Algorithm 2 in this case has converged to the optimal solution $P_{\eta}^{\star}$ with machine accuracy in $N_{1}+N_{2}=25$ iterations. We note that $P_{\eta}^{\star}$ is obtained by using Newton’s method, i.e., running Algorithm 2 without sparsification, and the stationarity condition has been verified to be satisfied under machine-epsilon accuracy. In Appendix G, we show that Algorithm 1 and Algorithm 2 both converge rapidly to the entropic optimal solution $P_{\eta}^{\star}$.

#### Ranking under constraints

In the second numerical example, we consider cost matrices associated with evaluation metrics in ranking problems (Liu et al., ; Manning, ). The problem of ranking can be naturally perceived as a transport problem where the transport matrix is analogous to the permutations of items (Mena et al., ). All linear additive ranking metrics, such as precision, recall, and discounted cumulative gain (DCG), are linear functions of the permutation matrix. For example, the DCG metric measures the quality of a ranking $\sigma:[n]\to[n]$. For a relevance vector $g\in\mathbb{R}^{n}$ and a discount vector $v\in\mathbb{R}^{n}$, a DCG score is defined by the following:

$$ $\mathrm{DCG}_{v,g}(\sigma)=\sum_{i=1}^{n}g_{\sigma(i)}v_{i}=(gv^{\top})\cdot P$ $$

where $P$ is the permutation matrix associated with $\sigma$. Thus, optimizing the DCG metric corresponds to an OT instance through the relaxation of the permutation matrix within the Birkhoff polytope.

In practical problems such as e-commerce ranking, there are typically multiple relevance labels in the form of different attributes of an item. Motivated by this, we consider the following constrained OT problem:

$$ $\displaystyle\max_{\begin{subarray}{c}P:P\bm{1}=\mathbf{1},\,P^{\top}\bm{1}= \mathbf{1}\end{subarray}}$ $\displaystyle(g_{c}v^{\top})\cdot P,$ (17) subject to $\displaystyle(g_{I}v^{\top})\cdot P\geq t_{I},$ $\displaystyle(g_{E}v^{\top})\cdot P=t_{E},$ $$

where $g_{c},g_{I},g_{E}$ are three relevance vectors associated with the cost, the inequality constraint, and the equality constraint, respectively. In accordance with the information retrieval literature, we consider the discount vector $v$ with $v_{i}=\frac{1}{\log_{2}(i+1)}$. The entries of $g_{c},g_{I},g_{E}$ are i.i.d. entries simulated from the Rademacher distribution. In equation 17, we choose the two threshold variables as $t_{I}=\frac{1}{n}D_{I}\cdot\mathbf{1}_{n\times n},t_{E}=\frac{1}{n}D_{E}\cdot
\mathbf{1}_{n\times n}$, so that the feasibility set is guaranteed to be non-empty.

In accordance with equation 16, we consider a problem size of $n=500$ and an entropy regularization of $\eta=\frac{1200}{500}=2.4$.
The performance is plotted in Figure 3, which shows that both Algorithm 1 and Algorithm 2 can quickly converge to an approximately optimal solution that approximately satisfies the additional constraints. Furthermore, Algorithm 2 in this case has converged to the optimal solution $P_{\eta}^{\star}$ with machine-epsilon accuracy in under $25$ iterations. Appendix G shows rapid convergence to $P_{\eta}^{\star}$ for both algorithms.

#### Pareto front for geometric transport problem under Euclidean distance and Manhattan distance

In the third numerical example, we consider the trade-off between two transport costs of a geometric nature. In this experiment, we run the task of Pareto front profiling between the Manhattan distance cost $c_{1}(x,y)=\lVert x-y\rVert_{1}$ and the Euclidean distance transport cost $c_{2}(x,y)=\lVert x-y\rVert_{2}^{2}$. As a byproduct, one obtains an interpolation between the Wasserstein $W_{1}$ transport plan under the $l_{1}$ distance and the $W_{2}$ transport plan under the $l_{2}$ distance. Let $C_{1},C_{2}$ be the cost matrix associated with the Manhattan distance and the Euclidean distance (Villani et al., ). To do so, we consider the following optimization task:

$$ $\displaystyle\min_{\begin{subarray}{c}P:P\bm{1}=r,\,P^{\top}\bm{1}=c\end{ subarray}}$ $\displaystyle C_{1}\cdot P,$ (18) subject to $\displaystyle C_{2}\cdot P\leq t^{2},$ $$

where $t$ is set so that the feasibility set is not empty. We let $t_{\mathrm{min}}$ be the $W_{2}$ distance between the source and target vectors, and let $t_{\mathrm{max}}$ be the Euclidean transport cost of a transport plan which minimizes Manhattan distance transport cost. We remark that $t_{\mathrm{min}}$ and $t_{\mathrm{max}}$ can be obtained through conventional OT algorithms. By tracing the value of $t$ in $[t_{\mathrm{min}},t_{\mathrm{max}}]$ and solving equation 18 through entropic regularization, one can effectively obtain a Pareto front between Euclidean distance transport cost and Manhattan distance transport cost.

Similar to previous cases, Algorithm 2 is able to reach optimal solution $P_{\eta}^{\star}$ within machine-epsilon accuracy. Thus, we run multiple instances of Algorithm 2 for $t\in[t_{\mathrm{min}},t_{\mathrm{max}}]$ and $\eta=10,100,1000$. Then, we plot the Pareto front formed by $P_{\eta}^{\star}$ for every choice of $\eta$. The true Pareto front can be obtained through running Algorithm 3 in Appendix F for $\eta_{\mathrm{target}}=8*10^{4}$, and convergence is verified through checking the first order derivative reaches machine accuracy.

Similar to Cuturi (), we illustrate the procedure through optimal transport on the MNIST dataset. We pick two images, which are converted to a vector of intensities on the $28\times 28$ pixel grid and normalized to sum to 1. The entry corresponding to the $(i_{1},i_{2})$-th pixel is conceptualized as the point $(i_{1}/28,i_{2}/28)\in\mathbb{R}^{2}$. In Figure 4, we see that the Pareto front formed by increasing $\eta$ indeed converges to the true Pareto front.

Figure: Figure 4: Pareto front profiling of Euclidean distance versus Manhattan distance.
Refer to caption: x4.png

## 6 Conclusion

We introduce an entropic formulation of optimal transport with a combination of additional equality and inequality constraints. We propose a Sinkhorn-type algorithm that has a novel constraint dual variable update step. We provide preliminary results on the approximation error of the entropic formulation and the convergence of the Sinkhorn-type algorithm, while more work remains to be done for the convergence property of Algorithm 1 in terms of convergence to the optimal entropic solution $P_{\eta}^{\star}$. Moreover, we provide a detailed discussion on improving the proposed Sinkhorn-type algorithm with acceleration techniques, especially Sinkhorn-Newton-Sparse and entropy regularization scheduling. The proposed work enables one to obtain approximately optimal solutions in more complicated OT instances efficiently. We contend that this work has the potential to be a vital tool in the field of optimal transport under constraint and in attracting the use of constrained optimal transport in machine learning.

## Appendix A Proof of Theorem 1

###### Definition 1 .

Let $\mathcal{S}$ be the constraint set defined in equation 4. Define $\mathcal{P}$ as the polyhedron formed by the transport matrix, i.e.

$$ $\mathcal{P}:=\{P\mid P\bm{1}=r,P^{\top}\bm{1}=c,P\geq 0,P\in\mathcal{S}\}.$ $$

The symbol $\mathcal{V}$ denotes the set of vertices of $\mathcal{P}$. The symbol $\mathcal{O}$
stands for the set of optimal vertex solutions, i.e.

$$ $\mathcal{O}:=\operatorname*{arg\,min}_{P\in\mathcal{V}}C\cdot P.$ (19) $$

The symbol $\Delta$ denotes the vertex optimality gap

$$ $\Delta=\min_{Q\in\mathcal{V}-\mathcal{O}}Q\cdot C-\min_{P\in\mathcal{O}}P\cdot C.$ $$

We can now finish the proof.

###### Proof.

This convergence result is mainly due to the application of Corollary 9 in Weed [] to this case. We define another polyhedron $\mathcal{Q}$ as follows:

$$ $\mathcal{Q}:=\{(P,s)\mid P\in\mathcal{P},\forall k\in[K],s_{k}=D_{k}\cdot P\}.$ $$

Let $P_{\eta}^{\star}$ be as defined in the statement, and for $k=1,\ldots,K$ we define $s_{\eta;k}^{\star}=D_{k}\cdot P_{\eta}^{\star}$. We use $R_{1}$ and $R_{H}$ to denote the $l_{1}$ and entropic radius of $\mathcal{Q}$ in the sense defined in Weed [].
It is easy to see that for $R_{1}$ one has

$$ $1\leq R_{1}=1+\max_{P\in\mathcal{P}}\sum_{k=1}^{K}P\cdot D_{k}\leq 1+K,$ $$

where the second inequality used Holder’s inequality and the assumption that $\lVert D_{k}\rVert_{\infty}\leq 1$

For $R_{H}$, one similarly has

$$ $\displaystyle R_{H}=$ $\displaystyle\max_{(s,P),(s^{\prime},P^{\prime})\in\mathcal{Q}}\sum_{ij}\left( p_{ij}\log(p_{ij})-p^{\prime}_{ij}\log(p^{\prime}_{ij})\right)+\sum_{k}\left(s _{k}\log(s_{k})-s^{\prime}_{k}\log(s^{\prime}_{k})\right)$ $\displaystyle\leq$ $\displaystyle\left(\max_{P,P^{\prime}\in\mathcal{P}}\sum_{ij}\left(p_{ij}\log( p_{ij})-p^{\prime}_{ij}\log(p^{\prime}_{ij})\right)\right)+\left(\max_{P,P^{ \prime}\in\mathcal{P}}\sum_{k}(P\cdot D_{k})\log(P\cdot D_{k})-(P^{\prime} \cdot D_{k})\log(P^{\prime}\cdot D_{k})\right)$ $\displaystyle\leq$ $\displaystyle\log(n^{2})+K/e,$ $$

where the second equality holds because $H(P)\in[0,\log{(n^{2})}]$ and $P\cdot D_{k}\log(P\cdot D_{k})\in[-1/e,0]$.

For $\eta\geq\frac{(K+1)(1+\ln(4n^{2}(K+1))}{\Delta}>\frac{R_{1}+R_{H}}{\Delta}$, one has

$$ $\displaystyle\lVert P^{\star}-P^{\star}_{\eta}\rVert_{1}\leq$ $\displaystyle\lVert(P^{\star},s^{\star})-(P^{\star}_{\eta},s_{\eta}^{\star}) \rVert_{1}$ $\displaystyle\leq$ $\displaystyle 2R_{1}\exp\left(-\eta\frac{\Delta}{R_{1}}+1+\frac{R_{H}}{R_{1}}\right)$ $\displaystyle=$ $\displaystyle 2R_{1}\exp\left(\frac{R_{H}-\eta\Delta}{R_{1}}+1\right)$ $\displaystyle\leq$ $\displaystyle 2(K+1)\exp\left(\frac{R_{H}-\eta\Delta}{K+1}+1\right)$ $\displaystyle=$ $\displaystyle 2(K+1)\exp\left(\frac{2\log(n)+K/e-\eta\Delta}{K+1}+1\right)$ $\displaystyle\leq$ $\displaystyle 8n^{\frac{2}{K+1}}(K+1)\exp\left(-\eta\frac{\Delta}{K+1}\right),$ $$

where the third inequality is because $R_{H}-\eta\Delta\leq 0$, and the last inequality holds because $\exp(\frac{K/e}{K+1}+1)\leq\exp(1+1/e)\leq 4$.
∎

###### Definition 1 .

###### Proof.

## Appendix B Proof of Theorem 2

For legibility, in what follows, we introduce the symbols $\hat{x},\hat{y},\hat{a}$ used exclusively for dummy variables.

###### Proof.

For the reader’s convenience, we list the explicit form of the Lyapunov function $f$ here:

$$ $f(\hat{x},\hat{y},\hat{a})=-\frac{1}{\eta}\sum_{ij}\exp{\left(\eta(-C_{ij}+ \sum_{m=1}^{L+K}\hat{a}_{m}\left(D_{m}\right)_{ij}+\hat{x}_{i}+\hat{y}_{j})-1 \right)}+\sum_{i}\hat{x}_{i}r_{i}+\sum_{j}\hat{y}_{j}c_{j}-\frac{1}{\eta}\sum_ {k=1}^{K}{\exp(-\eta\hat{a}_{k}-1)}.$ $$

First we prove $\sum_{ij}P_{ij}=1$. By assumption, the dual variable $(x,y,a)$ has gone through at least one update step. If the $x$ update step has been last performed, then one has $\nabla_{x}f=0$, which implies $\sum_{ij}P_{ij}=\sum_{i}r_{i}=1$. If the $y$ update step has been last performed, then one likewise has $\nabla_{y}f=0$ and $\sum_{ij}P_{ij}=\sum_{j}c_{j}=1$. If the $a$ update step has been performed, then the optimality in the $t$ variable as shown in equation 8 implies $\sum_{ij}P_{ij}=1$. Thus, one has $\sum_{ij}P_{ij}=1$ in all the three possible cases, as claimed.

The proof for equation 10 and equation 11 then largely follows from Lemma 2 of [Altschuler et al., ].
Suppose that an $x$ update step is performed and $P^{\prime}$ is the matrix formed by the dual variable $(x^{\prime},y,a)$. Due to the optimality of the $x^{\prime}$ variable, one has $\sum_{ij}P^{\prime}_{ij}=1$. Thus one has

$$ $f(x^{\prime},y,a)-f(x,y,a)=\frac{1}{\eta}(\sum_{ij}P_{ij}-\sum_{ij}P^{\prime}_ {ij})+\sum_{i}r_{i}(x^{\prime}_{i}-x_{i})=\sum_{i}r_{i}(x^{\prime}_{i}-x_{i})= \mathrm{KL}\left(r||P\mathbf{1}\right),$ $$

where the last equality is due to $r_{i}(x^{\prime}_{i}-x_{i})=r_{i}\log\frac{r_{i}}{\left(P\mathbf{1}\right)_{i}}$. The proof for equation 11 follows likewise:

$$ $f(x,y^{\prime},a)-f(x,y,a)=\frac{1}{\eta}(\sum_{ij}P_{ij}-\sum_{ij}P^{\prime}_ {ij})+\sum_{j}c_{j}(y^{\prime}_{j}-y_{j})=\sum_{j}c_{j}(y^{\prime}_{j}-y_{j})= \mathrm{KL}\left(c||P^{\top}\mathbf{1}\right).$ $$

We proceed with the proof for equation 13. First, we introduce an augmented Lyapunov function by maximizing over the $t$ variable, which provides a smoother optimization landscape for subsequent analysis. Define the augmented Lyapunov function $f_{\mathrm{aug}}$ as

$$ $f_{\mathrm{aug}}(\hat{x},\hat{y},\hat{a}):=\max_{\tilde{t}}f(\hat{x}+\tilde{t} \mathbf{1},\hat{y},\hat{a}).$ $$

As $\sum_{ij}P_{ij}=1$ implies optimality over the $t$ variable by equation 8, one has $f(x,y,a)=f_{\mathrm{aug}}(x,y,a)$ for the dual variable $(x,y,a)$. After the augmented $a$ update step, one has by definition that $f(x+t^{\prime}\mathbf{1},y,a^{\prime})=f_{\mathrm{aug}}(x,y,a^{\prime})$. Thus, one can simplify the left-hand side of equation 13 by the following equation:

$$ $f(x+t^{\prime}\mathbf{1},y,a^{\prime})-f(x,y,a)=f_{\mathrm{aug}}(x,y,a^{\prime })-f_{\mathrm{aug}}(x,y,a)=\max_{\tilde{a}}f_{\mathrm{aug}}(x,y,\tilde{a})-f_{ \mathrm{aug}}(x,y,a).$ $$

We then derive the formula for $f_{\mathrm{aug}}$. By direct calculation, one has

$$ $\displaystyle f(\hat{x}+t\mathbf{1},\hat{y},\hat{a})$ $\displaystyle=$ $\displaystyle-\frac{1}{\eta}\sum_{ij}\exp{\left(\eta(-C_{ij}+\sum_{m=1}^{L+K} \hat{a}_{m}\left(D_{m}\right)_{ij}+\hat{x}_{i}+\hat{y}_{j}+t)-1\right)}+\sum_{ i}\hat{x}_{i}r_{i}+\sum_{j}\hat{y}_{j}c_{j}+t-\frac{1}{\eta}\sum_{k=1}^{K}{ \exp(-\eta\hat{a}_{k}-1)}.$ $$

Let $t^{\star}=\operatorname*{arg\,max}_{\tilde{t}}f(\hat{x}+\tilde{t}\mathbf{1},
\hat{y},\hat{a})$.
Under the optimality condition $\partial_{t}f(\hat{x}+t^{\star}\mathbf{1},\hat{y},\hat{a})=0$, one has

$$ $\sum_{ij}\exp{\left(\eta(-C_{ij}+\sum_{m=1}^{L+K}\hat{a}_{m}\left(D_{m}\right) _{ij}+\hat{x}_{i}+\hat{y}_{j}+t^{\star})-1\right)}=1.$ $$

Taking log over the above equation, one has the following result for $t^{\star}$:

$$ $1-\eta t^{\star}=\mathrm{LSE}\left(\eta(-C+\sum_{m}\hat{a}_{m}D_{m}+\hat{x} \mathbf{1}^{\top}+\mathbf{1}\hat{y}^{\top})\right),$ (20) $$

where for a matrix $M$, the function $\mathrm{LSE}$ defines the log-sum-exponential function $\mathrm{LSE}\left(M\right)=\log\left(\sum_{ij}\exp{\left(m_{ij}\right)}\right)$. Thus one has

$$ $\displaystyle f_{\mathrm{aug}}(\hat{x},\hat{y},\hat{a})=$ $\displaystyle f(\hat{x}+t^{\star}\mathbf{1},\hat{y},\hat{a})$ $\displaystyle=$ $\displaystyle-\frac{1}{\eta}+\sum_{i}\hat{x}_{i}r_{i}+\sum_{j}\hat{y}_{j}c_{j} +t^{\star}-\frac{1}{\eta}\sum_{k=1}^{K}{\exp(-\eta\hat{a}_{k}-1)}$ $\displaystyle=$ $\displaystyle-\frac{1}{\eta}\mathrm{LSE}\left(\eta(-C+\sum_{m}\hat{a}_{m}D_{m} +\hat{x}\mathbf{1}^{\top}+\mathbf{1}\hat{y}^{\top})\right)+\sum_{i}\hat{x}_{i} r_{i}+\sum_{j}\hat{y}_{j}c_{j}-\frac{1}{\eta}\sum_{k=1}^{K}{\exp(-\eta\hat{a}_ {k}-1)}.$ $$

The rest of the proof follows from a perturbational argument on $f_{\mathrm{aug}}$ around the point $(x,y,a)$, and thus we provide a formula for the derivatives of $f_{\mathrm{aug}}$. Let $\hat{P}=\exp{\left(\eta(-C+\sum_{m}\hat{a}_{m}D_{m}+\hat{x}\bm{1}^{\top}+\bm{1
}\hat{y}^{\top})-1\right)}$ be the transport plan formed by the dual variable $(\hat{x},\hat{y},\hat{a})$. For $k=1,\ldots,K$, one has

$$ $\partial_{a_{k}}f_{\mathrm{aug}}(\hat{x},\hat{y},\hat{a})={\exp(-\eta\hat{a}_{ k}-1)}-D_{k}\cdot\frac{\hat{P}}{\sum_{ij}\hat{P}_{ij}},$ $$

and likewise for $l=1,\ldots,L$ one has

$$ $\partial_{a_{l+K}}f_{\mathrm{aug}}(\hat{x},\hat{y},\hat{a})=-D_{l+K}\cdot\frac {\hat{P}}{\sum_{ij}\hat{P}_{ij}}.$ $$

In terms of the second-order information, by direct calculation, one has

$$ $\partial_{a_{m}}\partial_{a_{m^{\prime}}}f_{\mathrm{aug}}(\hat{x},\hat{y},\hat {a})=-\eta(D_{m}\odot D_{m^{\prime}})\cdot\frac{\hat{P}}{\sum_{ij}\hat{P}_{ij} }+\eta\frac{\left(\hat{P}\cdot D_{m}\right)\left(\hat{P}\cdot D_{m^{\prime}} \right)}{\left(\sum_{ij}\hat{P}_{ij}\right)^{2}}+s_{mm^{\prime}},$ (21) $$

where $\odot$ is the Hadamard element-wise product, and $s_{mm^{\prime}}$ is zero except when $m=m^{\prime}=k$ for $k=1,\ldots,K$, in which case one has $s_{mm^{\prime}}=-\eta\exp(-\eta\hat{a}_{k}-1)$.

We then bound the spectrum of the Hessian matrix. Define

$$ $r_{mm^{\prime}}=-\eta(D_{m}\odot D_{m^{\prime}})\cdot\frac{\hat{P}}{\sum_{ij} \hat{P}_{ij}}+\eta\frac{\left(\hat{P}\cdot D_{m}\right)\left(\hat{P}\cdot D_{m ^{\prime}}\right)}{\left(\sum_{ij}\hat{P}_{ij}\right)^{2}},$ (22) $$

and $R=[r_{mm^{\prime}}]_{m,m^{\prime}=1}^{K+L}$. Let $(I,J)$ be a random vector supported on $[n]\times[n]$, and moreover let $(I,J)$ follow the multinomial distribution with $\mathbb{P}\left[(I,J)=(i,j)\right]=\frac{\hat{P}_{ij}}{\sum_{ij}\hat{P}_{ij}}$. Moreover, for $m=1,\ldots,K+L$, let $Y_{m}$ be the random variable defined by $Y_{m}=\left(D_{m}\right)_{IJ}$. Then, one can directly verify that $R=-\eta\mathrm{Cov}\left(Y_{1},\ldots,Y_{K+L}\right)$, and therefore $R$ is a negative semi-definite matrix. Let $S=[s_{mm^{\prime}}]_{m,m^{\prime}=1}^{K+L}$. The matrix $S$ is a diagonal matrix with non-positive entries and therefore is likewise negative semi-definite. By the trace bound one has

$$ $\displaystyle\lVert R\rVert_{2}\leq$ $\displaystyle-\sum_{m=1}^{K+L}r_{mm}$ $\displaystyle=$ $\displaystyle\sum_{m=1}^{K+L}\eta\left((D_{m}\odot D_{m})\cdot\frac{\hat{P}}{ \sum_{ij}\hat{P}_{ij}}-\left(\frac{\hat{P}\cdot D_{m}}{\sum_{ij}\hat{P}_{ij}} \right)^{2}\right)$ $\displaystyle\leq$ $\displaystyle\sum_{m=1}^{K+L}\eta\left((D_{m}\odot D_{m})\cdot\frac{\hat{P}}{ \sum_{ij}\hat{P}_{ij}}\right)$ $\displaystyle\leq$ $\displaystyle\eta(K+L)c_{d}^{2},$ $$

and therefore one has the spectral bound $S-\eta(K+L)c_{d}^{2}I\preceq\nabla_{a}^{2}f_{\mathrm{aug}}\preceq S$.

For simplicity, for the current variable $(x,y,a)$, define a function $g$ by $g(\hat{a})=f_{\mathrm{aug}}(x,y,\hat{a})$.
By Taylor’s remainder theorem, for any $\delta a$, one has

$$ $f_{\mathrm{aug}}(x,y,a+\delta a)-f_{\mathrm{aug}}(x,y,a)=g(a+\delta a)-g(a)= \delta a^{\top}\nabla_{a}g(a)+\frac{1}{2}\delta a^{\top}\nabla^{2}_{a}g(a+\xi \delta a)\delta a,$ $$

where $\xi\in[0,1]$ is an unknown quantity. Thus, the spectral bound on $\nabla_{a}^{2}f_{\mathrm{aug}}$ leads to the following inequality,

$$ $\displaystyle f_{\mathrm{aug}}(x,y,a+\delta a)-f_{\mathrm{aug}}(x,y,a)$ (23) $\displaystyle\geq$ $\displaystyle\delta a^{\top}\nabla_{a}g(a)-\frac{1}{2}\eta(K+L)c_{d}^{2}\lVert \delta a\rVert_{2}^{2}-\frac{1}{2}\sum_{k=1}^{K}\eta\exp(-\eta a_{k}-1)\exp(- \eta\xi\delta a_{k})\delta a_{k}^{2}$ $\displaystyle=$ $\displaystyle\sum_{m=1}^{K+L}\left(\delta a_{m}d_{m}-\frac{1}{2}\eta(K+L)c_{d} ^{2}\delta a_{m}^{2}\right)-\frac{1}{2}\sum_{k=1}^{K}\eta\exp(-\eta a_{k}-1) \exp(-\eta\xi\delta a_{k})\delta a_{k}^{2}$ $$

where $d_{m}=\partial_{a_{m}}g(a)$ is due to $\sum_{ij}P_{ij}=1$.

Importantly, the lower bound in equation 23 is fully separable in terms of the $\delta a_{m}$ terms. In what follows, we use equation 23 to give a construction for $\delta a$, which would then give a lower bound for the improvement in the Lyapunov function in the $a$ update step.
For $l=1,\ldots,L$, one sets $\delta a_{l+K}=\frac{d_{l+K}}{\eta(K+L)c_{d}^{2}}$. One then has

$$ $\delta a_{l+K}d_{l+K}-\frac{1}{2}\eta(K+L)c_{d}^{2}\delta a_{l+K}^{2}=\frac{d_ {l+K}^{2}}{2\eta(K+L)c_{d}^{2}}.$ $$

For $k=1,\ldots K$, one sets

$$ $\delta a_{k}=\max\left(-\frac{\log(2)}{\eta},\frac{d_{k}}{\eta\left(2\exp(- \eta a_{k}-1)+(K+L)c_{d}^{2}\right)}\right).$ (24) $$

We shall prove that the construction of $\delta a_{k}$ leads to the following bound:

$$ $\displaystyle\delta a_{k}d_{k}-\frac{1}{2}\eta(K+L)c_{d}^{2}\delta a_{k}^{2}- \frac{1}{2}\eta\exp(-\eta a_{k}-1)\exp(-\eta\xi\delta a_{k})\delta a_{k}^{2} \geq\lvert d_{k}\rvert\min\left(\frac{1}{8\eta},\frac{\lvert d_{k}\rvert}{8 \eta c_{d}+4\eta(K+L)c_{d}^{2}}\right).$ $$

By construction, one has $\exp(-\eta\xi\delta a_{k})\leq 2$, and thus

$$ $-\frac{1}{2}\eta\exp(-\eta a_{k}-1)\exp(-\eta\xi\delta a_{k})\delta a_{k}^{2} \geq-\frac{1}{2}\left(2\eta\exp(-\eta a_{k}-1)\delta\right)a_{k}^{2}.$ $$

If $-\frac{\log(2)}{\eta}\geq\frac{d_{k}}{\eta\left(2\exp(-\eta a_{k}-1)+(K+L)c_{d
}^{2}\right)}$, then one has $\delta a_{k}=-\frac{\log(2)}{\eta}\leq 0$ and $d_{k}\leq 0$. By equation 24 one has

$$ $\delta a_{k}\eta\left(2\exp(-\eta a_{k}-1)+(K+L)c_{d}^{2}\right)\geq d_{k},$ $$

and thus multiplying both sides by $-\frac{1}{2}\delta a_{k}$, the inequality becomes

$$ $-\frac{1}{2}\eta\left(2\exp(-\eta a_{k}-1)+(K+L)c_{d}^{2}\right)\delta a_{k}^{ 2}\geq-\frac{1}{2}\delta a_{k}d_{k}.$ $$

Thus

$$ $\displaystyle\delta a_{k}d_{k}-\frac{1}{2}\eta(K+L)c_{d}^{2}\delta a_{k}^{2}- \frac{1}{2}\eta\exp(-\eta a_{k}-1)\exp(-\eta\xi\delta a_{k})\delta a_{k}^{2}$ $\displaystyle\geq$ $\displaystyle\delta a_{k}d_{k}-\frac{1}{2}\eta\left(2\exp(-\eta a_{k}-1)+(K+L) c_{d}^{2}\right)\delta a_{k}^{2}$ $\displaystyle\geq$ $\displaystyle\delta a_{k}d_{k}-\frac{1}{2}\delta a_{k}d_{k}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\delta a_{k}d_{k}$ $\displaystyle=$ $\displaystyle\lvert d_{k}\rvert\frac{\log(2)}{2\eta},$ $$

which in particular implies the claimed bound as $\frac{\log(2)}{2}\geq\frac{1}{8}$.

Otherwise, if $-\frac{\log(2)}{\eta}\leq\frac{d_{k}}{\eta\left(2\exp(-\eta a_{k}-1)+(K+L)c_{d
}^{2}\right)}$, then one has $\delta a_{k}=\frac{d_{k}}{\eta\left(2\exp(-\eta a_{k}-1)+(K+L)c_{d}^{2}\right)}$ by construction. Then, one has

$$ $\displaystyle\delta a_{k}d_{k}-\frac{1}{2}\eta(K+L)c_{d}^{2}\delta a_{k}^{2}- \frac{1}{2}\eta\exp(-\eta a_{k}-1)\exp(-\eta\xi\delta a_{k})\delta a_{k}^{2}$ $\displaystyle\geq$ $\displaystyle\delta a_{k}d_{k}-\frac{1}{2}\eta\left(2\exp(-\eta a_{k}-1)+(K+L) c_{d}^{2}\right)\delta a_{k}^{2}$ $\displaystyle=$ $\displaystyle\frac{d_{k}^{2}}{2\eta\left(2\exp(-\eta a_{k}-1)+(K+L)c_{d}^{2} \right)}$ $\displaystyle=$ $\displaystyle\frac{d_{k}^{2}}{4\eta\exp(-\eta a_{k}-1)+2\eta(K+L)c_{d}^{2}}$ $\displaystyle\geq$ $\displaystyle\lvert d_{k}\rvert\frac{\lvert d_{k}\rvert}{4\eta\lvert d_{k} \rvert+4\eta c_{d}+2\eta(K+L)c_{d}^{2}}$ $\displaystyle\geq$ $\displaystyle\lvert d_{k}\rvert\min\left(\frac{1}{8\eta},\frac{\lvert d_{k} \rvert}{8\eta c_{d}+4\eta(K+L)c_{d}^{2}}\right),$ $$

where the third inequality is obtained by applying the mediant inequality, and the second inequality is because

$$ $\exp(-\eta a_{k}-1)\leq\lvert\exp(-\eta a_{k}-1)-P\cdot D_{k}\rvert+\lvert P \cdot D_{k}\rvert=\lvert d_{k}\rvert+\lvert P\cdot D_{k}\rvert\leq\lvert d_{k} \rvert+c_{d}.$ $$

The proof for equation 13 is by organizing the arranged results:

$$ $\displaystyle f(x+t^{\prime}\mathbf{1},y,a^{\prime})-f(x,y,a)$ $\displaystyle=$ $\displaystyle f_{\mathrm{aug}}(x,y,a^{\prime})-f_{\mathrm{aug}}(x,y,a)$ $\displaystyle=$ $\displaystyle\max_{\tilde{a}}f_{\mathrm{aug}}(x,y,\tilde{a})-f_{\mathrm{aug}}( x,y,a)$ $\displaystyle\geq$ $\displaystyle f_{\mathrm{aug}}(x,y,a+\delta a)-f_{\mathrm{aug}}(x,y,a)$ $\displaystyle=$ $\displaystyle g(a+\delta a)-g(a)$ $\displaystyle=$ $\displaystyle\delta a^{\top}\nabla_{a}g(a)+\frac{1}{2}\delta a^{\top}\nabla^{2 }_{a}g(a+\xi\delta a)\delta a$ $\displaystyle\geq$ $\displaystyle\sum_{k=1}^{K}\left(\delta a_{k}d_{k}-\frac{1}{2}\eta(K+L)c_{d}^{ 2}\delta a_{k}^{2}-\frac{1}{2}\eta\exp(-\eta a_{k}-1)\exp(-\eta\xi\delta a_{k} )\delta a_{k}^{2}\right)$ $\displaystyle\qquad+\sum_{l=1}^{L}\left(\delta a_{l+K}d_{l+K}-\frac{1}{2}\eta( K+L)c_{d}^{2}\delta a_{l+K}^{2}\right)$ $\displaystyle\geq$ $\displaystyle\sum_{k=1}^{K}\lvert d_{k}\rvert\min\left(\frac{1}{8\eta},\frac{ \lvert d_{k}\rvert}{8\eta c_{d}+4\eta(K+L)c_{d}^{2}}\right)+\sum_{l=1}^{L} \frac{d_{l+K}^{2}}{2\eta(K+L)c_{d}^{2}}.$ $$

∎

###### Remark 1 .

In the proof for Theorem 2, it might be advantageous to use an alternative spectral bound when the matrices $D_{m}$ have special sparsity structures. In such special cases, it is better to use the Gershgorin circle theorem instead of the trace bound in the numerical treatment for the $a$ step.

One defines

$$ $\displaystyle r_{1;mm^{\prime}}=-\eta(D_{m}\odot D_{m^{\prime}})\cdot\frac{ \hat{P}}{\sum_{ij}\hat{P}_{ij}},$ (25) $\displaystyle r_{2;mm^{\prime}}=\eta\frac{\left(\hat{P}\cdot D_{m}\right)\left (\hat{P}\cdot D_{m^{\prime}}\right)}{\left(\sum_{ij}\hat{P}_{ij}\right)^{2}},$ $$

and let $R_{1}=[r_{1;mm^{\prime}}]_{m,m^{\prime}=1}^{K+L},R_{2}=[r_{2;mm^{\prime}}]_{m,
m^{\prime}=1}^{K+L}$. Let $Y_{m}$ for $m=1,\ldots,K+L$ be as in the proof of Theorem 2. One can directly verify that $R_{2}$ is positive semi-definite as it is the outer product of the vector $\sqrt{\eta}\mathbb{E}\left[Y_{1},\ldots,Y_{K+L}\right]$ with itself, and thus $S+R_{1}\preceq\nabla_{a}^{2}f_{\mathrm{aug}}$. By the Gershgorin circle theorem, it follows that $\lVert R_{1}\rVert_{2}$ is bounded by the matrix 1-norm of $R_{1}$. In particular, define $E\in\mathbb{R}^{n\times n}$ as the sum of the constraint matrices up to absolute value, i.e.

$$ $\left(E\right)_{ij}=\sum_{m=1}^{K+L}\lvert\left(D_{m}\right)_{ij}\rvert.$ $$

Then, one has

$$ $\displaystyle\lVert R_{1}\rVert_{2}\leq$ $\displaystyle\max_{m}\sum_{m^{\prime}=1}^{K+L}\lvert r_{1;mm^{\prime}}\rvert$ $\displaystyle=$ $\displaystyle\eta\max_{m}\sum_{m^{\prime}=1}^{K+L}\left|(D_{m}\odot D_{m^{ \prime}})\cdot\frac{\hat{P}}{\sum_{ij}\hat{P}_{ij}}\right|$ $\displaystyle=$ $\displaystyle\eta\max_{m}\sum_{m^{\prime}=1}^{K+L}\left|\left(\frac{\hat{P}}{ \sum_{ij}\hat{P}_{ij}}\odot D_{m^{\prime}}\right)\cdot D_{m}\right|$ $\displaystyle\leq$ $\displaystyle\eta c_{d}\sum_{m^{\prime}=1}^{K+L}\lVert\frac{\hat{P}}{\sum_{ij} \hat{P}_{ij}}\odot D_{m^{\prime}}\rVert_{1}$ $\displaystyle=$ $\displaystyle\eta c_{d}\sum_{m^{\prime}=1}^{K+L}\left|\frac{\hat{P}}{\sum_{ij} \hat{P}_{ij}}\cdot D_{m^{\prime}}\right|$ $\displaystyle\leq$ $\displaystyle\eta c_{d}\left(\frac{\hat{P}}{\sum_{ij}\hat{P}_{ij}}\cdot E\right)$ $\displaystyle\leq$ $\displaystyle\eta c_{d}\lVert E\rVert_{\infty},$ $$

where the second and third equality is by the property of the Hadamard product, the second and fourth inequality is by Holder’s inequality, and the third inequality is by the definition of $E$. We define $c_{e}=\lVert E\rVert_{\infty}$, which gives one the following bound in contrast to equation 23:

$$ $\displaystyle f_{\mathrm{aug}}(x,y,a+\delta a)-f_{\mathrm{aug}}(x,y,a)$ (26) $\displaystyle\geq$ $\displaystyle\delta a^{\top}\nabla_{a}g(a)-\frac{1}{2}\eta c_{e}c_{d}\lVert \delta a\rVert_{2}^{2}-\frac{1}{2}\sum_{k=1}^{K}\eta\exp(-\eta a_{k}-1)\exp(- \eta\xi\delta a_{k})\delta a_{k}^{2}$ $\displaystyle=$ $\displaystyle\sum_{m=1}^{K+L}\left(\delta a_{m}d_{m}-\frac{1}{2}\eta c_{e}c_{d }\delta a_{m}^{2}\right)-\frac{1}{2}\sum_{k=1}^{K}\eta\exp(-\eta a_{k}-1)\exp( -\eta\xi\delta a_{k})\delta a_{k}^{2}.$ $$

Using equation 26, one can derive a tighter bound in the conclusion part of Theorem 2 by applying a similar argument in the main statement. Overall, the bound in equation 26 is more advantageous when $c_{e}$ is significantly smaller than $(K+L)c_{d}$, and indeed one can see that $c_{e}\leq(K+L)c_{d}$. One situation in which an advantage exists is when the constraints $D_{m}$ are themselves sparse, which is the case in capacity constrained OT.

###### Proof.

###### Remark 1 .

## Appendix C Proof of Theorem 3

###### Proof.

For any iteration index $u$, let $Q^{u}_{x},Q^{u}_{y},Q^{u}_{a}$ be the quantity as defined by the Theorem for $P_{u}$ in place of $P$. Similarly we define $d_{m}^{u}$ by equation 12 with $P_{u}$ in place of $P$. One then has by construction that

$$ $f(x^{u+1},y^{u+1},a^{u+1})-f(x^{u},y^{u},a^{u})\geq\frac{1}{3}\left(Q_{x}^{u}+ Q_{y}^{u}+Q_{a}^{u}\right).$ (27) $$

In particular, for any $k=1,\ldots,K$, one has

$$ $f(x^{u+1},y^{u+1},a^{u+1})-f(x^{u},y^{u},a^{u})\geq\frac{1}{3}\lvert d^{u}_{k} \rvert\min\left(\frac{1}{8\eta},\frac{\lvert d^{u}_{k}\rvert}{8\eta c_{d}+4 \eta(K+L)c_{d}^{2}}\right).$ $$

Suppose that there are $N_{I}$ iterations for which $\frac{1}{8\eta}\leq\frac{\lvert d^{u}_{k}\rvert}{8\eta c_{d}+4\eta(K+L)c_{d}^{
2}}$ for some $k=1,\ldots,K$. For such $u$, one has $\lvert d^{u}_{k}\rvert=\Omega(1)$ and therefore one has $f(x^{u+1},y^{u+1},a^{u+1})-f(x^{u},y^{u},a^{u})=\Omega(1)$. Thus $N_{I}=O(c_{g})$, and therefore for $k=1,\ldots,K$, one has $\frac{1}{8\eta}\geq\frac{\lvert d^{u}_{k}\rvert}{8\eta c_{d}+4\eta(K+L)c_{d}^{
2}}$ except for $O(c_{g})$ iterations. Let $C=\frac{1}{8\eta c_{d}+4\eta(K+L)c_{d}^{2}}$, and then the following condition holds except for $N_{I}$ iterations:

$$ $f(x^{u+1},y^{u+1},a^{u+1})-f(x^{u},y^{u},a^{u})\geq\frac{1}{3}\left(\mathrm{KL }\left(P_{u}\mathbf{1}||r\right)+\mathrm{KL}\left(P_{u}^{\top}\mathbf{1}||c \right)+C\sum_{m=1}^{K+L}(d_{m}^{u})^{2}\right).$ (28) $$

Let $N_{II}$ be the number of iterations for which the following condition is satisfied

$$ $\lVert\nabla f\rVert_{1}=\lVert P_{u}\mathbf{1}-r\rVert_{1}+\lVert P_{u}^{\top }\mathbf{1}-c\rVert_{1}+\sum_{m=1}^{K+L}\lvert d_{m}\rvert>\epsilon.$ (29) $$

Then, if $u$ is among the said $N_{II}$ iterations, one has

$$ $\displaystyle\epsilon^{2}<$ $\displaystyle\left(\lVert P_{u}\mathbf{1}-r\rVert_{1}+\lVert P_{u}^{\top} \mathbf{1}-c\rVert_{1}+\sum_{m=1}^{K+L}\lvert d^{u}_{m}\rvert\right)^{2}$ $\displaystyle\leq$ $\displaystyle(K+L+2)^{2}\left(\lVert P_{u}\mathbf{1}-r\rVert_{1}^{2}+\lVert P_ {u}^{\top}\mathbf{1}-c\rVert_{1}^{2}+\sum_{m=1}^{K+L}(d_{m}^{u})^{2}\right)$ $\displaystyle\leq$ $\displaystyle(K+L+2)^{2}\left(2\mathrm{KL}\left(P_{u}\mathbf{1}||r\right)+2 \mathrm{KL}\left(P_{u}^{\top}\mathbf{1}||c\right)+\sum_{m=1}^{K+L}(d_{m}^{u})^ {2}\right),$ $$

where the second inequality is by Cauchy-Schwartz, the third inequality is equation 12, and the third inequality is by Pinsker’s inequality. Thus, for at least $N_{II}-N_{I}$ iterations, both equation 28 and equation 29 are satisfied, under which there exists a constant $C^{\prime}$ for which one has

$$ $\displaystyle\epsilon^{2}\leq$ $\displaystyle(K+L+2)^{2}\left(2\mathrm{KL}\left(P_{u}\mathbf{1}||r\right)+2 \mathrm{KL}\left(P_{u}^{\top}\mathbf{1}||c\right)+\sum_{m=1}^{K+L}(d_{m}^{u})^ {2}\right)$ $\displaystyle\leq$ $\displaystyle C^{\prime}\left(f(x^{u+1},y^{u+1},a^{u+1})-f(x^{u},y^{u},a^{u})\right)$ $$

As $\sum_{u=0}^{\infty}f(x^{u+1},y^{u+1},a^{u+1})-f(x^{u},y^{u},a^{u})\leq c_{g}$, one must have $N_{II}-N_{I}=O\left((\epsilon)^{-2}c_{g}\right)$. One then has $N_{II}=N_{I}+O\left((\epsilon)^{-2}c_{g}\right)=O\left((\epsilon)^{-2}c_{g}\right)$, as desired, and thus the first $u$ for which equation 14 holds occurs after $O\left((\epsilon)^{-2}c_{g}\right)$ iterations.

Lastly, equation 15 holds as a result of equation 14 because one has $\lvert d^{u}_{K+l}\rvert=P^{u}\cdot D_{K+l}$ and $\lvert d^{u}_{k}\rvert\geq\lvert\min(0,P^{u}\cdot D_{k})\rvert$.

∎

As a result of Theorem 3, one has the following bound if one applies a rounding operation to the output of $P$ in Theorem 3:

###### Proposition 1 .

Let $P$ be the same as in Theorem 3, let $\mathcal{U}_{r,c}$ denote the set of transport matrix from $r$ to $c$, and let $\mathrm{Round}(P,\mathcal{U}_{r,c})$ be the result of applying the rounding algorithm in Altschuler et al. [] to obtain projection of $P$ into $\mathcal{U}_{r,c}$. Define constraint violation of $P$ by

$$ $\displaystyle\mathrm{Violation}(P)=$ $\displaystyle\sum_{k=1}^{K}\lvert\min\left(\mathrm{Round}(P,\mathcal{U}_{r,c}) \cdot D_{k},0\right)\rvert$ $\displaystyle+$ $\displaystyle\sum_{l=1}^{L}\lvert\mathrm{Round}(P,\mathcal{U}_{r,c})\cdot D_{l +K}\rvert.$ $$

One has

$$ $\mathrm{Violation}(P)\leq\epsilon\left(1+2(K+L)c_{d}\right)$ $$

where $c_{d}:=\max_{m\in[K+L]}\lVert D_{m}\rVert_{\infty}$.

###### Proof.

This proof is a simple consequence of existing results. By Lemma 7 in Altschuler et al. [], one has

$$ $\lVert P-\mathrm{Round}(P,\mathcal{U}_{r,c})\rVert_{1}\leq 2\epsilon.$ $$

Thus one has

$$ $\displaystyle\mathrm{Violation}(P)\leq$ $\displaystyle\sum_{k=1}^{K}\lvert\min\left(P\cdot D_{k},0\right)\rvert+\sum_{l =1}^{L}\lvert P\cdot D_{l+K}\rvert+\sum_{m=1}^{K+L}\lvert\left(P-\mathrm{Round }(P,\mathcal{U}_{r,c})\right)\cdot D_{m}\rvert$ $\displaystyle\leq$ $\displaystyle\epsilon+2(K+L)c_{d}\epsilon,$ $$

where the last inequality uses Holder’s inequality and equation 15.
∎

###### Proof.

###### Proposition 1 .

###### Proof.

## Appendix D Equivalence of primal and primal-dual form

We now show that the primal form in equation 5 can be obtained from the primal-dual form by eliminating the dual variables.

###### Proposition 2 .

Define

$$ $\displaystyle L(P,s,x,y,a)=$ $\displaystyle\frac{1}{\eta}P\cdot\log P+C\cdot P-x\cdot(P\bm{1}-r)-y\cdot(P^{ \top}\bm{1}-c)$ $\displaystyle+$ $\displaystyle\frac{1}{\eta}\sum_{k=1}^{K}s_{k}\log{s_{k}}+\sum_{k=1}^{K}a_{k}s _{k}-\sum_{m=1}^{K+L}a_{m}(D_{m}\cdot P),$ $$

and then, for $\mathcal{S}$ as in equation 4, the following equation holds:

$$ $\max_{x,y,a}\min_{P,s}L(P,s,x,y,a,b)=\min_{P:P\bm{1}=r,P^{\top}\bm{1}=c,P\in \mathcal{S}}\frac{1}{\eta}P\cdot\log P+\sum_{k=1}^{K}\frac{1}{\eta}(D_{k}\cdot P )\cdot\log(D_{k}\cdot P)+C\cdot P$ (30) $$

Moreover, for the Lyapunov potential function $f$ in equation 6, one has

$$ $f(x,y,a)=\min_{P,s}L(P,s,x,y,a).$ (31) $$

###### Proof.

As $L$ is concave in $x,y$ and convex in $P$, one can invoke the minimax theorem to interchange the operations of maximization and minimization. Therefore:

$$ $\displaystyle\min_{P,s}\max_{x,y,a}L(P,s,x,y,a)=$ $\displaystyle\min_{P,s:P\bm{1}=r,P^{\top}\bm{1}=c,P\in\mathcal{S},s_{k}=D_{k} \cdot P\,\forall k\in[K]}\frac{1}{\eta}P\cdot\log P+\frac{1}{\eta}s\cdot\log s +C\cdot P$ $\displaystyle=$ $\displaystyle\min_{P:P\bm{1}=r,P^{\top}\bm{1}=c,P\in\mathcal{S}}\frac{1}{\eta} P\cdot\log P+\sum_{k=1}^{K}\frac{1}{\eta}(D_{k}\cdot P)\cdot\log(D_{k}\cdot P) +C\cdot P.$ $$

In terms of entries, one writes $L(P,s,x,y,a)$ as follows:

$$ $\displaystyle\max_{x_{i},y_{j}}\min_{p_{ij}}L(P,s,x,y,a)=$ $\displaystyle\frac{1}{\eta}\sum_{ij}p_{ij}\log p_{ij}+\sum_{ij}C_{ij}p_{ij}- \sum_{i}x_{i}(\sum_{j}p_{ij}-r_{i})-\sum_{j}y_{j}(\sum_{i}p_{ij}-c_{j})$ (32) $\displaystyle+$ $\displaystyle\frac{1}{\eta}\sum_{k=1}^{K}s_{k}\log{s_{k}}+\sum_{k=1}^{K}a_{k}s _{k}-\sum_{m=1}^{K+L}\sum_{ij}a_{m}(D_{m})_{ij}p_{ij}$ $$

We then solve the inner min problem explicitly by taking the derivative of $p_{ij},s_{k}$ to zero, from which one obtains

$$ $p_{ij}=\exp(\eta(-C_{ij}+\sum_{m=1}^{K+L}a_{m}(D_{m})_{ij}+x_{i}+y_{j})-1).$ $$

and

$$ $s_{k}=\exp(-\eta a_{k}-1).$ $$

Plugging in the formula for $p_{ij}$ and $s_{k}$, one has

$$ $\displaystyle\min_{P,s}L(P,s,x,y,a)$ $\displaystyle=-\frac{1}{\eta}\sum_{ij}\exp(\eta(-C_{ij}+\sum_{m=1}^{K+L}a_{m}( D_{m})_{ij}+x_{i}+y_{j})-1)+\sum_{i}r_{i}x_{i}+\sum_{j}c_{j}y_{j}-\frac{1}{ \eta}\sum_{k=1}^{K}{\exp(-\eta a_{k}-1)},$ $$

which is equal to $f(x,y,a)$.
∎

###### Proposition 2 .

###### Proof.

## Appendix E Relation of Algorithm 1 to Bregman projection

In this section, we show that the proposed $a$ update step is equivalent to a Bregman projection for equality constraints. The Bregman projection step introduced in Benamou et al. [] introduces an iterative projection-based approach. As in the setting of Section 3, let $(x,y,a)$ be the current dual variable and let $P$ be the intermediate matrix corresponding to $(x,y,a)$.

Suppose that $K=0$ and thus there are no inequality constraints. As in main text, we define $\mathcal{E}$ by the following space

$$ $\mathcal{E}:=\bigcap_{l=1,\ldots,L}\{M\mid D_{l}\cdot M=0\}.$ $$

Moreover, assume $\sum_{i}r_{i}=\sum_{j}c_{j}=1$ and define $\Delta_{n\times n}$ as the $n^{2}$-dimensional simplex.
Let $a^{\prime},t^{\prime}=\operatorname*{arg\,max}_{\tilde{a},\tilde{t}}f(x+\tilde
{t}\mathbf{1},y,\tilde{a})$, and let $P^{\prime}$ be the intermediate matrix corresponding to $(x+t^{\prime}\mathbf{1},y,a^{\prime})$. Then, the claimed equivalence can be seen by proving the following equation:

$$ $P^{\prime}=\operatorname*{arg\,min}_{M\in\Delta_{n\times n}\cap\mathcal{E}} \mathrm{KL}\left(M||P\right)$ (33) $$

where for two entry-wise non-negative matrices $M,N$, the term $\mathrm{KL}\left(M||N\right)$ is defined by $\mathrm{KL}\left(M||N\right)=\sum_{ij}m_{ij}\left(\log(\frac{m_{ij}}{n_{ij}})-
1\right)$. Suppose $M=[m_{ij}]_{i,j=1}^{n}$ achieves the optimality condition set in equation 33. As it is a constrained optimization problem, the necessary condition for optimality is that there exists $\mu$ and $\lambda_{l}$ for $l=1,\ldots,L$, so that the following holds:

$$ $\forall{i,j\in[n]},\quad\partial_{m_{ij}}\mathrm{KL}\left(M||P\right)=\eta \left(\sum_{l=1}^{L}\lambda_{l}(D_{l})_{ij}+\mu\right).$ $$

Utilizing the definition of $\mathrm{KL}\left(\cdot||\cdot\right)$, one rewrites the above equation as below:

$$ $\log(m_{ij})-\log(p_{ij})=\eta\left(\sum_{l=1}^{L}\lambda_{l}(D_{l})_{ij}+\mu \right).$ $$

Thus, there exists $\mu$ and $\lambda_{l}$ for $l=1,\ldots,L$, for which one has $M=\exp(\sum_{l}\lambda_{l}D_{m}+\mu\mathbf{1}\mathbf{1}^{\top})\odot P$, where $\odot$ is the Hadamard product. Thus, one has

$$ $M_{ij}=\exp{\left(\eta(-C_{ij}+\sum_{l=1}^{L}(a_{m}+\lambda_{l})\left(D_{l} \right)_{ij}+x_{i}+y_{j}+\mu)-1\right)}.$ (34) $$

Furthermore, $M$ satisfies the following normalization condition

$$ $\sum_{ij}M_{ij}=1.$ (35) $$

Moreover, for $l\in[L]$, the following holds:

$$ $M\cdot D_{l}=0.$ (36) $$

As seen in the main text, $P^{\prime}$ satisfies equation 35, equation 36. Moreover, by setting $\lambda_{l}=a^{\prime}_{m}-a_{m}$ and $\mu=t^{\prime}$, one can show that $P^{\prime}$ also satisfies equation 34. As $f$ is concave and the equality constraints are affine, this shows that $P^{\prime}$ satisfies the optimality condition in that of the right-hand side of equation 33.

## Appendix F Practical implementation of accelerated Sinkhorn-type algorithm under constraint

In this section, we detail the procedure to combine Algorithm 1 with entropy regularization scheduling and Sinkhorn-Newton-Sparse (SNS). As the detail of entropic regularization scheduling is presented in Section 3, we shall give implementation details of SNS in the constrained case.
Similar to the construction in Tang et al. [], for a matrix $M\in\mathbb{R}^{n\times n}_{\geq 0}$, we use $\mathrm{Sparisfy}(M,\rho)$ to denote entry-wise truncation with a threshold $\rho$. For $\tilde{M}:=\mathrm{Sparisfy}(M,\rho)$, one has

$$ $\tilde{M}_{ij}=\begin{cases}M_{ij}&\text{if }M_{ij}\geq\rho,\\ 0&\text{otherwise}.\end{cases}$ $$

Let $P$ be the intermediate transport matrix formed by the current dual variable $(x,y,a)$. One can write down the Hessian term $\nabla^{2}f$ as follows:

$$ $\nabla^{2}f(x,y,a)=-\eta\begin{bmatrix}\mathrm{diag}(P\bm{1})&P&\nabla_{x} \nabla_{a}f\\ P^{\top}&\mathrm{diag}(P^{\top}\bm{1})&\nabla_{y}\nabla_{a}f\\ \nabla_{a}\nabla_{x}f&\nabla_{a}\nabla_{y}f&\nabla_{a}^{2}f\end{bmatrix}.$ (37) $$

As we assume $K+L=O(1)$ in this work, the terms such as $\nabla_{xy}\nabla_{a}f,\nabla_{a}^{2}f$ can be kept without posing significant challenges to the Newton step, which allows us to use the approximation below:

$$ $\nabla^{2}f(x,y,a)\approx H=-\eta\begin{bmatrix}\mathrm{diag}(P\bm{1})&\mathrm {Sparisfy}(P,\rho)&\nabla_{x}\nabla_{a}f\\ \mathrm{Sparisfy}(P^{\top},\rho)&\mathrm{diag}(P^{\top}\bm{1})&\nabla_{y} \nabla_{a}f\\ \nabla_{a}\nabla_{x}f&\nabla_{a}\nabla_{y}f&\nabla_{a}^{2}f\end{bmatrix}.$ (38) $$

The value of $\rho$ is a tunable parameter, and one sets $\rho$ so that $\mathrm{Sparisfy}(P,\rho)$ contains only $O(n)$ nonzero entries. As a result, applying $H^{-1}$ to a vector can be done with a $O((n+K+L)^{2})=O(n^{2})$ complexity through the conjugate-gradient algorithm [Golub and Van Loan, ]. Additionally, to ensure stability in the degenerate direction $(\delta x,\delta y,\delta a)=v:=(\mathbf{1}_{n},-\mathbf{1}_{n},\mathbf{0}_{K+
L})$, we in practice use the following modified version of Lyapunov function:

$$ $\tilde{f}(x,y,a)=f(x,y,a)-\frac{1}{2}(\sum_{i}x_{i}-\sum_{j}y_{j})^{2}.$ $$

Same as Tang et al. [], one can see that the maximizer of $\tilde{f}$ is also a maximizer of $f$. The final Sinkhorn-Newton-Sparse algorithm used is in Algorithm 2, where we include a Sinkhorn stage which uses Algorithm 1 as initialization, and a subsequent Newton stage which uses sparse Newton iteration to accelerate convergence.

Figure: Algorithm 2 Sinkhorn-Newton-Sparse (SNS) for OT under linear constraint

The combination of entropy regularization scheduling and SNS is described in Algorithm 3. We remark that the iteration count $N_{1,i},N_{2,i}$ within Algorithm 3 are typically set to be much smaller than $N_{1},N_{2}$ in Algorithm 2, which is possible because the optimization task at an entropy regularization parameter $\eta_{i}$ is initialized by the dual variables obtained for regularization parameter $\eta_{i-1}=\eta_{i}/2$.

As a special case, Algorithm 3 also provides substantial acceleration to the unconstrained entropic OT problem [Cuturi, ].

Figure: Algorithm 3 Sinkhorn-Newton-Sparse with entropy regularization scheduling for OT under linear constraint

## Appendix G Further details on numerical experiments

Despite the theoretical analysis in this work, it is an open problem as to whether the intermediate matrix $P$ formed by the dual variable $(x,y,a)$ converges to $P_{\eta}^{\star}$, which is the entropic optimal solution to equation 5. While the theoretical result is left for future work, we show that the numerical performance of Algorithm 1 lends strong support for rapid practical convergence. In Figure 6 and Figure 6, we show that both Algorithm 1 and Algorithm 2 can converge to $P_{\eta}^{\star}$ in the total variation (TV) distance. Furthermore, Algorithm 2 provides machine-epsilon accurate solutions rapidly during the Newton stage.

We further remark that the ranking example seems to suggest an exponential convergence result exists in practice for Algorithm 1 in this setting. After conducting repeated experiments in the ranking setting under different problem parameters, the numerical result obtained shows that exponential convergence is not peculiar to the experiment but is instead a persistent feature in the ranking case. The analysis in this special case might be of independent research interest for the field of information retrieval.

Figure: Figure 5: Random assignment problem. Plot of the proposed Sinkhorn-type algorithm in terms of TV distance to $P_{\eta}^{\star}$.
Refer to caption: x5.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 4
questions: Questions:
- For the Random assignment problem, how would the algorithm behave if cost matrix is not generated through uniform sampling and instead has some other distribution? Do you have any observations or thoughts?
- Sinkhorn supports backpropagation in neural networks owing to matrix/vector operations and hence can be used as part of the loss function or during intermediate steps in computation. Can you comment on whether your approach can be used directly in a similar fashion? Please explain your reaasoning.

Minor suggestions that won't affect the score:
- Line 166/167: “the following” repeated twice in the statement "We summarize the general form of constrained optimal transport by the following the following linear program (LP)"
- Consider citing combinatorial OT solvers that have been proposed that have successfully shown applications in ML to compute partial optimal transports
soundness: 3
strengths: - This paper combines ideas from several previous works and extends them meaningfully in novel ways to solve constrained OT problems under both equality and inequality constraints.
- The novelty lies in the use of Lyapunov function to characterize the optimization procedure to perform the constraint update dual step
- Authors have presented convergence analysis of the proposed Sinkhorn-type optimization procedure with acceleration mechanisms
- Authors have provided decent survey of related literature relating to OT, ML and efficient solvers for constrained and unconstrained OT.
- Authors have provided detailed proofs for theorems and propositions wherever necessary.
summary: This paper deals with constrained OT problems and the corresponding entropic regularization formulation which could potentially help researchers and practitioners of OT and Machine Learning to arrive at transport plans that have more complex structure than in the unconstrained case. Theoretically, the work is well grounded in existing literature and the authors have made use of multiple recent advancements to design a Sinkhorn-type algorithm to solve constrained OT problems and also have proposed accelerated convergence methods with corresponding bounds. The authors also provide proofs for the various theorems and propositions under certain assumptions. For methods that rely on kernel approximation, prior work do not account for dynamically evolving kernels in the constrained case such as in case of $K=exp(-C\eta)$ but with the authors approach of using a Lyapunov function to characterize the optimization in the variational formulation allows for solving transport problems with evolving kernels.
weaknesses: - The numerical experiments are based on (weak) assumptions such as the cost matrix entries being sampled from uniform distribution in case of random assignment problem or the Rademacher distribution in case of Ranking under constraints (appendix A). It is not clear how the proposed algorithm performs when the cost matrix may not conform to simple distributions.
-  Authors could present experiments that would be more relevant to the target ML community.
- Solving large scale problems would help ascertain the claims made by authors about the usability and efficiency of their proposed approach.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: I generally find the paper interesting, but have few concerns:

1- My main concern is that I am not sure if the choice of an entropy regularization makes sense for the additional constraints. For the standard OT problem, this choice is justified as the dual problem can be solved by exact block coordinate descent, but if one needs to employ the Newton's method, then, why should not one use e.g. a logarithmic barrier (instead of entropic regularization), which has the self-concordance property and is potentially better for the doubling strategy.

2- Most of the analysis seems like a straightforward extensions of previous works Weed (2018) and Altschuler et al. (2017). A general concern about Theorem 1 is that although it is formulated as an exponential decay, it really shows the requirement that the regularization parameter grows proportionally with the inverse of the duality gap. In practice, the gap can be extremely small for large problems, leading to extremely large regularization parameters. On the other hand, the result of Theorem 2 seems to be dimension-free, but it actually depends on the dimension through the regularization parameter. As such, I think that Theorem 2 should clearly reflect the dependency on eta (which is hopefully linear as in Altschuler et al. (2017)).

3- Although the authors mention some applications of their problem in their literature review, I think that the paper generally does not well motivate the study. The experiments are on toy scenarios and do not reflect scalability as they consider relatively small problems. For the MNIST case, for example, a more realistic domain adaptation scenario could be considered with at least few thousand points in each domain.
soundness: 4
strengths: Constrained optimal transport may have various applications. IN ML it can be used e.g. for various domain adaptation scenarios with structured data. However, OT is a difficult problem to solve at scale. The Sinkhorn-type algorithms are widely believed to be a suitable way to address the difficulties with the exact OT problem. The paper pursues such a solution and the discussion is supported by rigorous mathematical results.  The numerical results also reflect the superior convergence properties of the proposed algorithm.
summary: This paper addresses scalable solutions to the problem of optimal transport by Kantorovich relaxation, under arbitrary additional linear constraints. By introducing a set of slack variables and adopting an entropic regularization, it arrives at a dual formulation, which is then solved by a block coordinate descent algorithm with three blocks. Two of them can be exactly optimized leading to a Sinkhorn-type procedure, while the third block, corresponding to the additional inequality constraints is proposed to be solved by the Newton's method. Some discussions  about accelerating the Newton's method is provided on account of the sparsity of the Hessian.

On the theoretical side, the paper shows two main results: first the fact that entropic regularization has an exponentially decreasing effect on the solution, w.r.t. the regularization parameter and second a guarantee on a sublinear rate of convergence of the algorithm.

The results are finally verified in some numerical experiments.
weaknesses: I generally find the paper interesting, but have few concerns regarding the motivation of the problem, the justifications of the algorithmic choices and the implications of the theoretical results. These are mentioned as questions in the next part.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: - What is $n$ in Theorem 1?
- It is not clear why $f$ is a Lyapunov function. Can the authors explain this part in more detail?
soundness: 4
strengths: - The paper takes initial steps in the natural direction of OT problem research. 
- The theoretical results (such as the convergence rates and bounds) presented herein will likely be referenced in the foreseeable future. So, the results themselves are significant.
summary: The paper presents Entropy regularization formulation and a Sinkhorn algorithm to solve a more generalized class of Optimal Transport (OT) problems --- specifically, OT problems with inequality and equality constraints. 

The writing and presentation are clear and unambiguous for the most part. There is a significant theoretical contribution since the paper considers a larger class of problems than the ones typically considered in OT-related works. However, the motivation of the problem is lacking. Acceleration techniques that are presented in Section 3.1 are not necessarily original since they are an obvious extension of existing techniques. I think there is enough reason to accept the paper since it addresses an OT formulation that is rarely explored. More specifically, this solution is the first of its kind, and it will interest the ICLR audience. Having said that, I do have some reservations about giving a strong acceptance recommendation. 

One of the listed contributions (on Partial Optimal Transport) is only mentioned briefly and not included in the main paper. This seems like an attempt to bypass the page limit constraints, which, in my opinion, should not be encouraged. However, I will leave that to the discretion of the editors.
weaknesses: Firstly, the problem itself is not motivated well in the paper. The authors should consider working on the introduction to establish the relevance of this work better. 

Secondly, the only novelty I see is in Algorithm 1. However, I consider the theoretical contribution itself significant enough to overlook this shortcoming.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 2
questions: 1. Line 187, what does depending only on the LP mean? It is not clear to me that a constant can depend on the linear optimization problem.

2. Line 244, parameter setting can simply your analysis instead of ensuring the efficiency of the algorithm. Can the author improve their justifications of K=O(1), K+L=O(1)?

3. Line 227, in terms of computational complexity, APDAGD can not outperform.
soundness: 2
strengths: The paper extends the Sinkhorn algorithm to unconstrained optimal transport problem to constrained optimal transport problem successfully.
summary: This paper studies the algorithm for solving the optimal transport problem under equality and inequality constraints. The authors derive a Sinkhorn-type algorithm to solve the constrained optimal transport problem.
weaknesses: 1. The paper is hard to follow. The paper lists technical details and equations without further explanations.  For example, which is function L in Line 202? There are no further comments or interpretations of Theorem 2. Sometimes, the author refers to the equation in the appendix, for example, Line 297. The message is not direct. For example, what is Section 3.1 for? Merge some points in Section 1.1 to make the contributions clearer. It is highly recommended that the authors reorganize and polish the paper to make it accessible...

2. What is the motivation for focusing on the constrained optimal transport? Although it could be theoretically interesting to extend Sinkhorn to this kind of optimization problem, the author should make the significance of studying the problem clear.

3. What is the resulting computational complexity of the proposed Sinkhorn-type algorithm?  It seems that only polynomial time can be proved in the paper, which is not attractive enough. The complexity of the Sinkhorn-type algorithm for solving unconstrained optimal transport problems is O(n^2/\epsilon^2); it seems that a comparable complexity should be derived for constrained optimal transport problems.

Meta-review: The paper presents a novel Sinkhorn-type algorithm for constrained optimal transport problems, extending existing OT formulations to include both equality and inequality constraints. While the theoretical foundations are sound, reviewers highlighted weaknesses. A major concern is the lack of motivation and clear exemplification of the method’s relevance to machine learning, particularly its ability to provide new insights or practical improvements for ML applications involving OT. The experimental setup relies on toy scenarios and fails to demonstrate scalability or applicability to real-world ML tasks, such as domain adaptation or large-scale data alignment. Furthermore, the theoretical analysis of computational complexity is incomplete, and the algorithm’s scaling properties in terms of accuracy remain unclear. Despite revisions, these issues persist, with reviewers noting that the contributions—while technically valid—are incremental and primarily extend previous works without providing sufficient broader impact or practical utility for the ICLR community.

Decision: Reject