When $\alpha = 1$, we can split the subdifferential of $F_{\pi}^{(\alpha)}$ as $\partial F_{\pi}^{(1)}(\theta_{*}) = \pi(\Gamma) + \partial A(\theta_{*}) + \partial r(\theta_{*})$. This comes from the decomposition (23), Assumption 3 and the convexity and properness of $\theta \mapsto -\langle \theta, \pi(\Gamma) \rangle$, A and r, and [10, Corollary 16.38]. By the same arguments, when $\alpha > 1$, $\partial F_{\pi}^{(\alpha)}(\theta_{*}) = \partial \left( \frac{1}{\alpha-1} \log \circ h_{\pi}^{(\alpha)} \right) (\theta_{*}) + \partial A(\theta_{*}) + \partial r(\theta_{*})$

Assume by contradiction that $\theta_{*}$ belongs to the boundary of $\Theta$. Then $\partial A(\theta_{*}) = \emptyset$, because of Proposition 2, so Eq. (26) implies that $0 \in \emptyset$. This shows that $\theta_{*} \in \text{int} \Theta$.

Finally, since A is strictly convex on int $\Theta$ (Proposition 4), so is $F_{\pi}^{(1)}$, so such $\theta_{*}$ is unique. $\square$

## A.7 Proof of Proposition 11

*Proof of Proposition 11.* Consider the family of one-dimensional centered Gaussian distributions with variance $\sigma^2$, that we denote by $G_0^1$ in the following. It is an exponential family, with parameter $\theta = -\frac{1}{2\sigma^2}$, sufficient statistics $\Gamma(x) = x^2$ and log-partition function $A(\theta) = \frac{1}{2}\log(2\pi) - \frac{1}{2}\log(-2\theta)$, whose domain is $\Theta = \mathbb{R}_{-}$. We show that $(f_{\pi}^{(\alpha)})'$ is not Lipschitz for $\alpha > 0$ by showing that $(f_{\pi}^{(\alpha)})''$ is unbounded on $\Theta$.

Consider first the case $\alpha = 1$. From Proposition 6, $(f_{\pi}^{(\alpha)})''$ is independent of the choice of the target $\pi$, and is equal to

$$ (f_{\pi}^{(\alpha)})''(\theta) = A''(\theta) = \frac{1}{2\theta^2}. \qquad (27) $$

Now, for $\alpha \neq 1$, consider a target $\pi \in G_0^1$, meaning that there exists $\theta_{\pi} \in \Theta$ such that $\pi = q_{\theta_{\pi}}$. We can compute that $\pi_{\theta}^{(\alpha)} = q_{\alpha\theta_{\pi}+(1-\alpha)\theta}$, assuming that $\theta$ is such that $\alpha\theta_{\pi} + (1-\alpha)\theta \in \Theta$. This condition is always satisfied when $\alpha \le 1$, but when $\alpha > 1$, it is equivalent to having $\theta > \frac{\alpha}{\alpha-1}\theta_{\pi}$. In the case $\alpha > 1$, $f_{\pi}^{(\alpha)}$ is not even defined outside of $(\frac{\alpha}{\alpha-1}\theta_{\pi}, 0)$, showing that dom $f_{\pi}^{(\alpha)} = (0, \frac{\alpha}{\alpha-1}\theta_{\pi})$ for $\alpha > 1$. In the following, we consider $\theta \in \text{dom } f_{\pi}^{(\alpha)}$. Then we compute using the result of Proposition 6 that

$$ (f_{\pi}^{(\alpha)})''(\theta) = \frac{1}{2\theta^2} + (\alpha - 1) \left( \int x^4 q_{\alpha\theta_{\pi}+(1-\alpha)\theta}(x)dx - \left( \int x^2 q_{\alpha\theta_{\pi}+(1-\alpha)\theta}(x)dx \right)^2 \right). $$

To do so, we recall the following formulas

$$ \int x^4 \exp(-bx^2)dx = \frac{3\sqrt{\pi}}{4b^{5/2}}, \qquad \int x^2 \exp(-bx^2)dx = \frac{\sqrt{\pi}}{2b^{3/2}}, $$

and we note that $A(\theta) = \log(\sqrt{-\frac{\pi}{\theta}})$. We first compute

$$ \begin{aligned} \int x^4 q_{\alpha\theta_{\pi}+(1-\alpha)\theta}(x)dx &= \exp(A(\alpha\theta_{\pi} + (1-\alpha)\theta))^{-1} \int x^4 \exp(-((\alpha-1)\theta - \alpha\theta_{\pi})x^2)dx \\ &= \left( \frac{\pi}{(\alpha-1)\theta - \alpha\theta_{\pi}} \right)^{-1/2} \frac{3\sqrt{\pi}}{4((\alpha-1)\theta - \alpha\theta_{\pi})^{5/2}} \\ &= \frac{3}{4((\alpha-1)\theta - \alpha\theta_{\pi})^2}, \end{aligned} $$

and then by similar means $\int x^2 q_{\alpha\theta_{\pi}+(1-\alpha)\theta}(x)dx = \frac{1}{2((\alpha-1)\theta - \alpha\theta_{\pi})}$. These calculations yield

$$ (f_{\pi}^{(\alpha)})''(\theta) = \frac{1}{2\theta^2} + \frac{\alpha - 1}{2((\alpha - 1)\theta - \alpha\theta_{\pi})^2}. \qquad (28) $$

Equations (27) and (28) show that the absolute value of $\nabla^2 f_{\pi}^{(\alpha)}$ goes to $+\infty$ when $\theta$ approaches 0 or $\frac{\alpha}{\alpha-1}\theta_{\pi}$, which is in $\Theta$ if and only if $\alpha > 1$. $\square$
## A.8 Proof of Proposition 12

*Proof of Proposition 12.* Since $\pi = q_{\theta_{\pi}}$, we can compute that $\pi_{\theta}^{(\alpha)} = q_{\alpha\theta_{\pi}+(1-\alpha)\theta}$. Since $\Theta$ is convex (from Proposition 4) and $\theta_{\pi}, \theta \in \text{int} \Theta$, $\alpha\theta_{\pi} + (1-\alpha)\theta \in \text{int} \Theta$. Therefore, $\pi_{\theta}^{(\alpha)}(\Gamma) = \nabla A(\alpha\theta_{\pi} + (1-\alpha)\theta) \in \text{int dom } A^*$, since $A$ is Legendre and $\alpha\theta_{\pi} + (1-\alpha)\theta \in \text{int dom } A$, showing that Assumption 2 is satisfied.

Recall that $RD_{\alpha}(\pi, q_{\theta}) \ge 0$, and it is equal to zero if and only if $\pi = q_{\theta}$. In our case, this means that $f_{\pi}^{(\alpha)}(\theta) = 0$ if and only if $\theta = \theta_{\pi}$. Since $Q$ is minimal by assumption, this shows the existence and unicity of the minimizer of $f_{\pi}^{(\alpha)}$.

Consider now a stationary point of $f_{\pi}^{(\alpha)}$ i.e., $\theta \in \text{int} \Theta$ such that $\nabla f_{\pi}^{(\alpha)}(\theta) = 0$. This implies that

$$ q_{\theta}(\Gamma) = \pi_{\theta}^{(\alpha)}(\Gamma), \qquad (29) $$

due to the characterization given in Proposition 6. Under our assumptions, Eq. (29) now reads

$$ \nabla A(\theta) = \nabla A(\alpha\theta_{\pi} + (1-\alpha)\theta), $$

which is equivalent to having $\theta = \theta_{\pi}$ by inverting $\nabla A$ on both sides. Hence we have shown that $\nabla f_{\pi}^{(\alpha)}(\theta) = 0$ if and only $\theta = \theta_{\pi}$, showing the existence and unicity of the stationary point of $f_{\pi}^{(\alpha)}$. $\square$

## A.9 Proof of Proposition 13

*Proof of Proposition 13.* Under our Assumptions on $\pi$, we can compute

$$ f_{\pi}^{(\alpha)}(\theta) = \frac{1}{1-\alpha} (\alpha A(\theta_{\pi}) + (1-\alpha)A(\theta) - A(\alpha\theta_{\pi} + (1-\alpha)\theta)), \forall \theta \in \text{int} \Theta \qquad (30) $$

Consider in the following $\theta, \theta' \in \text{int} \Theta$, with $\theta, \theta' \in B(\theta_{\pi}, \nu)$.

(i) We can compute

$$ \begin{aligned} A(\alpha\theta_{\pi} + (1-\alpha)\theta) &= A(\theta_{\pi} + (1-\alpha)(\theta - \theta_{\pi})) \\ &= A(\theta_{\pi}) + (1-\alpha)\langle\nabla A(\theta_{\pi}), \theta - \theta_{\pi}\rangle + \frac{(1-\alpha)^2}{2}\|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})} + o(v^2), \\ A(\theta) &= A(\theta_{\pi}) + \langle\nabla A(\theta_{\pi}), \theta - \theta_{\pi}\rangle + \frac{1}{2}\|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})} + o(v^2). \end{aligned} $$

Using Eq. (30), these two equalities imply in particular that

$$ f_{\pi}^{(\alpha)}(\theta) = \frac{\alpha}{2} \|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})}^2 + o(v^2). \qquad (31) $$

(ii) We now turn to the Polyak-Łojasiewicz inequality. We begin by showing that $d_{f_{\pi}^{(\alpha)}}(\theta, \theta') = \frac{\alpha}{2}\|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2$ up to higher order terms. We now further compute that

$$ \begin{aligned} \nabla f_{\pi}^{(\alpha)}(\theta) &= \nabla A(\theta) - \nabla A(\alpha\theta_{\pi} + (1-\alpha)\theta_{\pi}) \\ &= \nabla A(\theta_{\pi} + (\theta - \theta_{\pi})) - \nabla A(\theta_{\pi} + (1-\alpha)(\theta - \theta_{\pi})) \\ &= \nabla A(\theta_{\pi}) + \nabla^2 A(\theta_{\pi})(\theta - \theta_{\pi}) - \nabla A(\theta_{\pi}) - (1-\alpha)\nabla^2 A(\theta_{\pi})(\theta - \theta_{\pi}) + o(v) \\ &= \alpha\nabla^2 A(\theta_{\pi})(\theta - \theta_{\pi}) + o(v), \end{aligned} $$
## A.8 Proof of Proposition 12

*Proof of Proposition 12.* Since $\pi = q_{\theta_{\pi}}$, we can compute that $\pi_{\theta}^{(\alpha)} = q_{\alpha\theta_{\pi}+(1-\alpha)\theta}$. Since $\Theta$ is convex (from Proposition 4) and $\theta_{\pi}, \theta \in \text{int} \Theta$, $\alpha\theta_{\pi} + (1-\alpha)\theta \in \text{int} \Theta$. Therefore, $\pi_{\theta}^{(\alpha)}(\Gamma) = \nabla A(\alpha\theta_{\pi} + (1-\alpha)\theta) \in \text{int dom } A^*$, since $A$ is Legendre and $\alpha\theta_{\pi} + (1-\alpha)\theta \in \text{int dom } A$, showing that Assumption 2 is satisfied.

Recall that $RD_{\alpha}(\pi, q_{\theta}) \ge 0$, and it is equal to zero if and only if $\pi = q_{\theta}$. In our case, this means that $f_{\pi}^{(\alpha)}(\theta) = 0$ if and only if $\theta = \theta_{\pi}$. Since $Q$ is minimal by assumption, this shows the existence and unicity of the minimizer of $f_{\pi}^{(\alpha)}$.

Consider now a stationary point of $f_{\pi}^{(\alpha)}$ i.e., $\theta \in \text{int} \Theta$ such that $\nabla f_{\pi}^{(\alpha)}(\theta) = 0$. This implies that

$$q_{\theta}(\Gamma) = \pi_{\theta}^{(\alpha)}(\Gamma), \qquad (29)$$

due to the characterization given in Proposition 6. Under our assumptions, Eq. (29) now reads

$$\nabla A(\theta) = \nabla A(\alpha\theta_{\pi} + (1-\alpha)\theta),$$

which is equivalent to having $\theta = \theta_{\pi}$ by inverting $\nabla A$ on both sides. Hence we have shown that $\nabla f_{\pi}^{(\alpha)}(\theta) = 0$ if and only $\theta = \theta_{\pi}$, showing the existence and unicity of the stationary point of $f_{\pi}^{(\alpha)}$. $\square$

## A.9 Proof of Proposition 13

*Proof of Proposition 13.* Under our Assumptions on $\pi$, we can compute

$$f_{\pi}^{(\alpha)}(\theta) = \frac{1}{1-\alpha}(\alpha A(\theta_{\pi}) + (1-\alpha)A(\theta) - A(\alpha\theta_{\pi} + (1-\alpha)\theta)), \forall \theta \in \text{int} \Theta \qquad (30)$$

Consider in the following $\theta, \theta' \in \text{int} \Theta$, with $\theta, \theta' \in B(\theta_{\pi}, \nu)$.

(i) We can compute

$$\begin{aligned} A(\alpha\theta_{\pi} + (1-\alpha)\theta) &= A(\theta_{\pi} + (1-\alpha)(\theta - \theta_{\pi})) \\ &= A(\theta_{\pi}) + (1-\alpha)\langle\nabla A(\theta_{\pi}), \theta - \theta_{\pi}\rangle + \frac{(1-\alpha)^2}{2}\|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})} + o(v^2), \\ A(\theta) &= A(\theta_{\pi}) + \langle\nabla A(\theta_{\pi}), \theta - \theta_{\pi}\rangle + \frac{1}{2}\|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})} + o(v^2). \end{aligned}$$

Using Eq. (30), these two equalities imply in particular that

$$f_{\pi}^{(\alpha)}(\theta) = \frac{\alpha}{2} \|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})}^2 + o(v^2). \qquad (31)$$

