

\subsection{\texorpdfstring{$f$}{f}-divergence GEMINI}

We detail here the derivation for 3 $f$-divergences that we previously chose: the \gls*{kl} divergence, the \gls*{tv} distance and the squared Hellinger distance, as well as the generic scenario for any function $f$.

\subsubsection{Generic scenario}

First, we recall that the definition of an $f$-divergence involves a convex function:

\begin{align*}
    f:\mathbb{R}^+ &\rightarrow \mathbb{R}\\
    x&\rightarrow f(x)\\
    \text{s.t.}\quad f(1)&=0,
\end{align*}

between two distributions $p$ and $q$ as described:

\begin{equation*
    D_\text{f-div}(p,q) = \mathbb{E}_{\vec{x} \sim q} \left[ f\left(\frac{p(\vec{x})}{q(\vec{x})}\right)\right].
\end{equation*}

We simply inject this definition in the \gls*{ova}-\gls*{gemini} and directly obtain both an expectation on the cluster assignment $y$ and on the data variable $\vec{x}$. We then merge the writing of the two expectations for the sake of clarity.

\begin{align*}
    \mathcal{I}^\text{ova}_\text{f-div} &= \mathbb{E}_{\p(y)} \left[ D_\text{f-div}(\p(\vec{x}|y) || p_\textup{data}(\vec{x}))\right],\\
    &= \mathbb{E}_{\p(y)} \left[ \mathbb{E}_{p_\textup{data}(\vec{x})}\left[ f\left(\frac{\p(\vec{x}|y)}{p_\textup{data}(\vec{x})}\right)\right]\right],\\
    &= \mathbb{E}_{\p(y),p_\textup{data}(\vec{x})} \left[ f\left(\frac{\p(y|\vec{x})}{\p(y)}\right)\right].\\
\end{align*}

Injecting the $f$-divergence in the \gls*{ovo}-\gls*{gemini} first yields:

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{f-div} &= \mathbb{E}_{\p(\ya),\p(\yb)} \left[ D_\text{f-div}(\p(\vec{x}|\ya) || \p(\vec{x}|\yb))\right],\\
    &= \mathbb{E}_{\p(\ya),\p(\yb)} \left[ \mathbb{E}_{\p(\vec{x}|\yb)}\left[ f\left(\frac{\p(\vec{x}|\ya)}{\p(\vec{x}|\yb)}\right)\right]\right].\\
\end{align*}

Now, by using Bayes theorem, we can perform the inner expectation over the data distribution. We then merge the expectations for the sake of clarity.

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{f-div} &= \mathbb{E}_{\p(\ya),\p(\yb)} \left[ \mathbb{E}_{p_\textup{data}(\vec{x})} \left[ \frac{\p(\yb|\vec{x})}{\p(\yb)}f\left(\frac{\p(\vec{x}|\ya)}{\p(\vec{x}|\yb)}\right)\right]\right],\\
    &= \mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\yb|\vec{x})}{\p(\yb)}f\left(\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\yb|\vec{x})\p(\ya)}\right)\right].\\
\end{align*}

Notice that we also changed the ratio of conditional distributions inside the function by a ratio of posteriors through Bayes' theorem, weighted by the relative cluster proportions.

Now, we can derive into details these equations for the 3 $f$-divergences we focused on: the \gls*{kl} divergence, the \gls*{tv} distance and the squared Hellinger distance.

\subsubsection{Kullback-Leibler divergence}

The function for Kullback-Leibler is $f(t) = t\log t$. We do not need to write the \gls*{ova} equation: it is straightforwardly the usual \gls*{mi}. For the \gls*{ovo}, we inject the function definition by replacing:

\begin{equation*}
    t=\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\yb|\vec{x})\p(\ya)},
\end{equation*}

in order to get:

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{KL} &=  \mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\yb|\vec{x})}{\p(\yb)}\times\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\yb|\vec{x})\p(\ya)}\log{\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\yb|\vec{x})\p(\ya)}}\right].
\end{align*}

We can first simplify the factors outside of the logs:

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{KL}&=\mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\ya|\vec{x})}{\p(\ya)}\log{\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\yb|\vec{x})\p(\ya)}}\right].
\end{align*}

If we use the properties of the log, we can separate the inner term in two sub-expressions:

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{KL} = \mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\ya|\vec{x})}{\p(\ya)}\log{\frac{\p(\ya|\vec{x})}{\p(\ya)}} + \frac{\p(\ya|\vec{x})}{\p(\ya)}\log{\frac{\p(\yb)}{\p(\yb|\vec{x})}}\right].
\end{align*}

Hence, we can use the linearity of the expectation to separate the two terms above. The first term is constant w.r.t. $y_b$, so we can remove this variable from the expectation among the subscripts:

\begin{multline*}
    \mathcal{I}^\text{ovo}_\text{KL} = \mathbb{E}_{\p(\ya),p_\textup{data}(\vec{x})} \left[ \frac{\p(\ya|\vec{x})}{\p(\ya)}\log{\frac{\p(\ya|\vec{x})}{\p(\ya)}} \right] + \mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\ya|\vec{x})}{\p(\ya)}\log{\frac{\p(\yb)}{\p(\yb|\vec{x})}} \right].
\end{multline*}

Since the variables $y_a$ and $y_b$ are independent, we can use the fact that:

\begin{equation*}
    \mathbb{E}_{\p(\ya)}\left[\frac{\p(\ya|\vec{x})}{\p(\ya)}\right] = \int \p(\ya) \frac{\p(\ya|\vec{x})}{\p(\ya)} dy_a = 1,
\end{equation*}

inside the second term to reveal the final form of the equation:

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{KL} = \mathbb{E}_{p_\textup{data}(\vec{x}),\p(y)} \left[ \frac{\p(y|\vec{x})}{\p(y)}\log{\frac{\p(y|\vec{x})}{\p(y)}} \right] + \mathbb{E}_{p_\textup{data}(\vec{x}),\p(y)} \left[ \log{\frac{\p(y)}{\p(y|\vec{x})}}\right].
\end{align*}

Notice that since both terms did not compare one cluster assignment $y_a$ against another $y_b$, we can switch to the same common variable $y$. Both terms are in fact \gls*{kl} divergences depending on the cluster assignment $y$. The first is the reverse of the second. This sum of \gls*{kl} divergences is sometimes called the \emph{symmetric} \gls*{kl}, and so can we write in two ways the \gls*{kl}-\gls*{ovo} \gls*{gemini}:

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{KL} &= \mathbb{E}_{p_\textup{data}(\vec{x})} \left[ D_\text{KL} (\p(y|\vec{x}) || \p(y))\right] + \mathbb{E}_{p_\textup{data}(\vec{x})} \left[ D_\text{KL} (\p(y) || \p(y|\vec{x}))\right],\\
    &= \mathbb{E}_{p_\textup{data}(\vec{x})} \left[ D_\text{KL-sym}(\p(y|\vec{x}) || \p(y))\right].\\
\end{align*}

We can also think of this equation as the usual \gls*{mi} with an additional term based on the reversed \gls*{kl} divergence.

\subsubsection{Total Variation distance}

For the total variation, the function is $f(t)=\frac{1}{2} |t-1|$. Thus, the \gls*{ova} \gls*{gemini} is:

\begin{equation*}
    \mathcal{I}^\text{ova}_\text{TV} = \frac{1}{2}\mathbb{E}_{\p(y),p_\textup{data}(\vec{x})} \left[ |\frac{\p(y|\vec{x})}{\p(y)}-1|\right].
\end{equation*}

And the \gls*{ovo} is:

