Considering LP duality, the following equality is established \(\textbf {d}(a, b, M) = \textbf {p}(a, b, M)\)  [1]}. Computing the exact Wasserstein distance in (6) is time consuming. To alleviate this problem, [2]} has introduced an interesting method that regularizes (6) using the entropy of the solution matrix \(H(T)\) , (i.e., \({\mathrm {min}} \big \langle T, M \big \rangle +\gamma H(T)\) ). It has been shown that if \(T^{\prime }_{\gamma }\)  is the solution of the regularized version of (6) and \(\alpha ^{\prime }_{\gamma }\)  is its dual solution in (7), then \(\exists ! u\in \mathbb {R}_+^n\) , \(v\in \mathbb {R}_+^m\)  such that the solution matrix is \(T^{\prime }_{\gamma }=\text{diag}(u)K\text{diag}(v)\)  and \(\alpha ^{\prime }_\gamma =-\log (u)/\gamma +(\log (u)^\top \textbf {1}_n)/({\gamma n})) \textbf {1}_n\)  where, \(K=exp(-M/\gamma )\) . The vectors \(u\)  and \(v\)  are updated iteratively between step 1 and 2 by using the well-known Sinkhorn algorithm as follows: step \(1) u=a/{Kv}\)  and step \(2) v=b/{K^{\top }u}\) , where\(/\)  denotes element-wise division operator [2]}.