(ii) We now turn to the Polyak-Łojasiewicz inequality. We begin by showing that $d_{f_{\pi}^{(\alpha)}}(\theta, \theta') = \frac{\alpha}{2}\|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2$ up to higher order terms. We now further compute that

$$\begin{aligned} \nabla f_{\pi}^{(\alpha)}(\theta) &= \nabla A(\theta) - \nabla A(\alpha\theta_{\pi} + (1-\alpha)\theta_{\pi}) \\ &= \nabla A(\theta_{\pi} + (\theta - \theta_{\pi})) - \nabla A(\theta_{\pi} + (1-\alpha)(\theta - \theta_{\pi})) \\ &= \nabla A(\theta_{\pi}) + \nabla^2 A(\theta_{\pi})(\theta - \theta_{\pi}) - \nabla A(\theta_{\pi}) - (1-\alpha)\nabla^2 A(\theta_{\pi})(\theta - \theta_{\pi}) + o(v) \\ &= \alpha\nabla^2 A(\theta_{\pi})(\theta - \theta_{\pi}) + o(v), \end{aligned}$$
which yields with Eq. (31) that

$$
\begin{align*}
d_{f_{\pi}^{(\alpha)}}(\theta, \theta') &= f_{\pi}^{(\alpha)}(\theta) - f_{\pi}^{(\alpha)}(\theta') - \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \theta - \theta' \rangle \\
&= \frac{\alpha}{2} \|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})}^2 - \frac{\alpha}{2} \|\theta' - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})}^2 - \alpha \langle \nabla^2 A(\theta_{\pi})(\theta' - \theta_{\pi}), \theta - \theta' \rangle + o(v^2) \\
&= \frac{\alpha}{2} \|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2 + o(v^2).
\end{align*}
$$

We thus have that

$$
f_{\pi}^{(\alpha)}(\theta) = f_{\pi}^{(\alpha)}(\theta') + \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \theta - \theta' \rangle + \frac{\alpha}{2} \|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2 + o(v^2). \quad (32)
$$

When $\theta'$ is fixed, the quantity $\theta \mapsto f_{\pi}^{(\alpha)}(\theta') + \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \theta - \theta' \rangle + \frac{\alpha}{2} \|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2$ is a quadratic form that is minimized when the following optimality condition is satisfied:

$$
\nabla f_{\pi}(\theta') + \alpha \nabla^2 A(\theta_{\pi})(\theta - \theta') = 0. \quad (33)
$$

We now compute the inverse of $\nabla^2 A(\theta_\pi)$. Consider $\eta \in \text{int dom } A^*$. Since $A$ is Legendre, we have that $\nabla A(\nabla A^*(\eta)) = \eta$, therefore differentiating this expression with respect to $\eta$ yields

$$
\nabla^2 A^*(\eta) \nabla^2 A(\nabla A^*(\eta)) = \mathrm{Id}.
$$

Now take $\eta = \nabla A(\theta_\pi)$ which belongs to int dom $A^*$ since $\theta_\pi \in \text{int dom } A$. We obtain that $\nabla^2 A^*(\nabla A(\theta_\pi))\nabla^2 A(\theta_\pi) = \text{Id}$. This shows that the optimality condition of Eq. (33) is equivalent to having

$$
\theta - \theta' = -\frac{1}{\alpha} \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta').
$$

This shows that the right-hand side of Eq. (32) can be minorized to obtain:

$$
\begin{align*}
f_{\pi}^{(\alpha)}(\theta) \ge f_{\pi}^{(\alpha)}(\theta') & - \frac{1}{\alpha} \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta') \rangle \\
& + \frac{1}{2\alpha} \langle \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta'), \nabla^2 A(\theta_{\pi}) \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta') \rangle + o(v^2) \\
& = f_{\pi}^{(\alpha)}(\theta') - \frac{1}{2\alpha} \|\nabla f_{\pi}^{(\alpha)}(\theta')\|_{\nabla^2 A^*(\nabla A(\theta_{\pi}))}^2 + o(v^2).
\end{align*}
$$

This is true in particular for $\theta = \theta_{\pi}$, which yields the result.

<figure><img src="image_8.png" /></figure>

B Convergence analysis of Algorithm 1

In order to prove Propositions 14 and 15, we start with a sufficient decrease lemma that reads as follows.

**Lemma 1.** Under Assumptions 1, 2, and 3, for $\tau > 0$ and $\alpha \in (0, 1]$, we have that for every $\theta \in \text{int } \Theta$,

$$
\tau \left( F_{\pi}^{(\alpha)}(T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta)) - F_{\pi}^{(\alpha)}(\theta) \right) \le -d_{A}(\theta, T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta)) + (\tau - 1)d_{A}(T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta), \theta). \quad (34)
$$

In the particular case where $\alpha = 1$, we further have

$$
\tau \left( F_{\pi}^{(1)}(T_{\tau F_{\pi}^{(1)}}^{A}(\theta)) - F_{\pi}^{(1)}(\theta') \right) \le (1-\tau)d_{A}(\theta', \theta) - (1-\tau)d_{A}(T_{\tau F_{\pi}^{(1)}}^{A}(\theta), \theta) \\
\qquad - d_{A}(\theta', T_{\tau F_{\pi}^{(1)}}^{A}(\theta)), \quad \forall \theta' \in \text{int } \Theta. \quad (35)
$$
which yields with Eq. (31) that

$$
\begin{align*}
d_{f_{\pi}^{(\alpha)}}(\theta, \theta') &= f_{\pi}^{(\alpha)}(\theta) - f_{\pi}^{(\alpha)}(\theta') - \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \theta - \theta' \rangle \\
&= \frac{\alpha}{2} \|\theta - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})}^2 - \frac{\alpha}{2} \|\theta' - \theta_{\pi}\|_{\nabla^2 A(\theta_{\pi})}^2 - \alpha \langle \nabla^2 A(\theta_{\pi})(\theta' - \theta_{\pi}), \theta - \theta' \rangle + o(v^2) \\
&= \frac{\alpha}{2} \|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2 + o(v^2).
\end{align*}
$$

We thus have that

$$
f_{\pi}^{(\alpha)}(\theta) = f_{\pi}^{(\alpha)}(\theta') + \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \theta - \theta' \rangle + \frac{\alpha}{2} \|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2 + o(v^2). \quad (32)
$$

When $\theta'$ is fixed, the quantity $\theta \mapsto f_{\pi}^{(\alpha)}(\theta') + \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \theta - \theta' \rangle + \frac{\alpha}{2} \|\theta - \theta'\|_{\nabla^2 A(\theta_{\pi})}^2$ is a quadratic form that is minimized when the following optimality condition is satisfied:

$$
\nabla f_{\pi}(\theta') + \alpha \nabla^2 A(\theta_{\pi})(\theta - \theta') = 0. \quad (33)
$$

We now compute the inverse of $\nabla^2 A(\theta_\pi)$. Consider $\eta \in \text{int dom } A^*$. Since $A$ is Legendre, we have that $\nabla A(\nabla A^*(\eta)) = \eta$, therefore differentiating this expression with respect to $\eta$ yields

$$
\nabla^2 A^*(\eta) \nabla^2 A(\nabla A^*(\eta)) = \mathrm{Id}.
$$

Now take $\eta = \nabla A(\theta_\pi)$ which belongs to int dom $A^*$ since $\theta_\pi \in \text{int dom } A$. We obtain that $\nabla^2 A^*(\nabla A(\theta_\pi))\nabla^2 A(\theta_\pi) = \text{Id}$. This shows that the optimality condition of Eq. (33) is equivalent to having

$$
\theta - \theta' = -\frac{1}{\alpha} \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta').
$$

This shows that the right-hand side of Eq. (32) can be minorized to obtain:

$$
\begin{align*}
f_{\pi}^{(\alpha)}(\theta) \ge f_{\pi}^{(\alpha)}(\theta') & - \frac{1}{\alpha} \langle \nabla f_{\pi}^{(\alpha)}(\theta'), \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta') \rangle \\
& + \frac{1}{2\alpha} \langle \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta'), \nabla^2 A(\theta_{\pi}) \nabla^2 A^*(\nabla A(\theta_{\pi})) \nabla f_{\pi}^{(\alpha)}(\theta') \rangle + o(v^2) \\
& = f_{\pi}^{(\alpha)}(\theta') - \frac{1}{2\alpha} \|\nabla f_{\pi}^{(\alpha)}(\theta')\|_{\nabla^2 A^*(\nabla A(\theta_{\pi}))}^2 + o(v^2).
\end{align*}
$$

This is true in particular for $\theta = \theta_{\pi}$, which yields the result.

<figure><img src="image_6.png" /></figure>

B Convergence analysis of Algorithm 1

In order to prove Propositions 14 and 15, we start with a sufficient decrease lemma that reads as follows.

**Lemma 1.** Under Assumptions 1, 2, and 3, for $\tau > 0$ and $\alpha \in (0, 1]$, we have that for every $\theta \in \text{int } \Theta$,

$$
\tau \left( F_{\pi}^{(\alpha)}(T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta)) - F_{\pi}^{(\alpha)}(\theta) \right) \le -d_{A}(\theta, T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta)) + (\tau - 1)d_{A}(T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta), \theta). \quad (34)
$$

In the particular case where $\alpha = 1$, we further have

