<figure><img src="image_9.png" /><figcaption>Figure 7: Rényi divergence between $q_{\mu_k, \Sigma_k} \in \mathcal{T}_{\nu}^d$ and $\pi$ in dimension $d = 5$ with $\kappa_{\pi} = 1000$ at every iteration $k$. The iterates $q_{\mu_k, \Sigma_k} \in \mathcal{T}_{\nu}^d$ are obtained using Algorithm 6. The line is the median Rényi divergence per iteration and the shaded area is the interval between the first and third quartiles. The quartiles are obtained by running the algorithm for 100 runs of 1000 iterations.</figcaption></figure>

**Synthesis of the results** Table 1 summarizes the results, for the three algorithms, in terms of final value of the Rényi divergence, averaged over 100 runs, after $10^3$ iterations. The results span the two scenarios, the first with a high-dimensional target with good conditioning, and the second with a poorly-conditioned target in lower dimension. The relative performance of Algorithms 5 and 6 depends on the scenario. On targets that are well-conditioned but high-dimensional, Algorithm 5 seems to behave better than Algorithm 6. On the contrary, Algorithm 6 yields here the best performance when the target is poorly conditioned, showing that our proposed scale adaptation mechanism is able to capture the geometry of the target distribution. We have also implemented the algorithms in the limit $\nu \to +\infty$, corresponding to the exponential family. In this case, the target is not properly captured by the proposals, showing the interest of our new result about the $\lambda$-exponential family. Finally, as expected, the idealized Algorithm 4 reaches the best results in most cases, confirming the validity of our optimality conditions.

<table><caption>Table 1: Median of the Rényi divergence $RD_{\alpha}(\pi, q_{\mu_K, \Sigma_K})$ over 100 runs of $K = 10^3$ iteration. "High d" corresponds to $d = 20$, $\kappa_{\pi} = 10$ and "High $\kappa_{\pi}$" to $d = 5$, $\kappa_{\pi} = 10^3$. The symbol $\times$ denotes situations when Equation (5.1) is not satisfied. For each target and each approximating family $\mathcal{T}_{\nu}^d$, we highlighted in bold font the algorithm achieving the lowest value between Algorithm 5 and 6. The values obtained with the idealized Algorithm 4 are indicated as a reference.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">ν<sub>π</sub> = 1</th><th colspan="2">ν<sub>π</sub> = 3</th><th colspan="2">ν<sub>π</sub> = 10</th></tr><tr><th>High d</th><th>High κ<sub>π</sub></th><th>High d</th><th>High κ<sub>π</sub></th><th>High d</th><th>High κ<sub>π</sub></th></tr></thead><tbody><tr><td rowspan="3">ν = 1</td><td>Alg. 4</td><td>2.25 · 10<sup>-2</sup></td><td>3.84 · 10<sup>-3</sup></td><td>2.61 · 10<sup>-1</sup></td><td>2.15 · 10<sup>-1</sup></td><td>6.81 · 10<sup>-1</sup></td><td>4.65 · 10<sup>-1</sup></td></tr><tr><td>Alg. 5</td><td><b>3.01</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>4.00</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>3.80</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>5.45</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>7.30</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>7.48</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>1.13 · 10<sup>0</sup></td><td><b>6.22</b> · <b>10<sup>-2</sup></b></td><td>1.25 · 10<sup>0</sup></td><td><b>2.29</b> · <b>10<sup>-1</sup></b></td><td>1.64 · 10<sup>0</sup></td><td><b>4.72</b> · <b>10<sup>-1</sup></b></td></tr><tr><td rowspan="3">ν = 3</td><td>Alg. 4</td><td>7.02 · 10<sup>-1</sup></td><td>7.78 · 10<sup>-1</sup></td><td>1.49 · 10<sup>-3</sup></td><td>4.50 · 10<sup>-4</sup></td><td>1.99 · 10<sup>-1</sup></td><td>1.03 · 10<sup>-1</sup></td></tr><tr><td>Alg. 5</td><td><b>1.08</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>1.83</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>1.46</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>4.18</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>2.49</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>4.76</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>2.08 · 10<sup>0</sup></td><td><b>1.02</b> · <b>10<sup>0</sup></b></td><td>1.03 · 10<sup>0</sup></td><td><b>1.78</b> · <b>10<sup>-2</sup></b></td><td>1.15 · 10<sup>0</sup></td><td><b>1.10</b> · <b>10<sup>-1</sup></b></td></tr><tr><td rowspan="3">ν = 10</td><td>Alg. 4</td><td>5.83 · 10<sup>0</sup></td><td>×</td><td>4.69 · 10<sup>-1</sup></td><td>2.41 · 10<sup>-1</sup></td><td>7.08 · 10<sup>-4</sup></td><td>2.24 · 10<sup>-4</sup></td></tr><tr><td>Alg. 5</td><td><b>7.20</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td>×</td><td><b>6.68</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>1.07</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>5.50</b> · <b>10<sup>-2</sup></b></td><td><b>5.25</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>9.36 · 10<sup>0</sup></td><td>×</td><td>1.71 · 10<sup>0</sup></td><td><b>2.74</b> · <b>10<sup>-1</sup></b></td><td>9.67 · 10<sup>-1</sup></td><td><b>8.63</b> · <b>10<sup>-3</sup></b></td></tr><tr><td rowspan="3">ν = +∞</td><td>Alg. 4</td><td>×</td><td>×</td><td>4.91 · 10<sup>1</sup></td><td>1.67 · 10<sup>1</sup></td><td>4.06 · 10<sup>1</sup></td><td>1.48 · 10<sup>1</sup></td></tr><tr><td>Alg. 5</td><td>×</td><td>×</td><td><b>5.01</b> · <b>10<sup>1</sup></b><sup>-1</sup></td><td><b>1.96</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>4.07</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>1.56</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>×</td><td>×</td><td>5.74 · 10<sup>1</sup></td><td><b>1.69</b> · <b>10<sup>1</sup></b></td><td>4.18 · 10<sup>1</sup></td><td><b>1.48</b> · <b>10<sup>1</sup></b></td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 7: Rényi divergence between $q_{\mu_k, \Sigma_k} \in \mathcal{T}_{\nu}^d$ and $\pi$ in dimension $d = 5$ with $\kappa_{\pi} = 1000$ at every iteration $k$. The iterates $q_{\mu_k, \Sigma_k} \in \mathcal{T}_{\nu}^d$ are obtained using Algorithm 6. The line is the median Rényi divergence per iteration and the shaded area is the interval between the first and third quartiles. The quartiles are obtained by running the algorithm for 100 runs of 1000 iterations.</figcaption></figure>

**Synthesis of the results** Table 1 summarizes the results, for the three algorithms, in terms of final value of the Rényi divergence, averaged over 100 runs, after $10^3$ iterations. The results span the two scenarios, the first with a high-dimensional target with good conditioning, and the second with a poorly-conditioned target in lower dimension. The relative performance of Algorithms 5 and 6 depends on the scenario. On targets that are well-conditioned but high-dimensional, Algorithm 5 seems to behave better than Algorithm 6. On the contrary, Algorithm 6 yields here the best performance when the target is poorly conditioned, showing that our proposed scale adaptation mechanism is able to capture the geometry of the target distribution. We have also implemented the algorithms in the limit $\nu \to +\infty$, corresponding to the exponential family. In this case, the target is not properly captured by the proposals, showing the interest of our new result about the $\lambda$-exponential family. Finally, as expected, the idealized Algorithm 4 reaches the best results in most cases, confirming the validity of our optimality conditions.

<table><caption>Table 1: Median of the Rényi divergence $RD_{\alpha}(\pi, q_{\mu_K, \Sigma_K})$ over 100 runs of $K = 10^3$ iteration. "High d" corresponds to $d = 20$, $\kappa_{\pi} = 10$ and "High $\kappa_{\pi}$" to $d = 5$, $\kappa_{\pi} = 10^3$. The symbol $\times$ denotes situations when Equation (5.1) is not satisfied. For each target and each approximating family $\mathcal{T}_{\nu}^d$, we highlighted in bold font the algorithm achieving the lowest value between Algorithm 5 and 6. The values obtained with the idealized Algorithm 4 are indicated as a reference.</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">ν<sub>π</sub> = 1</th><th colspan="2">ν<sub>π</sub> = 3</th><th colspan="2">ν<sub>π</sub> = 10</th></tr><tr><th>High d</th><th>High κ<sub>π</sub></th><th>High d</th><th>High κ<sub>π</sub></th><th>High d</th><th>High κ<sub>π</sub></th></tr></thead><tbody><tr><td rowspan="3">ν = 1</td><td>Alg. 4</td><td>2.25 · 10<sup>-2</sup></td><td>3.84 · 10<sup>-3</sup></td><td>2.61 · 10<sup>-1</sup></td><td>2.15 · 10<sup>-1</sup></td><td>6.81 · 10<sup>-1</sup></td><td>4.65 · 10<sup>-1</sup></td></tr><tr><td>Alg. 5</td><td><b>3.01</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>4.00</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>3.80</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>5.45</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>7.30</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>7.48</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>1.13 · 10<sup>0</sup></td><td><b>6.22</b> · <b>10<sup>-2</sup></b></td><td>1.25 · 10<sup>0</sup></td><td><b>2.29</b> · <b>10<sup>-1</sup></b></td><td>1.64 · 10<sup>0</sup></td><td><b>4.72</b> · <b>10<sup>-1</sup></b></td></tr><tr><td rowspan="3">ν = 3</td><td>Alg. 4</td><td>7.02 · 10<sup>-1</sup></td><td>7.78 · 10<sup>-1</sup></td><td>1.49 · 10<sup>-3</sup></td><td>4.50 · 10<sup>-4</sup></td><td>1.99 · 10<sup>-1</sup></td><td>1.03 · 10<sup>-1</sup></td></tr><tr><td>Alg. 5</td><td><b>-1.08</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-1.83</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-1.46</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-4.18</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-2.49</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-4.76</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>2.08 · 10<sup>0</sup></td><td><b>1.02</b> · <b>10<sup>0</sup></b></td><td>1.03 · 10<sup>0</sup></td><td><b>1.78</b> · <b>10<sup>-2</sup></b></td><td>1.15 · 10<sup>0</sup></td><td><b>1.10</b> · <b>10<sup>-1</sup></b></td></tr><tr><td rowspan="3">ν = 10</td><td>Alg. 4</td><td>5.83 · 10<sup>0</sup></td><td>×</td><td>4.69 · 10<sup>-1</sup></td><td>2.41 · 10<sup>-1</sup></td><td>7.08 · 10<sup>-4</sup></td><td>2.24 · 10<sup>-4</sup></td></tr><tr><td>Alg. 5</td><td><b>-7.20</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td>×</td><td><b>-6.68</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-1.07</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-5.50</b> · <b>10<sup>-2</sup></b><sup>-1</sup></td><td><b>-5.25</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>9.36 · 10<sup>0</sup></td><td>×</td><td>1.71 · 10<sup>0</sup></td><td><b>2.74</b> · <b>10<sup>-1</sup></b></td><td>9.67 · 10<sup>-1</sup></td><td><b>8.63</b> · <b>10<sup>-3</sup></b></td></tr><tr><td rowspan="3">ν = +∞</td><td>Alg. 4</td><td>×</td><td>×</td><td>4.91 · 10<sup>1</sup></td><td>1.67 · 10<sup>1</sup></td><td>4.06 · 10<sup>1</sup></td><td>1.48 · 10<sup>1</sup></td></tr><tr><td>Alg. 5</td><td>×</td><td>×</td><td><b>-5.01</b> · <b>10<sup>1</sup></b><sup>-1</sup></td><td><b>-1.96</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-4.07</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td><td><b>-1.56</b> · <b>10<sup>-1</sup></b><sup>-1</sup></td></tr><tr><td>Alg. 6</td><td>×</td><td>×</td><td>5.74 · 10<sup>1</sup></td><td><b>1.69</b> · <b>10<sup>1</sup></b></td><td>4.18 · 10<sup>1</sup></td><td><b>1.48</b> · <b>10<sup>1</sup></b></td></tr></tbody></table>
We can observe on Table 1 that Algorithms 5 and 6 yield lower performance than Algorithm 4. However, implementing this last algorithm is unrealistic in practice, as it needs samples from the escort of the target. However, we can notice that, when $\nu_{\pi} \neq \nu$, i.e., the approximating family does not match with the target, the algorithms based on MALA are able to reach similar performance than Algorithm 4.

