In the book Computational Optimal Transport, by Peyré and Cuturi, the authors describe an algorithm called Sinkhorn to solve the Entropic Regularized Optimal Transport problem.
The algorithm is quite straight forward, here is a summary.
Let $C_{n \times n}$ be a cost matrix, and $K_{n\times n} = e^{-C}$ the Gibbs kernel matrix. For two vectors $a_{n\times 1}$ and $b_{n\times 1}$, the algorithm is:
- Set $v^{(0)} = [1,...,1]^T \in \mathbb R^n$, and $u^{(0)} = [1,...,1]^T \in \mathbb R^n$;
- Then for $N$ iterations, make $$u^{N+1}_{i}= \frac{a_i}{Kv^{(N)}}\quad \text{and} \quad v^{N+1}_{i}= \frac{b_i}{Ku^{(N+1)}}, \forall i=1,...,n$$
- Finally, calculate $P = diag(u) K diag(v)$.
So, the algorithm consists of some simple calculations. The authors then claim that there is some numerical instability, and that the numbers might get too small. Hence, they recommend doing the calculations in the logarithmic domain.
My question is, how would one modify this algorithm to perform the same calculations in the logarithmic domain? I mean, I could use $K_{modified }:= \log(K) = -C$. But then, how to perform the matrix multiplications and then transform the problem back to the original values?