\begin{align*}
    \mathcal{I}^\text{ovo}_\text{TV} &=\frac{1}{2}\mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\yb|\vec{x})}{\p(\yb)}|\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\yb|\vec{x})\p(\ya)}-1|\right],\\
    &=\frac{1}{2}\mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ |\frac{\p(\ya|\vec{x})}{\p(\ya)} - \frac{\p(\yb|\vec{x})}{\p(\yb)} | \right].
\end{align*}

We did not find any further simplification of these equations.

\subsubsection{Squared Hellinger distance}

Finally, the squared Hellinger distance is based on $f(t)=2(1-\sqrt{t})$. Hence the \gls*{ova} unfolds as:

\begin{align*}
    \mathcal{I}^\text{ova}_{\text{H}^2}&= \mathbb{E}_{\p(y),p_\textup{data}(\vec{x})} \left[ 2\left(1-\sqrt{\frac{\p(y|\vec{x})}{\p(y)}}\right)\right],\\
    &= 2-2\mathbb{E}_{p_\textup{data}(\vec{x}),\p(y)} \left[\sqrt{\frac{\p(y|\vec{x})}{\p(y)}}\right].
\end{align*}

The idea of the squared Hellinger-\gls*{ova} \gls*{gemini} is therefore to minimise the expected square root of the relative certainty between the posterior and cluster proportion.

For the \gls*{ovo} setting, the definition yields:

\begin{equation*}
    \mathcal{I}^\text{ovo}_{\text{H}^2}=\mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\yb|\vec{x})}{\p(\yb)}\times 2 \times\left(1-\sqrt{\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\yb|\vec{x})\p(\ya)}}\right)\right],
\end{equation*}

which we can already simplify by putting the constant 2 outside of the expectation, and by inserting all factors inside the square root before simplifying and separating the expectation:

\begin{align*}
    \mathcal{I}^\text{ovo}_{\text{H}^2} &= 2\mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\yb|\vec{x})}{\p(\yb)} - \frac{\p(\yb|\vec{x})}{\p(\yb)}\sqrt{\frac{\p(\ya|\vec{x})\p(\yb)}{\p(\ya)\p(\yb|\vec{x})}}\right],\\
    &= 2\mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \frac{\p(\yb|\vec{x})}{\p(\yb)} \right] - 2\mathbb{E}_{\p(\ya),\p(\yb),p_\textup{data}(\vec{x})} \left[ \sqrt{\frac{\p(\ya|\vec{x}) \p(\yb|\vec{x})}{\p(\ya) \p(\yb)}}\right].
\end{align*}

We can replace the first term by the constant 1, as shown for the \gls*{kl}-\gls*{ovo} derivation. Since we can split the square root into the product of two square roots, we can apply twice the expectation over $y_a$ and $y_b$ because these variables are independent:

\begin{equation*}
    \mathcal{I}^\text{ovo}_{\text{H}^2} = 2-2\mathbb{E}_{p_\textup{data}(\vec{x})} \left[ \mathbb{E}_{\p(y)} \left[\sqrt{\frac{\p(y|\vec{x})}{\p(y)}}\right]^2\right].
\end{equation*}

To avoid computing this squared expectation, we use the equation of the variance $\mathbb{V}$ to replace it. Thus:

\begin{align*}
    \mathcal{I}^\text{ovo}_{\text{H}^2} &= 2-2\mathbb{E}_{p_\textup{data}(\vec{x})} \left[ \mathbb{E}_{\p(y)}\left[\frac{\p(y|\vec{x})}{\p(y)}\right] - \mathbb{V}_{\p(y)}\left[\sqrt{\frac{\p(y|\vec{x})}{\p(y)}}\right]\right],\\
    &= 2 - 2\mathbb{E}_{p_\textup{data}(\vec{x})}\left[ \mathbb{E}_{\p(y)} \left[ \frac{\p(y|\vec{x})}{\p(y)}\right]\right] + 2 \mathbb{E}_{p_\textup{data}(\vec{x})} \left[ \mathbb{V}_{\p(y)} \left[ \sqrt{\frac{\p(y|\vec{x})}{\p(y)}}\right]\right].
\end{align*}

Then, for the same reason as before, the second term is worth 1, which cancels the first constant. We therefore end up with:

\begin{equation*}
    \mathcal{I}^\text{ovo}_{\text{H}^2}= 2\mathbb{E}_{p_\textup{data}(\vec{x})} \left[ \mathbb{V}_{\p(y)}\left[\sqrt{\frac{\p(y|\vec{x})}{\p(y)}}\right]\right].
\end{equation*}

Similar to the \gls*{kl}-\gls*{ovo} case, the squared Hellinger \gls*{ovo} converges to an \gls*{ova} setting, i.e. we only need information about the cluster distribution itself without comparing it to another. Furthermore, the idea of maximising the variance of the cluster assignments is straightforward for clustering.

\subsection{Maximum Mean Discrepancy}

When using an IPM with a family of functions that project an input of $\mathcal{X}$ to the unit ball of an RKHS $\mathcal{H}$, the \gls*{ipm} becomes the \gls*{mmd} distance.

\begin{align*}
\text{MMD}(p,q) &= \sup_{f: ||f||_\mathcal{H}\leq 1} \mathbb{E}_{\vec{x}_a\sim p}[f(\vec{x}_a)] - \mathbb{E}_{\vec{x}_b \sim q} [f(\vec{x}_b)],\\
&= \| \mathbb{E}_{\vec{x}_a\sim p} [\varphi(\vec{x}_a)] - \mathbb{E}_{\vec{x}_b \sim q}[\varphi(\vec{x}_b)]\|_{\mathcal{H}},\\
\end{align*}

where $\varphi$ is a embedding function of the RKHS.

By using a kernel function $\kappa(\vec{x}_a,\vec{x}_b) = <\varphi(\vec{x}_a), \varphi(\vec{x}_b)>$, we can express the square of this distance thanks to inner product space properties \cite{gretton_kernel_2012}:

\begin{align*}
\text{MMD}^2 (p,q) &= \mathbb{E}_{\vec{x}_a, \vec{x}_a^\prime \sim p}[\kappa(\vec{x}_a, \vec{x}_a^\prime)] + \mathbb{E}_{\vec{x}_b,\vec{x}_b^\prime \sim q}[\kappa(\vec{x}_b, \vec{x}_b^\prime)] - 2 \mathbb{E}_{\vec{x}_a\sim p, \vec{x}_b\sim q}[\kappa(\vec{x}_a, \vec{x}_b)].
\end{align*}


Now, we can derive each term of this equation using our distributions $p\equiv\p(\vec{x}|y)$ and $q\equiv p_\textup{data}(\vec{x})$ for the \gls*{ova} case, and $p\equiv\p(\vec{x}|\ya), q\equiv\p(\vec{x}|\yb)$ for the \gls*{ovo} case. In both scenarios, we aim at finding an expectation over the data variable $x$ using only the respectively known and estimable terms $\p(y|\vec{x})$ and $\p(y)$.

\paragraph{OvA scenario}

For the first term, we use Bayes' theorem twice to get an expectation over two variables $\vec{x}_a$ and $\vec{x}_b$ drawn from the data distribution.

\begin{equation*}
    \begin{split}
        \mathbb{E}_{\vec{x}_a, \vec{x}_a^\prime \sim p}&= \mathbb{E}_{\vec{x}_a, \vec{x}_a^\prime \sim \p(\vec{x}|y)} \left[ \kappa(\vec{x}_a,\vec{x}_a^\prime)\right],\\
        &= \mathbb{E}_{\vec{x}_a,\vec{x}_a^\prime \sim p_\textup{data}(\vec{x})} \left[ \frac{p_\theta(y|\vec{x}_a)p_\theta(y|\vec{x}_a^\prime)}{\p(y)^2} \kappa(\vec{x}_a,\vec{x}_a^\prime)\right].
    \end{split}