We see in Table 1 that Algorithm 6 outperforms Algorithm 5 on the high $\kappa_{\pi}$ scenario, sometimes by one or two orders of magnitude. This gain can be explained by the fact that Algorithm 6 better handles the shape of the target. This indicates that as soon as the target may be poorly conditioned, it is best to turn to Algorithm 6 instead of Algorithm 5.

On the other hand, the situation is reversed on the high $d$ scenario, where the performance of Algorithm 6 decreases. This indicates that on high-dimensional and well-conditioned targets, it may be beneficial to use Algorithm 5 instead of its scaled version, in Algorithm 6.

Finally, let us mention that when the algorithm matches the scenario, that is Algorithm 5 is used for high $d$ or Algorithm 6 is used for high $\kappa_{\pi}$, it is especially important to choose $\nu = \nu_{\pi}$. Indeed, this is when we observe the biggest degradation if $\nu \neq \nu_{\pi}$.

## 5.2 Maximum likelihood estimation with Student distributions

We consider now maximum likelihood problems of the form ($P_{MLE}$) and ($P_{MLE-Mixt}$) over the $\lambda$-exponential family $Q_{\lambda}$. We will work in the case where $Q_{\lambda}$ is the Student family $T_{\nu}^{d}$.

### 5.2.1 Online maximum likelihood with approximate proximal updates

We now consider a maximum likelihood estimation problem of the form ($P_{MLE}$). The approximating family is hereagain the Student family, $T_{\nu}^{d}$. The samples processed for the maximum likelihood estimation are also distributed following a Student distribution $\pi \in T_{\nu}^{d}$. Following [24], we consider an online setting, where one sample is delivered at each iteration of the algorithm. We implement Algorithm 2 in this setting and study how they approach the true maximum likelihood estimator, depending on the value of $\lambda$.

We assume that at every iteration $k \in \mathbb{N}$ one point $x_k \sim \pi$ is sampled. We implement Algorithm 2 and apply, at each iteration, the operator $P_{\tau_k}^{\{x_k\}}$, with a single data point, namely $x_k$, and we set $\tau_k = \frac{1}{k}$, ensuring an averaging effect. In our setting, this leads to Algorithm 7.

**Algorithm 7:** Online algorithm to solve Problem ($P_{MLE}$) on Student families.

Choose an approximating family $T_{\nu}^{d}$ and initialize $\mu_0$ and $\Sigma_0$.

for $k = 0, \dots, \text{do}$

Using the new sample $x_k$, update $\mu_{k+1}, \Sigma_{k+1}$ following

