The Sinkhorn updates can be formulated as matrix-vector products for computations between discrete measures, hence the speed-up on computations with GPUs (see [1]}).
Setting \(K \triangleq (e^{-x_i,y_i) /\epsilon })_{ij}\) , \(u \triangleq (e^{x_i) / \epsilon })_i\)  and \(v \triangleq (e^{(y_j) / \epsilon })_j\) .
It reads in the case of balanced OT
\(u\leftarrow \frac{{\phi }(-_())}{K(v\odot )},\qquad v\leftarrow \frac{{\phi }(-_()}{K^\top (u\odot )},\) 