\end{equation*}

For the second term, we do not need to perform anything particular as we directy get an expectation over the data variabes $\vec{x}_a$ and $\vec{x}_b$.

\begin{equation*}
    \mathbb{E}_{\vec{x}_b, \vec{x}_b^\prime \sim q}= \mathbb{E}_{\vec{x}_b, \vec{x}_b^\prime \sim p_\textup{data}(\vec{x})} \left[ \kappa(\vec{x}_b,\vec{x}_b^\prime)\right].
\end{equation*}

The last term only needs Bayes theorem once, for the distribution $q$ is directly replaced by the data distribution $p_\textup{data}(\vec{x})$:

\begin{align*}
    \mathbb{E}_{\vec{x}_a\sim p, \vec{x}_b \sim q}&=\mathbb{E}_{\vec{x}_a \sim \p(\vec{x}|y), \vec{x}_b \sim p_\textup{data}(\vec{x})}\left[ \kappa(\vec{x}_a,\vec{x}_b)\right],\\
    &= \mathbb{E}_{\vec{x}_a, \vec{x}_b \sim px} \left[ \frac{p_\theta(y|\vec{x}_a)}{\p(y)} \kappa(\vec{x}_a,\vec{x}_b)\right].
\end{align*}

Note that for the last term, we could replace $p_\theta(y|\vec{x}_a)$ by $p_\theta(y|\vec{x}_b)$; that would not affect the result since $\vec{x}_a$ and $\vec{x}_b$ are independently drawn from $p_\textup{data}(\vec{x})$. We thus replace all terms, and do not forget to put a square root on the entire sum since we have computed so far the squared \gls*{mmd}:

\begin{align*}
        \mathcal{I}_\text{MMD}^\text{ova} &= \mathbb{E}_{\p(y)} \left[ \text{MMD}(\p(\vec{x}|y),p_\textup{data}(\vec{x}))\right],\\
        &=\mathbb{E}_{\p(y)} \left[ \left(\mathbb{E}_{\vec{x}_a,\vec{x}_a^\prime \sim p_\textup{data}(\vec{x})} \left[ \frac{p_\theta(y|\vec{x}_a)p_\theta(y|\vec{x}_a^\prime)}{\p(y)^2} \kappa(\vec{x}_a,\vec{x}_a^\prime)\right] \right.\right.\\
        &\left.\left.\qquad+ \mathbb{E}_{\vec{x}_b, \vec{x}_b^\prime \sim p_\textup{data}(\vec{x})} \left[ \kappa(\vec{x}_b,\vec{x}_b^\prime)\right] -2 \mathbb{E}_{\vec{x}_a, \vec{x}_b \sim p_\textup{data}(\vec{x})} \left[ \frac{p_\theta(y|\vec{x}_a)}{\p(y)} \kappa(\vec{x}_a,\vec{x}_b)\right] \right)^{\frac{1}{2}}  \right].\\
\end{align*}

Since all variables $\vec{x}_a$, $\vec{x}_a^\prime$, $\vec{x}_b$ and $\vec{x}_b^\prime$ are independently drawn from the same distribution $p_\textup{data}(\vec{x})$, we can replace all of them by the variables $\vec{x}$ and $\vec{x}^\prime$. We then use the linearity of the expectation and factorise by the kernel $\kappa(\vec{x},\vec{x}^\prime)$:

\begin{equation*}
        \mathcal{I}_\text{MMD}^\text{ova} = \mathbb{E}_{\p(y)} \left[ \mathbb{E}_{\vec{x},\vec{x}^\prime \sim p_\textup{data}(\vec{x})}\left[ \kappa(\vec{x},\vec{x}^\prime) \left( \frac{p_\theta(y|\vec{x})p_\theta(y|\vec{x}^\prime)}{\p(y)^2} + 1 - 2\frac{p_\theta(y|\vec{x})}{\p(y)}\right) \right]^{\frac{1}{2}}\right].
\end{equation*}

\paragraph{OvO scenario}

The two first terms of the \gls*{ovo} \gls*{mmd} are the same as the first term of the \gls*{ova} setting, with a simple subscript $a$ or $b$ at the appropriate place. Only the negative term changes. We once again use Bayes' theorem twice:

\begin{align*}
        \mathbb{E}_{\vec{x}_a \sim p, \vec{x}_b \sim q} [\kappa(\vec{x}_a, \vec{x}_b)]&= \mathbb{E}_{\vec{x}_a \sim \p(\vec{x}|\ya), \vec{x}_b \p(\vec{x}|\yb)} \left[ \kappa(\vec{x}_a,\vec{x}_b) \right],\\
        &= \mathbb{E}_{\vec{x}_a, \vec{x}_b \sim p_\textup{data}(\vec{x})} \left[\frac{p_\theta(y_a|\vec{x}_a)}{\p(\ya)}\frac{p_\theta(y_b|\vec{x}_b)}{\p(\yb)} \kappa(\vec{x}_a,\vec{x}_b) \right].
\end{align*}

The final sum is hence similar to the \gls*{ova}:

\begin{align*}
        \mathcal{I}_\text{MMD}^\text{ovo} &= \mathbb{E}_{\p(\ya), \p(\yb)} \left[ \text{MMD} (\p(\vec{x}|\ya), \p(\vec{x}|\yb))\right],\\
        &= \mathbb{E}_{\p(\ya), \p(\yb)} \left[ \mathbb{E}_{\vec{x},\vec{x}^\prime \sim p_\textup{data}(\vec{x})} \left[ \kappa(\vec{x},\vec{x}^\prime) \left( \frac{\p(\ya|\vec{x}) p_\theta(y_a|\vec{x}^\prime)}{\p(\ya)^2} + \frac{\p(\yb|\vec{x}) p_\theta(y_b|\vec{x}^\prime)}{\p(\yb)^2}\right.\right.\right.\\&\quad\left.\left.\left.-2 \frac{\p(\ya|\vec{x}) p_\theta(y_b|\vec{x}^\prime) }{\p(\ya) \p(\yb)} \right) \right]^{\frac{1}{2}}\right].
\end{align*}

\subsection{Wasserstein distance}
\label{app:wasserstein_convergence}

To compute the Wasserstein distance between the distributions $p_\theta(\vec{x}|y=k)$, we estimate it using approximate distributions. We replace $p_\theta(\vec{x}|y=k)$ by a weighted sum of Dirac measures on specific samples $\vec{x}_i$: $p_N^k$:

\begin{equation*}
    p_\theta(\vec{x}|y=k) \approx \sum_{i=1}^N m_i^k\delta_{\vec{x}_i} = p_N^k,
\end{equation*}

where $\{m_i^k\}_{i=1}^N$ is the set of weights. We now show that computing the Wasserstein distance between these approximates converges to the correct distance. We first need to show that $p_N^k$ weakly converges to $p$. To that end, let $f$ be any bounded and continuous function. Computing the expectation of such through $p_\theta$ is:

\begin{equation*}
\mathbb{E}_{\vec{x}\sim p_\theta(\vec{x}|y=k)}[f(\vec{x})] = \int_{\mathcal{X}}f(\vec{x})p_\theta(\vec{x}|y=k)d\vec{x},
\end{equation*}

which can be estimated using self-normalised importance sampling~\cite[Chapter 9]{owen_monte_2009}. The proposal distribution we take for sampling is $p_\textup{data}(\vec{x})$. Although we cannot evaluate both $\p(\vec{x}|y)$ and $p_\textup{data}(\vec{x})$ up to a constant, we can evaluate their ratio up to a constant which is sufficient:

