Hi Cédric.

Thanks for the clarifications and for posting such a detailed document. In fact, the comment in ot.gromov._util.gwggrad regarding the typo in Peyré et al prompted me to check the calculations on my own. I need to understand these functions to adapt them to another project, so I decided to open this report after my calculations didn't agree with either the paper or POT.

Regarding your computations, I now agree that the current implementation of the exact line-search step is correct. However, I'm still not sure about the gradient. My issue is that $\mathbf{p}$ and $\mathbf{q}$ are inputs to the optimization problem, so they are constants with respect to the coupling $\mathbf{T}$, not to mention that $\mathbf{T}$ is also constrained by $\mathbf{T} \cdot \mathbf{1}_m = \mathbf{p}$ and $\mathbf{T}^\top \cdot \mathbf{1}_n = \mathbf{q}$. Considering that:

• $c_{\mathbf{C}, \overline{\mathbf{C}}} = f_1(\mathbf{C}) \mathbf{p} \mathbf{1}_m^\top + \mathbf{1}_n \mathbf{q}^\top f_2(\overline{\mathbf{C}})^\top$ is also constant w.r.t. $\mathbf{T}$,
• $\mathcal{L}(\mathbf{C}, \overline{\mathbf{C}}) \otimes \mathbf{T} = c_{\mathbf{C}, \overline{\mathbf{C}}} - h_1(\mathbf{C}) \mathbf{T} h_2(\overline{\mathbf{C}})^\top$, and
• $\mathcal{E}_L^{GW}(\mathbf{C}, \overline{\mathbf{C}}, \mathbf{T}) = \langle \mathcal{L}(\mathbf{C}, \overline{\mathbf{C}}) \otimes \mathbf{T}, \mathbf{T} \rangle = \langle c, \mathbf{T} \rangle - \langle h_1(\mathbf{C}) \mathbf{T} h_2(\overline{\mathbf{C}})^\top, \mathbf{T} \rangle$,

I find $\nabla_{\mathbf{T}} \mathcal{E}_L^{GW}(\mathbf{C}, \overline{\mathbf{C}}, \mathbf{T})$ by computing the gradient of the two summands above. In the following, I'll write $\mathbf{D}$ instead of $\overline{\mathbf{C}}$ (I have trouble rendering \overline consistently). We get:

• $\nabla_{\mathbf{T}} \langle c_{\mathbf{C},\mathbf{D}}, \mathbf{T} \rangle = c_{\mathbf{C},\mathbf{D}}$ because $c_{\mathbf{C},\mathbf{D}}$ is constant w.r.t. $\mathbf{T}$.
• $\frac{\partial}{\partial T_{pq}} \langle h_1(\mathbf{C}) \mathbf{T} h_2(\mathbf{D})^\top, \mathbf{T} \rangle = \frac{\partial}{\partial T_{pq}} \sum_{ijkl} h_1(C_{ik}) h_2(D_{jl}) T_{ij} T_{kl} = \sum_{kl} h_1(C_{pk}) h_2(D_{ql}) T_{kl} + \sum_{ij} h_1(C_{ip}) h_2(D_{jq}) T_{ij}$. These sums are the $(p,q)$ entries of $h_1(\mathbf{C}) \mathbf{T} h_2(\mathbf{D})^\top$ and $h_1(\mathbf{C})^\top \mathbf{T} h_2(\mathbf{D})$, respectively. Then $\nabla \langle h_1(\mathbf{C}) \mathbf{T} h_2(\mathbf{D})^\top, \mathbf{T} \rangle = h_1(\mathbf{C}) \mathbf{T} h_2(\mathbf{D})^\top + h_1(\mathbf{C})^\top \mathbf{T} h_2(\mathbf{D})$.

Thus, I get that the gradient of $\mathcal{E}_L^{GW}(\mathbf{C}, \mathbf{D}, \mathbf{T})$ is

$$c_{\mathbf{C},\mathbf{D}} - h_1(\mathbf{C}) \mathbf{T} h_2(\mathbf{D})^\top - h_1(\mathbf{C})^\top \mathbf{T} h_2(\mathbf{D})$$

in the general case, and

$$c_{\mathbf{C},\mathbf{D}} - 2 h_1(\mathbf{C}) \mathbf{T} h_2(\mathbf{D})^\top$$

in the symmetric case. In both cases, I'm missing a $c_{\mathbf{C},\mathbf{D}}$ summand compared to $\mathcal{L}(\mathbf{C}, \mathbf{D}) \otimes \mathbf{T} + \mathcal{L}(\mathbf{C}^\top, \mathbf{D}^\top) \otimes \mathbf{T}$ and 2ドル \mathcal{L}(\mathbf{C}, \mathbf{D}) \otimes \mathbf{T}$ in the symmetric case.

