$$

\begin{aligned}

\min_{\mu \in \mathbb{Z}_+^{X \times Y}}& \sum_{(x,y) \in X \times Y} \mu_{xy}|x-y|^{1/\zeta} \\

\text{s.t. }& \sum_{x \in X} \mu_{xy} = n_x \\

& \sum_{y \in Y} \mu_{xy} = m_y

\text{s.t. }& \sum_{y \in Y} \mu_{xy} = n_x \\

& \sum_{x \in X} \mu_{xy} = m_y

\end{aligned}

$$

$$

\begin{aligned}

V_P = \min_{\mu \geq 0}& \sum_{(x,y) \in X \times Y} \mu_{xy}c_{xy} \\

\text{s.t. }& \sum_{x \in X} \mu_{xy} = n_x \\

& \sum_{y \in Y} \mu_{xy} = m_y

\text{s.t. }& \sum_{y \in Y} \mu_{xy} = n_x \\

& \sum_{x \in X} \mu_{xy} = m_y

\end{aligned}

$$

$$

\begin{aligned}

W_P = \max_{\mu \geq 0}& \sum_{(x,y) \in X \times Y} \mu_{xy}y_{xy} \\

\text{s.t. }& \sum_{x \in X} \mu_{xy} = n_x \\

& \sum_{y \in Y} \mu_{xy} = m_y

\text{s.t. }& \sum_{y \in Y} \mu_{xy} = n_x \\

& \sum_{x \in X} \mu_{xy} = m_y

\end{aligned}

$$

We proceed to describe and implement the algorithm to compute the dual solution.

As already mentioned, the algorithm starts from the matched pairs $(x_0,y_0)$ with no subpairs and assigns the (temporary) values $\psi_{x_0} = c_{x_0 y_0}$ and $\psi_{y_0} = 0,$ i.e. the $x$ type sustains the whole cost of matching.

As already mentioned, the algorithm starts from the matched pairs $(x_0,y_0)$ with no subpairs and assigns the (temporary) values $\phi_{x_0} = c_{x_0 y_0}$ and $\psi_{y_0} = 0,$ i.e. the $x$ type sustains the whole cost of matching.

\leq \min (c_{x_0 y_j} + c_{x_i y_0} - c_{x_0 y_0} , c_{x_i y_j}) - c_{x_j y_j} , \quad \text{for all } 1 \leq i < j \leq p.

$$

Then for all $i \in [p]$ compute the adjustment $ \Delta_i = \sum_{k = i+1}^p \beta_k + \phi_{x_p} - \phi_{x_1}$ and modify the dual variables

Then for all $i \in [p]$ compute the adjustment $ \Delta_i = \sum_{k = i+1}^p \beta_k + \phi_{x_p} - \phi_{x_i}$ and modify the dual variables

$$

\begin{aligned}