$$
\tau \left( F_{\pi}^{(1)}(T_{\tau F_{\pi}^{(1)}}^{A}(\theta)) - F_{\pi}^{(1)}(\theta') \right) \le (1-\tau)d_{A}(\theta', \theta) - (1-\tau)d_{A}(T_{\tau F_{\pi}^{(1)}}^{A}(\theta), \theta) \\
\qquad - d_{A}(\theta', T_{\tau F_{\pi}^{(1)}}^{A}(\theta)), \quad \forall \theta' \in \text{int } \Theta. \quad (35)
$$
Proof. Using [74, Lemma 4.1], which still holds in our finite-dimensional Hilbert setting, we get that

$$ \tau(F_{\pi}^{(\alpha)}(T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta)) - F_{\pi}^{(\alpha)}(\theta')) \le d_{A}(\theta', \theta) - (1 - \tau)d_{A}(T_{\tau F_{\pi}^{(1)}}^{A}(\theta), \theta) \\ - d_{A}(\theta', T_{\tau F_{\pi}^{(1)}}^{A}(\theta)) - \tau d_{f_{\pi}^{(\alpha)}}(\theta', \theta), \forall \theta' \in \text{int } \Theta, $$

where $d_{f_{\pi}^{(\alpha)}}(\theta', \theta) = f_{\pi}^{(\alpha)}(\theta') - f_{\pi}^{(\alpha)}(\theta) - \langle \nabla f_{\pi}^{(\alpha)}(\theta), \theta' - \theta \rangle$.

Equation (34) comes by evaluating the above at $\theta' = \theta$. To get Eq. (35), the strong convexity of $f_{\pi}^{(1)}$ relatively to A yields

$$ d_{f_{\pi}^{(1)}}(\theta', \theta) \ge d_{A}(\theta', \theta), \forall \theta', \theta \in \text{int } \Theta, $$

showing the result.

<figure><img src="image_9.png" /></figure>

We also give a sequential consistency lemma, that links the Bregman divergence $d_A$ with the Euclidean distance.

**Lemma 2.** Consider two sequences $\{\theta_k\}_{k \in \mathbb{N}}$ and $\{\theta'_k\}_{k \in \mathbb{N}}$ and assume that there exists a compact set $C \subset \text{int } \Theta$ such that $\theta_k, \theta'_k \in C$ for every $k \in \mathbb{N}$. In this case, if $d_A(\theta_k, \theta'_k) \xrightarrow{k \to +\infty} 0$, then $\|\theta_k - \theta'_k\| \xrightarrow{k \to +\infty} 0$.

Proof. We introduce the convex hull of $C$, denoted by $\text{conv } C$ which is the intersection of every convex set containing $C$. Therefore $\text{conv } C \subset \text{int } \Theta$. Since we are in finite dimension, we also have that $\text{conv } C$ is compact. Thus, $\text{conv } C$ is a convex compact included in $\text{int } \Theta$.

A is proper, strictly convex, and continuous on $\text{conv } C \subset \text{int } \Theta$, therefore, A is uniformly convex, following the definition of [10, Definition 10.5] on $\text{conv } C$ [10, Proposition 10.15]. This means that there exists an increasing function $\psi: \mathbb{R}_+ \to [0, +\infty]$ that vanishes only at 0, such that for every $\theta, \theta' \in \text{conv } C$,

$$ \psi(\|\theta - \theta'\|) \le \frac{1}{2} A(\theta) + \frac{1}{2} A(\theta') - A\left(\frac{1}{2}\theta + \frac{1}{2}\theta'\right). $$

Because A is convex on $\text{conv } C$, we prove following the proof of [19, Eq. (1.32)] that for every $\theta, \theta' \in \text{conv } C$,

$$ d_A(\theta, \theta') \ge \psi(\|\theta - \theta'\|). $$

Suppose now by contradiction that $d_A(\theta_k, \theta'_k) \xrightarrow{k \to +\infty} 0$ while there exists some $\epsilon > 0$ such that $\|\theta_k - \theta'_k\| \ge \epsilon$ for every $k \in \mathbb{N}$. Then we have that

$$ d_A(\theta_k, \theta'_k) \ge \psi(\epsilon) > 0, $$

which is a contradiction, hence showing the result.

<figure><img src="image_10.png" /></figure>

## B.1 Proof of Proposition 14

*Proof of Proposition 14.* We first give an intermediate result that will be used several times. Using Lemma 1, we can get for any $k \in \mathbb{N}$ that

$$ F_{\pi}^{(\alpha)}(\theta_{k+1}) - F_{\pi}^{(\alpha)}(\theta_k) \le -\frac{1}{\tau_{k+1}} d_A(\theta_k, \theta_{k+1}) - \left(\frac{1}{\tau_{k+1}} - 1\right) d_A(\theta_{k+1}, \theta_k). \quad (36) $$

(i) Due to the non-negativity of the Bregman divergences and the hypothesis on the step-size $\tau_{k+1}$, the right-hand side of Eq. (36) is non-negative, yielding the decrease property $F_{\pi}^{(\alpha)}(\theta_{k+1}) \le F_{\pi}^{(\alpha)}(\theta_k)$.
Proof. Using [74, Lemma 4.1], which still holds in our finite-dimensional Hilbert setting, we get that

$$ \tau(F_{\pi}^{(\alpha)}(T_{\tau F_{\pi}^{(\alpha)}}^{A}(\theta)) - F_{\pi}^{(\alpha)}(\theta')) \le d_{A}(\theta', \theta) - (1-\tau)d_{A}(T_{\tau F_{\pi}^{(1)}}^{A}(\theta), \theta) \\ - d_{A}(\theta', T_{\tau F_{\pi}^{(1)}}^{A}(\theta)) - \tau d_{f_{\pi}^{(\alpha)}}(\theta', \theta), \forall \theta' \in \text{int } \Theta, $$

where $d_{f_{\pi}^{(\alpha)}}(\theta', \theta) = f_{\pi}^{(\alpha)}(\theta') - f_{\pi}^{(\alpha)}(\theta) - \langle \nabla f_{\pi}^{(\alpha)}(\theta), \theta' - \theta \rangle$.

Equation (34) comes by evaluating the above at $\theta' = \theta$. To get Eq. (35), the strong convexity of $f_{\pi}^{(1)}$ relatively to A yields

$$ d_{f_{\pi}^{(1)}}(\theta', \theta) \ge d_{A}(\theta', \theta), \forall \theta', \theta \in \text{int } \Theta, $$

showing the result.

<figure><img src="image_7.png" /></figure>

We also give a sequential consistency lemma, that links the Bregman divergence $d_A$ with the Euclidean distance.

**Lemma 2.** Consider two sequences $\{\theta_k\}_{k \in \mathbb{N}}$ and $\{\theta'_k\}_{k \in \mathbb{N}}$ and assume that there exists a compact set $C \subset \text{int } \Theta$ such that $\theta_k, \theta'_k \in C$ for every $k \in \mathbb{N}$. In this case, if $d_A(\theta_k, \theta'_k) \xrightarrow{k \to +\infty} 0$, then $\|\theta_k - \theta'_k\| \xrightarrow{k \to +\infty} 0$.

Proof. We introduce the convex hull of $C$, denoted by $\text{conv } C$ which is the intersection of every convex set containing $C$. Therefore $\text{conv } C \subset \text{int } \Theta$. Since we are in finite dimension, we also have that $\text{conv } C$ is compact. Thus, $\text{conv } C$ is a convex compact included in $\text{int } \Theta$.

A is proper, strictly convex, and continuous on $\text{conv } C \subset \text{int } \Theta$, therefore, A is uniformly convex, following the definition of [10, Definition 10.5] on $\text{conv } C$ [10, Proposition 10.15]. This means that there exists an increasing function $\psi: \mathbb{R}_+ \to [0, +\infty]$ that vanishes only at 0, such that for every $\theta, \theta' \in \text{conv } C$,

$$ \psi(\|\theta - \theta'\|) \le \frac{1}{2}A(\theta) + \frac{1}{2}A(\theta') - A\left(\frac{1}{2}\theta + \frac{1}{2}\theta'\right). $$

Because A is convex on $\text{conv } C$, we prove following the proof of [19, Eq. (1.32)] that for every $\theta, \theta' \in \text{conv } C$,

$$ d_{A}(\theta, \theta') \ge \psi(\|\theta - \theta'\|). $$

Suppose now by contradiction that $d_A(\theta_k, \theta'_k) \xrightarrow{k \to +\infty} 0$ while there exists some $\epsilon > 0$ such that $\|\theta_k - \theta'_k\| \ge \epsilon$ for every $k \in \mathbb{N}$. Then we have that

$$ d_{A}(\theta_{k}, \theta'_{k}) \geq \psi(\epsilon) > 0, $$

which is a contradiction, hence showing the result.

<figure><img src="image_8.png" /></figure>

## B.1 Proof of Proposition 14

*Proof of Proposition 14.* We first give an intermediate result that will be used several times. Using Lemma 1, we can get for any $k \in \mathbb{N}$ that

$$ F_{\pi}^{(\alpha)}(\theta_{k+1}) - F_{\pi}^{(\alpha)}(\theta_k) \le -\frac{1}{\tau_{k+1}} d_A(\theta_k, \theta_{k+1}) - \left(\frac{1}{\tau_{k+1}} - 1\right) d_A(\theta_{k+1}, \theta_k). \quad (36) $$

(i) Due to the non-negativity of the Bregman divergences and the hypothesis on the step-size $\tau_{k+1}$, the right-hand side of Eq. (36) is non-negative, yielding the decrease property $F_{\pi}^{(\alpha)}(\theta_{k+1}) \le F_{\pi}^{(\alpha)}(\theta_k)$.
(ii) If $F_{\pi}^{(\alpha)}(\theta_{K+1}) = F_{\pi}^{(\alpha)}(\theta_{K})$, then, using Lemma 1 and $\tau_{K+1} \le 1$, $d_{A}(\theta_{K}, \theta_{K+1}) \le 0$. By Proposition 3, this shows that $\theta_{K+1} = \theta_{K}$. Since $\theta_{K+1} = T^{A}_{\tau_{K+1}F_{\pi}^{(\alpha)}(\theta_{K})}$, $\theta_{K}$ is a fixed point of $T^{A}_{\tau_{K+1}F_{\pi}^{(\alpha)}(\theta_{K})}$. From Proposition 7, it is a stationary point of $F_{\pi}^{(\alpha)}$.

(iii) By summing Eq. (36) over $k \in [0, K]$, one gets

$$ \sum_{k=0}^{K} \left( \frac{1}{\tau_{k+1}} d_{A}(\theta_{k}, \theta_{k+1}) + \left( \frac{1}{\tau_{k+1}} - 1 \right) d_{A}(\theta_{k+1}, \theta_{k}) \right) \leq F_{\pi}^{(\alpha)}(\theta_{0}) - F_{\pi}^{(\alpha)}(\theta_{K+1}). $$

Because of the hypothesis on the step-sizes, one has for any $k \in \mathbb{N}$ that $1 \le \frac{1}{\tau_{k+1}}$ and $0 \le \frac{1}{\tau_{k+1}} - 1$, giving the inequality

$$ \sum_{k=0}^{K} d_{A}(\theta_{k}, \theta_{k+1}) \leq F_{\pi}^{(\alpha)}(\theta_{0}) - \vartheta_{\pi}^{(\alpha)}. $$

Since the right-hand side of the above is uniform in $K$, this implies the result.

(iv) Suppose that $K$ is the first iterate such that $d_{A}(\theta_{K}, \theta_{K+1}) \le \varepsilon$. Then for any $k \in [0, K - 1]$, $d_{A}(\theta_{k}, \theta_{k+1}) > \varepsilon$. Consider Eq. (B.1) where the sum goes only from $k = 0$ to $k = K - 1$, then one can show the desired result with

$$ \varepsilon K \le \sum_{k=0}^{K-1} d_{A}(\theta_{k}, \theta_{k+1}) \le F_{\pi}^{(\alpha)}(\theta_{0}) - \vartheta_{\pi}^{(\alpha)}. $$

(v) This proof relies on two notions of subdifferentials: the limiting subdifferential $\partial_{L}$ [70, Chapter 6] and the Fréchet subdifferential $\partial_{F}$ [70, Chapter 4]. Our working space $\mathcal{H}$ is a finite-dimensional Hilbert space, which is included in the setting of [70].

Set $k \in \mathbb{N}$. Under Assumptions 1, 2, and 3, since $\theta_{0} \in \text{int} \Theta$, Proposition 8 applies and thus $\theta_{k+1} = T^{A}_{\tau_{k+1}F_{\pi}^{(\alpha)}(\theta_{k})}$. This implies that there exists $g_{k+1} \in \partial_{r}(\theta_{k+1})$ such that

$$ \frac{1}{\tau_{k+1}}(\nabla A(\theta_{k+1}) - \nabla A(\theta_{k})) + \nabla f_{\pi}^{(\alpha)}(\theta_{k}) + g_{k+1} = 0. \quad (37) $$

According to [70, Corollary 4.35],

$$ \nabla f_{\pi}^{(\alpha)}(\theta_{k+1}) + g_{k+1} \in \partial_{F} F_{\pi}^{(\alpha)}(\theta_{k+1}). \quad (38) $$

Using Eq. (37) and the assumptions on $\tau_{k+1}$,

$$ \| \nabla f_{\pi}^{(\alpha)}(\theta_{k+1}) + g_{k+1} \| \le \| \nabla f_{\pi}^{(\alpha)}(\theta_{k+1}) - \nabla f_{\pi}^{(\alpha)}(\theta_{k}) \| + \frac{1}{\epsilon} \| \nabla A(\theta_{k+1}) - \nabla A(\theta_{k}) \|. $$

The additional hypothesis introduced in (iv) ensures that both $\theta_{k+1}$ and $\theta_{k}$ belong to $C$, a compact set included in $\text{int} \Theta$. Since $\nabla^{2} f_{\pi}^{(\alpha)}$ is continuous on $C$ (by Proposition 6) and $C$ is bounded, $\nabla f_{\pi}^{(\alpha)}$ is Lipschitz on $C$. The same reasoning applies for $\nabla A$. This shows that there exists a scalar $s > 0$ such that, for every $k \in \mathbb{N}$, there exists $\varrho_{k+1} \in \partial_{F} F_{\pi}^{(\alpha)}(\theta_{k+1})$ satisfying

$$ \| \varrho_{k+1} \| \le s \| \theta_{k+1} - \theta_{k} \|. \quad (39) $$

Now, we deduce from (iii) that $d_{A}(\theta_{k+1}, \theta_{k}) \xrightarrow{k \to +\infty} 0$. Using Lemma 2, this yields $\| \theta_{k+1} - \theta_{k} \| \xrightarrow{k \to +\infty} 0$, showing that the sequence $\{\varrho_{k}\}_{k \in \mathbb{N}}$ is such that

$$ \varrho_{k} \in \partial_{F} F_{\pi}^{(\alpha)}(\theta_{k}), \forall k \in \mathbb{N}, \text{ and } \varrho_{k} \xrightarrow{k \to +\infty} 0. \quad (40) $$
(ii) If $F_{\pi}^{(\alpha)}(\theta_{K+1}) = F_{\pi}^{(\alpha)}(\theta_{K})$, then, using Lemma 1 and $\tau_{K+1} \le 1$, $d_{A}(\theta_{K}, \theta_{K+1}) \le 0$. By Proposition 3, this shows that $\theta_{K+1} = \theta_{K}$. Since $\theta_{K+1} = T^{A}_{\tau_{K+1}F_{\pi}^{(\alpha)}(\theta_{K})}$, $\theta_{K}$ is a fixed point of $T^{A}_{\tau_{K+1}F_{\pi}^{(\alpha)}(\theta_{K})}$. From Proposition 7, it is a stationary point of $F_{\pi}^{(\alpha)}$.

(iii) By summing Eq. (36) over $k \in [0, K]$, one gets

$$ \sum_{k=0}^{K} \left( \frac{1}{\tau_{k+1}} d_{A}(\theta_{k}, \theta_{k+1}) + \left( \frac{1}{\tau_{k+1}} - 1 \right) d_{A}(\theta_{k+1}, \theta_{k}) \right) \leq F_{\pi}^{(\alpha)}(\theta_{0}) - F_{\pi}^{(\alpha)}(\theta_{K+1}). $$

Because of the hypothesis on the step-sizes, one has for any $k \in \mathbb{N}$ that $1 \le \frac{1}{\tau_{k+1}}$ and $0 \le \frac{1}{\tau_{k+1}} - 1$, giving the inequality

$$ \sum_{k=0}^{K} d_{A}(\theta_{k}, \theta_{k+1}) \leq F_{\pi}^{(\alpha)}(\theta_{0}) - \vartheta_{\pi}^{(\alpha)}. $$

Since the right-hand side of the above is uniform in $K$, this implies the result.

(iv) Suppose that $K$ is the first iterate such that $d_{A}(\theta_{K}, \theta_{K+1}) \le \varepsilon$. Then for any $k \in [0, K - 1]$, $d_{A}(\theta_{k}, \theta_{k+1}) > \varepsilon$. Consider Eq. (B.1) where the sum goes only from $k = 0$ to $k = K - 1$, then one can show the desired result with

$$ \varepsilon K \le \sum_{k=0}^{K-1} d_{A}(\theta_{k}, \theta_{k+1}) \le F_{\pi}^{(\alpha)}(\theta_{0}) - \vartheta_{\pi}^{(\alpha)}. $$

(v) This proof relies on two notions of subdifferentials: the limiting subdifferential $\partial_{L}$ [70, Chapter 6] and the Fréchet subdifferential $\partial_{F}$ [70, Chapter 4]. Our working space $\mathcal{H}$ is a finite-dimensional Hilbert space, which is included in the setting of [70].

Set $k \in \mathbb{N}$. Under Assumptions 1, 2, and 3, since $\theta_{0} \in \text{int} \Theta$, Proposition 8 applies and thus $\theta_{k+1} = T^{A}_{\tau_{k+1}F_{\pi}^{(\alpha)}(\theta_{k})}$. This implies that there exists $g_{k+1} \in \partial_{r}(\theta_{k+1})$ such that

$$ \frac{1}{\tau_{k+1}}(\nabla A(\theta_{k+1}) - \nabla A(\theta_{k})) + \nabla f_{\pi}^{(\alpha)}(\theta_{k}) + g_{k+1} = 0. \quad (37) $$

According to [70, Corollary 4.35],

$$ \nabla f_{\pi}^{(\alpha)}(\theta_{k+1}) + g_{k+1} \in \partial_{F} F_{\pi}^{(\alpha)}(\theta_{k+1}). \quad (38) $$

Using Eq. (37) and the assumptions on $\tau_{k+1}$,

$$ \| \nabla f_{\pi}^{(\alpha)}(\theta_{k+1}) + g_{k+1} \| \le \| \nabla f_{\pi}^{(\alpha)}(\theta_{k+1}) - \nabla f_{\pi}^{(\alpha)}(\theta_{k}) \| + \frac{1}{\epsilon} \| \nabla A(\theta_{k+1}) - \nabla A(\theta_{k}) \|. $$

The additional hypothesis introduced in (iv) ensures that both $\theta_{k+1}$ and $\theta_{k}$ belong to $C$, a compact set included in $\text{int} \Theta$. Since $\nabla^{2} f_{\pi}^{(\alpha)}$ is continuous on $C$ (by Proposition 6) and $C$ is bounded, $\nabla f_{\pi}^{(\alpha)}$ is Lipschitz on $C$. The same reasoning applies for $\nabla A$. This shows that there exists a scalar $s > 0$ such that, for every $k \in \mathbb{N}$, there exists $\varrho_{k+1} \in \partial_{F} F_{\pi}^{(\alpha)}(\theta_{k+1})$ satisfying

$$ \|\varrho_{k+1}\| \le s \|\theta_{k+1} - \theta_{k}\|. \quad (39) $$

Now, we deduce from (iii) that $d_{A}(\theta_{k+1}, \theta_{k}) \underset{k \to +\infty}{\longrightarrow} 0$. Using Lemma 2, this yields $\|\theta_{k+1} - \theta_{k}\| \underset{k \to +\infty}{\longrightarrow} 0$, showing that the sequence $\{\varrho_{k}\}_{k \in \mathbb{N}}$ is such that

$$ \varrho_{k} \in \partial_{F} F_{\pi}^{(\alpha)}(\theta_{k}), \forall k \in \mathbb{N}, \text{ and } \varrho_{k} \underset{k \to +\infty}{\longrightarrow} 0. \quad (40) $$
On the other hand, the sequence {$\theta_k$}$_{k \in \mathbb{N}}$ is contained in the compact set C by assumption. Hence, there exists $\theta_{\lim} \in C$, and a strictly increasing function $\varphi: \mathbb{N} \to \mathbb{N}$ such that $\theta_{\varphi(k)} \xrightarrow{k \to +\infty} \theta_{\lim}$. The regularizing term r is continuous on C as assumed in (iv), implying

$$
\left\{
\begin{array}{l}
\theta_{\varphi(k)} \xrightarrow{k \to +\infty} \theta_{\lim}, \\
F_{\pi}^{(\alpha)}(\theta_{\varphi(k)}) \xrightarrow{k \to +\infty} F_{\pi}^{(\alpha)}(\theta_{\lim}), \\
\varrho_{\varphi(k)} \in \partial_F F_{\pi}^{(\alpha)}(\theta_{\varphi(k)}), \quad \varrho_{\varphi(k)} \xrightarrow{k \to +\infty} 0.
\end{array}
\right.
$$

By definition of the limiting subdifferential $\partial_L F_{\pi}^{(\alpha)}$ [70, Definition 6.1], this shows that $0 \in \partial_L F_{\pi}^{(\alpha)}(\theta_{\lim})$.
Hence $\theta_{\lim}$ is a stationary point of $F_{\pi}^{(\alpha)}$ concluding the proof. $\square$

## B.2 Proof of Proposition 15

*Proof of Proposition 15.* We first give an inequality to prove (i)-(ii). Consider iteration $k$ of Algorithm 1, and evaluate Eq. (35) from Lemma 1 at $\theta' = \theta_*$, yielding

$$
\tau_{k+1} \left( F_{\pi}^{(\alpha)}(\theta_{k+1}) - F_{\pi}^{(\alpha)}(\theta_{*}) \right) \le (1 - \tau_{k+1}) d_A(\theta_{*}, \theta_k) \\ \qquad - (1 - \tau_{k+1}) d_A(\theta_{k+1}, \theta_k) - d_A(\theta_{*}, \theta_{k+1}). \quad (41)
$$

(i) Since $\tau_{k+1} \in [\epsilon, 1]$, $F_{\pi}^{(1)}(\theta_{k+1}) \ge F_{\pi}^{(1)}(\theta_*)$, and $d_A$ takes non-negative values (from Proposition 3), Eq. (41) gives

$$
d_A(\theta_*, \theta_{k+1}) \le (1 - \tau_{k+1}) d_A(\theta_*, \theta_k), \quad (42)
$$

from which we deduce the results since $\tau_{k+1} \in [\epsilon, 1]$.

(ii) Since $\tau_{k+1} \in [\epsilon, 1]$ and $d_A$ takes non-negative values, we get from Eq. (41) that

$$
\tau_{k+1} \left( F_{\pi}^{(1)}(\theta_{k+1}) - F_{\pi}^{(1)}(\theta_{*}) \right) \le (1 - \tau_{k+1}) d_A(\theta_{*}, \theta_k).
$$

With Eq. (42) and the condition on $\tau_{k+1}$, we obtain

$$
(F_{\pi}^{(1)}(\theta_{k+1}) - F_{\pi}^{(1)}(\theta_{*})) \le \frac{1}{\epsilon} d_{A}(\theta_{*}, \theta_{k+1}),
$$

from which we conclude using point (i) and Proposition 10.

(iii) Using Proposition 14 (i), we obtain that for every $k \in \mathbb{N}$, $F_{\pi}^{(1)}(\theta_k) \le F_{\pi}^{(1)}(\theta_0)$, meaning that the sequence {$\theta_k$}$_{k \in \mathbb{N}}$ is contained in a sub-level set of $F_{\pi}^{(1)}$. $F_{\pi}^{(1)}$ is coercive under our assumptions (see the proof of Proposition 10), and it is lower semicontinuous from Proposition 4, so its sub-level sets are compact. This means that we can extract converging subsequences from {$\theta_k$}$_{k \in \mathbb{N}}$.

Consider now such a subsequence {$\theta_{\varphi(k)}$}$_{k \in \mathbb{N}}$, with $\theta_{\varphi(k)} \xrightarrow{k \to +\infty} \theta_{\lim}$. $F_{\pi}^{(1)}$ is lower semicontinuous, so

$$
\liminf F_{\pi}^{(1)}(\theta_{\varphi(k)}) \geq F_{\pi}^{(1)}(\theta_{\lim}).
$$

However, because of (ii), $\liminf F_{\pi}^{(1)}(\theta_{\varphi(k)}) = F_{\pi}^{(1)}(\theta_*)$, so we obtain that $F_{\pi}^{(1)}(\theta_{\lim}) = F_{\pi}^{(1)}(\theta_*)$. Using Proposition 10, this shows that $\theta_{\lim} = \theta_*$.

We have shown that {$\theta_k$}$_{k \in \mathbb{N}}$ is contained in a compact set and that each of its converging subsequences converges to $\theta_*$, which implies the result. $\square$
On the other hand, the sequence {$\theta_k$}$_{k \in \mathbb{N}}$ is contained in the compact set C by assumption. Hence, there exists $\theta_{\lim} \in C$, and a strictly increasing function $\varphi: \mathbb{N} \to \mathbb{N}$ such that $\theta_{\varphi(k)} \xrightarrow{k \to +\infty} \theta_{\lim}$. The regularizing term r is continuous on C as assumed in (iv), implying

$$
\left\{
\begin{array}{l}
\theta_{\varphi(k)} \xrightarrow{k \to +\infty} \theta_{\lim}, \\
F_{\pi}^{(\alpha)}(\theta_{\varphi(k)}) \xrightarrow{k \to +\infty} F_{\pi}^{(\alpha)}(\theta_{\lim}), \\
\varrho_{\varphi(k)} \in \partial_F F_{\pi}^{(\alpha)}(\theta_{\varphi(k)}), \quad \varrho_{\varphi(k)} \xrightarrow{k \to +\infty} 0.
\end{array}
\right.
$$

By definition of the limiting subdifferential $\partial_L F_{\pi}^{(\alpha)}$ [70, Definition 6.1], this shows that $0 \in \partial_L F_{\pi}^{(\alpha)}(\theta_{\lim})$.
Hence $\theta_{\lim}$ is a stationary point of $F_{\pi}^{(\alpha)}$ concluding the proof. $\square$

## B.2 Proof of Proposition 15

*Proof of Proposition 15.* We first give an inequality to prove (i)-(ii). Consider iteration $k$ of Algorithm 1, and evaluate Eq. (35) from Lemma 1 at $\theta' = \theta_*$, yielding

$$
\tau_{k+1} \left( F_{\pi}^{(\alpha)}(\theta_{k+1}) - F_{\pi}^{(\alpha)}(\theta_{*}) \right) \le (1 - \tau_{k+1}) d_A(\theta_{*}, \theta_{k}) \\ - (1 - \tau_{k+1}) d_A(\theta_{k+1}, \theta_{k}) - d_A(\theta_{*}, \theta_{k+1}). \quad (41)
$$

(i) Since $\tau_{k+1} \in [\epsilon, 1]$, $F_{\pi}^{(1)}(\theta_{k+1}) \ge F_{\pi}^{(1)}(\theta_*)$, and $d_A$ takes non-negative values (from Proposition 3), Eq. (41) gives

$$
d_A(\theta_*, \theta_{k+1}) \le (1 - \tau_{k+1}) d_A(\theta_*, \theta_k), \quad (42)
$$

from which we deduce the results since $\tau_{k+1} \in [\epsilon, 1]$.

(ii) Since $\tau_{k+1} \in [\epsilon, 1]$ and $d_A$ takes non-negative values, we get from Eq. (41) that

$$
\tau_{k+1} \left( F_{\pi}^{(1)}(\theta_{k+1}) - F_{\pi}^{(1)}(\theta_{*}) \right) \le (1 - \tau_{k+1}) d_A(\theta_{*}, \theta_{k}).
$$

With Eq. (42) and the condition on $\tau_{k+1}$, we obtain

$$
(F_{\pi}^{(1)}(\theta_{k+1}) - F_{\pi}^{(1)}(\theta_{*})) \le \frac{1}{\epsilon} d_{A}(\theta_{*}, \theta_{k+1}),
$$

from which we conclude using point (i) and Proposition 10.

(iii) Using Proposition 14 (i), we obtain that for every $k \in \mathbb{N}$, $F_{\pi}^{(1)}(\theta_k) \le F_{\pi}^{(1)}(\theta_0)$, meaning that the sequence {$\theta_k$}$_{k \in \mathbb{N}}$ is contained in a sub-level set of $F_{\pi}^{(1)}$. $F_{\pi}^{(1)}$ is coercive under our assumptions (see the proof of Proposition 10), and it is lower semicontinuous from Proposition 4, so its sub-level sets are compact. This means that we can extract converging subsequences from {$\theta_k$}$_{k \in \mathbb{N}}$.

Consider now such a subsequence {$\theta_{\varphi(k)}$}$_{k \in \mathbb{N}}$, with $\theta_{\varphi(k)} \xrightarrow{k \to +\infty} \theta_{\lim}$. $F_{\pi}^{(1)}$ is lower semicontinuous, so

$$
\liminf F_{\pi}^{(1)}(\theta_{\varphi(k)}) \geq F_{\pi}^{(1)}(\theta_{\lim}).
$$

However, because of (ii), $\liminf F_{\pi}^{(1)}(\theta_{\varphi(k)}) = F_{\pi}^{(1)}(\theta_*)$, so we obtain that $F_{\pi}^{(1)}(\theta_{\lim}) = F_{\pi}^{(1)}(\theta_*)$. Using Proposition 10, this shows that $\theta_{\lim} = \theta_*$.

We have shown that {$\theta_k$}$_{k \in \mathbb{N}}$ is contained in a compact set and that each of its converging subsequences converges to $\theta_*$, which implies the result. $\square$
## B.3 Proof of Proposition 16

*Proof of Proposition 16.* We first prove that we can apply the results of Proposition 14. Assumption 2 holds because of the result of Proposition 12. Since $r \equiv 0$, Assumption 3 holds. Therefore, all the hypotheses of Proposition 14 hold, showing the result.

(iii) Because of Proposition 14, every converging subsequence of $\{\theta_k\}_{k \in \mathbb{N}}$ converges to $S_{\pi}^{(\alpha)}$. However, $S_{\pi}^{(\alpha)} = \{\theta_{\pi}\}$ in our case (see Proposition 12). Therefore, all the converging subsequences of $\{\theta_k\}_{k \in \mathbb{N}}$ converge to $\theta_{\pi}$, showing that the sequence of iterates converges to $\theta_{\pi}$.

(ii) Due to the smoothness of $f_{\pi}^{(\alpha)}$ relatively to A shown in Proposition 9, we have for any $k \in \mathbb{N}$ that

$$f_{\pi}^{(\alpha)}(\theta_{k+1}) \le f_{\pi}^{(\alpha)}(\theta_k) + \langle \nabla f_{\pi}^{(\alpha)}(\theta_k), \theta_{k+1} - \theta_k \rangle + d_A(\theta_{k+1}, \theta_k). \quad (43)$$

Since $r \equiv 0$, we have the relation $\nabla A(\theta_{k+1}) = \nabla A(\theta_k) - \tau_{k+1} \nabla f_{\pi}^{(\alpha)}(\theta_k)$, therefore, Eq. (43) reads

$$f_{\pi}^{(\alpha)}(\theta_{k+1}) \le f_{\pi}^{(\alpha)}(\theta_k) - \frac{1}{\tau_{k+1}} \langle \nabla A(\theta_{k+1}) - \nabla A(\theta_k), \theta_{k+1} - \theta_k \rangle + d_A(\theta_{k+1}, \theta_k). \quad (44)$$

Since $\tau_{k+1} \le 1$ and $\langle \nabla A(\theta_{k+1}) - \nabla A(\theta_k), \theta_{k+1} - \theta_k \rangle = d_A(\theta_{k+1}, \theta_k) + d_A(\theta_k, \theta_{k+1})$, we get

$$\begin{aligned} f_{\pi}^{(\alpha)}(\theta_{k+1}) &\le f_{\pi}^{(\alpha)}(\theta_k) - \frac{1}{\tau_{k+1}} (d_A(\theta_k, \theta_{k+1}) + d_A(\theta_{k+1}, \theta_k)) + \frac{1}{\tau_{k+1}} d_A(\theta_{k+1}, \theta_k) \\ &= f_{\pi}^{(\alpha)}(\theta_k) - \frac{1}{\tau_{k+1}} d_{A^*}(\nabla A(\theta_{k+1}), \nabla A(\theta_k)). \end{aligned}$$

Now, we consider some $v > 0$. Since $\theta_k \to \theta_{\pi}$, there exists $K \in \mathbb{N}$ such that for any $k \ge K$, $\theta_k \in B(\theta_{\pi}, v)$ and $\nabla A(\theta_k) \in B(\nabla A(\theta_{\pi}), v)$ (recall that $\nabla A$ is continuous on int $\Theta$). Thus, we can write following the same steps as in the proof of Proposition 13 that for any $k \ge K$,

$$d_{A^*}(\nabla A(\theta_{k+1}), \nabla A(\theta_k)) = \frac{1}{2} \|\nabla A(\theta_{k+1}) - \nabla A(\theta_k)\|_{\nabla^2 A^*(\nabla A(\theta_{\pi}))}^2 + o(v^2).$$

From there, we obtain that

$$\begin{aligned} d_{A^*}(\nabla A(\theta_{k+1}), \nabla A(\theta_k)) &= \frac{\tau_{k+1}^2}{2} \|\nabla f_{\pi}^{(\alpha)}(\theta_k)\|_{\nabla^2 A^*(\nabla A(\theta_{\pi}))}^2 + o(v^2) \\ &\ge \tau_{k+1}^2 \alpha f_{\pi}^{(\alpha)}(\theta_k) + o(v^2), \end{aligned}$$

where the last inequality comes from Proposition 13. With our previous point, this yields

$$\begin{aligned} f_{\pi}^{(\alpha)}(\theta_{k+1}) &\le (1 - \alpha \tau_{k+1}) f_{\pi}^{(\alpha)}(\theta_k) + o(v^2) \\ &< (1 - \alpha \delta \epsilon) f_{\pi}^{(\alpha)}(\theta_k) + o(v^2) \end{aligned}$$

for any constant $\delta \in (0, 1)$. By choosing $\delta$ or $v$ small enough, we finally obtain

$$f_{\pi}^{(\alpha)}(\theta_{k+1}) \le (1 - \alpha \delta \epsilon) f_{\pi}^{(\alpha)}(\theta_k).$$

This means that for any $k \ge K$, we have $f_{\pi}^{(\alpha)}(\theta_k) \le (1 - \alpha \delta \epsilon)^{k-K} f_{\pi}^{(\alpha)}(\theta_K)$. Since the sequence $\{f_{\pi}^{(\alpha)}(\theta_k)\}_{k \in \mathbb{N}}$ is decreasing, $f_{\pi}^{(\alpha)}(\theta_K) \le f_{\pi}^{(\alpha)}(\theta_0)$, which shows the result with $C = (1 - \alpha \delta \epsilon)^{-K}$. $\square$

<figure><img src="image_11.png" /></figure>
## B.3 Proof of Proposition 16

*Proof of Proposition 16.* We first prove that we can apply the results of Proposition 14. Assumption 2 holds because of the result of Proposition 12. Since $r \equiv 0$, Assumption 3 holds. Therefore, all the hypotheses of Proposition 14 hold, showing the result.

(iii) Because of Proposition 14, every converging subsequence of $\{\theta_k\}_{k \in \mathbb{N}}$ converges to $S_{\pi}^{(\alpha)}$. However, $S_{\pi}^{(\alpha)} = \{\theta_{\pi}\}$ in our case (see Proposition 12). Therefore, all the converging subsequences of $\{\theta_k\}_{k \in \mathbb{N}}$ converge to $\theta_{\pi}$, showing that the sequence of iterates converges to $\theta_{\pi}$.

(ii) Due to the smoothness of $f_{\pi}^{(\alpha)}$ relatively to A shown in Proposition 9, we have for any $k \in \mathbb{N}$ that

$$f_{\pi}^{(\alpha)}(\theta_{k+1}) \le f_{\pi}^{(\alpha)}(\theta_k) + \langle \nabla f_{\pi}^{(\alpha)}(\theta_k), \theta_{k+1} - \theta_k \rangle + d_A(\theta_{k+1}, \theta_k). \quad (43)$$

Since $r \equiv 0$, we have the relation $\nabla A(\theta_{k+1}) = \nabla A(\theta_k) - \tau_{k+1} \nabla f_{\pi}^{(\alpha)}(\theta_k)$, therefore, Eq. (43) reads

$$f_{\pi}^{(\alpha)}(\theta_{k+1}) \le f_{\pi}^{(\alpha)}(\theta_k) - \frac{1}{\tau_{k+1}} \langle \nabla A(\theta_{k+1}) - \nabla A(\theta_k), \theta_{k+1} - \theta_k \rangle + d_A(\theta_{k+1}, \theta_k). \quad (44)$$

Since $\tau_{k+1} \le 1$ and $\langle \nabla A(\theta_{k+1}) - \nabla A(\theta_k), \theta_{k+1} - \theta_k \rangle = d_A(\theta_{k+1}, \theta_k) + d_A(\theta_k, \theta_{k+1})$, we get

$$\begin{aligned} f_{\pi}^{(\alpha)}(\theta_{k+1}) &\le f_{\pi}^{(\alpha)}(\theta_k) - \frac{1}{\tau_{k+1}} (d_A(\theta_k, \theta_{k+1}) + d_A(\theta_{k+1}, \theta_k)) + \frac{1}{\tau_{k+1}} d_A(\theta_{k+1}, \theta_k) \\ &= f_{\pi}^{(\alpha)}(\theta_k) - \frac{1}{\tau_{k+1}} d_{A^*}(\nabla A(\theta_{k+1}), \nabla A(\theta_k)). \end{aligned}$$

Now, we consider some $v > 0$. Since $\theta_k \to \theta_{\pi}$, there exists $K \in \mathbb{N}$ such that for any $k \ge K$, $\theta_k \in B(\theta_{\pi}, v)$ and $\nabla A(\theta_k) \in B(\nabla A(\theta_{\pi}), v)$ (recall that $\nabla A$ is continuous on int $\Theta$). Thus, we can write following the same steps as in the proof of Proposition 13 that for any $k \ge K$,

$$d_{A^*}(\nabla A(\theta_{k+1}), \nabla A(\theta_k)) = \frac{1}{2} \| \nabla A(\theta_{k+1}) - \nabla A(\theta_k) \|^2_{\nabla^2 A^*(\nabla A(\theta_{\pi}))} + o(v^2).$$

From there, we obtain that

$$\begin{aligned} d_{A^*}(\nabla A(\theta_{k+1}), \nabla A(\theta_k)) &= \frac{\tau_{k+1}^2}{2} \| \nabla f_{\pi}^{(\alpha)}(\theta_k) \|^2_{\nabla^2 A^*(\nabla A(\theta_{\pi}))} + o(v^2) \\ &\ge \tau_{k+1}^2 \alpha f_{\pi}^{(\alpha)}(\theta_k) + o(v^2), \end{aligned}$$

where the last inequality comes from Proposition 13. With our previous point, this yields

$$\begin{aligned} f_{\pi}^{(\alpha)}(\theta_{k+1}) &\le (1 - \alpha \tau_{k+1}) f_{\pi}^{(\alpha)}(\theta_k) + o(v^2) \\ &< (1 - \alpha \delta \epsilon) f_{\pi}^{(\alpha)}(\theta_k) + o(v^2) \end{aligned}$$

for any constant $\delta \in (0, 1)$. By choosing $\delta$ or $v$ small enough, we finally obtain

$$f_{\pi}^{(\alpha)}(\theta_{k+1}) \le (1 - \alpha \delta \epsilon) f_{\pi}^{(\alpha)}(\theta_k).$$

This means that for any $k \ge K$, we have $f_{\pi}^{(\alpha)}(\theta_k) \le (1 - \alpha \delta \epsilon)^{k-K} f_{\pi}^{(\alpha)}(\theta_K)$. Since the sequence $\{f_{\pi}^{(\alpha)}(\theta_k)\}_{k \in \mathbb{N}}$ is decreasing, $f_{\pi}^{(\alpha)}(\theta_K) \le f_{\pi}^{(\alpha)}(\theta_0)$, which shows the result with $C = (1 - \alpha \delta \epsilon)^{-K}$. $\square$

<figure><img src="image_9.png" /></figure>
C Computation of a Bregman proximal operator

Consider an orthonormal matrix $Q$ and the family of Gaussian distribution with covariance of the form $\Sigma = Q \text{diag}(\sigma_1^2, ..., \sigma_d^2)Q^\top$ and mean $\mu \in \mathbb{R}^d$. It is an exponential family with parameters $\theta = (\theta_1, \theta_2)^\top$, with $\theta_1 = \text{diag}(\frac{1}{\sigma_1^2}, ..., \frac{1}{\sigma_d^2})Q^\top\mu$ and $\theta_2 = -(\frac{1}{2\sigma_1^2}, ..., \frac{1}{2\sigma_d^2})^\top$. Its sufficient statistics is $\Gamma(x) = (Q^\top x, (Q^\top x_1)^2, ..., (Q^\top x_d)^2)$. Its log-partition function is $A(\theta) = -\frac{1}{4}\theta_1^\top(\text{diag}(\theta_2))^{-1}\theta_1 + \frac{d}{2}\log(2\pi) - \frac{1}{2}\sum_{i=1}^d \log(-2(\theta_2)_i)$, and its natural parameters $\nabla A(\theta)$ are $Q^\top\mu$ and $((Q^\top\mu)_1^2 + \sigma_1^2, ..., (Q^\top\mu)_d^2 + \sigma_d^2)^\top$.

We consider a regularizer that enforces sparsity on some components of the mean. We propose to this end

$$r(\theta) = \sum_{i=1}^{d} \eta_i |(\theta_1)_i|, \qquad (45)$$

where $\eta_i \ge 0$ for $i \in [1, d]$.

Since $\sigma_i^2 > 0$ for all $i \in [1, d]$, having a null component in $\theta_1$ means that $Q^\top\mu$ has a null component, promoting sparsity in $Q^\top\mu$. We aim at computing $\tilde{\theta} = \text{prox}_{\tau r}^A(\theta)$.

**Lemma 3.** Consider the Gaussian family defined above. Consider $q_\theta$ in this family, with $\theta \in \text{int} \Theta$ and whose mean and covariance are respectively $\mu$ and $Q \text{diag}(\sigma_1^2, ..., \sigma_d^2)Q^\top$. If we consider the regularizing function defined in Eq. (45), then $\tilde{\theta} = \text{prox}_{\tau r}^A(\theta)$ is such that the mean $\tilde{\mu}$ and covariance $Q \text{diag}(\tilde{\sigma}_1^2, ..., \tilde{\sigma}_d^2)Q^\top$ of $q_\theta$ satisfy for any $i \in [1, d]$

$$ (Q^\top \tilde{\mu})_i = \begin{cases} 0 & \text{if } (Q^\top \mu)_i \in [-\tau \eta_i, \tau \eta_i], \\ -\tau \eta_i + (Q^\top \mu)_i & \text{if } (Q^\top \mu)_i > \tau \eta_i, \\ \tau \eta_i + (Q^\top \mu)_i & \text{if } (Q^\top \mu)_i < -\tau \eta_i. \end{cases} \\ \\ \tilde{\sigma}_i^2 = (\sigma_i)^2 + ((Q^\top \mu)_i^2 - (Q^\top \tilde{\mu})_i^2). $$

Consider $i \in [1, d]$. In the particular case where $\eta_i = 0$, then $\mu_i^* = \mu_i$ and $\tilde{\sigma}_i^2 = (\sigma_i)^2$. We can also remark that we always have $\tilde{\sigma}_i^2 \ge \sigma_i^2$, with equality if and only if $(Q^\top \mu)_i = 0$. Therefore, the operator $\text{prox}_{\tau r}^A$ modifies $q_\theta$ by shrinking certain values of the mean to zero, but it increases the variance. In particular, the bigger the $(Q^\top \mu)_i$, the bigger the variance increase.

When $Q = I$, the exponential family is the family of Gaussian distributions with diagonal covariance. The above results can thus be applied to this family too.

*Proof.* The regularizing function $r$ is separable, so we study the optimality condition for every $i \in [1, d]$. This is justified by [10, Proposition 16.8], which shows that $\partial r(\theta)$ is the Cartesian product of its subdifferentials with respect to each of its variable. Therefore, for $i \in [1, d]$, we have

$$ \begin{cases} \frac{1}{\tau}((Q^\top \mu)_i - (Q^\top \tilde{\mu})_i) & \in \eta_i \partial |((\tilde{\theta}_1)_i), \\ \frac{1}{\tau}((Q^\top \mu)_i^2 + \sigma_i^2 - ((Q^\top \tilde{\mu})_i^2 + \tilde{\sigma}_i^2)) & = 0, \end{cases} $$

from which we already deduce the result about the standard deviation.

Because $(\tilde{\Sigma}_i)^2 > 0$, the sign of $(\tilde{\theta}_1)_i = \frac{1}{(\tilde{\Sigma}_i)^2}(Q^\top \tilde{\mu})_i$ is the sign of $(Q^\top \tilde{\mu})_i$ and we get that

$$ (Q^\top \mu)_i - (Q^\top \tilde{\mu})_i \in \begin{cases} [-\tau \eta_i, \tau \eta_i] & \text{if } (Q^\top \tilde{\mu})_i = 0, \\ \{\tau \eta_i\} & \text{if } (Q^\top \tilde{\mu})_i > 0, \\ \{-\tau \eta_i\} & \text{if } (Q^\top \tilde{\mu})_i < 0, \end{cases} $$

from which we can obtain the result.

□
C Computation of a Bregman proximal operator

Consider an orthonormal matrix $Q$ and the family of Gaussian distribution with covariance of the form $\Sigma = Q \text{diag}(\sigma_1^2, ..., \sigma_d^2)Q^\top$ and mean $\mu \in \mathbb{R}^d$. It is an exponential family with parameters $\theta = (\theta_1, \theta_2)^\top$, with $\theta_1 = \text{diag}(\frac{1}{\sigma_1^2}, ..., \frac{1}{\sigma_d^2})Q^\top\mu$ and $\theta_2 = -(\frac{1}{2\sigma_1^2}, ..., \frac{1}{2\sigma_d^2})^\top$. Its sufficient statistics is $\Gamma(x) = (Q^\top x, (Q^\top x_1)^2, ..., (Q^\top x_d)^2)$. Its log-partition function is $A(\theta) = -\frac{1}{4}\theta_1^\top(\text{diag}(\theta_2))^{-1}\theta_1 + \frac{d}{2}\log(2\pi) - \frac{1}{2}\sum_{i=1}^d \log(-2(\theta_2)_i)$, and its natural parameters $\nabla A(\theta)$ are $Q^\top\mu$ and $((Q^\top\mu)_1^2 + \sigma_1^2, ..., (Q^\top\mu)_d^2 + \sigma_d^2)^\top$.

We consider a regularizer that enforces sparsity on some components of the mean. We propose to this end

$$r(\theta) = \sum_{i=1}^{d} \eta_i |(\theta_1)_i|, \qquad (45)$$

where $\eta_i \ge 0$ for $i \in [1, d]$.

Since $\sigma_i^2 > 0$ for all $i \in [1, d]$, having a null component in $\theta_1$ means that $Q^\top\mu$ has a null component, promoting sparsity in $Q^\top\mu$. We aim at computing $\tilde{\theta} = \text{prox}_{\tau r}^A(\theta)$.

**Lemma 3.** Consider the Gaussian family defined above. Consider $q_\theta$ in this family, with $\theta \in \text{int} \Theta$ and whose mean and covariance are respectively $\mu$ and $Q \text{diag}(\sigma_1^2, ..., \sigma_d^2)Q^\top$. If we consider the regularizing function defined in Eq. (45), then $\tilde{\theta} = \text{prox}_{\tau r}^A(\theta)$ is such that the mean $\tilde{\mu}$ and covariance $Q \text{diag}(\tilde{\sigma}_1^2, ..., \tilde{\sigma}_d^2)Q^\top$ of $q_\theta$ satisfy for any $i \in [1, d]$

$$ (Q^\top \tilde{\mu})_i = \begin{cases} 0 & \text{if } (Q^\top \mu)_i \in [-\tau \eta_i, \tau \eta_i], \\ -\tau \eta_i + (Q^\top \mu)_i & \text{if } (Q^\top \mu)_i > \tau \eta_i, \\ \tau \eta_i + (Q^\top \mu)_i & \text{if } (Q^\top \mu)_i < -\tau \eta_i. \end{cases} \\ \tilde{\sigma}_i^2 = (\sigma_i)^2 + ((Q^\top \mu)_i^2 - (Q^\top \tilde{\mu})_i^2). $$

Consider $i \in [1, d]$. In the particular case where $\eta_i = 0$, then $\mu_i^* = \mu_i$ and $\tilde{\sigma}_i^2 = (\sigma_i)^2$. We can also remark that we always have $\tilde{\sigma}_i^2 \ge \sigma_i^2$, with equality if and only if $(Q^\top \mu)_i = 0$. Therefore, the operator $\text{prox}_{\tau r}^A$ modifies $q_\theta$ by shrinking certain values of the mean to zero, but it increases the variance. In particular, the bigger the $(Q^\top \mu)_i$, the bigger the variance increase.

When $Q = I$, the exponential family is the family of Gaussian distributions with diagonal covariance. The above results can thus be applied to this family too.

**Proof.** The regularizing function $r$ is separable, so we study the optimality condition for every $i \in [1, d]$. This is justified by [10, Proposition 16.8], which shows that $\partial r(\theta)$ is the Cartesian product of its subdifferentials with respect to each of its variable. Therefore, for $i \in [1, d]$, we have

$$ \begin{cases} \frac{1}{\tau}((Q^\top \mu)_i - (Q^\top \tilde{\mu})_i) & \in \eta_i \partial |((\tilde{\theta}_1)_i), \\ \frac{1}{\tau}((Q^\top \mu)_i^2 + \sigma_i^2 - ((Q^\top \tilde{\mu})_i^2 + \tilde{\sigma}_i^2)) & = 0, \end{cases} $$

from which we already deduce the result about the standard deviation.

Because $(\tilde{\Sigma}_i)^2 > 0$, the sign of $(\tilde{\theta}_1)_i = \frac{1}{(\tilde{\Sigma}_i)^2}(Q^\top \tilde{\mu})_i$ is the sign of $(Q^\top \tilde{\mu})_i$ and we get that

$$ (Q^\top \mu)_i - (Q^\top \tilde{\mu})_i \in \begin{cases} [-\tau \eta_i, \tau \eta_i] & \text{if } (Q^\top \tilde{\mu})_i = 0, \\ \{\tau \eta_i\} & \text{if } (Q^\top \tilde{\mu})_i > 0, \\ \{-\tau \eta_i\} & \text{if } (Q^\top \tilde{\mu})_i < 0, \end{cases} $$

from which we can obtain the result.

□
# D Supplementary numerical experiments

## D.1 Understanding the influence of the parameters

To this end, we use Gaussian targets in various dimensions $d$, with unnormalized density of the form

$$ \tilde{\pi}(x) = \exp \left( -\frac{1}{2} (x - \bar{\mu})^{\top} \bar{\Sigma}_{\kappa}^{-1} (x - \bar{\mu}) \right), \forall x \in \mathbb{R}^{d}. \quad (46) $$

Their means $\bar{\mu}$ are chosen uniformly in $[-0.5, 0.5]^d$ and their covariance matrices $\bar{\Sigma}_{\kappa}$ are chosen with a condition number equal to $\kappa$, following the procedure in [66, Section 5].

We now discuss the influence of $\alpha, \tau$ on the practical speed and robustness of Algorithm 2, in its non-regularized version RMM. We recall that this algorithm resorts to importance sampling to approximate the integrals involved in the computation of $\pi_{\theta}^{(\alpha)}(\Gamma)$, which creates an approximation error linked with the sample size, $N$. The influence of $\tau$ can be understood through the theory on stochastic Bregman gradient descent with fixed step-size. In particular, [40, Theorem 5.3] states that such methods converge to a neighborhood of the optimum, whose size decreases with $\tau$. On the other hand, low values of $\alpha$ amount to a concave transformation of the importance weights, which is known in the importance sampling field to lead to a higher effective sample size [54].

In order to highlight this compromise between speed and robustness, we use the RMM algorithm to approximate the target described in Eq. (46) with $\kappa = 10$. We use a constant number of samples per iteration $N = 500$, for $d \in \{5, 10, 20, 40\}$. It is recommended for importance sampling procedures that the sample size grows as $\exp(d)$ to avoid weight degeneracy [12]. In our setting, $d$ increases while $N$ remains constant, thus creating approximation errors that increase with $d$.

For each dimension, we test $\alpha \in \{0.5, 1.0\}$ and $\tau \in \{0.25, 0.5, 1.0\}$. We track the square errors $\|\bar{\mu} - \mu_k\|^2$ and $\|\bar{\Sigma}_{\kappa} - \Sigma_k\|^2_F$, that are averaged over $10^3$ independent runs.

<figure><img src="image_12.png" /><figcaption>Figure 5: MSE averaged over $10^3$ runs, in dimension $d = 5$.</figcaption></figure>

In dimension $d = 5$, all the choices of parameters lead to convergence, as shown in Fig. 5. We can notice that the lowest values of $\tau$ lead to the slowest convergence, but the values reached are lower. On the contrary, when $\tau = 1.0$, the algorithm stops early at higher values.

Finally, for $d = 40$, only the lowest values of $\alpha$ and $\tau$ yield a significant decrease of the MSE as shown in Fig. 6. This shows that low values of $\alpha$ and $\tau$ can counteract high approximation errors. As expected, the convergence is slower and the final MSE values are higher than in lower dimensions.
# D Supplementary numerical experiments

## D.1 Understanding the influence of the parameters

To this end, we use Gaussian targets in various dimensions $d$, with unnormalized density of the form

$$ \tilde{\pi}(x) = \exp \left( -\frac{1}{2} (x - \bar{\mu})^{\top} \bar{\Sigma}_{\kappa}^{-1} (x - \bar{\mu}) \right), \forall x \in \mathbb{R}^{d}. \quad (46) $$

Their means $\bar{\mu}$ are chosen uniformly in $[-0.5, 0.5]^d$ and their covariance matrices $\bar{\Sigma}_{\kappa}$ are chosen with a condition number equal to $\kappa$, following the procedure in [66, Section 5].

We now discuss the influence of $\alpha, \tau$ on the practical speed and robustness of Algorithm 2, in its non-regularized version RMM. We recall that this algorithm resorts to importance sampling to approximate the integrals involved in the computation of $\pi_{\theta}^{(\alpha)}(\Gamma)$, which creates an approximation error linked with the sample size, $N$. The influence of $\tau$ can be understood through the theory on stochastic Bregman gradient descent with fixed step-size. In particular, [40, Theorem 5.3] states that such methods converge to a neighborhood of the optimum, whose size decreases with $\tau$. On the other hand, low values of $\alpha$ amount to a concave transformation of the importance weights, which is known in the importance sampling field to lead to a higher effective sample size [54].

In order to highlight this compromise between speed and robustness, we use the RMM algorithm to approximate the target described in Eq. (46) with $\kappa = 10$. We use a constant number of samples per iteration $N = 500$, for $d \in \{5, 10, 20, 40\}$. It is recommended for importance sampling procedures that the sample size grows as $\exp(d)$ to avoid weight degeneracy [12]. In our setting, $d$ increases while $N$ remains constant, thus creating approximation errors that increase with $d$.

For each dimension, we test $\alpha \in \{0.5, 1.0\}$ and $\tau \in \{0.25, 0.5, 1.0\}$. We track the square errors $\|\bar{\mu} - \mu_k\|^2$ and $\|\bar{\Sigma}_{\kappa} - \Sigma_k\|^2_F$, that are averaged over $10^3$ independent runs.

<figure><img src="image_10.png" /><figcaption>Figure 5: MSE averaged over $10^3$ runs, in dimension $d = 5$.</figcaption></figure>

In dimension $d = 5$, all the choices of parameters lead to convergence, as shown in Fig. 5. We can notice that the lowest values of $\tau$ lead to the slowest convergence, but the values reached are lower. On the contrary, when $\tau = 1.0$, the algorithm stops early at higher values.

Finally, for $d = 40$, only the lowest values of $\alpha$ and $\tau$ yield a significant decrease of the MSE as shown in Fig. 6. This shows that low values of $\alpha$ and $\tau$ can counteract high approximation errors. As expected, the convergence is slower and the final MSE values are higher than in lower dimensions.
<figure><img src="image_13.png" /><figcaption>Figure 6: MSE averaged over $10^3$ runs, in dimension $d = 40$.</figcaption></figure>

This study shows that the parameters $\alpha$ and $\tau$ should be lowered to compensate for high approximation errors possibly arising in Algorithm 2. On the contrary, when these errors are low, one can increase the values of $\tau$ to create faster algorithms.

## D.2 Comparison with the variational Rényi bound on a Gaussian target

Our theoretical analysis provides guidelines to choose the step-size $\tau$ for our RMM algorithm (Propositions 14 and 15) but also shows that there is no equivalent guarantees for the VRB algorithm (see Proposition 11). In particular, poorly chosen step-sizes could create unstable behaviors. We thus investigate these effects in the following by comparing our novel RMM algorithm with the VRB algorithm on Gaussian targets.

We use Gaussian target from Eq. (46), with $\kappa = 10$, and $d = 5$. Each algorithm is run with constant number of samples $N = 500$, and constant values of the step-size $\tau$. We test values of $\alpha$ corresponding to the Hellinger distance ($\alpha = 0.5$) and the KL divergence ($\alpha = 1.0$). We test two different exponential families: Gaussian with full covariance, and Gaussian with diagonal covariance. For each tested value of $\tau$, $10^3$ runs are performed.

Figure 7 shows that the VRB algorithm used with diagonal covariance in the approximation family exhibits two distinct regimes. For sufficiently low values of $\tau$, it is able to improve the estimates compared to initialization, but once $\tau$ crosses a certain threshold, the MSE reaches very high values, showing a degradation from the initialization. The VRB algorithm with full covariance in the approximation family is not able to create covariance matrices that are positive definite, hence it stops after initialization. On the contrary, our RMM algorithm does not degrade the values reached at initialization even for the worst settings of $\tau$, and reaches the lowest MSE values for properly chosen step-sizes.

This confirms that the lack of Euclidean smoothness of $f_{\pi}^{(\alpha)}$ translates numerically into a high level of instability of VRB with respect to the choice of the step-size. On the contrary, the RMM algorithm has a more stable behavior even for poorly chosen step-sizes, confirming the theoretical study of Section 5.

## References

[1] O. Akyildiz and J. Míguez. Convergence rates for optimised adaptive importance samplers. *Statistics and Computing*, 31(12), 2021.
<figure><img src="image_11.png" /><figcaption>Figure 6: MSE averaged over $10^3$ runs, in dimension $d = 40$.</figcaption></figure>

This study shows that the parameters $\alpha$ and $\tau$ should be lowered to compensate for high approximation errors possibly arising in Algorithm 2. On the contrary, when these errors are low, one can increase the values of $\tau$ to create faster algorithms.

## D.2 Comparison with the variational Rényi bound on a Gaussian target

Our theoretical analysis provides guidelines to choose the step-size $\tau$ for our RMM algorithm (Propositions 14 and 15) but also shows that there is no equivalent guarantees for the VRB algorithm (see Proposition 11). In particular, poorly chosen step-sizes could create unstable behaviors. We thus investigate these effects in the following by comparing our novel RMM algorithm with the VRB algorithm on Gaussian targets.

We use Gaussian target from Eq. (46), with $\kappa = 10$, and $d = 5$. Each algorithm is run with constant number of samples $N = 500$, and constant values of the step-size $\tau$. We test values of $\alpha$ corresponding to the Hellinger distance ($\alpha = 0.5$) and the KL divergence ($\alpha = 1.0$). We test two different exponential families: Gaussian with full covariance, and Gaussian with diagonal covariance. For each tested value of $\tau$, $10^3$ runs are performed.

Figure 7 shows that the VRB algorithm used with diagonal covariance in the approximation family exhibits two distinct regimes. For sufficiently low values of $\tau$, it is able to improve the estimates compared to initialization, but once $\tau$ crosses a certain threshold, the MSE reaches very high values, showing a degradation from the initialization. The VRB algorithm with full covariance in the approximation family is not able to create covariance matrices that are positive definite, hence it stops after initialization. On the contrary, our RMM algorithm does not degrade the values reached at initialization even for the worst settings of $\tau$, and reaches the lowest MSE values for properly chosen step-sizes.

This confirms that the lack of Euclidean smoothness of $f_{\pi}^{(\alpha)}$ translates numerically into a high level of instability of VRB with respect to the choice of the step-size. On the contrary, the RMM algorithm has a more stable behavior even for poorly chosen step-sizes, confirming the theoretical study of Section 5.

## References

[1] O. Akyildiz and J. Míguez. Convergence rates for optimised adaptive importance samplers. *Statistics and Computing*, 31(12), 2021.
<figure><img src="image_14.png" /><figcaption>Figure 7: MSE in the estimation of $\bar{\mu}$ and $\bar{\Sigma}_{\kappa}$ ($d = 5$) after 100 iterations, against values of $\tau$. For each value of $\tau$, $10^3$ runs with 500 samples per iteration are conducted. The dotted black lines represent the MSE at initialization. The prefix dG refer to the family of diagonal Gaussians, while the prefix G refers to Gaussians with full covariance.</figcaption></figure>

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

[3] S. Amari. Natural gradient works efficiently in learning. *Neural Computation*, 10(2):251–276, 1998.

[4] Y. Atchadé, G. Fort, and E. Moulines. On perturbed proximal gradient algorithms. *Journal of Machine Learning Research*, 18(10):1–33, 2017.

[5] O. Banerjee, L. E. Ghaoui, and A. d’Aspremont. Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data. *Journal of Machine Learning Research*, 9(15):485–516, 2008.

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

[7] H. Bauschke, J. Bolte, and M. Teboulle. A descent lemma beyond Lipschitz gradient continuity: revisited and applications. *Mathematics of Operations Research*, 42(2):330–348, 2017.
<figure><img src="image_12.png" /><figcaption>Figure 7: MSE in the estimation of $\bar{\mu}$ and $\bar{\Sigma}_{\kappa}$ ($d = 5$) after 100 iterations, against values of $\tau$. For each value of $\tau$, $10^3$ runs with 500 samples per iteration are conducted. The dotted black lines represent the MSE at initialization. The prefix dG refer to the family of diagonal Gaussians, while the prefix G refers to Gaussians with full covariance.</figcaption></figure>

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

[3] S. Amari. Natural gradient works efficiently in learning. *Neural Computation*, 10(2):251–276, 1998.

[4] Y. Atchadé, G. Fort, and E. Moulines. On perturbed proximal gradient algorithms. *Journal of Machine Learning Research*, 18(10):1–33, 2017.

[5] O. Banerjee, L. E. Ghaoui, and A. d’Aspremont. Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data. *Journal of Machine Learning Research*, 9(15):485–516, 2008.

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

[7] H. Bauschke, J. Bolte, and M. Teboulle. A descent lemma beyond Lipschitz gradient continuity: revisited and applications. *Mathematics of Operations Research*, 42(2):330–348, 2017.
[8] H. Bauschke and J. Borwein. Legendre functions and the method of random Bregman projections. *Journal of Convex Analysis*, 4:27-67, 1997.

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

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

[11] A. Benfenati, E. Chouzenoux, and J.-C. Pesquet. Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation. *Signal Processing*, 169, 2020.

[12] T. Bengsston, P. Bickel, and B. Li. *Curse-of-dimensionality revisited: Collapse of the particle filter in very large scale systems*, pages 316–334. Institute of Mathematical Statistics, 2008.

[13] D. Blei, A. Kucukelbir, and J. McAuliffe. Variational inference: A review for the statistician. *Journal of the American Statistical Association*, 112(518):859-877, 2017.

[14] J. Bolte, S. Sabach, M. Teboulle, and Y. Vaisbourd. First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems. *SIAM Journal on Optimization*, 28(3):2131-2151, 2018.

[15] L. D. Brown. *Fundamentals of Statistical Exponential Families with Applications in Statistical Decision Theory*. Institute of Mathematical Statistics, 1986.

[16] M. F. Bugallo, V. Elvira, and L. Martino. A new strategy for effective learning in population Monte Carlo sampling. In *Asilomar Conference on Signals, Systems and Computers*, pages 1540–1544, 2016.

[17] T. Bui. Connecting the thermodynamic variational objective and annealed importance sampling. Technical report, 2020.

[18] L. Bungert, T. Roith, D. Tenbrinck, and M. Burger. A Bregman learning framework for sparse neural networks. *Journal of Machine Learning Research*, 23(192):1-43, 2022.

[19] D. Butnariu and A. Iusem. *Totally convex functions*, chapter 1. Kluwer Academic Publisher, 2000.

[20] T. Campbell and B. Beronov. Sparse variational inference: Bayesian coresets from scratch. In *Advances in Neural Information Processing Systems (NeurPIS)*, volume 32, 2019.

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

[22] N. L. Carothers. *Real Analysis*. Cambridge University Press, 2000.

[23] P. Combettes and J.-C. Pesquet. *Proximal Splitting Methods in Signal Processing*, page 185-212. Springer-Verlag, New York, 2010.

[24] P. L. Combettes and C. L. Müller. Perspective functions: Proximal calculus and applications in high-dimensional statistics. *Journal of Mathematical Analysis and Applications*, 457(2):1283-1306, 2018.

[25] J. M. Cornuet, J. M. Marin, A. Mira, and C. P. Robert. Adaptive multiple importance sampling. *Scandinavian Journal of Statistics*, 39(4):798-812, December 2012.

[26] I. Csiszár and P. Shields. *Information Theory and Statistics: A Tutorial*. Now Foundations and Trends, 2004.
[8] H. Bauschke and J. Borwein. Legendre functions and the method of random Bregman projections. *Journal of Convex Analysis*, 4:27-67, 1997.

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

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

[11] A. Benfenati, E. Chouzenoux, and J.-C. Pesquet. Proximal approaches for matrix optimization problems: Application to robust precision matrix estimation. *Signal Processing*, 169, 2020.

[12] T. Bengsston, P. Bickel, and B. Li. *Curse-of-dimensionality revisited: Collapse of the particle filter in very large scale systems*, pages 316–334. Institute of Mathematical Statistics, 2008.

[13] D. Blei, A. Kucukelbir, and J. McAuliffe. Variational inference: A review for the statistician. *Journal of the American Statistical Association*, 112(518):859–877, 2017.

[14] J. Bolte, S. Sabach, M. Teboulle, and Y. Vaisbourd. First order methods beyond convexity and Lipschitz gradient continuity with applications to quadratic inverse problems. *SIAM Journal on Optimization*, 28(3):2131–2151, 2018.

[15] L. D. Brown. *Fundamentals of Statistical Exponential Families with Applications in Statistical Decision Theory*. Institute of Mathematical Statistics, 1986.

[16] M. F. Bugallo, V. Elvira, and L. Martino. A new strategy for effective learning in population Monte Carlo sampling. In *Asilomar Conference on Signals, Systems and Computers*, pages 1540–1544, 2016.

[17] T. Bui. Connecting the thermodynamic variational objective and annealed importance sampling. Technical report, 2020.

[18] L. Bungert, T. Roith, D. Tenbrinck, and M. Burger. A Bregman learning framework for sparse neural networks. *Journal of Machine Learning Research*, 23(192):1–43, 2022.

[19] D. Butnariu and A. Iusem. *Totally convex functions*, chapter 1. Kluwer Academic Publisher, 2000.

[20] T. Campbell and B. Beronov. Sparse variational inference: Bayesian coresets from scratch. In *Advances in Neural Information Processing Systems (NeurPIS)*, volume 32, 2019.

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

[22] N. L. Carothers. *Real Analysis*. Cambridge University Press, 2000.

[23] P. Combettes and J.-C. Pesquet. *Proximal Splitting Methods in Signal Processing*, page 185-212. Springer-Verlag, New York, 2010.

[24] P. L. Combettes and C. L. Müller. Perspective functions: Proximal calculus and applications in high-dimensional statistics. *Journal of Mathematical Analysis and Applications*, 457(2):1283-1306, 2018.

[25] J. M. Cornuet, J. M. Marin, A. Mira, and C. P. Robert. Adaptive multiple importance sampling. *Scandinavian Journal of Statistics*, 39(4):798-812, December 2012.

[26] I. Csiszár and P. Shields. *Information Theory and Statistics: A Tutorial*. Now Foundations and Trends, 2004.
[27] K. Daudel, R. Douc, and F. Portier. Infinite-dimensional gradient-based descent for alpha-divergence minimisation. *The Annals of Statistics*, 49(4):2250–2270, 2021.

[28] K. Daudel, R. Douc, and F. Roueff. Monotonic Alpha-divergence minimization. *Journal of Machine Learning Research*, 62(24):1–76, 2023.

[29] A. B. Dieng, D. Tran, R. Ranganath, J. Paisley, and D. Blei. Variational inference via $\chi$ upper bound minimization. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 30, 2017.

[30] J. Domke. Provable smoothness guarantees for black-box variational inference. In *International Conference on Machine Learning (ICML)*, pages 2587–2596, 2020.

[31] R. Douc, A. Guillou, J. Marin, and C. Robert. Convergence of adaptive mixtures of importance sampling schemes. *Annals of Statistics*, 35:420–448, 2007.

[32] A. Durmus, S. Majewski, and B. Miasojedow. Analysis of Langevin Monte Carlo via convex optimization. *Journal of Machine Learning Research*, 20(73):1–49, 2019.

[33] R. Dykstra. An iterative procedure for obtaining *I*-projections onto the intersection of convex sets. *The Annals of Probability*, 13(3):975–984, 1985.

[34] Y. El-Laham, V. Elvira, and M. F. Bugallo. Recursive shrinkage covariance learning in adaptive importance sampling. In *IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)*, pages 624–628, 2019.

[35] Y. El-Laham, L. Martino, V. Elvira, and M. Bugallo. Efficient adaptive multiple importance sampling. In *European Signal Processing Conference (EUSIPCO)*, pages 1–5, 2019.

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

[37] T. Gao, S. Lu, J. Liu, and C. Chu. Randomized Bregman coordinate descent methods for non-Lipschitz optimization. https://arxiv.org/pdf/2001.05202, 2020.

[38] M. E. Gheche, G. Chierchia, and J.-C. Pesquet. Proximity operators of discrete information divergences. *IEEE Transactions on Information Theory*, 64(2):1092–1104, 2018.

[39] R. B. Grosse, C. J. Maddison, and R. R. Salakhutdinov. Annealing between distributions by averaging moments. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 26, 2013.

[40] F. Hanzely and P. Richtárik. Fastest rates for stochastic mirror descent methods. *Computational Optimization and Applications*, 79(3):717–766, 2021.

[41] T. Hastie, R. Tibshirani, and J. Firedman. *The Elements of Statistical Learning*. Springer, 2009.

[42] J. Hensman, M. Rattray, and N. Lawrence. Fast variational inference in the conjugate exponential family. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 25, 2012.

[43] J. Hernandez-Lobato, Y. Li, M. Rowland, T. Bui, D. Hernandez-Lobato, and R. Turner. Black-box alpha divergence minimization. In *International Conference on Machine Learning (ICML)*, volume 48, pages 1511–1520, 2016.

[44] J.-B. Hiriart-Urruty and C. Lemaréchal. *Abstract Duality for Practitioners*, pages 137–193. Springer, 1993.
[27] K. Daudel, R. Douc, and F. Portier. Infinite-dimensional gradient-based descent for alpha-divergence minimisation. *The Annals of Statistics*, 49(4):2250–2270, 2021.

[28] K. Daudel, R. Douc, and F. Roueff. Monotonic Alpha-divergence minimization. *Journal of Machine Learning Research*, 62(24):1–76, 2023.

[29] A. B. Dieng, D. Tran, R. Ranganath, J. Paisley, and D. Blei. Variational inference via $\chi$ upper bound minimization. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 30, 2017.

[30] J. Domke. Provable smoothness guarantees for black-box variational inference. In *International Conference on Machine Learning (ICML)*, pages 2587–2596, 2020.

[31] R. Douc, A. Guillou, J. Marin, and C. Robert. Convergence of adaptive mixtures of importance sampling schemes. *Annals of Statistics*, 35:420–448, 2007.

[32] A. Durmus, S. Majewski, and B. Miasojedow. Analysis of Langevin Monte Carlo via convex optimization. *Journal of Machine Learning Research*, 20(73):1–49, 2019.

[33] R. Dykstra. An iterative procedure for obtaining *I*-projections onto the intersection of convex sets. *The Annals of Probability*, 13(3):975–984, 1985.

[34] Y. El-Laham, V. Elvira, and M. F. Bugallo. Recursive shrinkage covariance learning in adaptive importance sampling. In *IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)*, pages 624–628, 2019.

[35] Y. El-Laham, L. Martino, V. Elvira, and M. Bugallo. Efficient adaptive multiple importance sampling. In *European Signal Processing Conference (EUSIPCO)*, pages 1–5, 2019.

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

[37] T. Gao, S. Lu, J. Liu, and C. Chu. Randomized Bregman coordinate descent methods for non-Lipschitz optimization. https://arxiv.org/pdf/2001.05202, 2020.

[38] M. E. Gheche, G. Chierchia, and J.-C. Pesquet. Proximity operators of discrete information divergences. *IEEE Transactions on Information Theory*, 64(2):1092–1104, 2018.

[39] R. B. Grosse, C. J. Maddison, and R. R. Salakhutdinov. Annealing between distributions by averaging moments. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 26, 2013.

[40] F. Hanzely and P. Richtárik. Fastest rates for stochastic mirror descent methods. *Computational Optimization and Applications*, 79(3):717–766, 2021.

[41] T. Hastie, R. Tibshirani, and J. Firedman. *The Elements of Statistical Learning*. Springer, 2009.

[42] J. Hensman, M. Rattray, and N. Lawrence. Fast variational inference in the conjugate exponential family. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 25, 2012.

[43] J. Hernandez-Lobato, Y. Li, M. Rowland, T. Bui, D. Hernandez-Lobato, and R. Turner. Black-box alpha divergence minimization. In *International Conference on Machine Learning (ICML)*, volume 48, pages 1511–1520, 2016.

[44] J.-B. Hiriart-Urruty and C. Lemaréchal. *Abstract Duality for Practitioners*, pages 137–193. Springer, 1993.
[45] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. *Journal of Machine Learning Research*, 14(4):1303-1347, 2013.

[46] A. Honkela, T. Raiko, M. Kuusela, M. Tornio, and J. Karhunen. Approximate Riemannian conjugate gradient learning for fixed-form variational Bayes. *Journal of Machine Learning Research*, 11(106):3235-3268, 2010.

[47] Y. Huang, E. Chouzenoux, and J.-C. Pesquet. Unrolled variational Bayesian algorithm for image blind deconvolution. *IEEE Transactions on Image Processing*, 32:430-445, 2022.

[48] E. Ionides. Truncated importance sampling. *Journal of Computational and Graphical Statistics*, 17(2):295-311, 2008.

[49] G. Ji, D. Sujono, and E. Sudderth. Marginalized stochastic natural gradients for black-box variational inference. In *International Conference on Machine Learning (ICML)*, volume 139, pages 4870-4881, 2021.

[50] H. Karimi, J. Nutini, and M. W. Schmidt. Linear convergence of gradient and proximal-gradient methods under the Polyak-Łojasiewicz condition. In *Joint European Conference on Machine Learning and Knowledge Discovery in Databases (ECML PKDD)*, 2016.

[51] M. Khan, R. Babanezhad, W. Lin, M. Schmidt, and M. Sugiyama. Faster stochastic variational inference using proximal-gradient methods with general divergence functions. In *Conference on Uncertainty in Artificial Intelligence (UAI)*, page 319-328, 2016.

[52] 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.

[53] J. Knoblauch, J. Jewson, and T. Damoulas. An optimization-centric view on Bayes' rule: Reviewing and generalizing variational inference. *Journal of Machine Learning Research*, 23(132):1-109, 2022.

[54] E. Koblents and J. Míguez. A population Monte Carlo scheme with transformed weights and its application to stochastic kinetic models. *Statistics and Computing*, 25(2):407-425, 2013.

[55] A. Korba and F. Portier. Adaptive importance sampling meets mirror descent: a bias-variance tradeoff. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, volume 151, pages 11503-11527, 2022.

[56] M. A. Kumar and I. Sason. Projection theorems for the Rényi divergence on $\alpha$-convex sets. *IEEE Transactions on Information Theory*, 62(9):4924-4935, 2016.

[57] M. Lambert, S. Chewi, F. Bach, S. Bonnabel, and P. Rigollet. Variational inference via Wasserstein gradient flows. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 35, 2022.

[58] Y. Li and R. Turner. Rényi divergence variational inference. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 29, 2016.

[59] W. Lin, M. Khan, and M. Schmidt. Fast and simple natural-gradient variational inference with mixture of exponential-family approximations. In *International Conference on Machine Learning (ICML)*, volume 97, pages 3992-4002, 2019.

[60] 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.
[45] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. *Journal of Machine Learning Research*, 14(4):1303–1347, 2013.

[46] A. Honkela, T. Raiko, M. Kuusela, M. Tornio, and J. Karhunen. Approximate Riemannian conjugate gradient learning for fixed-form variational Bayes. *Journal of Machine Learning Research*, 11(106):3235–3268, 2010.

[47] Y. Huang, E. Chouzenoux, and J.-C. Pesquet. Unrolled variational Bayesian algorithm for image blind deconvolution. *IEEE Transactions on Image Processing*, 32:430–445, 2022.

[48] E. Ionides. Truncated importance sampling. *Journal of Computational and Graphical Statistics*, 17(2):295–311, 2008.

[49] G. Ji, D. Sujono, and E. Sudderth. Marginalized stochastic natural gradients for black-box variational inference. In *International Conference on Machine Learning (ICML)*, volume 139, pages 4870–4881, 2021.

[50] H. Karimi, J. Nutini, and M. W. Schmidt. Linear convergence of gradient and proximal-gradient methods under the Polyak-Łojasiewicz condition. In *Joint European Conference on Machine Learning and Knowledge Discovery in Databases (ECML PKDD)*, 2016.

[51] M. Khan, R. Babanezhad, W. Lin, M. Schmidt, and M. Sugiyama. Faster stochastic variational inference using proximal-gradient methods with general divergence functions. In *Conference on Uncertainty in Artificial Intelligence (UAI)*, page 319–328, 2016.

[52] 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.

[53] J. Knoblauch, J. Jewson, and T. Damoulas. An optimization-centric view on Bayes' rule: Reviewing and generalizing variational inference. *Journal of Machine Learning Research*, 23(132):1–109, 2022.

[54] E. Koblents and J. Míguez. A population Monte Carlo scheme with transformed weights and its application to stochastic kinetic models. *Statistics and Computing*, 25(2):407–425, 2013.

[55] A. Korba and F. Portier. Adaptive importance sampling meets mirror descent: a bias-variance tradeoff. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, volume 151, pages 11503–11527, 2022.

[56] M. A. Kumar and I. Sason. Projection theorems for the Rényi divergence on $\alpha$-convex sets. *IEEE Transactions on Information Theory*, 62(9):4924–4935, 2016.

[57] M. Lambert, S. Chewi, F. Bach, S. Bonnabel, and P. Rigollet. Variational inference via Wasserstein gradient flows. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 35, 2022.

[58] Y. Li and R. Turner. Rényi divergence variational inference. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 29, 2016.

[59] W. Lin, M. Khan, and M. Schmidt. Fast and simple natural-gradient variational inference with mixture of exponential-family approximations. In *International Conference on Machine Learning (ICML)*, volume 97, pages 3992–4002, 2019.

[60] 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.
[61] L. Martino, V. Elvira, D. Luengo, and J. Corander. An adaptive population importance sampler. In *IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, pages 8038–8042, 2014.

[62] L. Martino, V. Elvira, J. Míguez, A. Artés-Rodríguez, and P. Djurić. A comparison of clipping strategies for importance sampling. In *IEEE Statistical Signal Processing Workshop (SSP)*, pages 558–562, 2018.

[63] 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.

[64] T. Minka. Divergence measures and message passing. Technical report, 2005.

[65] P. D. Moral, A. Doucet, and A. Jasra. Sequential Monte Carlo samplers. *Journal of the Royal Statistical Society: Series B*, 68(3):411–436, 2006.

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

[67] M. Mukkamala, P. Ochs, T. Pock, and S. Sabach. Convex-concave backtracking for inertial Bregman proximal gradient algorithms in nonconvex optimization. *SIAM Journal on Mathematics of Data Science*, 2(3):658–682, 2020.

[68] R. Neal. Annealed importance sampling. *Statistics and Computing*, 11:125–139, 2001.

[69] F. Nielsen and R. Nock. Entropies and cross-entropies of exponential families. In *IEEE International Conference on Image Processing (ICIP)*, pages 3621–3624, 2010.

[70] J.-P. Penot. *Calculus without Derivatives*. Springer, 2013.

[71] R. Ranganath, S. Gerrish, and D. Blei. Black box variational inference. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, volume 33, pages 814–822, 2014.

[72] A. Saha, K. Barath, and S. Kurtek. A geometric variational approach to Bayesian inference. *Journal of the American Statistical Association*, 115(530):822–835, 2020.

[73] Y. Shao, Y. Zhou, and D. Cai. Variational inference with graph regularization for image annotation. *ACM Transactions on Intelligent Systems and Technology*, 2(2):1–21, 2011.

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

[75] M. Titsias and M. Lázaro-Gredilla. Local expectation gradients for black box variational inference. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 28, 2015.

[76] T. van Erven and P. Harremoes. Rényi divergence and Kullback-Leibler divergence. *IEEE Transactions of Information Theory*, 60(7):3797–3820, 2014.

[77] A. Vehtari, D. Simpson, A. Gelman, Y. Yao, and J. Gabry. Pareto smoothed importance sampling. https://arxiv.org/abs/1507.02646, 2015.

[78] S. Vempala and A. Wibosono. Rapid convergence of the unadjusted Langevin algorithm: Isoperimetry suffices. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 32, 2019.
[61] L. Martino, V. Elvira, D. Luengo, and J. Corander. An adaptive population importance sampler. In *IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, pages 8038–8042, 2014.

[62] L. Martino, V. Elvira, J. Míguez, A. Artés-Rodríguez, and P. Djurić. A comparison of clipping strategies for importance sampling. In *IEEE Statistical Signal Processing Workshop (SSP)*, pages 558–562, 2018.

[63] 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.

[64] T. Minka. Divergence measures and message passing. Technical report, 2005.

[65] P. D. Moral, A. Doucet, and A. Jasra. Sequential Monte Carlo samplers. *Journal of the Royal Statistical Society: Series B*, 68(3):411–436, 2006.

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

[67] M. Mukkamala, P. Ochs, T. Pock, and S. Sabach. Convex-concave backtracking for inertial Bregman proximal gradient algorithms in nonconvex optimization. *SIAM Journal on Mathematics of Data Science*, 2(3):658–682, 2020.

[68] R. Neal. Annealed importance sampling. *Statistics and Computing*, 11:125–139, 2001.

[69] F. Nielsen and R. Nock. Entropies and cross-entropies of exponential families. In *IEEE International Conference on Image Processing (ICIP)*, pages 3621–3624, 2010.

[70] J.-P. Penot. *Calculus without Derivatives*. Springer, 2013.

[71] R. Ranganath, S. Gerrish, and D. Blei. Black box variational inference. In *International Conference on Artificial Intelligence and Statistics (AISTATS)*, volume 33, pages 814–822, 2014.

[72] A. Saha, K. Barath, and S. Kurtek. A geometric variational approach to Bayesian inference. *Journal of the American Statistical Association*, 115(530):822–835, 2020.

[73] Y. Shao, Y. Zhou, and D. Cai. Variational inference with graph regularization for image annotation. *ACM Transactions on Intelligent Systems and Technology*, 2(2):1–21, 2011.

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

[75] M. Titsias and M. Lázaro-Gredilla. Local expectation gradients for black box variational inference. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 28, 2015.

[76] T. van Erven and P. Harremoes. Rényi divergence and Kullback-Leibler divergence. *IEEE Transactions of Information Theory*, 60(7):3797–3820, 2014.

[77] A. Vehtari, D. Simpson, A. Gelman, Y. Yao, and J. Gabry. Pareto smoothed importance sampling. https://arxiv.org/abs/1507.02646, 2015.

[78] S. Vempala and A. Wibosono. Rapid convergence of the unadjusted Langevin algorithm: Isoperimetry suffices. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 32, 2019.
[79] X. Xiao. A unified convergence analysis of stochastic Bregman proximal gradient and extra-gradients methods. *Journal of Optimization Theory and Applications*, 188(3):605-627, 2021.

[80] R. Yao and Y. Yang. Mean field variational inference via Wasserstein gradient flow. https://arxiv.org/abs/2207.08074, 2022.

[81] C. Zhang, J. Bütepage, H. Kjellström, and S. Mandt. Advances in variational inference. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 41(8):2008-2026, 2019.

[82] Y. Zheng, A. Fraysse, and T. Rodet. Efficient unsupervised variational Bayesian image reconstruction using a sparse gradient prior. *Neurocomputing*, 359:449-465, 2019.
[79] X. Xiao. A unified convergence analysis of stochastic Bregman proximal gradient and extra-gradients methods. *Journal of Optimization Theory and Applications*, 188(3):605-627, 2021.

[80] R. Yao and Y. Yang. Mean field variational inference via Wasserstein gradient flow. https://arxiv.org/abs/2207.08074, 2022.

[81] C. Zhang, J. Bütepage, H. Kjellström, and S. Mandt. Advances in variational inference. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 41(8):2008-2026, 2019.

[82] Y. Zheng, A. Fraysse, and T. Rodet. Efficient unsupervised variational Bayesian image reconstruction using a sparse gradient prior. *Neurocomputing*, 359:449-465, 2019.