# Subspace Detours Meet Gromov-Wasserstein
## 4 Triangular Coupling as Limit of Optimal Transport Plans for Quadratic Cost
### Solving Hadamard-Wasserstein in the Discrete Setting  
In this part, we derive formulas to solve numerically \(\mathcal{H}\mathcal{W}\) (12). Let \(x_{1},\ldots,x_{n}\in\mathbb{R}^{d}\), \(y_{1},\ldots,y_{m}\in\mathbb{R}^{d}\), \(\alpha\in\Sigma_{m}\), \(\beta\in\Sigma_{m}\), \(p=\sum_{i=1}^{n}\alpha_{i}\delta_{x_{i}}\) and \(q=\sum_{j=1}^{m}\beta_{j}\delta_{y_{j}}\) two discrete measures in \(\mathbb{R}^{d}\). The Hadamard Wasserstein problem (12) becomes in the discrete setting:  
\[\mathcal{H}\mathcal{W}^{2}(p,q) =\inf_{\gamma\in\Pi(p,q)}\ \sum_{i,j}\sum_{k,\ell}\|x_{i}\odot x_{k}-y_{j} \odot y_{\ell}\|_{2}^{2}\ \gamma_{i,j}\gamma_{k,\ell}\] \[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]  
with \(\mathcal{E}(\gamma)=\sum_{i,j}\sum_{k,\ell}\|x_{i}\odot x_{k}-y_{j}\odot y_{ \ell}\|_{2}^{2}\ \gamma_{i,j}\gamma_{k,\ell}\). As denoted in [9], if we note:  
\[\mathcal{L}_{i,j,k,\ell}=\|x_{i}\odot x_{k}-y_{j}\odot y_{\ell}\|_{2}^{2},\]  
then we have:  
\[\mathcal{E}(\gamma)=\langle\mathcal{L}\otimes\gamma,\gamma\rangle,\]  
where \(\otimes\) is defined as:  
\[\mathcal{L}\otimes\gamma=\Big{(}\sum_{k,\ell}\mathcal{L}_{i,j,k,\ell}\gamma_{ k,\ell}\Big{)}_{i,j}\in\mathbb{R}^{n\times m}.\]  
We show in the next proposition a decomposition of \(\mathcal{L}\otimes\gamma\), which allows us to compute this tensor product more efficiently. Proposition 7: _Let \(\gamma\in\Pi(p,q)=\{M\in(\mathbb{R}_{+})^{n\times m},\ M\mathbbm{1}_{m}=p,\ M^{T} \mathbbm{1}_{n}=q\}\), where \(\mathbbm{1}_{n}=(1,\ldots,1)^{T}\in\mathbb{R}^{n}\). Let us note \(X=(x_{i}\odot x_{k})_{i,k}\in\mathbb{R}^{n\times n\times d}\), \(Y=(y_{j}\odot y_{\ell})_{j,\ell}\in\mathbb{R}^{m\times m\times d}\), \(X^{(2)}=(\|X_{i,k}\|_{2}^{2})_{i,k}\in\mathbb{R}^{n\times n}\), \(Y^{(2)}=(\|Y_{j,l}\|_{2}^{2})_{j,l}\in\mathbb{R}^{m\times m}\), and \(\forall t\in\{1,\ldots,d\}\), \(X_{t}=(X_{i,k,t})_{i,k}\in\mathbb{R}^{n\times n}\) and \(Y_{t}=(Y_{j,\ell,t})_{j,\ell}\in\mathbbProof. : First, we can start by writing:  
\[\mathcal{L}_{i,j,k,\ell} =\|x_{i}\odot x_{k}-y_{j}\odot y_{\ell}\|_{2}^{2}\] \[=\|X_{i,k}-Y_{j,\ell}\|_{2}^{2}\] \[=\|X_{i,k}\|_{2}^{2}+\|Y_{j,\ell}\|_{2}^{2}-2\langle X_{i,k},Y_{j, \ell}\rangle\] \[=[X^{(2)}]_{i,k}+[Y^{(2)}]_{j,\ell}-2\langle X_{i,k},Y_{j,\ell}\rangle.\]  
We cannot directly apply proposition 1 from [9] (as the third term is a scalar product), but by performing the same type of computation, we obtain:  
\[\mathcal{L}\otimes\gamma=A+B+C\]  
with  
\[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]0  
\[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]1  
and  
\[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]2 \[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]3 \[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]4 \[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]5  
Finally, we have:  
\[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]6  
From this decomposition, we can compute the tensor product \(\mathcal{L}\otimes\gamma\) with a complexity of \(O(d(n^{2}m+m^{2}n))\) using only multiplications of matrices (instead of \(O(dn^{2}m^{2})\) for a naive computation). Remark 1: _For the degenerated cost function (14), we just need to replace \(X\) and \(Y\) by \(\tilde{X}_{t}=A_{t}^{\frac{1}{2}}X\) and \(\tilde{Y}_{t}=A_{t}^{\frac{1}{2}}Y\) in the previous proposition._  
To solve this problem numerically, we can use the conditional gradient algorithm (Algorithm 2 in [41]). This algorithm only requires to compute the gradient:  
\[=\inf_{\gamma\in\Pi(p,q)}\ \mathcal{E}(\gamma)\]7  
at each step and a classical OT problem. This algorithm is more efficient than solving the quadratic problem directly. Moreover, while it is a non-convex problem, it actually converges to a local stationary point [42]. On Figure 3, we generated 30 points of 2 Gaussian distributions, and computed the optimal coupling of \(\mathcal{HW}_{t}\) for several \(t\). These points have the same uniform weight. We plot the couplings between the points on the second row, and between the projected points on their first coordinate on the first row.