Thus, am I right to treat $\mathbf{p}$, $\mathbf{q}$, and ultimately $c_{\mathbf{C}, \overline{\mathbf{C}}}$ as constants w.r.t $\mathbf{T}$? Alternatively, is there a theoretical or practical justification for treating $\mathbf{p}$ and $\mathbf{q}$ as functions of $\mathbf{T}$ even though $\mathbf{p}$ and $\mathbf{q}$ are part of the input data?

Best,
Mario

Hello Mario,

Sorry if I was not explicit enough and actually my previous answer is not exhaustive. We need to be more rigorous to talk about this, because I believe that our two solvers are correct but not really for the reason you mentioned ;)

1. Solving the original GW problem. The analysis in my note aims at solving for equation 1, considering the objective function $\mathcal{E}^{GW}_L$ as a whole. To do so, you cannot formally proceed as you did.
Considering my previous message, you can express the GW cost i.e our objective function as $f(\mathbf{X}) = g(\mathbf{X}) + h(\mathbf{X})$ where $g(\mathbf{X}) = < f_1(\mathbf{C}) \mathbf{X 1}_m \mathbf{1}_m^\top + \mathbf{1}_n \mathbf{1}_n^\top \mathbf{X} f_2(\overline{\mathbf{C}})^\top , \mathbf{X} >$. Indeed, we know that for $\mathbf{T} \in \mathcal{U}(\mathbf{p}, \mathbf{q})$, we have $g(\mathbf{T}) = < c_{\mathbf{C}, \mathbf{\overline{C}}} , \mathbf{T} > = < f_1(\mathbf{C}), \mathbf{p} \mathbf{p}^\top> + < f_2(\mathbf{\overline{C}}), \mathbf{q} \mathbf{q}^\top>$.
However, $g$ remains a quadratic function in its variable. You have to take that into account to express its partial derivatives whose expressions are independent of the set over which you apply them a posteriori.
Saying that $\nabla g( \mathbf{.} ) = \nabla < c_{\mathbf{C}, \mathbf{\overline{C}}} , \mathbf{.} >$ is wrong as you take into account the set over which you apply these matrix valued functions before computing their expressions. So equation 3 is the right way to go to avoid these mistakes.

I hope it is clear because I don't see how to put it more clearly. Maybe the following exaggerated analogy can help: consider $f(x) = x$. you seek for its derivative in $y = {2}$. So you compute $f'(x) = 1$ which is constant and implies that $f'(2) = 1$ too. You do not compute the derivative of $g(x) = f(2)$ i.e $g'(x) = 0$.

2. Solving an equivalent probem to the GW problem. I believe that the algorithm you implemented actually comes down to solving an equivalent problem to the original one that I detailed in the note, and which is implemented in POT. The true reason behind this is that we are not talking about the same objective functions.

We have the following equivalent formulations for the GW problem (1):
$$\min_{\mathbf{T} \in \mathcal{U}(\mathbf{p}, \mathbf{q})} < f_1(\mathbf{C}) \mathbf{T 1}_m \mathbf{1}_m^\top + \mathbf{1}_n \mathbf{1}_n^\top \mathbf{T} f_2(\overline{\mathbf{C}})^\top , \mathbf{T} > + h(\textbf{T})$$

As you mentioned, we have $f_1(\mathbf{C}) \mathbf{T 1}_m \mathbf{1}_m^\top + \mathbf{1}_n \mathbf{1}_n^\top \mathbf{T} f_2(\overline{\mathbf{C}})^\top = c_{\mathbf{C}, \mathbf{\overline{C}}}$ when $\mathbf{T} \in \mathcal{U}(\mathbf{p}, \mathbf{q})$. So we could equivalently consider the following equivalent optimization problem (2):

$$\min_{\mathbf{T} \in \mathcal{U}(\mathbf{p}, \mathbf{q})} < c_{\mathbf{C}, \mathbf{\overline{C}}} , \mathbf{T} > + h(\textbf{T})$$
This objective function does indeed match your computation. For which your gradient computation would be right but we were not talking about the same thing exactly.

Notice that we could also go further by considering the other equivalent problem (3)
$$< f_1(\mathbf{C}), \mathbf{p} \mathbf{p}^\top> + < f_2(\mathbf{\overline{C}}), \mathbf{q} \mathbf{q}^\top> + \min_{\mathbf{T} \in \mathcal{U}(\mathbf{p}, \mathbf{q})} h(\textbf{T})$$

Best,
Cédric

Hi Cédric,

I understand the difference now. Indeed, I thought that the gradient of $\mathcal{E}_{L}^{GW}$ had to be computed for $\mathbf{T} \in \mathcal{U}(\mathbf{p}, \mathbf{q})$ because we want the optimal $\mathbf{T}$ to lie on $\mathcal{U}(\mathbf{p}, \mathbf{q})$.

Thanks for your help!
Mario