$$ \left\{ \begin{array}{l} \mu_{k+1} = \frac{1}{k+1}x_k + \frac{k}{k+1}\mu_k, \\ \Sigma_{k+1} = \frac{1}{k+1}x_k x_k^\top + \frac{k}{k+1}(\Sigma_k + \mu_k \mu_k^\top) - \mu_{k+1} \mu_{k+1}^\top. \end{array} \right. \qquad (5.9) $$

As discussed in Section 4.3, Algorithm 7 cannot exactly recover the parameters $(\mu_{\pi}, \Sigma_{\pi})$ of the distribution of the data points, even when $k \to +\infty$. From Propositions 11 and 4, the sequence $\{(\mu_k, \Sigma_k)\}_{k \in \mathbb{N}}$ converges to $(\mu_*, \Sigma_*)$ satisfying $\mu_* = \mu_{\pi}$ and $\Sigma_* = \frac{\nu}{\nu-2}\Sigma_{\pi}$, provided that $\nu > 2$.

We illustrate the behavior of Algorithm 7 by showing several runs of it, in dimension $d = 1$, with $\nu \in \{3, 10\}$. This yields trajectories in the plane $(\mu, \sigma^2)$. In the Gaussian case, recovered when $\nu \to +\infty$, we have from Corollary 2 that $q_* = \pi$. Trying different values of $\nu$ allows to explore situations that are far from the Gaussian setting when $\nu = 3$, or closer to it when $\nu = 10$. In the latter case, we expect a lower mismatch between $\pi$ and $q_*$.
We can observe on Table 1 that Algorithms 5 and 6 yield lower performance than Algorithm 4. However, implementing this last algorithm is unrealistic in practice, as it needs samples from the escort of the target. However, we can notice that, when $\nu_{\pi} \neq \nu$, i.e., the approximating family does not match with the target, the algorithms based on MALA are able to reach similar performance than Algorithm 4.

We see in Table 1 that Algorithm 6 outperforms Algorithm 5 on the high $\kappa_{\pi}$ scenario, sometimes by one or two orders of magnitude. This gain can be explained by the fact that Algorithm 6 better handles the shape of the target. This indicates that as soon as the target may be poorly conditioned, it is best to turn to Algorithm 6 instead of Algorithm 5.

On the other hand, the situation is reversed on the high $d$ scenario, where the performance of Algorithm 6 decreases. This indicates that on high-dimensional and well-conditioned targets, it may be beneficial to use Algorithm 5 instead of its scaled version, in Algorithm 6.

Finally, let us mention that when the algorithm matches the scenario, that is Algorithm 5 is used for high $d$ or Algorithm 6 is used for high $\kappa_{\pi}$, it is especially important to choose $\nu = \nu_{\pi}$. Indeed, this is when we observe the biggest degradation if $\nu \neq \nu_{\pi}$.

## 5.2 Maximum likelihood estimation with Student distributions

We consider now maximum likelihood problems of the form ($P_{MLE}$) and ($P_{MLE-Mixt}$) over the $\lambda$-exponential family $Q_{\lambda}$. We will work in the case where $Q_{\lambda}$ is the Student family $T_{\nu}^{d}$.

### 5.2.1 Online maximum likelihood with approximate proximal updates

We now consider a maximum likelihood estimation problem of the form ($P_{MLE}$). The approximating family is hereagain the Student family, $T_{\nu}^{d}$. The samples processed for the maximum likelihood estimation are also distributed following a Student distribution $\pi \in T_{\nu}^{d}$. Following [24], we consider an online setting, where one sample is delivered at each iteration of the algorithm. We implement Algorithm 2 in this setting and study how they approach the true maximum likelihood estimator, depending on the value of $\lambda$.

We assume that at every iteration $k \in \mathbb{N}$ one point $x_k \sim \pi$ is sampled. We implement Algorithm 2 and apply, at each iteration, the operator $P_{\tau_k}^{\{x_k\}}$, with a single data point, namely $x_k$, and we set $\tau_k = \frac{1}{k}$, ensuring an averaging effect. In our setting, this leads to Algorithm 7.

<caption>**Algorithm 7:** Online algorithm to solve Problem ($P_{MLE}$) on Student families.</caption>

Choose an approximating family $T_{\nu}^{d}$ and initialize $\mu_0$ and $\Sigma_0$.

for $k = 0, \dots, \text{do}$

Using the new sample $x_k$, update $\mu_{k+1}, \Sigma_{k+1}$ following

$$ \left\{ \begin{array}{l} \mu_{k+1} = \frac{1}{k+1}x_k + \frac{k}{k+1}\mu_k, \\ \Sigma_{k+1} = \frac{1}{k+1}x_k x_k^\top + \frac{k}{k+1}(\Sigma_k + \mu_k \mu_k^\top) - \mu_{k+1} \mu_{k+1}^\top. \end{array} \right. \qquad (5.9) $$

As discussed in Section 4.3, Algorithm 7 cannot exactly recover the parameters $(\mu_{\pi}, \Sigma_{\pi})$ of the distribution of the data points, even when $k \to +\infty$. From Propositions 11 and 4, the sequence $\{(\mu_k, \Sigma_k)\}_{k \in \mathbb{N}}$ converges to $(\mu_*, \Sigma_*)$ satisfying $\mu_* = \mu_{\pi}$ and $\Sigma_* = \frac{\nu}{\nu-2}\Sigma_{\pi}$, provided that $\nu > 2$.

We illustrate the behavior of Algorithm 7 by showing several runs of it, in dimension $d = 1$, with $\nu \in \{3, 10\}$. This yields trajectories in the plane $(\mu, \sigma^2)$. In the Gaussian case, recovered when $\nu \to +\infty$, we have from Corollary 2 that $q_* = \pi$. Trying different values of $\nu$ allows to explore situations that are far from the Gaussian setting when $\nu = 3$, or closer to it when $\nu = 10$. In the latter case, we expect a lower mismatch between $\pi$ and $q_*$.
<figure><img src="image_10.png" /><figcaption>(a) Plot of the log-likelihood of one trajectory of Algorithm 7, with the log-likelihood of the true parameters in orange and the bound of Corollary 2 in red.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Plot of 10 trajectories of Algorithm 7, with the point $(\mu_*, \sigma_*^2)$ to which the trajectories converge and the point $(\mu_\pi, \sigma_\pi^2)$ encoding the distribution of the samples.</figcaption></figure>

Figure 8: Plots of trajectories of Algorithm 7, initialized at $\mu_0 = -2$ and $\sigma_0^2 = 10$, in dimension $d = 1$, with samples generated following $\pi \in \mathcal{T}_\nu^d$, $\nu = 3$, with parameters $(\mu_\pi, \sigma_\pi^2)$.

Figure 8 shows that the iterates $\{((\mu_k, \sigma_k^2))\}_{k \in \mathbb{N}}$ generated by Algorithm 7 converge to the point $(\mu_*, \sigma_*^2)$, which is different from the true parameters $(\mu_\pi, \sigma_\pi^2)$. We can also observe in this figure that the log-likelihood of the iterates gets very close to the one of $\pi$. The bound on the sub-optimal log-likelihood, predicted by Proposition 10 and Corollary 2, is satisfied by the iterates $\{((\mu_k, \sigma_k^2))\}_{k \in \mathbb{N}}$ after a small number of iterations.

<figure><img src="image_12.png" /><figcaption>(a) Plot of the log-likelihood of one trajectory of Algorithm 7, with the log-likelihood of the true parameters in orange and the bound of Corollary 2 in red.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Plot of 10 trajectories of Algorithm 7, with the point $(\mu_*, \sigma_*^2)$ to which the trajectories converge and the point $(\mu_\pi, \sigma_\pi^2)$ encoding the distribution of the samples.</figcaption></figure>

Figure 9: Plots of trajectories of Algorithm 7 initialized at $\mu_0 = -2$ and $\sigma_0^2 = 10$, in dimension $d = 1$, with samples generated following $\pi \in \mathcal{T}_\nu^d$, $\nu = 10$, with parameters $(\mu_\pi, \sigma_\pi^2)$.

Figure 9 considers a higher value of $\nu$. This setting is closer to the Gaussian case, reached in the limit $\nu \to +\infty$, for which our algorithm reaches the true distribution of the samples. We thus observe that in Figure 9, the log-likelihood converges to the value of the log-likelihood of $\pi$. This is in contrast with Figure
<figure><img src="image_14.png" /><figcaption>(a) Plot of the log-likelihood of one trajectory of Algorithm 7, with the log-likelihood of the true parameters in orange and the bound of Corollary 2 in red.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Plot of 10 trajectories of Algorithm 7, with the point $(\mu_*, \sigma_*^2)$ to which the trajectories converge and the point $(\mu_\pi, \sigma_\pi^2)$ encoding the distribution of the samples.</figcaption></figure>

Figure 8: Plots of trajectories of Algorithm 7, initialized at $\mu_0 = -2$ and $\sigma_0^2 = 10$, in dimension $d = 1$, with samples generated following $\pi \in \mathcal{T}_\nu^d$, $\nu = 3$, with parameters $(\mu_\pi, \sigma_\pi^2)$.

Figure 8 shows that the iterates $\{((\mu_k, \sigma_k^2))\}_{k \in \mathbb{N}}$ generated by Algorithm 7 converge to the point $(\mu_*, \sigma_*^2)$, which is different from the true parameters $(\mu_\pi, \sigma_\pi^2)$. We can also observe in this figure that the log-likelihood of the iterates gets very close to the one of $\pi$. The bound on the sub-optimal log-likelihood, predicted by Proposition 10 and Corollary 2, is satisfied by the iterates $\{((\mu_k, \sigma_k^2))\}_{k \in \mathbb{N}}$ after a small number of iterations.

<figure><img src="image_16.png" /><figcaption>(a) Plot of the log-likelihood of one trajectory of Algorithm 7, with the log-likelihood of the true parameters in orange and the bound of Corollary 2 in red.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Plot of 10 trajectories of Algorithm 7, with the point $(\mu_*, \sigma_*^2)$ to which the trajectories converge and the point $(\mu_\pi, \sigma_\pi^2)$ encoding the distribution of the samples.</figcaption></figure>

Figure 9: Plots of trajectories of Algorithm 7 initialized at $\mu_0 = -2$ and $\sigma_0^2 = 10$, in dimension $d = 1$, with samples generated following $\pi \in \mathcal{T}_\nu^d$, $\nu = 10$, with parameters $(\mu_\pi, \sigma_\pi^2)$.

Figure 9 considers a higher value of $\nu$. This setting is closer to the Gaussian case, reached in the limit $\nu \to +\infty$, for which our algorithm reaches the true distribution of the samples. We thus observe that in Figure 9, the log-likelihood converges to the value of the log-likelihood of $\pi$. This is in contrast with Figure
8, in which we can observe gap. We again observe that the iterates converge to the point $(\mu_*, \sigma_*^2)$, which is very close this time to the true parameters $(\mu_\pi, \sigma_\pi^2)$. Compared to Figure 8 in the case $\nu = 3$, we see that the bound predicted by Corollary 2 is reached from the first iterates, meaning that it is not a tight bound for the log-likelihood of $q_*$.

According to our theoretical results, Algorithm 7 converges to a sub-optimal solution of Problem (P_MLE). Such solution is very easy to implement and could be used to initialize a more complex but exact maximum likelihood estimation algorithm [21, 4]. Moreover, the obtained sub-optimal solution has links with the probability distribution that generated the data, as discussed in Section 4.3 and thus remains relevant for computing exact maximum likelihood estimators.

### 5.2.2 Maximum likelihood estimation with mixtures using relaxed EM

We consider here a maximum likelihood estimation problem over a mixture of Student distributions, that is, Problem (P_MLE-Mixt) where $Q_\lambda = \mathcal{T}_\nu^d$. The samples are also considered to be distributed from a mixture of Student distributions from $\mathcal{T}_\nu^d$, denoted by $\pi$ such that $\pi = \sum_{j=1}^J \xi_{*,j}q_{\mu_{*,j},\Sigma_{*,j}}$. We implement the relaxed EM method described in Algorithm 3 in this particular case. The resulting scheme is summarized in Algorithm 8. Algorithm 3 only requires to work with a $\lambda$-exponential family satisfying Assumptions 1 and 2, so Algorithm 8 is a particular instance for a specific choice of family. We notice that Student distributions benefit from specific properties that would also allow the design of exact EM algorithms [21], so we aim here at illustrating as a proof of concept the use of our mixture-based algorithm.

**Algorithm 8:** A sub-optimal EM algorithm to solve Problem (P_MLE-Mixt) on Student families.

Choose an approximating family $\mathcal{T}_\nu^d$. Initialize Let $\mu_{0,j} \in \mathbb{R}^d$, $\Sigma_{0,j} \in \mathcal{S}_{++}^d$, and $\xi_{0,j} \ge 0$ for $j = 1, \dots, J$ such that $\sum_{j=1}^J \xi_{0,j} = 1$. 

**for** $k = 0, \dots, \text{do}$

For every $j = 1, \dots, J$, define the function $\gamma_{k,j}$ following Equation (4.20), and update the parameters $\xi_{k+1,j}$ and $\mu_{k+1,j}$, $\Sigma_{k+1,j}$ such that they satisfy

$$ \xi_{j,k+1} = \frac{1}{N} \sum_{i=1}^{N} \gamma_{k,j}(x_i), \qquad (5.10) $$

$$ \left\{ \begin{array}{l} \mu_{k+1,j} = \sum_{i=1}^{N} \frac{\gamma_{k,j}(x_i)}{\sum_{i'=1}^{N} \gamma_{k,j}(x_{i'})} x_i, \\ \Sigma_{k+1,j} = \sum_{i=1}^{N} \frac{\gamma_{k,j}(x_i)}{\sum_{i'=1}^{N} \gamma_{k,j}(x_{i'})} x_i x_i^\top - \mu_{k+1,j} \mu_{k+1,j}^\top. \end{array} \right. \qquad (5.11) $$

We illustrate the behavior of Algorithm 8 in dimension $d = 2$, with $\nu \in \{3, 10\}$. Note that a greater value of $\nu$ corresponds to a value of $\lambda$ closer to 0, in which case the approximate M-steps are closer to being optimal (they are optimal for $\lambda = 0$). We use $N = 200$ samples, from $\pi = \sum_{j=1}^J \xi_{*,j}q_{\mu_{*,j},\Sigma_{*,j}}$ with $J = 4$. We use $\{\xi_{*,j}\}_{j=1}^J = \{0.4, 0.1, 0.2, 0.3\}$, with locations parameters $\mu_{*,1} = (10, 10)^\top$, $\mu_{*,2} = (-10, 10)^\top$, $\mu_{*,3} = -\mu_{*,1}$, and $\mu_{*,4} = -\mu_{*,3}$. The shape matrices $\Sigma_{*,j}$, $j = 1, \dots, J$ are constructed in $\mathcal{S}_{++}^d$ with condition number $\kappa = 10$ following [34]. This is a controlled setting which allows to observe precisely the behavior of Algorithm 8 (i.e., an instance of Algorithm 3). Algorithm 8 is initialized with mixture weights satisfying $\xi_{j,0} = 1/J$ for $j = 1, \dots, J$, initial locations parameters $\mu_{j,0}$ sampled from a normal distribution with zero mean and covariance $10I$ and shape parameters being $\Sigma_{j,0} = 10I$ for $j = 1, \dots, J$.

Figure 10 shows the performance of Algorithm 8 when mixture components are from $\mathcal{T}_\nu^d$, with $\nu = 3$. The resulting mixture is able to identify the different components of the data-generating distribution $\pi$ and to achieve a significant increase in terms of log-likelihood from initialization. In this setting, the suboptimality in solving the M-step of the EM algorithm is more pronounced, as the iterates generated by the algorithm
8, in which we can observe gap. We again observe that the iterates converge to the point $(\mu_*, \sigma_*^2)$, which is very close this time to the true parameters $(\mu_\pi, \sigma_\pi^2)$. Compared to Figure 8 in the case $\nu = 3$, we see that the bound predicted by Corollary 2 is reached from the first iterates, meaning that it is not a tight bound for the log-likelihood of $q_*$.

According to our theoretical results, Algorithm 7 converges to a sub-optimal solution of Problem (P_MLE). Such solution is very easy to implement and could be used to initialize a more complex but exact maximum likelihood estimation algorithm [21, 4]. Moreover, the obtained sub-optimal solution has links with the probability distribution that generated the data, as discussed in Section 4.3 and thus remains relevant for computing exact maximum likelihood estimators.

### 5.2.2 Maximum likelihood estimation with mixtures using relaxed EM

We consider here a maximum likelihood estimation problem over a mixture of Student distributions, that is, Problem (P_MLE-Mixt) where $Q_\lambda = \mathcal{T}_\nu^d$. The samples are also considered to be distributed from a mixture of Student distributions from $\mathcal{T}_\nu^d$, denoted by $\pi$ such that $\pi = \sum_{j=1}^J \xi_{*,j} q_{\mu_{*,j}, \Sigma_{*,j}}$. We implement the relaxed EM method described in Algorithm 3 in this particular case. The resulting scheme is summarized in Algorithm 8. Algorithm 3 only requires to work with a $\lambda$-exponential family satisfying Assumptions 1 and 2, so Algorithm 8 is a particular instance for a specific choice of family. We notice that Student distributions benefit from specific properties that would also allow the design of exact EM algorithms [21], so we aim here at illustrating as a proof of concept the use of our mixture-based algorithm.

**Algorithm 8:** A sub-optimal EM algorithm to solve Problem (P_MLE-Mixt) on Student families.

Choose an approximating family $\mathcal{T}_\nu^d$. Initialize Let $\mu_{0,j} \in \mathbb{R}^d$, $\Sigma_{0,j} \in \mathcal{S}_{++}^d$, and $\xi_{0,j} \ge 0$ for $j = 1, \dots, J$ such that $\sum_{j=1}^J \xi_{0,j} = 1$. 

**for** $k = 0, \dots, \text{do}$

For every $j = 1, \dots, J$, define the function $\gamma_{k,j}$ following Equation (4.20), and update the parameters $\xi_{k+1,j}$ and $\mu_{k+1,j}$, $\Sigma_{k+1,j}$ such that they satisfy

$$ \xi_{j,k+1} = \frac{1}{N} \sum_{i=1}^{N} \gamma_{k,j}(x_i), \qquad (5.10) $$

$$ \left\{ \begin{aligned} \mu_{k+1,j} &= \sum_{i=1}^{N} \frac{\gamma_{k,j}(x_i)}{\sum_{i'=1}^{N} \gamma_{k,j}(x_{i'})} x_i, \\ \Sigma_{k+1,j} &= \sum_{i=1}^{N} \frac{\gamma_{k,j}(x_i)}{\sum_{i'=1}^{N} \gamma_{k,j}(x_{i'})} x_i x_i^\top - \mu_{k+1,j} \mu_{k+1,j}^\top. \end{aligned} \right. \qquad (5.11) $$

