Problematic (Debiased Sinkhorn barycenter, proposed by H.Janti et al.): Let $\alpha_1, \ldots, \alpha_K \in \Delta_n$ and $\mathbf{K}=e^{-\frac{\mathrm{C}}{\varepsilon}}$. Let $\pi$ denote a sequence $\pi_1, \ldots, \pi_K$ of transport plans in $\mathbb{R}_{+}^{n \times n}$ and the constraint sets $\mathcal{H}_1=\left\{\pi \mid \forall k, \pi_k \mathbb{1}=\alpha_k\right\}$, and $\mathcal{H}_2=\left\{\pi \mid \forall k \forall k^{\prime}, \pi_k^{\top} \mathbb{1}=\pi_{k^{\prime}} \mathbb{1}\right\}$. The barycenter problem $\min _{\alpha \in \Delta_n} \sum_{k=1}^K w_k \mathrm{~S}_{\varepsilon}\left(\alpha_k, \alpha\right)$ is equivalent to: \begin{align*} & \min _{\substack{\pi \in \mathcal{H}_1 \cap \mathcal{H}_2 \\ d \in \mathbb{R}_{+}^n}}\left[\varepsilon \sum_{k=1}^K w_k \mathrm{KL}\left(\pi_k \mid \mathbf{K} \operatorname{diag}(d)\right) +\frac{\varepsilon}{2}\langle d-\mathbb{1}, \mathbf{K}(d-\mathbb{1})\rangle\right]. \tag{3.} \end{align*} where $\mathrm{KL}(\mathbf{P} \mid \mathbf{Q})=\sum_{i, j} \mathbf{P}_{i j} \log \left(\frac{\mathbf{P}_{i j}}{\mathbf{Q}_{i j}}\right)+\mathbf{Q}_{i j}-\mathbf{P}_{i j}$.
To compute the barycenter numerically, we can use Algorithm 1 in the same paper of Janati.
Comment: According to the, we can observe that Algorithm 1 converges very quickly when $\varepsilon$ is small (in this case, 0.005), corresponding to the case with the smallest number of small eigenvalues. When $\varepsilon$ increases, the algorithm takes much more time to converge, and the kernel $\mathbf{K}$ has more small eigenvalues. Based on the common understanding, this issue might come from the ill-condition of the Gaussian kernel $\mathbf{K} = \exp(-C/\varepsilon)$ when the size of the cost matrix $C$ (the number of discretization) is too large. Then, a part of the optimization problem may be \emph{too flat} in some directions. For instance, the convexity of the map $d \mapsto \langle (d-\mathbb{1}), \mathbf{K}(d-\mathbb{1}) \rangle$ is directly related to the eigenvalues of $\mathbf{K}$ (since the Hessian matrix is $\nabla^2 = 2\mathbf{K}$).
Question: I wonder if there is any good way to minimize the objective function w.r.t the variable $d$ (so that I believe it will make the updating of $\pi$ more robust as a consequence) when the objective function is too flat.