\begin{align*}
\mathbb{E}_{\vec{x}\sim p_\theta(\vec{x}|y=k)}[f(\vec{x})]&= \int_{\mathcal{X}}f(\vec{x})\frac{p_\theta(\vec{x}|y=k)}{p_\textup{data}(\vec{x})}p_\textup{data}(\vec{x}) d\vec{x},\\
&= \int_{\mathcal{X}}f(\vec{x})\frac{p_\theta(y=k|\vec{x})}{p_\theta(y=k)}p_\textup{data}(\vec{x}) d\vec{x},\\
&\approx \sum_{i=1}^N f(\vec{x}_i) \frac{p_\theta(y=k|\vec{x}=\vec{x}_i)}{\sum_{j=1}^N p_\theta(y=k|\vec{x}=\vec{x}_j)}.
\end{align*}

Now, by noticing in the last line that the importance weights are self normalised and add up to 1, we can identify them as the point masses of our previous Dirac approximations:

\begin{equation*}
m_i^k = \frac{p_\theta(y=k|\vec{x}=\vec{x}_i)}{\sum_{j=1}^N p_\theta(y=k|\vec{x}=\vec{x}_j)}.
\end{equation*}

This allows to write that the Monte Carlo estimation through importance sampling of the expectation w.r.t $p_\theta(\vec{x}|y=k)$ is directly the expectation taken on the discrete approximation $p_N^k$. We can conclude that there is a convergence between the two expectations owing to the law of large numbers:

\begin{equation*}
\lim_{N\rightarrow +\infty}\mathbb{E}_{\vec{x}\sim p_N^k}[f(\vec{x})]  = \mathbb{E}_{\vec{x}\sim p_\theta(\vec{x}|y=k)}[f(\vec{x})].
\end{equation*}

Since $f$ is bounded and continuous, the portmanteau theorem~\cite{billingsley_convergence_2013} states that $p_N^k$ weakly converges to $p_\theta(\vec{x}|y=k)$ when defining the importance weights as the normalised predictions cluster-wise.

To conclude, when two series of measures $p_N$ and $q_N$ weakly converge respectively to $p$ and $q$, so does their Wasserstein distance ~\cite[Corollary 6.9]{villani_optimal_2009}, hence:

\begin{equation}
    \lim_{N\rightarrow+\infty}\mathcal{W}_c(p_N^{k_1},p_N^{k_2})= \mathcal{W}_c\left(p_\theta(\vec{x}|y=k_1)\|p_\theta(\vec{x}|y=k_2)\right).
\end{equation}













\subsection{Generative process of Gaussian and Student Mixture}

We describe here the generative protocol for the Gaussian and Student mixture dataset. Each cluster distribution is centered around a mean $\mu_i$ which proximity is controlled by a scalar $\alpha$. For simplicity, all covariance matrices are the identity scaled by a scalar $\sigma$. We define:

\begin{equation*}
    \mu_1 = [\alpha, \alpha],\quad\mu_2 = [\alpha, -\alpha],\quad\mu_3 = [-\alpha, \alpha],\quad\mu_4 = [-\alpha, -\alpha].
\end{equation*}

To sample from a multivariate Student-t distribution, we first draw samples $x$ from a centered multivariate Gaussian distribution. We then sample another variable $u$ from a $\chi^2$-distribution using the degrees of freedom $\rho$ as parameter. Finally, $x$ is multiplied by $\sqrt{\frac{\rho}{u}}$, yielding samples from the Student-t distribution.

\subsection{Extended experiment}

\input{figs_tex/tab_complete_gstm}

For the main experiment, we fixed $\sigma=1$, $\alpha=5$ and the degree of freedom $\rho=1$. We further tested our method when training a \gls*{mlp} for 4 or 8 clusters. For both cases, we also considered $\rho=1$ and $\rho=2$. We include the complete results in Table~\ref{tab:complete_gstm_exp}.

\subsection{Models definition}
Let us consider a mixture of two Gaussian distributions, both with different means $\mu_0$ and $\mu_1$, s.t. $\mu_0<\mu_1$ and of same standard deviation $\sigma$:

\begin{equation*}
p(x|y=0) = \mathcal{N}(x|\mu_0,\sigma^2), p(x|y=1) = \mathcal{N}(x|\mu_1,\sigma^2) ,
\end{equation*}

where $y$ is the cluster assignment. We take balanced clusters proportions, i.e. $p(y=0)=p(y=1)=\frac{1}{2}$. This first model is the basis that generated the complete dataset $p(x)$. When performing clustering with our discriminative model, we are not aware of the distribution. Consequently: we create other models. We want to compute the difference of mutual information between two decision boundaries that discriminative models $p_\theta(y|x)$ may yield.

We define two decision boundaries: one which splits evenly the data space called $p_A$ and another which splits it on a closed set $p_B$:

\begin{equation*}\label{eq:p_a}
p_A(y=1|x) = \left\{ \begin{array}{c c}
1-\epsilon&x>\frac{\mu_1-\mu_0}{2}\\
\epsilon&\text{otherwise}
\end{array}\right. ,
\end{equation*}

\begin{equation}\label{eq:p_b}
p_B(y=1|x) = \left\{ \begin{array}{c c}
1-\epsilon&x \in [\mu_0, \mu_1]\\
\epsilon&\text{otherwise}
\end{array}\right. .
\end{equation}

Our goal is to show that both models $p_A$ and $p_B$ will converge to the same value of mutual information as $\epsilon$ converges to 0.

\subsection{Computing cluster proportions}
\subsubsection{Cluster proportion of the correct decision boundary}

To compute the cluster proportions, we estimate with samples $x$ from the distribution $p_\text{data}(x)$. Since we are aware for this demonstration of the true nature of the data distribution, we can use $p(x)$ for sampling. Consequently, we can compute the two marginals:

\begin{align*}
p_A(y=1) &= \int_\mathcal{X} p(x) p_A(y=1|x)dx ,\\
&= \int_{-\infty}^{\frac{\mu_1-\mu_0}{2}}p(x)\epsilon dx + \int_{\frac{\mu_1-\mu_0}{2}}^{+\infty}p(x)(1-\epsilon)dx ,\\
&=\epsilon \left(\int_{-\infty}^{\frac{\mu_1-\mu_0}{2}}p(x)dx \right) + (1-\epsilon) \left( \int_{\frac{\mu_1-\mu_0}{2}}^{+\infty}p(x)dx\right) ,\\
&= \frac{1}{2} .
\end{align*}

\subsubsection{Cluster proportion of the misplaced decision boundary}
For the misplaced decision boundary, the marginal is different:

\begin{equation}\begin{split}
p_B(y=1) &= \int_\mathcal{X} p(x) p_B(y=1|x)dx ,\\
&= \epsilon\left(\int_{-\infty}^{\mu_0}p(x)dx + \int_{\mu_1}^{+\infty}p(x)dx\right) + (1-\epsilon) \int_{\mu_0}^{\mu_1}p(x)dx ,\\
&=\epsilon \left(1-\int_{\mu_0}^{\mu_1}p(x)dx\right) + (1-\epsilon)\int_{\mu_0}^{\mu_1}p(x)dx .
\end{split}\end{equation}

Here, we simply introduce a new variable named $\beta$ that will be a shortcut for noting the proportion of data between $\mu_0$ and $\mu_1$:

\begin{equation*
\beta = \int_{\mu_0}^{\mu_1}p(x)dx .
\end{equation*}

And so can we simply write the cluster proportion of decision boundary model B as:

\begin{equation*}\label{eq:py_b}
p_B(y=1) = \epsilon (1-\beta) + (1-\epsilon)\beta ,
\end{equation*}

Leading to the summary of proportions in Table~\ref{tab:proportions}. For convenience, we will write the proportions of model B using the shortcuts:

\begin{equation*}\label{eq:pi_b}
\pi_B = p_B(y=1) = \epsilon + \beta(1-2\epsilon) ,
\end{equation*}
\begin{equation*
\bar{\pi}_B = p_B(y=0) = 1-\epsilon - \beta(1-2\epsilon) .
\end{equation*}

\input{figs_tex/tab_clustering_proportions}

\subsection{Computing the KL divergences}
\subsubsection{Correct decision boundary}

We first start by computing the Kullback-Leibler divergence for some arbitrary value of $x\in\mathbb{R}$:

\begin{equation*}
D_\text{KL}(p_A(y|x)||p_A(y))= \sum_{i=0}^1 p_A(y=i|x) \log{\frac{p_A(y=i|x)}{p_A(y=i)}} .
\end{equation*}

We now need to detail the specific cases, for the value of $p(y=i|x)$ is dependent on $x$. We start $\forall x <\frac{\mu_1-\mu_0}{2}$:

\begin{align*}
D_\text{KL}(p_A(y|x)||p_A(y))&= p_A(y=0|x)\log{\frac{p_A(y=0|x)}{\frac{1}{2}}} + p_A(y=1|x)\log{\frac{p_A(y=1|x)}{\frac{1}{2}}} ,\\
&=(1-\epsilon)\log{2(1-\epsilon)} + \epsilon \log{2\epsilon} .
\end{align*}

The opposite case, $\forall x\geq \frac{\mu_1-\mu_0}{2}$ yields:
\begin{align*}
D_\text{KL}(p_A(y|x)||p_A(y))&= p_A(y=0|x)\log{\frac{p_A(y=0|x)}{\frac{1}{2}}} + p_A(y=1|x)\log{\frac{p_A(y=1|x)}{\frac{1}{2}}} ,\\
&=\epsilon\log{2\epsilon} + (1-\epsilon) \log{2(1-\epsilon)} .
\end{align*}

Since both cases are equal, we can write down:

\begin{equation}\label{eq:correct_kl_div}
D_\text{KL}(p_A(y|x)||p_A(y)) = \epsilon\log{2\epsilon} + (1-\epsilon)\log{2(1-\epsilon)} ,\forall x\in\mathbb{R} .
\end{equation}

\subsubsection{Misplaced boundary}

We proceed to the same detailing of the Kullback-Leibler divergence computation for the misplaced decision boundary. We start with the set $x\in [\mu_0,\mu_1]$:

\begin{align*}
D_\text{KL}(p_B(y|x)||p_B(y))&= p_B(y=0|x) \log{\frac{p_B(y=0|x)}{p_B(y=0)}} + p_B(y=1|x)\log{\frac{p_B(y=1|x)}{p_B(y=1)}} ,\\
&=\epsilon \log{\frac{\epsilon}{\bar{\pi}_B}} + (1-\epsilon) \log{\frac{1-\epsilon}{\pi_B}} .
\end{align*}

When $x$ is out of this set, the divergence becomes:

\begin{align*}
D_\text{KL}(p_B(y|x)||p_B(y))&=p_B(y=0|x) \log{\frac{p_B(y=0|x)}{p_B(y=0)}} + p_B(y=1|x)\log{\frac{p_B(y=1|x)}{p_B(y=1)}} ,\\
&=(1-\epsilon)\log{\frac{1-\epsilon}{\bar{\pi}_B}} + \epsilon \log{\frac{\epsilon}{\pi_B}} .
\end{align*}

To fuse the two results, we will write the KL divergence as such:

\begin{equation}\label{eq:odd_kl_div}
D_\text{KL}(p_B(y|x)||p_B(y)) = \epsilon\log{\epsilon}+ (1-\epsilon)\log(1-\epsilon) - C(x) ,\forall x\in\mathbb{R} ,
\end{equation}

where $C(x)$ is a constant term depending on $x$ defined by:

\begin{equation}\label{eq:odd_constant}
C(x)=\left\{\begin{array}{c c}
\epsilon\log{\bar{\pi}_B} + (1-\epsilon)\log{\pi_B}&x\in[\mu_0,\mu_1]\\
\epsilon\log{\pi_B} + (1-\epsilon)\log{\bar{\pi}_B}&x\in\mathbb{R}\setminus[\mu_0,\mu_1]\\
\end{array}  .\right.
\end{equation}

For simplicity of later writings, we will shorten the notations by:

\begin{equation*}
C(x)=\left\{ \begin{array}{c c}
\alpha_1&x\in[\mu_0,\mu_1]\\
\alpha_0&x\in\mathbb{R}\setminus[\mu_0,\mu_1]
\end{array}\right. .
\end{equation*}

\subsection{Evaluating the mutual information}

\subsubsection{Correct decision boundary}

We inject the value of the Kullback-Leibler divergence from Eq.~(\ref{eq:correct_kl_div}) inside an expectation performed over the data distribution $p_\text{data}(x)$:

\begin{align}
\mathcal{I}_A(X;Y) &= \mathbb{E}_{x\sim p_\text{data}(x)}\left[ D_\text{KL}(p_A(y|x)||p_A(y))\right] ,\\
&=\int_\mathcal{X}p(x) \left(\epsilon \log(2\epsilon) + (1-\epsilon)\log(2(1-\epsilon)) \right)dx ,\\
&=\epsilon \log(2\epsilon) + (1-\epsilon)\log(2(1-\epsilon)) .\label{eq:mi_a}
\end{align}

Since the KL divergence was independent of $x$, we could leave the constant outside of the integral which is equal to 1.

We can assess the coherence of Eq.~(\ref{eq:mi_a}) since its limit as $\epsilon$ approaches 0 is $\log 2$. In terms of bits, this is the same as saying that the information on $X$ directly gives us information on the $Y$ of the cluster.

\subsubsection{Odd decision boundary}

We inject the value of the KL divergence from Eq.~(\ref{eq:odd_kl_div}) inside the expectation of the mutual information:

\begin{align*}
\mathcal{I}_B(X;Y) &= \mathbb{E}_{x\sim p_\text{data}(x)}\left[ D_\text{KL}(p_B(y|x)||p_B(y))\right] ,\\
&=\int_{\mathcal{X}}p(x)\left(\epsilon\log\epsilon +(1-\epsilon)\log(1-\epsilon) -C(x)\right) ,dx\\
&=\epsilon\log\epsilon +(1-\epsilon)\log(1-\epsilon) - \int_{\mathcal{X}}p(x)C(x)dx .
\end{align*}

The first terms are constant with respect to $x$ and the integral of $p(x)$ over $\mathcal{X}$ adds up to 1. We finally need to detail the expectation of the constant $C(x)$ from Eq.~(\ref{eq:odd_constant}):

\begin{align*}
\mathbb{E}_x[C(x)] &= \int_{-\infty}^{\mu_0}C(x)p(x)dx + \int_{\mu_0}^{\mu_1}C(x)p(x)dx + \int_{\mu_1}^{+\infty}C(x)p(x)dx ,\\
&= \alpha_0\left(\int_{-\infty}^{\mu_0}p(x)dx + \int_{\mu_1}^{+\infty}p(x)dx \right) + \alpha_1\int_{\mu_0}^{\mu_1}p(x)dx ,\\
&=\alpha_0(1-\beta) + \alpha_1\beta .
\end{align*}

This can be further improved by unfolding the description of $\alpha_0$ and $\alpha_1$ from Eq.~(\ref{eq:odd_constant}):

\begin{align*}
\alpha_0(1-\beta)+\beta\alpha_1&= \alpha_0 +\beta(\alpha_1-\alpha_0),\\
&=\epsilon\log{\pi_B}+(1-\epsilon)\log{\bar{\pi}_B} + \beta\left[\epsilon\log{\bar{\pi}_B}+(1-\epsilon)\log{\pi_B}\right.\\&\quad\left.-\epsilon\log{\pi_B}- (1-\epsilon)\log{\bar{\pi}_B}\right],\\
&= \left[1-\epsilon+\beta\epsilon-\beta+\beta\epsilon\right]\log{\bar{\pi}_B} + \left[\epsilon+\beta-\beta\epsilon-\beta\epsilon\right]\log{\pi_B},\\
&=\log{\bar{\pi}_B} + \left[2\beta\epsilon-\beta-\epsilon\right]\log{\frac{\bar{\pi}_B}{\pi_B}} .
\end{align*}

We can finally write down the mutual information for the model B:

\begin{equation}\label{eq:mi_b}
\mathcal{I}_B(X;Y) = \epsilon\log{\epsilon} +(1-\epsilon)\log(1-\epsilon) -\log{\bar{\pi}_B} - \left[2\beta\epsilon-\beta-\epsilon\right]\log{\frac{\bar{\pi}_B}{\pi_B}} .
\end{equation}

\subsection{Differences of mutual information}

Now that we have the exact value of both mutual informations, we can compute their differences:

\begin{align*}
\Delta_\mathcal{I} &= \mathcal{I}_A(X;Y) - \mathcal{I}_B(X;Y),\\
&= \epsilon\log(2\epsilon) + (1-\epsilon)\log(2(1-\epsilon)) - \epsilon\log{\epsilon} -(1-\epsilon)\log(1-\epsilon)\\&\quad+\log{\bar{\pi}_B} + \left[2\beta\epsilon-\beta-\epsilon\right]\log{\frac{\bar{\pi}_B}{\pi_B}} ,\\
&=\epsilon\log{2}+(1-\epsilon)\log2 +\log{\bar{\pi}_B} + \left[2\beta\epsilon-\beta-\epsilon\right]\log{\frac{\bar{\pi}_B}{\pi_B}} .
\end{align*}

We then deduce how the difference of mutual information evolves as the decision boundary becomes sharper, i.e. when $\epsilon$ approaches 0:

\begin{equation*}
\lim_{\epsilon\rightarrow 0}\Delta_\mathcal{I} = \log{2} + \log{\bar{\pi}_B} -\beta\log{\frac{\bar{\pi}_B}{\pi_B}} .
\end{equation*}

However, the cluster proportions by B $\pi_B$ also take a different value as $\epsilon$ approaches 0. Recalling Eq.~(\ref{eq:p_b}):

\begin{equation*}
\lim_{\epsilon\rightarrow 0} p_B(y=1) = \beta , \lim_{\epsilon\rightarrow 0} p_B(y=0) = 1-\beta .
\end{equation*}

And finally can we write that:

\begin{equation}\begin{split}
\lim_{\epsilon\rightarrow 0} \Delta_\mathcal{I} &= \log{2} + \log(1-\beta) - \beta \log\frac{1-\beta}{\beta} ,\\
&=\log{2} + (1-\beta)\log(1-\beta) + \beta \log{\beta} ,\\
&=\log{2} -\mathcal{H}(\beta) .
\end{split}\end{equation}

\input{figs_tex/fig_convergence_differences}

To conclude, as the decision boundaries turn sharper, i.e. when $\epsilon$ approaches 0, the difference of mutual information between the two models is controlled by the entropy of proportion of data $\beta$ between the two means $\mu_0$ and $\mu_1$. We know that for binary entropies, the optimum is reached for $\beta=0.5$. In other words having $\mu_0$ and $\mu_1$ distant enough to ensure balance of proportions between the two clusters of model B leads to a difference of mutual information equal to 0. We experimentally highlight this convergence in Figure~\ref{fig:mi_convergence_differences} where we compute the mutual information of models A and B as the distance between the two means $\mu_0$ and $\mu_1$ increases in the Gaussian distribution mixture.

















\section{Introduction}
\label{sec:introduction}
\input{sections/introduction}

\section{Is mutual information a good clustering objective?}
\label{sec:about_mi}
\input{sections/mi_clustering_objective}

\section{Extending the mutual information to the generalised mutual information}
\label{sec:gemini}
\input{sections/mi_to_gemini}

\section{Experiments}
\label{sec:experiments}
\input{sections/experiments}

\section{Conclusion}
\label{sec:conclusions}
\input{sections/conclusion}

\section*{Acknowledgements}
This work has been supported by the French government, through the 3IA C\^ote d'Azur, Investment in the Future, project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. We would also like to thank the France Canada Research Fund (FFCR) for their contribution to the project. This work was partly supported by EU Horizon 2020 project AI4Media, under contract no. 951911. The authors are grateful to the OPAL infrastructure from Université Côte d'Azur for providing resources and support.




\subsection{When the MI fails because of the modelling}
\label{ssec:exp_categorical}
\input{sections/experiments_section/destroying_mi_once_and_for_all}

\subsection{Resistance to outliers}
\label{ssec:exp_gstm}
\input{sections/experiments_section/gaussian_student_mixture}

\subsection{Leveraging a manifold geometry}
\label{ssec:exp_moons}
\input{sections/experiments_section/moons}

\subsection{Fitting MNIST}
\label{ssec:exp_mnist}
\input{sections/experiments_section/mnist}

\subsection{Cifar10 clustering using a SIMCLR-derived kernel}
\label{ssec:exp_cifar10}
\input{sections/experiments_section/simclr_kernel}








\subsection{Representation learning}
\label{ssec:mi_representation}




Representation learning consists in finding high-level features $\vec{z}_i$ extracted from the data $\vec{x}_i$ in order to perform a \emph{downstream task}, e.g. clustering or classification. \gls*{mi} between $\vec{x}$ and $\vec{z}$ is a common choice for learning features\cite{hjelm_learning_2019}. However, estimating correctly \gls*{mi} between two random variables in continuous domains is often intractable when $p(\vec{x}|\vec{z})$ or $p(\vec{z}|\vec{x})$ is unknown, thus lower bounds are preferred, e.g. variational estimators such as MINE \cite{belghazi_mine_2018}, $\mathcal{I}_\text{NCE}$\cite{van_den_oord_representation_2018}. Another common choice of loss function to train features are contrastive losses such as NT-XENT \cite{chen_simple_2020} where the similarity between the features $\vec{z}_i$ from data $\vec{x}_i$ is maximised with the features $\tilde{\vec{z}}$ from a data-augmented $\tilde{\vec{x}}_i$ against any other features $\vec{z}_j$. Recently, \citet{do_clustering_2021} achieved excellent performances in single-stage methods by highlighting the link between the $\mathcal{I}_\text{NCE}$ estimator \cite{van_den_oord_representation_2018} and contrastive learning losses. Representation learning therefore comes at the cost of a complex lower bound estimator on \gls*{mi}, which often requires data augmentation.
Moreover, it was noticed that the \gls*{mi} is hardly predictive of downstream tasks \cite{tschannen_mutual_2019} when the variable $y$ is continuous, i.e. a high value of \gls*{mi} does not clarify whether the discovered representations are insightful with regards to the target of the downstream task.

\subsection{Discriminative clustering}
\label{sssec:mi_clustering}


The \gls*{mi} has been first used as an objective for learning discriminative clustering models \cite{bridle_unsupervised_1992}. Associated architectures went from simple logistic regression \cite{krause_discriminative_2010} to deeper architectures \cite{hu_learning_2017,ji_invariant_2019}. Beyond architecture improvement, the \gls*{mi} maximisation was also carried with several regularisations. These regularisations include penalty terms such as weight decay \cite{krause_discriminative_2010} or Virtual Adversarial Training (VAT, \citealp{hu_learning_2017,miyato_virtual_2018}). Data augmentation was further used to provide invariances in clustering, as well as specific architecture designs like auxiliary clustering heads \cite{ji_invariant_2019}. Rewriting the \gls*{mi} in terms of entropies:
\begin{equation}
\label{eq:mi_kl_entropies}
\mathcal{I} (\vec{x};y) = \mathcal{H}(y) - \mathcal{H}(y|\vec{x})
\end{equation}
highlights a requirement for balanced clusters, through the cluster entropy term $\mathcal{H}(y)$. Indeed, a uniform distribution maximises the entropy. This hints that an unregularised discrete mutual information for clustering can possibly produce uniformly distributed clusters among samples, regardless of how close they could be. We highlight this claim in section~\ref{ssec:local_maxima}. As an example of regularisation impact: maximising the \gls*{mi} with $\ell_2$ constraint can be equivalent to a soft and regularised K-Means in a feature space \cite{jabi_deep_2019}. In clustering, the number of clusters to find is usually not known in advance. Therefore, an interesting clustering algorithm should be able to find a relevant number of clusters, i.e. perform model selection. However, model selection for parametric deep clustering models is expensive \cite{ronen_deepdpm_2022}. Cluster selection through \gls*{mi} maximisation has been little studied in related works, since experiments usually tasked models to find the (supervised) classes of datasets. Furthermore, the literature diverged towards deep learning methods focusing mainly on images, yet rarely on other type of data such as tabular data \cite{min_survey_2018}.

\subsection{Maximising the MI can lead to bad decision boundaries}
\label{ssec:local_maxima}


Maximising the \gls*{mi} directly can be a poor objective: a high \gls*{mi} value is not necessarily predictive of the quality of the features regarding downstream tasks \cite{tschannen_mutual_2019} when $y$ is continuous. We support a similar argument for the case where the data $x$ is a continuous random variable and the cluster assignment $y$ a categorical variable. Indeed, the \gls*{mi} can be maximised by setting appropriately a sharp decision boundary which partitions evenly the data. This reasoning can be seen in the entropy-based formulation of the \gls*{mi} (Eq. \ref{eq:mi_kl_entropies}): any sharp decision boundary minimises the negative conditional entropy, while ensuring balanced clusters maximises the entropy of cluster proportions. Consider for example Figure ~\ref{fig:example_good_odd_mi}, where a mixture of Gaussian distributions with equal variances is separated by a sharp decision boundary. We highlight that both models will have the same mutual information on condition that the misplaced decision boundary of Figure~\ref{sfig:odd_decision_boundary} splits evenly the dataset (see Appendix~\ref{app:mi_convergence}).

\input{figs_tex/fig_example_bad_mi}

Globally, \gls*{mi} misses the idea in clustering that any two points close to one another may be in the same cluster according to some chosen metric. Hence regularisations are required to ensure this constraint. A sketch of these insights was mentionned by \citet{bridle_unsupervised_1992,corduneanu_information_2012}.

\subsection{The discriminative clustering framework for GEMINIs}
\label{ssec:discriminative_clustering}

We change our view on the mutual information by seeing it as a discriminative clustering objective that aims at separating the data distribution given cluster assignments $p(\vec{x}|y)$ from the data distribution $p(\vec{x})$ according to the \gls*{kl} divergence:
\begin{equation}\label{eq:base_mi}
    \mathcal{I}(\vec{x};y) = \mathbb{E}_{y\sim p(y)} \left[ D_\text{KL}(p(\vec{x}|y)\|p(\vec{x}))\right].
\end{equation}
To highlight the discriminative clustering design, we explicitly remove our hypotheses on the data distribution by writing $p_\textup{data}(\vec{x})$.  The only part of the model that we design is a conditional distribution $\p(y|\vec{x})$ that assigns a cluster $y$ to a sample $\vec{x}$ using the parameters $\theta$ \cite{minka_discriminative_2005}. This conditional distribution can typically be a neural network of adequate design regarding the data, e.g. a CNN, or a simple categorical distribution. Consequently, the cluster proportions are controlled by $\theta$ because $\p(y)=\mathbb{E}[\p(y|\vec{x})]$ and so is the conditional distribution $\p(\vec{x}|y)$ even though intractable. This questions how Eq. (\ref{eq:base_mi}) can be computed. Fortunately, well-known properties of \gls*{mi} can invert the distributions on which the KL divergence is computed \cite{bridle_unsupervised_1992,krause_discriminative_2010} via Bayes' theorem
\begin{equation} \label{eq:tractable_discriminative_mi}
    \mathcal{I}(\vec{x};y) = \mathbb{E}_{\vec{x} \sim p_\textup{data}(\vec{x})} \left[ D_\text{KL} (\p(y|\vec{x}) \| \p(y))\right],
\end{equation}
which is possible to estimate.
Since we highlighted earlier that the KL divergence in the \gls*{mi} can lead to inappropriate decision boundaries, we are interested in replacing it by other distances or divergences. However, changing it in Eq. (\ref{eq:tractable_discriminative_mi}) would focus on the separation of cluster assignments from the cluster proportions which may be irrelevant to the data distribution. We rather alter Eq. (\ref{eq:base_mi}) to clearly show that we separate data distributions given clusters from the entire data distribution because it allows us to take into account the data space geometry.

\subsection{The GEMINI}
\label{ssec:gemini}

\input{figs_tex/tab_geminis}

The goal of the \gls*{gemini} is to separate data distributions according to an arbitrary distance $D$, i.e. changing the \gls*{kl} divergence for another divergence or distance in the \gls*{mi}. Moreover, we question the evaluation of the distance between the distribution of the data given a cluster assumption $\p(\vec{x}|y)$ and the entire data distribution $p_\textup{data}(\vec{x})$. We argue that it is intuitive in clustering to compare the distribution of one cluster against the distribution of \emph{another cluster} rather than the data distribution. This raises the definition of two \glspl*{gemini}, one named \gls*{ova}:
\begin{equation}\label{eq:gemini_ova}
    \mathcal{I}^\text{OvA}_D(\vec{x};y) = \mathbb{E}_{y \sim \p(y)} \left[ D(\p(\vec{x}|y)\|p_\textup{data}(\vec{x}))\right],
\end{equation}
which compares the cluster distributions to the data distribution, and the \gls*{ovo} in which we independently draw cluster assignments $y_a$ and $y_b$ (see App.~\ref{app:ovo_and_ova} for an \gls*{ovo} justification):
\begin{equation}\label{eq:gemini_ovo}
    \mathcal{I}^\text{OvO}_D(\vec{x};y) = \mathbb{E}_{y_a,y_b \sim \p(y)} \left[ D(\p(\vec{x}|\ya) \| \p(\vec{x}|\yb))\right],
\end{equation}
There exists other distances than the KL to measure how far two distributions $p$ and $q$ are one from the other. We can make a clear distinction between two types of distances, Csiszar's $f$-divergences \cite{csiszar_information-type_1967} and \gls*{ipm} \cite{sriperumbudur_integral_2009}. However, unlike $f$-divergences, \gls*{ipm}-derived distances like the Wasserstein distance or the \gls*{mmd}\cite{gneiting_strictly_2007,gretton_kernel_2012} bring knowledge about the data throughout either a distance $c$ or a kernel $\kappa$: these distances are geometry-aware.

\paragraph{$f$-divergence GEMINIs:} These divergences involve a convex function $f:\mathbb{R}^+\rightarrow\mathbb{R}$ such that $f(1)=0$. This function is applied to evaluate the ratio between two distributions $p$ and $q$, as in Eq.~(\ref{eq:f_divergences_definition}):
\begin{equation}
\label{eq:f_divergences_definition}
D_\text{f-div}(p,q) = \mathbb{E}_{\vec{z} \sim q(\vec{z})} \left[ f\left(\frac{p(\vec{z})}{q(\vec{z})}\right)\right].
\end{equation}
We will focus on three $f$-divergences: the \gls*{kl} divergence, the \gls*{tv} distance and the squared Hellinger distance. While the \gls*{kl} divergence is the usual divergence for the \gls*{mi}, the \gls*{tv} and the squared Hellinger distance present different advantages among $f$-divergences. First of all, both of them are bounded between 0 and 1. It is consequently easy to check when any \gls*{gemini} using those is maximised contrarily to the \gls*{mi} that is bounded by the minimum of the entropies of $\vec{x}$ and $y$ \cite{gray_maximum_1977}. When used as distance between data conditional distribution $\p(\vec{x}|y)$ and data distribution $p_\textup{data}(\vec{x})$, we can apply Bayes' theorem in order to get an estimable equation to maximise, which only involves cluster assignment $\p(y|\vec{x})$ and marginals $\p(y)$ (see Table~\ref{tab:all_geminis}).

\paragraph{MMD GEMINIs:} The \gls*{mmd} corresponds to the distance between the respective expected embedding of samples from distribution $p$ and distribution $q$ in a reproducing kernel hilbert space (RKHS) $\mathcal{H}$:
\begin{equation}\label{eq:mmd_definition}
    \text{MMD}(p\|q) = \| \mathbb{E}_{\vec{z} \sim p(\vec{z})} [\varphi(\vec{z})] - \mathbb{E}_{\vec{z}\sim q(\vec{z})} [\varphi(\vec{z})]\|_\mathcal{H},
\end{equation}
where $\varphi$ is the RKHS embedding. To compute this distance we can use the kernel trick \cite{gretton_kernel_2012} by involving the kernel function $\kappa(\vec{a},\vec{b})=\langle\varphi(\vec{a}),\varphi(\vec{b})\rangle$. We use Bayes' theorem to uncover a version of the \gls*{mmd} that can be estimated through Monte Carlo using only the predictions $\p(y|\vec{x})$ (see Table~\ref{tab:all_geminis}).
\looseness=-1

\paragraph{Wasserstein GEMINIs:} This distance is an optimal transport distance, defined as:
\begin{equation}\label{eq:wasserstein_definition}
    \mathcal{W}_c(p,q) = \left(\inf_{\gamma \in \Gamma(p,q)} \mathbb{E}_{\vec{x}, \vec{z} \sim \gamma(\vec{x},\vec{z})}\left[c(\vec{x},\vec{z}) \right]\right),
\end{equation}
where $\Gamma(p,q)$ is the set of all couplings between $p$ and $q$ and $c$ a distance function in $\mathcal{X}$. Computing the Wasserstein distance between two distributions $p_\theta(\vec{x}|y=k_1)$ and $p_\theta(\vec{x}|y=k_2)$ is difficult in our discriminative context because we only have access to a finite set of samples $N$. To achieve the Wasserstein-GEMINI, we instead use approximations of the distributions with weighted sums of Diracs:
\looseness=-1
\begin{equation}\label{eq:dirac_approximation}
p_\theta(\vec{x}|y=k) \approx \sum_{i=1}^N m_i^k \delta_{\vec{x}_i} = p_N^k,\quad\text{with}\quad m_i^k = \frac{p_\theta(y=k|\vec{x}_i)}{\sum_{j=1}^Np_\theta(y=k|\vec{x}_j)},
\end{equation}
where $\delta_{\vec{x}_i}$ is a Dirac located on sample location $\vec{x}_i\in\mathcal{X}$. The Wasserstein-\gls{ova} and -\gls{ovo} applied to Dirac sums are compatible with the \verb+emd2+ function of the Python optimal transport package~\cite{flamary_pot_2021} which gracefully supports automatic differentation (see Appendix~\ref{app:wasserstein_convergence} for convergence to the expectation). All \glspl*{gemini} are summarised in Table~\ref{tab:all_geminis}, (see Appendix~\ref{app:deriving_geminis} for derivations).
\looseness=-1




\section*{Checklist}

\begin{enumerate}


\item For all authors...
\begin{enumerate}
  \item Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
    \answerYes{We clearly state that we focus on mutual information for clustering and how we improve it.}
  \item Did you describe the limitations of your work?
    \answerYes{We describe shortcomings in Sec.~\ref{ssec:exp_cifar10},\ref{sec:conclusions} and provide more in App.~\ref{app:exp_complexity}.}
  \item Did you discuss any potential negative societal impacts of your work?
    \answerNA{}
  \item Have you read the ethics review guidelines and ensured that your paper conforms to them?
    \answerYes{Yes, we confirm having read the ethics review guidelines and conform to it.}
\end{enumerate}


\item If you are including theoretical results...
\begin{enumerate}
  \item Did you state the full set of assumptions of all theoretical results?
    \answerYes{We clearly highlighted that we compute the GEMINI between a continuous variable, namely the data, and a discrete variable, namely the cluster assignment. The framework was detailed precisely in Sec.~\ref{ssec:discriminative_clustering} and all distances have been introduced.}
        \item Did you include complete proofs of all theoretical results?
    \answerYes{Please refer to appendices~\ref{app:mi_convergence}, \ref{app:wasserstein_convergence} or \ref{app:deriving_geminis} for proofs of our claims.}
\end{enumerate}


\item If you ran experiments...
\begin{enumerate}
  \item Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)?
    \answerYes{Link to github will be available upon release of paper and a sample is zipped aside in the supplementary material.}
  \item Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)?
    \answerYes{Data generation processes are reported in supplementary materials and main parameters are described at the head of Sec.~\ref{sec:experiments}.}
        \item Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)?
    \answerYes{When applicable due to multiple runs in table form, e.g. Table~\ref{tab:mnist_experiment}.}
        \item Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)?
    \answerYes{We mention it at the beginning of Sec.~\ref{sec:experiments}}
\end{enumerate}


\item If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
\begin{enumerate}
  \item If your work uses existing assets, did you cite the creators?
    \answerYes{All packages for the experiments are mentionned in Appendix~\ref{app:requirements}}
  \item Did you mention the license of the assets?
    \answerYes{We are under licenses of MIT (NumPy, POT,), BSD (torch,pandas, scikit-learn) and the python software foundation (matplotlib, python).}
  \item Did you include any new assets either in the supplemental material or as a URL?
    \answerYes{The code once released will contain the specific file \verb+losses.py+ that will contain only the definitions of GEMINIs, alongside all configurations to reproduce the experiments.}
  \item Did you discuss whether and how consent was obtained from people whose data you're using/curating?
    \answerNA{}
  \item Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content?
    \answerNA{}
\end{enumerate}


\item If you used crowdsourcing or conducted research with human subjects...
\begin{enumerate}
  \item Did you include the full text of instructions given to participants and screenshots, if applicable?
    \answerNA{}
  \item Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable?
    \answerNA{}
  \item Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation?
    \answerNA{}
\end{enumerate}


\end{enumerate}