We illustrate the behavior of Algorithm 8 in dimension $d = 2$, with $\nu \in \{3, 10\}$. Note that a greater value of $\nu$ corresponds to a value of $\lambda$ closer to 0, in which case the approximate M-steps are closer to being optimal (they are optimal for $\lambda = 0$). We use $N = 200$ samples, from $\pi = \sum_{j=1}^J \xi_{*,j} q_{\mu_{*,j}, \Sigma_{*,j}}$ with $J = 4$. We use $\{\xi_{*,j}\}_{j=1}^J = \{0.4, 0.1, 0.2, 0.3\}$, with locations parameters $\mu_{*,1} = (10, 10)^\top$, $\mu_{*,2} = (-10, 10)^\top$, $\mu_{*,3} = -\mu_{*,1}$, and $\mu_{*,4} = -\mu_{*,3}$. The shape matrices $\Sigma_{*,j}$, $j = 1, \dots, J$ are constructed in $\mathcal{S}_{++}^d$ with condition number $\kappa = 10$ following [34]. This is a controlled setting which allows to observe precisely the behavior of Algorithm 8 (i.e., an instance of Algorithm 3). Algorithm 8 is initialized with mixture weights satisfying $\xi_{j,0} = 1/J$ for $j = 1, \dots, J$, initial locations parameters $\mu_{j,0}$ sampled from a normal distribution with zero mean and covariance $10I$ and shape parameters being $\Sigma_{j,0} = 10I$ for $j = 1, \dots, J$.

Figure 10 shows the performance of Algorithm 8 when mixture components are from $\mathcal{T}_\nu^d$, with $\nu = 3$. The resulting mixture is able to identify the different components of the data-generating distribution $\pi$ and to achieve a significant increase in terms of log-likelihood from initialization. In this setting, the suboptimality in solving the M-step of the EM algorithm is more pronounced, as the iterates generated by the algorithm
cannot reach the log-likelihood achieved by the data-generating distribution. This is to be expected, as the corresponding value of $\lambda$ is far from $\lambda = 0$ where the M-step is optimal.

<figure><img src="image_14.png" /><figcaption>Figure 10: Plot of the performance achieved by Algorithm 8 after $K = 100$ iterations in terms of log-likelihood and graphical representation in the sample space, for mixtures of distributions in $T_{\nu}^d$, with $\nu = 3$ and $d = 2$.</figcaption></figure>

Figure 11 shows the performance of Algorithm 8 when mixture components are from $T_{\nu}^d$, with $\nu = 10$. We observe that the proposed algorithm generates iterates whose log-likelihood matches the one of the data-generating distribution. Indeed, this setting is closer to the case $\lambda = 0$ where our algorithm solves the M-step in the EM algorithm exactly, showing that the sub-optimality has no severe effect in this case.
cannot reach the log-likelihood achieved by the data-generating distribution. This is to be expected, as the
corresponding value of $\lambda$ is far from $\lambda = 0$ where the M-step is optimal.

<figure><img src="image_18.png" /><figcaption>Figure 10: Plot of the performance achieved by Algorithm 8 after $K = 100$ iterations in terms of log-likelihood and graphical representation in the sample space, for mixtures of distributions in $T_\nu^d$, with $\nu = 3$ and $d = 2$.</figcaption></figure>

Figure 11 shows the performance of Algorithm 8 when mixture components are from $T_\nu^d$, with $\nu = 10$. We observe that the proposed algorithm generates iterates whose log-likelihood matches the one of the data-generating distribution. Indeed, this setting is closer to the case $\lambda = 0$ where our algorithm solves the M-step in the EM algorithm exactly, showing that the sub-optimality has no severe effect in this case.
<figure><img src="image_15.png" /><figcaption>Figure 11: Plot of the performance achieved by Algorithm 8 after $K = 100$ iterations in terms of log-likelihood and graphical representation in the sample space, for mixtures of distributions in $T_\nu^d$, with $\nu = 10$ and $d = 2$.</figcaption></figure>

# 6 Conclusion

In this work, we have studied variational inference and maximum likelihood estimation problems over the $\lambda$-exponential family, and we have proposed algorithms to solve these problems. Several known results on the standard exponential family are retrieved as special cases.

First, we have shown that variational inference problems over the $\lambda$-exponential family can be solved by satisfying a generalized moment-matching condition that extends the existing one for the standard exponential family. We have also proposed an iterative algorithm to solve this problem, which identifies with a Bregman proximal algorithm in the particular case of the exponential family. The usefulness of our optimality conditions and our algorithm is confirmed by numerical experiments on heavy-tailed targets. These experiments show that the $\lambda$-exponential family can be used to capture phenomenon that the standard exponential family fails to represent.

Second, in maximum likelihood estimation problems, we exhibited sub-optimal solutions with a novel algorithm converging to it. In the case of the exponential family, the solutions become optimal and the algorithm reads again as a Bregman proximal algorithm. In the general case, our algorithm is quick and easy to implement, as demonstrated through numerical experiments. For problems with mixtures, we also proposed a relaxed EM algorithm that recovers the standard EM algorithm in the case of the exponential family. An interesting line of research would be the combination of our algorithms, which leads to sub-optimal solutions, with exact methods.

We achieved these results by extending convex analysis notions to a more general setting, replacing the scalar product by a well-chosen non-linear coupling. By leveraging the specific structure of the problems we consider, we have been able to exhibit optimality conditions and proximal-like algorithm using such tools, which is one of the main novelties of our work. Extending our results and techniques to more general problems, including other divergences and distances over probabilities, or more general couplings, related for instance with elliptical distributions, seems to be an exciting area of research.
<figure><img src="image_19.png" /><figcaption>(a) Plot of the log-likelihood achieved by the iterates of Algorithm 8, with the log-likelihood of data-generating distribution in orange.</figcaption></figure>

<caption>(b) Plot of the samples generated by $\pi$, the different colors denoting the component of the mixture from which they have been drawn. The level lines of the final distribution generated by Algorithm 8 are shown in grey.</caption>

Figure 11: Plot of the performance achieved by Algorithm 8 after $K = 100$ iterations in terms of log-likelihood and graphical representation in the sample space, for mixtures of distributions in $T_\nu^d$, with $\nu = 10$ and $d = 2$.

# 6 Conclusion

In this work, we have studied variational inference and maximum likelihood estimation problems over the $\lambda$-exponential family, and we have proposed algorithms to solve these problems. Several known results on the standard exponential family are retrieved as special cases.

First, we have shown that variational inference problems over the $\lambda$-exponential family can be solved by satisfying a generalized moment-matching condition that extends the existing one for the standard exponential family. We have also proposed an iterative algorithm to solve this problem, which identifies with a Bregman proximal algorithm in the particular case of the exponential family. The usefulness of our optimality conditions and our algorithm is confirmed by numerical experiments on heavy-tailed targets. These experiments show that the $\lambda$-exponential family can be used to capture phenomenon that the standard exponential family fails to represent.

Second, in maximum likelihood estimation problems, we exhibited sub-optimal solutions with a novel algorithm converging to it. In the case of the exponential family, the solutions become optimal and the algorithm reads again as a Bregman proximal algorithm. In the general case, our algorithm is quick and easy to implement, as demonstrated through numerical experiments. For problems with mixtures, we also proposed a relaxed EM algorithm that recovers the standard EM algorithm in the case of the exponential family. An interesting line of research would be the combination of our algorithms, which leads to sub-optimal solutions, with exact methods.

We achieved these results by extending convex analysis notions to a more general setting, replacing the scalar product by a well-chosen non-linear coupling. By leveraging the specific structure of the problems we consider, we have been able to exhibit optimality conditions and proximal-like algorithm using such tools, which is one of the main novelties of our work. Extending our results and techniques to more general problems, including other divergences and distances over probabilities, or more general couplings, related for instance with elliptical distributions, seems to be an exciting area of research.
A Proof of Proposition 3

Proof of Proposition 3. (i) Consider a distribution in $\mathcal{T}_\nu^d$ with location parameter $\mu$ and scale matrix $\Sigma$. Then we compute for any $x \in \mathbb{R}^d$ the following.

$$ 
\begin{aligned}
q_{\mu, \Sigma}(x) &\propto \left(1 + \frac{1}{\nu}(x - \mu)^{\top}\Sigma^{-1}(x - \mu)\right)^{-\frac{\nu+d}{2}} \\
&\propto \left(1 + \frac{1}{\nu}\mu^{\top}\Sigma^{-1}\mu - \frac{2}{\nu}\mu^{\top}x + \frac{1}{\nu}x^{\top}\Sigma^{-1}x\right)^{-\frac{\nu+d}{2}} \\
&\propto \left(1 + \frac{1}{\nu}\mu^{\top}\Sigma^{-1}\mu\right)^{-\frac{\nu+d}{2}} \left(1 + \left(-\frac{2}{\nu+d}\right) \left(-\frac{\nu+d}{2\nu(1 + \frac{1}{\nu}\mu^{\top}\Sigma^{-1}\mu)}(-2\mu^{\top}x + x^{\top}\Sigma^{-1}x)\right)\right)^{-\frac{\nu+d}{2}}
\end{aligned}
$$

and since $\mu^\top x = \langle \mu, x \rangle$ and $x^\top \Sigma^{-1} x = \langle \Sigma^{-1}, xx^\top \rangle$, we can identify that $q_{\mu, \Sigma} = q_\vartheta$.

We can identify from the above that $\mathcal{T}_\nu^d$ is an instance of the $\lambda$-exponential family with $\lambda = -\frac{2}{\nu+d}$. Its parameters are

$$ \vartheta_1 = \frac{\nu + d}{\nu + \mu^{\top}\Sigma^{-1}\mu}\Sigma^{-1}\mu, \quad \vartheta_2 = -\frac{\nu + d}{2(\nu + \mu^{\top}\Sigma^{-1}\mu)}\Sigma^{-1}. \quad (A.1) $$

In order to compute $\varphi_\lambda$, let us inverse the mapping $\mu, \Sigma \mapsto \vartheta_1, \vartheta_2$. First, we can easily compute that $\mu = -\frac{1}{2}\vartheta_2^{-1}\vartheta_1$. Now, we compute the intermediate quantity $\mu^\top \Sigma^{-1} \mu$. Remark that

$$ \vartheta_1^\top \vartheta_2^{-1} \vartheta_1 = - \frac{2(\nu + d)}{\nu + \mu^\top \Sigma^{-1} \mu} \mu^\top \Sigma^{-1} \mu. \quad (A.2) $$

Hence we deduce that

$$ \nu + \mu^{\top}\Sigma^{-1}\mu = \frac{2\nu(\nu + d)}{2(\nu + d) + \vartheta_1^{\top}\vartheta_2^{-1}\vartheta_1}. \quad (A.3) $$

From Equations (A.1) and (A.3), it comes that $\Sigma^{-1} = -\frac{4\nu}{2(\nu+d)+\vartheta_1^{\top}\vartheta_2^{-1}\vartheta_1}\vartheta_2$. Summarizing our results, we thus obtained

$$ \mu = -\frac{1}{2}\vartheta_2^{-1}\vartheta_1, \quad \Sigma = -\frac{2(\nu + d) + \vartheta_1^{\top}\vartheta_2^{-1}\vartheta_1}{4\nu}\vartheta_2^{-1}. \quad (A.4) $$

Finally, we turn to the computation of $\varphi_\lambda$. We can identify

$$ 
\begin{aligned}
\varphi_\lambda(\vartheta) &= -\log\left(\frac{\det(\Sigma)^{-\frac{1}{2}}}{Z_\nu}\left(1 + \frac{1}{\nu}\mu\Sigma^{-1}\mu^\top\right)^{-\frac{\nu+d}{2}}\right) \\
&= \frac{1}{2}\log\det(\Sigma) + \frac{\nu+d}{2}\log\left(1 + \frac{1}{\nu}\mu\Sigma^{-1}\mu^\top\right) + \log(Z_\nu) \\
&= \frac{d}{2}\log\left(\frac{2(\nu+d) + \vartheta_1^\top\vartheta_2\vartheta_1}{4\nu}\right) + \frac{1}{2}\log\det(-\vartheta_2^{-1}) + \frac{\nu+d}{2}\log\left(\frac{2(\nu+d)}{2(\nu+d) + \vartheta_1^\top\vartheta_2\vartheta_1}\right) + \log Z_\nu \\
&= -\frac{d}{2}\log(4\nu) + \frac{1}{2}\log\det(-\vartheta_2^{-1}) + \frac{\nu+d}{2}\log(2(\nu+d)) - \frac{\nu}{2}\log(2(\nu+d) + \vartheta_1^\top\vartheta_2^{-1}\vartheta_1) + \log Z_\nu.
\end{aligned}
$$

This shows in particular that $\det \varphi_\lambda(\vartheta) = \{\vartheta \in \mathbb{R}^d \times S_{-\cdot}^d, 2(\nu+d) + \vartheta_1^\top \vartheta_2^{-1} \vartheta_1 > 0\}$, which is non-empty. This shows that $\mathcal{T}_\nu^d$ satisfies Assumption 1.
# A Proof of Proposition 3

*Proof of Proposition 3.* (i) Consider a distribution in $\mathcal{T}_\nu^d$ with location parameter $\mu$ and scale matrix $\Sigma$. Then we compute for any $x \in \mathbb{R}^d$ the following.

$$ 
\begin{aligned}
q_{\mu, \Sigma}(x) &\propto \left(1 + \frac{1}{\nu}(x - \mu)^{\top}\Sigma^{-1}(x - \mu)\right)^{-\frac{\nu+d}{2}} \\
&\propto \left(1 + \frac{1}{\nu}\mu^{\top}\Sigma^{-1}\mu - \frac{2}{\nu}\mu^{\top}x + \frac{1}{\nu}x^{\top}\Sigma^{-1}x\right)^{-\frac{\nu+d}{2}} \\
&\propto \left(1 + \frac{1}{\nu}\mu^{\top}\Sigma^{-1}\mu\right)^{-\frac{\nu+d}{2}} \left(1 + \left(-\frac{2}{\nu+d}\right) \left(-\frac{\nu+d}{2\nu(1 + \frac{1}{\nu}\mu^{\top}\Sigma^{-1}\mu)}(-2\mu^{\top}x + x^{\top}\Sigma^{-1}x)\right)\right)^{-\frac{\nu+d}{2}}
\end{aligned}
$$

and since $\mu^\top x = \langle \mu, x \rangle$ and $x^\top \Sigma^{-1} x = \langle \Sigma^{-1}, xx^\top \rangle$, we can identify that $q_{\mu, \Sigma} = q_\vartheta$.

We can identify from the above that $\mathcal{T}_\nu^d$ is an instance of the $\lambda$-exponential family with $\lambda = -\frac{2}{\nu+d}$. Its parameters are

$$ \vartheta_1 = \frac{\nu + d}{\nu + \mu^{\top}\Sigma^{-1}\mu}\Sigma^{-1}\mu, \quad \vartheta_2 = -\frac{\nu + d}{2(\nu + \mu^{\top}\Sigma^{-1}\mu)}\Sigma^{-1}. \quad (A.1) $$

In order to compute $\varphi_\lambda$, let us inverse the mapping $\mu, \Sigma \mapsto \vartheta_1, \vartheta_2$. First, we can easily compute that $\mu = -\frac{1}{2}\vartheta_2^{-1}\vartheta_1$. Now, we compute the intermediate quantity $\mu^\top \Sigma^{-1} \mu$. Remark that

$$ \vartheta_1^\top \vartheta_2^{-1} \vartheta_1 = - \frac{2(\nu + d)}{\nu + \mu^\top \Sigma^{-1} \mu} \mu^\top \Sigma^{-1} \mu. \quad (A.2) $$

Hence we deduce that

$$ \nu + \mu^{\top}\Sigma^{-1}\mu = \frac{2\nu(\nu + d)}{2(\nu + d) + \vartheta_1^{\top}\vartheta_2^{-1}\vartheta_1}. \quad (A.3) $$

From Equations (A.1) and (A.3), it comes that $\Sigma^{-1} = -\frac{4\nu}{2(\nu+d)+\vartheta_1^{\top}\vartheta_2^{-1}\vartheta_1}\vartheta_2$. Summarizing our results, we thus obtained

$$ \mu = -\frac{1}{2}\vartheta_2^{-1}\vartheta_1, \quad \Sigma = -\frac{2(\nu + d) + \vartheta_1^{\top}\vartheta_2^{-1}\vartheta_1}{4\nu}\vartheta_2^{-1}. \quad (A.4) $$

Finally, we turn to the computation of $\varphi_\lambda$. We can identify

$$ 
\begin{aligned}
\varphi_\lambda(\vartheta) &= -\log\left(\frac{\det(\Sigma)^{-\frac{1}{2}}}{Z_\nu}\left(1 + \frac{1}{\nu}\mu\Sigma^{-1}\mu^\top\right)^{-\frac{\nu+d}{2}}\right) \\
&= \frac{1}{2}\log\det(\Sigma) + \frac{\nu+d}{2}\log\left(1 + \frac{1}{\nu}\mu\Sigma^{-1}\mu^\top\right) + \log(Z_\nu) \\
&= \frac{d}{2}\log\left(\frac{2(\nu+d) + \vartheta_1^\top\vartheta_2\vartheta_1}{4\nu}\right) + \frac{1}{2}\log\det(-\vartheta_2^{-1}) + \frac{\nu+d}{2}\log\left(\frac{2(\nu+d)}{2(\nu+d) + \vartheta_1^\top\vartheta_2\vartheta_1}\right) + \log Z_\nu \\
&= -\frac{d}{2}\log(4\nu) + \frac{1}{2}\log\det(-\vartheta_2^{-1}) + \frac{\nu+d}{2}\log(2(\nu+d)) - \frac{\nu}{2}\log(2(\nu+d) + \vartheta_1^\top\vartheta_2^{-1}\vartheta_1) + \log Z_\nu.
\end{aligned}
$$

This shows in particular that $\text{dom } \varphi_\lambda(\vartheta) = \{\vartheta \in \mathbb{R}^d \times S_{-\cdot}^d, 2(\nu+d) + \vartheta_1^\top\vartheta_2^{-1}\vartheta_1 > 0\}$, which is non-empty. This shows that $\mathcal{T}_\nu^d$ satisfies Assumption 1.
(ii) We now turn to the study of the escort probabilities. We can compute for $x \in \mathbb{R}^d$ the following:

$$ 
\begin{aligned}
q_{\mu, \Sigma}^{(\alpha)}(x) &= \frac{1}{Z^{(\alpha)}} \left( 1 + \frac{1}{\nu}(x - \mu)\Sigma^{-1}(x - \mu)^{\top} \right)^{-\alpha \frac{\nu+d}{2}} \\
&= \frac{1}{Z^{(\alpha)}} \left( 1 + \frac{1}{\nu+2}(x - \mu) \left( \frac{\nu}{\nu+2}\Sigma \right)^{-1} (x - \mu)^{\top} \right)^{-\frac{(\nu+2)+d}{2}}.
\end{aligned}
$$

We recognize that $q_{\mu, \Sigma}^{(\alpha)}$ is a Student distribution with $\nu + 2 > 2$ degrees of freedom, location parameter $\mu$ and scale matrix $\frac{\nu}{\nu+2}\Sigma$. Hence, we obtain that

$$ 
\left\{ 
\begin{aligned}
q_{\mu, \Sigma}^{(\alpha)}(x) &= \mu, \\
q_{\mu, \Sigma}^{(\alpha)}((x - \mu)(x - \mu)^{\top}) &= \frac{\nu+2}{(\nu+2)-2} \left( \frac{\nu}{\nu+2}\Sigma \right) = \Sigma.
\end{aligned}
\right. 
\qquad (A.5) 
$$

To show the bijection result, we show that the map $(\mu, \Sigma) \mapsto (\vartheta_1, \vartheta_2)$ is a bijection between $\mathbb{R}^d \times S_{++}^d$ and dom $\varphi_\lambda$. Consider $\mu \in \mathbb{R}^d$, $\Sigma \in S_{++}^d$ and $\vartheta_1, \vartheta_2$ defined as in Equation (A.1). We can first remark that $\vartheta_1 \in \mathbb{R}^d$ and that $\vartheta_2 \in S_{--}^d$. Using the result of Equation (A.2), we now compute

$$ 
\begin{aligned}
2(\nu + d) + \vartheta_1^{\top} \vartheta_2^{-1} \vartheta_1 &= 2(\nu + d) - \frac{2(\nu + d)}{\nu + \mu^{\top} \Sigma^{-1} \mu} \mu^{\top} \Sigma^{-1} \mu \\
&= \frac{2\nu(\nu + d)}{\nu + \mu^{\top} \Sigma^{-1} \mu} > 0,
\end{aligned}
$$

showing that $\vartheta_1, \vartheta_2 \in \text{dom } \varphi_\lambda$. Consider now $\vartheta_1, \vartheta_2 \in \text{dom } \varphi_\lambda$, and $\mu, \Sigma$ as given by Equation (A.4). By definition of $\text{dom } \varphi_\lambda$, $\mu \in \mathbb{R}^d$ and $\Sigma \in S_{++}^d$, showing the result.

We now compute the Rényi entropy of $q_{\mu, \Sigma} \in \mathcal{T}_\nu^d$ for $\alpha = 1 - \lambda$ with $\lambda = -\frac{2}{\nu+d}$. By using similar steps as above, we obtain

$$ H_{\alpha}(q_{\mu, \Sigma}) = \frac{1}{1-\alpha} \log \left( \frac{1}{Z_{\nu}^{\alpha}(\det \Sigma)^{\frac{\alpha}{2}}} Z_{\nu+2} \det \left( \frac{\nu}{\nu+2} \Sigma \right)^{\frac{1}{2}} \right). \qquad (A.6) $$

(iii) Consider $\vartheta \in \text{dom } \varphi_\lambda$, and $p \in \mathcal{P}(\mathcal{X}, dx)$. Consider $\mu, \Sigma \in \mathbb{R}^d \times S_{++}^d$ given by Equation (A.4). We can then compute

$$ 1 + \lambda \langle \vartheta, p^{(\alpha)}(T) \rangle = 1 - \frac{2}{\nu + \mu^{\top} \Sigma^{-1} \mu} p^{(\alpha)}(x)^{\top} \Sigma^{-1} \mu + \frac{1}{\nu + \mu^{\top} \Sigma^{-1} \mu} \text{tr}(\Sigma^{-1} p^{(\alpha)}(x x^{\top})), \qquad (A.7) $$

which is defined if $p^{(\alpha)}$ has finite first and second order moments. Introducing the quantity $V := p^{(\alpha)}((x - p^{(\alpha)}(x))(x - p^{(\alpha)}(x))^{\top}) = p^{(\alpha)}(xx^{\top}) - p^{(\alpha)}(x)p^{(\alpha)}(x)^{\top} \in S_{++}^d$, we get for any $\vartheta \in \text{dom } \varphi_\lambda$ that

$$ 
\begin{aligned}
1 + \lambda \langle \vartheta, p^{(\alpha)}(T) \rangle &= 1 - \frac{2}{\nu + \mu^{\top} \Sigma^{-1} \mu} p^{(\alpha)}(x)^{\top} \Sigma^{-1} \mu + \frac{1}{\nu + \mu^{\top} \Sigma^{-1} \mu} \text{tr}(\Sigma^{-1} V) \\
&\quad + \frac{1}{\nu + \mu^{\top} \Sigma^{-1} \mu} p^{(\alpha)}(x)^{\top} \Sigma^{-1} p^{(\alpha)}(x) \\
&= \frac{1}{\nu + \mu \Sigma^{-1} \mu} \left( \nu + (\mu - p^{(\alpha)}(x))^{\top} \Sigma^{-1} (\mu - p^{(\alpha)}(x))^{\top} + \text{tr}(\Sigma^{-1} V) \right) \\
&\geq \frac{\nu}{\nu + \mu \Sigma^{-1} \mu}.
\end{aligned}
$$

This shows that for any $\vartheta \in \text{dom } \varphi_\lambda$, and $p^{(\alpha)} \in \mathcal{P}(\mathcal{X}, dx)$ with finite first and second order moments, the quantity $c_\lambda(\vartheta, p^{(\alpha)}(T))$ is in $\mathbb{R}$. With the result of (ii), this shows that $\mathcal{T}_\nu^d$, seen as an instance of the $\lambda$-exponential family, satisfies Assumptions 2. $\square$
(ii) We now turn to the study of the escort probabilities. We can compute for $x \in \mathbb{R}^d$ the following:

$$ 
\begin{aligned}
q_{\mu, \Sigma}^{(\alpha)}(x) &= \frac{1}{Z^{(\alpha)}} \left( 1 + \frac{1}{\nu}(x - \mu)\Sigma^{-1}(x - \mu)^{\top} \right)^{-\alpha \frac{\nu+d}{2}} \\
&= \frac{1}{Z^{(\alpha)}} \left( 1 + \frac{1}{\nu+2}(x - \mu) \left( \frac{\nu}{\nu+2}\Sigma \right)^{-1} (x - \mu)^{\top} \right)^{-\frac{(\nu+2)+d}{2}}.
\end{aligned}
$$

We recognize that $q_{\mu, \Sigma}^{(\alpha)}$ is a Student distribution with $\nu + 2 > 2$ degrees of freedom, location parameter $\mu$ and scale matrix $\frac{\nu}{\nu+2}\Sigma$. Hence, we obtain that

$$ 
\left\{ 
\begin{aligned}
q_{\mu, \Sigma}^{(\alpha)}(x) &= \mu, \\
q_{\mu, \Sigma}^{(\alpha)}((x-\mu)(x-\mu)^{\top}) &= \frac{\nu+2}{(\nu+2)-2} \left( \frac{\nu}{\nu+2}\Sigma \right) = \Sigma.
\end{aligned}
\right. 
\qquad (A.5) 
$$

To show the bijection result, we show that the map $(\mu, \Sigma) \mapsto (\vartheta_1, \vartheta_2)$ is a bijection between $\mathbb{R}^d \times S_{++}^d$ and dom $\varphi_\lambda$. Consider $\mu \in \mathbb{R}^d$, $\Sigma \in S_{++}^d$ and $\vartheta_1, \vartheta_2$ defined as in Equation (A.1). We can first remark that $\vartheta_1 \in \mathbb{R}^d$ and that $\vartheta_2 \in S_{--}^d$. Using the result of Equation (A.2), we now compute

$$ 
\begin{aligned}
2(\nu + d) + \vartheta_1^{\top} \vartheta_2^{-1} \vartheta_1 &= 2(\nu + d) - \frac{2(\nu + d)}{\nu + \mu^{\top} \Sigma^{-1} \mu} \mu^{\top} \Sigma^{-1} \mu \\
&= \frac{2\nu(\nu + d)}{\nu + \mu^{\top} \Sigma^{-1} \mu} > 0,
\end{aligned}
$$

showing that $\vartheta_1, \vartheta_2 \in \text{dom } \varphi_\lambda$. Consider now $\vartheta_1, \vartheta_2 \in \text{dom } \varphi_\lambda$, and $\mu, \Sigma$ as given by Equation (A.4). By definition of $\text{dom } \varphi_\lambda$, $\mu \in \mathbb{R}^d$ and $\Sigma \in S_{++}^d$, showing the result.

We now compute the Rényi entropy of $q_{\mu, \Sigma} \in \mathcal{T}_{\nu}^d$ for $\alpha = 1 - \lambda$ with $\lambda = -\frac{2}{\nu+d}$. By using similar steps as above, we obtain

$$ H_{\alpha}(q_{\mu, \Sigma}) = \frac{1}{1-\alpha} \log \left( \frac{1}{Z_{\nu}^{\alpha} (\det \Sigma)^{\frac{\alpha}{2}}} Z_{\nu+2} \det \left( \frac{\nu}{\nu+2} \Sigma \right)^{\frac{1}{2}} \right). \qquad (A.6) $$

(iii) Consider $\vartheta \in \text{dom } \varphi_\lambda$, and $p \in \mathcal{P}(\mathcal{X}, dx)$. Consider $\mu, \Sigma \in \mathbb{R}^d \times S_{++}^d$ given by Equation (A.4). We can then compute

$$ 1 + \lambda \langle \vartheta, p^{(\alpha)}(T) \rangle = 1 - \frac{2}{\nu + \mu^{\top} \Sigma^{-1} \mu} p^{(\alpha)}(x)^{\top} \Sigma^{-1} \mu + \frac{1}{\nu + \mu^{\top} \Sigma^{-1} \mu} \text{tr}(\Sigma^{-1} p^{(\alpha)}(x x^{\top})), \qquad (A.7) $$

which is defined if $p^{(\alpha)}$ has finite first and second order moments. Introducing the quantity $V := p^{(\alpha)}((x - p^{(\alpha)}(x))(x - p^{(\alpha)}(x))^{\top}) = p^{(\alpha)}(xx^{\top}) - p^{(\alpha)}(x)p^{(\alpha)}(x)^{\top} \in S_{++}^d$, we get for any $\vartheta \in \text{dom } \varphi_\lambda$ that

$$ 
\begin{aligned}
1 + \lambda \langle \vartheta, p^{(\alpha)}(T) \rangle &= 1 - \frac{2}{\nu + \mu^{\top} \Sigma^{-1} \mu} p^{(\alpha)}(x)^{\top} \Sigma^{-1} \mu + \frac{1}{\nu + \mu^{\top} \Sigma^{-1} \mu} \text{tr}(\Sigma^{-1} V) \\
&\quad + \frac{1}{\nu + \mu^{\top} \Sigma^{-1} \mu} p^{(\alpha)}(x)^{\top} \Sigma^{-1} p^{(\alpha)}(x) \\
&= \frac{1}{\nu + \mu^{\top} \Sigma^{-1} \mu} \left( \nu + (\mu - p^{(\alpha)}(x))^{\top} \Sigma^{-1} (\mu - p^{(\alpha)}(x))^{\top} + \text{tr}(\Sigma^{-1} V) \right) \\
&\geq \frac{\nu}{\nu + \mu^{\top} \Sigma^{-1} \mu}.
\end{aligned}
$$

This shows that for any $\vartheta \in \text{dom } \varphi_\lambda$, and $p^{(\alpha)} \in \mathcal{P}(\mathcal{X}, dx)$ with finite first and second order moments, the quantity $c_\lambda(\vartheta, p^{(\alpha)}(T))$ is in $\mathbb{R}$. With the result of (ii), this shows that $\mathcal{T}_{\nu}^d$, seen as an instance of the $\lambda$-exponential family, satisfies Assumptions 2. $\square$
## References

[1] O. Akyildiz and J. Míguez. Convergence rates for optimised adaptive importance samplers. *Statistics and Computing*, 31(12), 2021.

[2] S.-I. Amari. *Differential-Geometrical Methods in Statistics*. Springer New York, 1985.

[3] S.-I. Amari and A. Ohara. Geometry of q-exponential family of probability distributions. *Entropy*, 13(6):1170–1185, 2011.

[4] I. Ayadi, F. Bouchard, and F. Pascal. Elliptical Wishart distribution: Maximum likelihood estimator from information geometry. In *IEEE International Conference on Speech, Acoustics and Signal Processing (ICASSP)*, 2023.

[5] A. Banerjee, S. Merugu, I. S. Dhillon, and J. Ghosh. Clustering with Bregman divergences. *Journal of Machine Learning Research*, 6(58):1705–1749, 2005.

[6] O. Barndorff-Nielsen. *Information and Exponential Families in Statistical Theory*. John Wiley & Sons, Ltd, 2014.

[7] H. Bauschke, J. Borwein, and P. Combettes. Bregman monotone optimization algorithms. *SIAM Journal on Control and Optimization*, 42(2):596–636, 2003.

[8] H. Bauschke and P. Combettes. *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*. Springer, 2011.

[9] E. M. Bednarczuk and M. Syga. On duality for nonconvex minimization problems within the framework of abstract convexity. *Optimization*, 71(4):949–971, 2022.

[10] C. Bishop. *Pattern Recognition and Machine Learning*. Springer, 2006.

[11] M. F. Bugallo, V. Elvira, L. Martino, D. Luengo, J. Míguez, and P. M. Djuric. Adaptive importance sampling: The past, the present, and the future. *IEEE Signal Processing Magazine*, 34(4):60–79, 2017.

[12] L. L. Campbell. Equivalence of Gauss's principle and minimum discrimination information estimation of probabilities. *Annals of Mathematical Statistics*, 41(3):1011–1015, 1970.

[13] O. Cappé, R. Douc, A. Guillin, J. M. Marin, and C. P. Robert. Adaptive importance sampling in general mixture classes. *Statistics and Computing*, 18:447–459, 2008.

[14] J.-P. Chancelier and M. De Lara. Constant along primal rays conjugacies and the $l_0$ pseudonorm. *Optimization*, 71(2):355–386, 2020.

[15] G. Chantas, N. Galatsanos, A. Likas, and M. Saunders. Variational Bayesian image restoration based on a product of t-distributions image prior. *IEEE Transactions on Image Processing*, 17(10):1795–1805, 2008.

[16] P. Douglas, S. Bergamini, and F. Renzoni. Tunable Tsallis distributions in dissipative optical lattices. *Physical Review Letters*, 96:110601, 2006.

[17] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo. Generalized multiple importance sampling. *Statistical Science*, 34(1):129–155, 2019.

[18] M. Fajardo and J. Vidal. On subdifferentials via a generalized conjugation scheme: an application to DC problems and optimality conditions. *Set-Valued and Variational Analysis*, 30:1313–1331, 2022.
References

[1] O. Akyildiz and J. Míguez. Convergence rates for optimised adaptive importance samplers. *Statistics and Computing*, 31(12), 2021.

[2] S.-I. Amari. *Differential-Geometrical Methods in Statistics*. Springer New York, 1985.

[3] S.-I. Amari and A. Ohara. Geometry of q-exponential family of probability distributions. *Entropy*, 13(6):1170–1185, 2011.

[4] I. Ayadi, F. Bouchard, and F. Pascal. Elliptical Wishart distribution: Maximum likelihood estimator from information geometry. In *IEEE International Conference on Speech, Acoustics and Signal Processing (ICASSP)*, 2023.

[5] A. Banerjee, S. Merugu, I. S. Dhillon, and J. Ghosh. Clustering with Bregman divergences. *Journal of Machine Learning Research*, 6(58):1705–1749, 2005.

[6] O. Barndorff-Nielsen. *Information and Exponential Families in Statistical Theory*. John Wiley & Sons, Ltd, 2014.

[7] H. Bauschke, J. Borwein, and P. Combettes. Bregman monotone optimization algorithms. *SIAM Journal on Control and Optimization*, 42(2):596–636, 2003.

[8] H. Bauschke and P. Combettes. *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*. Springer, 2011.

[9] E. M. Bednarczuk and M. Syga. On duality for nonconvex minimization problems within the framework of abstract convexity. *Optimization*, 71(4):949–971, 2022.

[10] C. Bishop. *Pattern Recognition and Machine Learning*. Springer, 2006.

[11] M. F. Bugallo, V. Elvira, L. Martino, D. Luengo, J. Míguez, and P. M. Djuric. Adaptive importance sampling: The past, the present, and the future. *IEEE Signal Processing Magazine*, 34(4):60–79, 2017.

[12] L. L. Campbell. Equivalence of Gauss's principle and minimum discrimination information estimation of probabilities. *Annals of Mathematical Statistics*, 41(3):1011–1015, 1970.

[13] O. Cappé, R. Douc, A. Guillin, J. M. Marin, and C. P. Robert. Adaptive importance sampling in general mixture classes. *Statistics and Computing*, 18:447–459, 2008.

[14] J.-P. Chancelier and M. De Lara. Constant along primal rays conjugacies and the $l_0$ pseudonorm. *Optimization*, 71(2):355–386, 2020.

[15] G. Chantas, N. Galatsanos, A. Likas, and M. Saunders. Variational Bayesian image restoration based on a product of t-distributions image prior. *IEEE Transactions on Image Processing*, 17(10):1795–1805, 2008.

[16] P. Douglas, S. Bergamini, and F. Renzoni. Tunable Tsallis distributions in dissipative optical lattices. *Physical Review Letters*, 96:110601, 2006.

[17] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo. Generalized multiple importance sampling. *Statistical Science*, 34(1):129–155, 2019.

[18] M. Fajardo and J. Vidal. On subdifferentials via a generalized conjugation scheme: an application to DC problems and optimality conditions. *Set-Valued and Variational Analysis*, 30:1313–1331, 2022.
[19] A. Gelman, A. Jakulin, M. G. Pittau, and Y.-S. Su. A weakly informative default prior distribution for logistic and other regression models. *Annals of Applied Statistics*, 2(4):1360–1383, 2008.

[20] T. Guilmeau, V. Elvira, and E. Chouzenoux. Regularized Rényi divergence minimization through Bregman proximal gradient algorithms. Preprint, https://arxiv.org/abs/2211.04776, 2022.

[21] M. Hasanasab, J. Hertrich, and G. Steidl. Alternatives to the EM algorithm for estimating the parameters of the Student t-distribution. *Numerical Algorithms*, 87:77–118, 2021.

[22] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. *Journal of Machine Learning Research*, 14(4):1303–1347, 2013.

[23] S. F. Jarner and G. O. Roberts. Convergence of heavy-tailed Monte Carlo Markov chain algorithms. *Scandinavian Journal of Statistics*, 34(4):781–815, 2007.

[24] A. S. Kainth, T.-K. L. Wong, and F. Rudzicz. Conformal mirror descent with logarithmic divergences. *Information Geometry*, 2022.

[25] M. Khan and W. Lin. Conjugate-computation variational inference: Converting variational inference in non-conjugate models to inferences in conjugate models. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, pages 878–887, 2017.

[26] M. Kot, M. A. Lewis, and P. van Den Driessche. Dispersal data and the spread of invading organisms. *Ecology*, 77(7):2027–2042, 1996.

[27] S. Kullback and R. A. Leibler. On information and sufficiency. *The Annals of Mathematical Statistics*, 22(1):79–86, 1951.

[28] F. Laus, F. Pierre, and G. Steidl. Nonlocal myriad filters for Cauchy noise removal. *Journal of Mathematical Imaging and Vision*, 60:1324–1354, 2018.

[29] A. Le Franc, J.-P. Chancelier, and M. De Lara. The Capra-subdifferential of the $l_0$ pseudonorm. *Optimization*, pages 1–23, 2022.

[30] Y. Marnissi, E. Chouzenoux, A. Benazza-Benyahia, and J.-C. Pesquet. Majorize-minimize adapted Metropolis-Hastings algorithm. *IEEE Transactions on Signal Processing*, 68:2356–2369, 2020.

[31] J. Martin, L. C. Wilcox, C. Burstedde, and O. Ghattas. A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion. *SIAM Journal of Scientific Computing*, 34(3):A1460–A1487, 2012.

[32] A. F. T. Martins, M. Treviso, A. Farinhas, P. M. Q. Aguiar, M. A. T. Figueiredo, M. Blondel, and V. Niculae. Sparse continuous distributions and Fenchel Young losses. *Journal of Machine Learning Research*, 23(257):1–74, 2022.

[33] V. Masrani, R. Brekelmans, T. Bui, F. Nielsen, A. Galstyan, G. V. Steeg, and F. Wood. q-paths: Generalizing the geometric annealing path using power means. In *Conference on Uncertainty in Artificial Intelligence (UAI)*, volume 161, pages 1938–1947, 2021.

[34] J. J. Moré and G. Toraldo. Algorithms for bound constrained quadratic programming problems. *Numerische Mathematik*, 55(4):377–400, 1989.

[35] F. Nielsen and R. Nock. Entropies and cross-entropies of exponential families. In *IEEE International Conference on Image Processing (ICIP)*, pages 3621–3624, 2010.
[19] A. Gelman, A. Jakulin, M. G. Pittau, and Y.-S. Su. A weakly informative default prior distribution for logistic and other regression models. *Annals of Applied Statistics*, 2(4):1360–1383, 2008.

[20] T. Guilmeau, V. Elvira, and E. Chouzenoux. Regularized Rényi divergence minimization through Bregman proximal gradient algorithms. Preprint, https://arxiv.org/abs/2211.04776, 2022.

[21] M. Hasanasab, J. Hertrich, and G. Steidl. Alternatives to the EM algorithm for estimating the parameters of the Student t-distribution. *Numerical Algorithms*, 87:77–118, 2021.

[22] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. *Journal of Machine Learning Research*, 14(4):1303–1347, 2013.

[23] S. F. Jarner and G. O. Roberts. Convergence of heavy-tailed Monte Carlo Markov chain algorithms. *Scandinavian Journal of Statistics*, 34(4):781–815, 2007.

[24] A. S. Kainth, T.-K. L. Wong, and F. Rudzicz. Conformal mirror descent with logarithmic divergences. *Information Geometry*, 2022.

[25] M. Khan and W. Lin. Conjugate-computation variational inference: Converting variational inference in non-conjugate models to inferences in conjugate models. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, pages 878–887, 2017.

[26] M. Kot, M. A. Lewis, and P. van Den Driessche. Dispersal data and the spread of invading organisms. *Ecology*, 77(7):2027–2042, 1996.

[27] S. Kullback and R. A. Leibler. On information and sufficiency. *The Annals of Mathematical Statistics*, 22(1):79–86, 1951.

[28] F. Laus, F. Pierre, and G. Steidl. Nonlocal myriad filters for Cauchy noise removal. *Journal of Mathematical Imaging and Vision*, 60:1324–1354, 2018.

[29] A. Le Franc, J.-P. Chancelier, and M. De Lara. The Capra-subdifferential of the $l_0$ pseudonorm. *Optimization*, pages 1–23, 2022.

[30] Y. Marnissi, E. Chouzenoux, A. Benazza-Benyahia, and J.-C. Pesquet. Majorize-minimize adapted Metropolis-Hastings algorithm. *IEEE Transactions on Signal Processing*, 68:2356–2369, 2020.

[31] J. Martin, L. C. Wilcox, C. Burstedde, and O. Ghattas. A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion. *SIAM Journal of Scientific Computing*, 34(3):A1460–A1487, 2012.

[32] A. F. T. Martins, M. Treviso, A. Farinhas, P. M. Q. Aguiar, M. A. T. Figueiredo, M. Blondel, and V. Niculae. Sparse continuous distributions and Fenchel Young losses. *Journal of Machine Learning Research*, 23(257):1–74, 2022.

[33] V. Masrani, R. Brekelmans, T. Bui, F. Nielsen, A. Galstyan, G. V. Steeg, and F. Wood. q-paths: Generalizing the geometric annealing path using power means. In *Conference on Uncertainty in Artificial Intelligence (UAI)*, volume 161, pages 1938–1947, 2021.

[34] J. J. Moré and G. Toraldo. Algorithms for bound constrained quadratic programming problems. *Numerische Mathematik*, 55(4):377–400, 1989.

[35] F. Nielsen and R. Nock. Entropies and cross-entropies of exponential families. In *IEEE International Conference on Image Processing (ICIP)*, pages 3621–3624, 2010.
[36] R. Nock, Z. Cranko, A. K. Menon, L. Qu, and R. C. Williamson. f-GANs in an information geometric nutshell. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 30, 2017.

[37] S. Nowozin, B. Cseke, and R. Tomioka. f-GAN: Training generative neural samplers using variational divergence minimization. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 29, 2016.

[38] D. Peel and G. J. McLachlan. Robust mixture modelling using the t-distribution. *Statistics and Computing*, 10:339-348, 2000.

[39] S. T. Rachev and L. Rüschendorf. *Mass Transportation Problems*. Springer-Verlag, 1998.

[40] G. O. Roberts and J. S. Rosenthal. Optimal scaling for various Metropolis-Hastings algorithms. *Statistical Science*, 16(4):351-367, 2001.

[41] G. O. Roberts and O. Stramer. Langevin diffusions and Metropolis-Hastings algorithms. *Methodology and Computing in Applied Probability*, 4(4):337-357, 2002.

[42] A. Rényi. On measures of entropy and information. In *Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability*, volume 1, pages 547-561, 1961.

[43] M. Teboulle. A simplified view of first order methods for optimization. *Mathematical Programming*, 170(1):67-96, 2018.

[44] Y. Tikochinsky, N. Z. Tishby, and R. D. Levine. Alternative approach to maximum-entropy inference. *Physical Review A*, 30:2638-2644, 1984.

[45] T. van Erven and P. Harremoës. Rényi divergence and Kullback-Leibler divergence. *IEEE Transactions on Information Theory*, 60(7):3797-3820, 2014.

[46] M. J. Wainwright and M. I. Jordan. Graphical models, exponential families, and variational inference. *Foundations and Trends® in Machine Learning*, 1(1-2):1-305, 2008.

[47] S. Wang and T. Swartz. Moment matching adaptive importance sampling with skew-Student proposals. *Monte Carlo Methods and Applications*, 28(2):149-162, 2022.

[48] H. White. Maximum likelihood estimation of misspecified models. *Econometrica*, 50(1):1-25, 1982.

[49] T.-K. L. Wong. Logarithmic divergences from optimal transport and Rényi geometry. *Information Geometry*, 1(1):39-78, 2018.

[50] T.-K. L. Wong and J. Zhang. Tsallis and Rényi deformations linked via a new λ-duality. *IEEE Transactions on Information Theory*, 68(8):5353-5373, 2022.
[36] R. Nock, Z. Cranko, A. K. Menon, L. Qu, and R. C. Williamson. f-GANs in an information geometric nutshell. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 30, 2017.

[37] S. Nowozin, B. Cseke, and R. Tomioka. f-GAN: Training generative neural samplers using variational divergence minimization. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 29, 2016.

[38] D. Peel and G. J. McLachlan. Robust mixture modelling using the t-distribution. *Statistics and Computing*, 10:339-348, 2000.

[39] S. T. Rachev and L. Rüschendorf. *Mass Transportation Problems*. Springer-Verlag, 1998.

[40] G. O. Roberts and J. S. Rosenthal. Optimal scaling for various Metropolis-Hastings algorithms. *Statistical Science*, 16(4):351-367, 2001.

[41] G. O. Roberts and O. Stramer. Langevin diffusions and Metropolis-Hastings algorithms. *Methodology and Computing in Applied Probability*, 4(4):337-357, 2002.

[42] A. Rényi. On measures of entropy and information. In *Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability*, volume 1, pages 547-561, 1961.

[43] M. Teboulle. A simplified view of first order methods for optimization. *Mathematical Programming*, 170(1):67-96, 2018.

[44] Y. Tikochinsky, N. Z. Tishby, and R. D. Levine. Alternative approach to maximum-entropy inference. *Physical Review A*, 30:2638-2644, 1984.

[45] T. van Erven and P. Harremoës. Rényi divergence and Kullback-Leibler divergence. *IEEE Transactions on Information Theory*, 60(7):3797-3820, 2014.

[46] M. J. Wainwright and M. I. Jordan. Graphical models, exponential families, and variational inference. *Foundations and Trends® in Machine Learning*, 1(1-2):1-305, 2008.

[47] S. Wang and T. Swartz. Moment matching adaptive importance sampling with skew-Student proposals. *Monte Carlo Methods and Applications*, 28(2):149-162, 2022.

[48] H. White. Maximum likelihood estimation of misspecified models. *Econometrica*, 50(1):1-25, 1982.

[49] T.-K. L. Wong. Logarithmic divergences from optimal transport and Rényi geometry. *Information Geometry*, 1(1):39-78, 2018.

[50] T.-K. L. Wong and J. Zhang. Tsallis and Rényi deformations linked via a new λ-duality. *IEEE Transactions on Information Theory*, 68(8):5353-5373, 2022.