arXiv:2604.18310v1 [stat.ML] 20 Apr 2026
Symmetry Guarantees Statistic Recovery in Variational Inference Daniel Marks Department of Computer Science University of Oxford [email protected]
Dario Paccagnan Department of Computing Imperial College London [email protected]
Mark van der Wilk Department of Computer Science University of Oxford [email protected]
Abstract Variational inference (VI) is a central tool in modern machine learning, used to approximate an intractable target density by optimising over a tractable family of distributions. As the variational family cannot typically represent the target exactly, guarantees on the quality of the resulting approximation are crucial for understanding which of its properties VI can faithfully capture. Recent work has identified instances in which symmetries of the target and the variational family enable the recovery of certain statistics, even under model misspecification. However, these guarantees are inherently problem-specific and offer little insight into the fundamental mechanism by which symmetry forces statistic recovery. In this paper, we overcome this limitation by developing a general theory of symmetryinduced statistic recovery in variational inference. First, we characterise when variational minimisers inherit the symmetries of the target and establish conditions under which these pin down identifiable statistics. Second, we unify existing results by showing that previously known statistic recovery guarantees in location–scale families arise as special cases of our theory. Third, we apply our framework to distributions on the sphere to obtain novel guarantees for directional statistics in von Mises–Fisher families. Together, these results provide a modular blueprint for deriving new recovery guarantees for VI in a broad range of symmetry settings.
1
Introduction
One of the central problems of modern statistics and machine learning is the accurate estimation of intractable probability distributions. It most prominently arises in the context of Bayesian statistics, where the intractability of the marginal likelihood renders exact computation of the posterior distribution infeasible in all but the simplest statistical models. This has led to the development of a range of approximate inference methods to approximate or sample from such intractable distributions. Variational inference (VI) has emerged as a popular approximate inference technique (Jordan et al., 1999; Wainwright and Jordan, 2008; Blei et al., 2017), balancing the scalability of simpler methods such as Laplace approximations (MacKay, 1992) with the fidelity of more computationally intensive approaches such as MCMC (Robert and Casella, 2004). At its core, VI recasts approximate inference as an optimisation problem: given an intractable target distribution, it seeks the best approximation to it within a tractable variational family by minimising a measure of dissimilarity between the two. The most commonly employed objective is the reverse Kullback–Leibler (KL) divergence, Preprint.
for which a readily computable optimisation surrogate is provided by the evidence lower bound. However, the principles of VI are agnostic to the particular choice of divergence, and alternative objectives have been successfully explored, including the forward KL divergence (Naesseth et al., 2020), α-divergences (Hernandez-Lobato et al., 2016; Dieng et al., 2017; Daudel et al., 2023), α–β divergences (Regli and Silva, 2018), and Rényi divergences (Li and Turner, 2016). A crucial drawback of variational methods lies in the fact that the variational family may not contain the target distribution, in which case the resulting approximation can be arbitrarily bad. This occurs, for example, when the variational family cannot match the target’s number of modes, or is restricted to factorised distributions, as in mean-field variational inference (Peterson and Anderson, 1987; Hinton and Van Camp, 1993), despite the target exhibiting strong dependencies. As a result, guarantees on the quality of VI’s estimates are crucial for assessing the reliability of the resulting approximations and understanding when and why they may be relied upon in practice. In this work, we derive formal guarantees for the statistic recovery of variational approximations, even under severe model misspecification, in the presence of symmetries of the target distribution. Our goal is to develop a general framework for symmetry propagation in variational inference, together with a precise correspondence between symmetries and recoverable statistics, that allows us to answer the following open problem first articulated by Margossian and Saul (2025a): For a given symmetry of the target distribution, which statistics does variational inference recover exactly? Main contributions 1. We prove a general statistic recovery theorem showing that if the target is invariant under a group of symmetries that the variational family respects, a unique variational minimiser inherits this invariance and thus recovers all symmetry-determined statistics. 2. We show that the recovery results of Margossian and Saul (2025a,b) in location–scale families under even and elliptical symmetry arise as special cases of our general framework. 3. We apply our theory to directional statistics on the sphere and derive novel guarantees for recovering the axis of symmetry in von Mises–Fisher families under rotational symmetry. Our results are formulated at the level of f -divergences, thereby encompassing many of the aforementioned variational objectives, and are agnostic to the particular variational inference scheme used to optimise them (Hoffman et al., 2013; Ranganath et al., 2014; Kucukelbir et al., 2017). Conceptually, our theory identifies the fundamental mechanism by which symmetry propagates through variational inference and provides a modular framework for deriving new statistic recovery guarantees in symmetry settings beyond those considered here.
2
Related work
A large body of work has concentrated on asymptotic guarantees for variational approximations. Wang and Blei (2019) proved a variational Bernstein–von Mises theorem, establishing frequentist consistency of variational estimators. Pati et al. (2018) and Yang et al. (2020) proved convergence rates for point estimates constructed from variational approximations with reverse KL and α-divergence objectives, while Zhang and Gao (2020) proved convergence rates for entire posteriors in nonparametric settings. Alquier and Ridgway (2020) similarly derived such rates for tempered posteriors. A parallel line of work studies non-asymptotic guarantees for variational methods. Han and Yang (2019) derived explicit bounds on the KL divergence between a mean-field variational posterior and a normal distribution centred at the maximum likelihood estimator, and showed that the variational mean matches the MLE up to higher-order terms. Motivated by empirical evidence of VI’s ability to recover the mean (MacKay, 2003; Giordano et al., 2018), subsequent work focused on the quality of variational estimates of the target’s statistics. In particular, Katsevich and Rigollet (2024) derived non-asymptotic bounds on the mean and covariance error in Gaussian VI in terms of total variation. Most closely related to our work is that of Margossian and Saul (2025a,b), who showed that, provided the variational minimiser is unique, VI in location–scale families exactly recovers the mean when the target and variational family exhibit even symmetry, and additionally recovers the covariance up to a 2
multiplicative constant when they exhibit elliptical symmetry. Our work generalises these results by moving beyond recovery guarantees tied to specific symmetry classes, distributional families, and statistics, and instead developing a unified symmetry-based framework for statistic recovery.
3
A symmetry-based framework for statistic recovery
Let (X , B) be a measurable space, and let P denote the set of associated probability measures. VI starts with a target distribution P ∈ P, and seeks the best approximation to it inside a tractable variational family Q ⊆ P by minimising a statistical divergence D : P × P → [0, ∞]: Q⋆ := argmin D(P∥Q).
(⋆)
Q∈Q
Throughout this paper, we consider f -divergences Df , a broad class of discrepancy measures that encompasses many of the objectives used in variational inference, including the usual reverse KL divergence.1 In most practical applications, P ∈ / Q, meaning that the variational family is not rich enough to contain the true target and the variational problem is misspecified. For example, P may be a complicated multimodal distribution and Q a simple family of Gaussian distributions. Our goal in this section is therefore not necessarily to recover the full target distribution. Instead, we wish to characterise when, despite misspecification, VI can faithfully recover different statistics of the target. We model a statistic as a partial map S : P ⇀ Y, where Y is an output space. This reflects the fact that many statistics are only well-defined on a natural subclass of P, which we call the domain of S and denote by dom(S) ⊆ P. For instance, for X = Rd , the mean is a statistic mapping distributions with finite first moment to vectors in Rd , while the covariance is a statistic mapping distributions with finite second moment to positive-semidefinite d × d matrices. Since recovery statements are only meaningful when the statistic is well-defined, we assume that P ∈ dom(S) and Q ⊆ dom(S). We are particularly interested in statistic recovery guarantees that arise from symmetries, such as coordinate permutations, rotations, or sign flips. Such symmetries emerge from transformations of the sample space itself, and in our setting it is therefore natural to represent them by a group G of measurable bijections g : X → X with measurable inverses — a measure-theoretic generalisation of a change of coordinates — so that they may be composed and inverted.2 Each g ∈ G acts on any probability measure π ∈ P via the pushforward map g# π(A) := π(g −1 (A)),
A ∈ B.
Equivalently, if X ∼ π, then g(X) ∼ g# π. We say that π ∈ P is G-invariant if g# π = π for all g ∈ G, that is, if applying any symmetry in the group leaves the distribution unchanged. For example, a centred isotropic Gaussian in Rd is invariant under the action of the orthogonal group O(d) of real d × d orthogonal matrices, since it is left unchanged by rotations and reflections about the origin. Similarly, we say that a collection F ⊆ P is G-stable if for all g ∈ G, we have g# F ⊆ F, meaning that the pushforward of any distribution in F also lies in F. For example, the Gaussian family in Rd is stable under the action of the general linear group GL(d, R) of real d × d invertible matrices, since every linear transformation of a Gaussian random variable is also a Gaussian random variable. It turns out that symmetries of the target distribution can constrain the values that its statistics may take, and these constraints can be inherited by a variational minimiser, thus guaranteeing statistic recovery even when P ∈ / Q. To make this possible, we impose three structural conditions describing how the target, the variational family, and the statistic must interact with G: Conditions (1) The target P ∈ dom(S) is G-invariant. (2) The variational family Q ⊆ dom(S) is G-stable. (3) The statistic S : P ⇀ Y has G-stable domain, and for all g ∈ G and π1 , π2 ∈ dom(S), S(π1 ) = S(π2 )
=⇒
S(g# π1 ) = S(g# π2 ).
1 The definition and further background are deferred to Appendix A.1. 2 Formally, G is a subgroup of the group of measurable automorphisms Aut(X , B).
3
The final condition requires that if two distributions have the same statistic, then their pushforwards by any transformation in G also have the same statistic. This allows us to study the effect of symmetries of a distribution directly at the level of its statistic values. To this end, fix y ∈ Im(S). By definition, there exists at least one distribution π ∈ dom(S) such that S(π) = y. Now, for any g ∈ G, consider the quantity S(g# π), which is well-defined since dom(S) is assumed to be G-stable. At first glance, this may seem to depend on the particular choice of π that generated y, since multiple distributions can share the same statistic. However, Condition (3) ensures precisely that if π1 and π2 both satisfy S(π1 ) = S(π2 ) = y then S(g# π1 ) = S(g# π2 ). We may thus introduce, for every g ∈ G, a well-defined map ρg : Im(S) → Im(S) given by ρg (y) := S(g# π),
for any π such that S(π) = y.
Since our ultimate goal is statistic recovery, we are naturally led to consider the statistic values pinned down by the symmetries in G. We collect these in the fixed set ΓS (G) := {y ∈ Im(S) | ρg (y) = y, ∀g ∈ G} . We are now ready to state our main statistic recovery result. At a high level, Conditions (1) and (2) ensure that pushing forward a variational minimiser by any element of G yields another minimiser. If the minimiser is unique, it must therefore inherit the G-invariance of the target. Condition (3) then lifts this symmetry from distributions to statistic values, yielding the desired recovery guarantee. Theorem 3.1 (Statistic recovery) Let S : P ⇀ Y be a statistic, P ∈ dom(S) a target distribution, and Q ⊆ dom(S) a variational family. Assume the variational problem (⋆) has a unique minimiser, denoted by Q⋆ . If Conditions (1)–(3) are satisfied for some group of symmetries G, then S(P), S(Q⋆ ) ∈ ΓS (G). A proof can be found in Appendix A.2. Theorem 3.1 shows that even when P ∈ / Q, variational inference in a G-stable family against a G-invariant target forces every symmetry-compatible statistic of the target and a unique variational minimiser to lie in the corresponding fixed set. When the fixed set is a singleton, exact recovery of the statistic is achieved. In general, however, the fixed set may instead be a low-dimensional subset of the codomain, in which case the recovery is only partial. We will see examples of both cases in Section 4. Remark 3.2 In practice, variational inference is carried out over a parameterised variational family Q = {ϕ(θ) | θ ∈ Θ}, where Θ is a set of parameters and ϕ : Θ → P is a parameterisation map, by minimising Df (P∥ϕ(θ)) over Θ. While uniqueness of the minimising distribution need not imply uniqueness of the minimising parameter due to non-identifiability of the parameterisation, the converse does hold. As parameter-level uniqueness is often easier to verify in concrete models, this is the form in which we will use Theorem 3.1 in later sections.
4
Location–scale families: guarantees under even and elliptical symmetry
In this section, we show that our framework generalises existing statistic recovery guarantees for VI in location–scale families. In particular, we demonstrate that the recovery results of Margossian and Saul (2025a,b) under even and elliptical symmetry arise as direct corollaries of Theorem 3.1, obtained by verifying conditions (1)–(3) and computing the corresponding fixed set in each setting.3 We begin by fixing the common setup and notation for this section. Throughout, we work on X = Rd with d ≥ 1, equipped with its Borel σ-algebra B, and denote by P the associated set of probability measures. For n ∈ N, we let Pn ⊆ P denote the set of probability measures on Rd with finite n-th d d moment. We also write S+ and S++ for the sets of real d × d positive-semidefinite and positivedefinite matrices, respectively. In all examples, we consider variational families of location–scale form and reason about the recovery of the mean statistic µ : P ⇀ Rd with domain dom(µ) = P1 , d and the covariance statistic Σ : P ⇀ S+ with domain dom(Σ) = P2 , in different symmetry settings. Definition 4.1 (Location–scale family) Fix a base distribution Q0 ∈ P. The location–scale family d generated by Q0 is the collection of distributions Q := {Qν,S | ν ∈ Rd , S ∈ S++ } where Qν,S := (Tν,S )# Q0 ,
with
Tν,S (x) := ν + S 1/2 x.
3 See Appendix B.1 for a comparison between our presentation and that of Margossian and Saul (2025a,b).
4
Equivalently, if X0 ∼ Q0 , then ν +S 1/2 X0 ∼ Qν,S . We note that since each Qν,S is the pushforward of Q0 under an invertible affine map, Q ⊆ Pn if and only if Q0 ∈ Pn , meaning that the mean and covariance are well-defined for the whole variational family whenever they are well-defined for the base distribution. 4.1
Mean recovery under even symmetry
We first consider the setting in which the target and variational family exhibit even symmetry, and show that a unique variational minimiser exactly recovers the mean statistic. Definition 4.2 (Even symmetry) A distribution π ∈ P is even symmetric about a ∈ Rd if π(A) = π(2a − A) for all A ∈ B, := where 2a − A {2a − x | x ∈ A}. To obtain the desired recovery guarantee, we specialise our framework to the even symmetric setting. Fix m ∈ Rd and let P ∈ P1 be even symmetric about m. Further, let Q ⊆ P1 be the location– scale family generated by a base distribution Q0 ∈ P1 that is even symmetric about 0.4 We apply Theorem 3.1 to the mean statistic µ under the symmetry group G := {e, rm }, where e denotes the identity map and ra (x) := 2a − x for all a ∈ Rd . We begin by verifying Conditions (1)–(3): −1 (1) Trivially, we know that e# P = P. Moreover, since rm = rm , we have for every A ∈ B, −1 (rm )# P(A) = P(rm (A)) = P(rm (A)) = P(A), where the last equality follows from the even symmetry of P about m. Thus P is G-invariant. d (2) For any (ν, S) ∈ Rd × S++ , we have e# Qν,S = Qν,S ∈ Q. Moreover, it is easy to verify that rm ◦ Tν,S = T2m−ν,S ◦ r0 , and therefore (rm )# Qν,S = (rm ◦ Tν,S )# Q0 = (T2m−ν,S ◦ r0 )# Q0 (1) = (T2m−ν,S )# (r0 )# Q0 = Q2m−ν,S ∈ Q, where the last equality follows from the even symmetry of Q0 about 0. Hence, Q is G-stable. (3) First, the domain dom(µ) = P1 is G-stable. Indeed, if π ∈ P1 , then e# π = π ∈ P1 and, moreover, (rm )# π ∈ P1 , since Z Z Z ∥x∥ d((rm )# π)(x) = ∥2m − x∥ dπ(x) ≤ 2∥m∥ + ∥x∥ dπ(x) < ∞. Rd
Rd
Rd
Further, for every π ∈ P1 , we can compute µ(g# π) for every g ∈ G as µ(π), g = e, µ(g# π) = 2m − µ(π), g = rm .
(2)
As the right-hand side of (2) depends only on µ(π), it follows that if µ(π1 ) = µ(π2 ), then µ(g# π1 ) = µ(g# π2 ) for all π1 , π2 ∈ P1 and all g ∈ G. We have thus verified the conditions of Theorem 3.1. Moreover, by (2), Condition (3) allows us to define the maps ρe (x) = x and ρrm (x) = 2m − x on Im(µ) = Rd and compute the corresponding fixed set Γµ (G) = {x ∈ Rd | ρg (x) = x, ∀g ∈ G} = {x ∈ Rd | 2m − x = x} = {m}. Hence, Theorem 3.1 guarantees exact recovery of the mean when specialised to the even-symmetric setting, provided the variational minimiser is unique. We record this result in the following corollary, which coincides with the mean recovery theorem of Margossian and Saul (2025b). Corollary 4.3 (cf. Theorem 10 in Margossian and Saul (2025b)) Fix m ∈ Rd and let P ∈ P1 be even symmetric about m. Let Q0 ∈ P1 be even symmetric about 0, and let Q be the location–scale d family it generates. If Df (P∥Qν,S ) has a unique minimiser (ν ⋆ , S ⋆ ) over Rd × S++ , then µ(Qν ⋆ ,S ⋆ ) = µ(P) = m. The left panel of Figure 1 illustrates this result. It is worth noting that in (1), the scale matrix remains unchanged under pushforward by elements of G. Corollary 4.3 therefore also applies mutatis mutandis to VI over a location family, i.e., a location–scale family with fixed scale matrix. 4 This is equivalent to requiring that each Q ν,S ∈ Q is even symmetric about its location parameter ν.
5
4.2
Mean and covariance recovery under elliptical symmetry
We next consider the setting where the target and variational family exhibit elliptical symmetry. Here, we show that a unique variational minimiser exactly recovers the mean and further recovers the covariance matrix up to a multiplicative constant. This, in turn, yields exact recovery of the correlation matrix. Definition 4.4 (Elliptical symmetry) A distribution π ∈ P is elliptically symmetric about γ ∈ Rd d with scale matrix C ∈ S++ if there exists an O(d)-invariant distribution π0 ∈ P such that Tγ,C (x) := γ + C 1/2 x.
π = (Tγ,C )# π0 ,
Here, O(d) is the orthogonal group, which consists of real d × d orthogonal matrices. Equivalently, if X0 ∼ π0 , then γ + C 1/2 X0 ∼ π. We proceed, as before, by specialising Theorem 3.1 to this new d setting. Let P ∈ P be elliptically symmetric about m ∈ Rd with scale matrix M ∈ S++ , and let Q ⊆ P be the location–scale family generated by an O(d)-invariant base distribution Q0 ∈ P.5 Mean recovery. Observe that a distribution that is elliptically symmetric about γ ∈ Rd is also even symmetric about γ. Indeed, if π = (Tγ,C )# π0 for some O(d)-invariant π0 , then for every A ∈ B −1 −1 −1 π(2γ − A) = π0 (Tγ,C (2γ − A)) = π0 (−Tγ,C (A)) = π0 (Tγ,C (A)) = π(A),
where the penultimate step follows from the fact that O(d)-invariance of π0 implies (−Id )# π0 = π0 . Therefore, if P, Q0 ∈ P1 , mean recovery is immediate by Corollary 4.3. Covariance recovery. To reason about covariance recovery, we assume that P, Q0 ∈ P2 and consider the symmetry group of ellipsoids centred at m with shape matrix M , defined as G := {gR | R ∈ O(d)},
AR := M 1/2 RM −1/2 .
gR (x) := m + AR (x − m),
We begin by verifying Conditions (1)–(3) under the action of this group: (1) By elliptical symmetry, there exists an O(d)-invariant P0 ∈ P2 such that P = (Tm,M )# P0 . It is easy to check that gR ◦ Tm,M = Tm,M ◦ R for every R ∈ O(d), and therefore (gR )# P = (gR ◦ Tm,M )# P0 = (Tm,M ◦ R)# P0 = (Tm,M )# (R# P0 ) = P. d (2) Let (ν, S) ∈ Rd × S++ and R ∈ O(d). Define
ν ′ := m + AR (ν − m),
S ′ := AR SA⊤ R,
U := (S ′ )−1/2 AR S 1/2 .
Then U ∈ O(d) because U U ⊤ = Id , and since Q0 is O(d)-invariant, we have U# Q0 = Q0 . Moreover, we have that gR ◦ Tν,S = Tν ′ ,S ′ ◦ U and therefore (gR )# Qν,S = (gR ◦ Tν,S )# Q0 = (Tν ′ ,S ′ ◦ U )# Q0 = (Tν ′ ,S ′ )# (U# Q0 ) = Qν ′ ,S ′ ∈ Q. (3) First, the domain P2 is G-stable. Indeed, for all x ∈ Rd ∥gR (x)∥ ≤ ∥m∥ + ∥AR ∥op (∥x∥ + ∥m∥), 2
2
(3)
2
and squaring (3), together with the identity (a + b) ≤ 2a + 2b , yields that for any π ∈ P2 Z Z ∥x∥2 d((gR )# π)(x) ≤ 2∥AR ∥2op ∥x∥2 dπ(x) + 2(1 + ∥AR ∥op )2 ∥m∥2 < ∞, Rd
Rd
since ∥AR ∥op < ∞ by the submultiplicativity of the operator norm. Moreover, for every π ∈ P2 and every R ∈ O(d), we have µ((gR )# π) = gR (µ(π)) and so Z ⊤ Σ((gR )# π) = gR (x) − gR (µ(π)) gR (x) − gR (µ(π)) dπ(x) d R Z (4) = AR (x − µ(π))(x − µ(π))⊤ dπ(x) A⊤ R Rd
= AR Σ(π) A⊤ R, meaning that if Σ(π1 ) = Σ(π2 ), then Σ((gR )# π1 ) = Σ((gR )# π2 ) for all R ∈ O(d). 5 This is equivalent to requiring that each Q ν,S ∈ Q is elliptically symmetric about ν with scale matrix S.
6
p
p
0.8
q mean
−1
0 marginal
mean
1
0
x2
density
0.6
q
2
0.4
−2 0.2 −4 0.0
−1
0
1
2
3
4
0
1
x
2
3
4
x1
Figure 1: Unique variational minimisers q from a Gaussian family exactly recover symmetrydetermined statistics of highly non-Gaussian targets p. Left: under even symmetry, the mean is recovered exactly as per Corollary 4.3. Right: under elliptical symmetry, the mean and correlation are recovered exactly, whereas the covariance is recovered only up to scale, as per Corollary 4.5. By (4), Condition (3) allows us to define for every gR ∈ G the map ρgR (V ) = AR V A⊤ R on d Im(Σ) = S+ and find the positive-semidefinite matrices satisfying ρgR (V ) = V for all gR ∈ G: d ΓΣ (G) = {V ∈ S+ | AR V A⊤ R = V, ∀R ∈ O(d)} = {λM | λ ≥ 0}. d The last equality follows from the fact that a matrix V ∈ S+ lies in the fixed set if and only if −1/2 −1/2 M VM is invariant under conjugation by every R ∈ O(d), since AR = M 1/2 RM −1/2 , which in turn implies that M −1/2 V M −1/2 = λId for some λ ≥ 0. It follows that if P is not a Dirac measure, then Σ(P) = αM for some α > 0. Similarly, if Q0 is not a Dirac measure and (ν ⋆ , S ⋆ ) d is a unique minimiser of the variational objective over Rd × S++ , then Σ(Qν ⋆ ,S ⋆ ) = bM for some b > 0. Therefore, outside the trivial point-mass case, a unique variational minimiser recovers the covariance up to a positive multiplicative constant.
Correlation recovery. We conclude this section by noting an immediate implication of covariance recovery for the correlation statistic. To that end, let P2◦ ⊆ P2 denote the set of distributions with finite second moment and positive marginal variances, and assume P, Q0 ∈ P2◦ . Since Q0 is O(d)-invariant, its covariance is equal to a positive scalar multiple of the identity, which guarantees that Q ⊆ P2◦ . From covariance recovery, we know that Σ(P) = aM for some a > 0, and Σ(Qν ⋆ ,S ⋆ ) = bM for some b > 0. To remove this scale ambiguity, we may define for every π ∈ P2◦ the correlation statistic ρ(π) := D(π)−1/2 Σ(π)D(π)−1/2 ,
D(π) := diag(Σ(π)).
Then D(P) = a diag(M ) and D(Qν ⋆ ,S ⋆ ) = b diag(M ), and thus ρ(P) = diag(M )−1/2 M diag(M )−1/2 = ρ(Qν ⋆ ,S ⋆ ). Therefore, the correlation is recovered exactly by a unique minimiser. We summarise our discussion on statistic recovery under elliptical symmetry in the following corollary, which coincides with Theorem 11 of Margossian and Saul (2025b). An example is provided in the right panel of Figure 1. d Corollary 4.5 (cf. Theorem 11 in Margossian and Saul (2025b)) Fix m ∈ Rd and M ∈ S++ , and let P be elliptically symmetric about m with scale matrix M . Let Q0 be an O(d)-invariant base distribution, and let Q be the location–scale family it generates. If Df (P∥Qν,S ) has a unique d minimiser (ν ⋆ , S ⋆ ) over Rd × S++ , then
(i) if P, Q0 ∈ P1 , then µ(Qν ⋆ ,S ⋆ ) = µ(P) = m, (ii) if P, Q0 ∈ P2 and are not Dirac measures then Σ(Qν ⋆ ,S ⋆ ), Σ(P) ∈ {λM | λ > 0}, (iii) if in addition P, Q0 ∈ P2◦ , then ρ(Qν ⋆ ,S ⋆ ) = ρ(P) = diag(M )−1/2 M diag(M )−1/2 . Remark 4.6 Note that although the conditions of Theorem 3.1 are sufficient for statistic recovery, they are not always necessary. Indeed, correlation is recovered in Corollary 4.5 through the covariance, even though it does not itself satisfy Condition (3) for d ≥ 2; see Appendix B.2 for a discussion. 7
5
von Mises–Fisher family: new guarantees under rotational symmetry
So far, we have demonstrated that previously known statistic recovery results arise as special instances of our general framework. We now use this framework to obtain a new recovery guarantee for variational inference on the sphere. Specifically, we show that when the target exhibits rotational symmetry, variational inference in a von Mises–Fisher family can exactly recover the axis of symmetry. Unlike the location–scale examples, where the relevant symmetry, variational family, and recovered statistic were known in advance, here we show how our framework can be applied more organically: we start from a symmetry of the target distribution, determine a statistic that is pinned down by it, and purposely select a stable variational family in order to guarantee its recovery. This section therefore not only establishes a new result, but also provides a blueprint illustrating how our framework can be applied to new symmetries and statistics beyond the examples considered here. Let X = S d−1 be the unit sphere in Rd , d ≥ 3, equipped with its Borel σ-algebra B, and let P denote the associated set of probability measures. To ease the presentation, we work with distributions that are absolutely continuous with respect to the uniform probability measure σ on the unit sphere, which we collect inside the set Pσ ⊆ P. For any π ∈ Pσ , we denote its density with respect to σ by fπ . Definition 5.1 (Rotational symmetry) A distribution π ∈ Pσ is rotationally symmetric about axis direction w ∈ S d−1 if there exists a measurable function ψ : [−1, 1] → R≥0 such that fπ (x) = ψ(w⊤ x),
for σ-a.e. x ∈ S d−1 .
We call ψ the axial profile of π. If ψ is constant, then π is the uniform distribution on S d−1 . Suppose we are given a non-uniform target P ∈ Pσ that is rotationally symmetric about u ∈ S d−1 . Observe that the level sets of the map x 7→ u⊤ x are precisely the parallels orthogonal to u. Thus, a rotationally symmetric distribution may vary only with the polar angle away from u. In particular, it is invariant under the group G := {gR | Ru = u, R ∈ SO(d)}, where SO(d) is the special orthogonal group of orthogonal d × d matrices with determinant 1, and gR (x) = Rx for all x ∈ S d−1 . Indeed, (gR )# P = P for all gR ∈ G, since by rotation invariance of σ we have for σ-a.e. x ∈ S d−1 , f(gR )# P (x) = fP (R−1 x) = ψ(u⊤ R−1 x) = ψ((Ru)⊤ x) = ψ(u⊤ x) = fP (x). Given this target symmetry, we now ask which statistic, if any, it pins down. A natural first candidate is the axis direction. However, if π ∈ Pσ is rotationally symmetric about w with axial profile ψ, it is also rotationally symmetric about −w with axial profile ψ̃(t) := ψ(−t). Thus, the most the symmetry structure can uniquely determine is a one-dimensional subspace of Rd . More precisely, if π ∈ Pσ is non-uniform and rotationally symmetric about both w1 , w2 ∈ S d−1 , then span(w1 ) = span(w2 ); proof in Appendix C.1. We therefore consider the space of lines through the origin in Rd denoted by RP d−1 , let PA ⊆ Pσ be the set of all non-uniform rotationally symmetric distributions, and define the axis statistic A : P ⇀ RP d−1 on dom(A) = PA by A(π) = span(w), where w is any axis direction of rotational symmetry of π. PA is G-stable, since for any π ∈ PA rotationally symmetric about w with axial profile ψ and any gR ∈ G, f(gR )# π (x) = fπ R−1 x = ψ(w⊤ R−1 x) = ψ (Rw)⊤ x , for σ-a.e. x ∈ S d−1 , meaning (gR )# π is rotationally symmetric about Rw with profile ψ and thus (gR )# π ∈ PA . Hence, A((gR )# π) = span(Rw) = R span(w) = RA(π),
(5)
which implies that if A(π1 ) = A(π2 ), then A((gR )# π1 ) = A((gR )# π2 ) for all gR ∈ G. By (5), we can define for every gR ∈ G the map ρgR (L) := RL on Im(A) = RP d−1 to obtain the fixed set ΓA (G) = {L ∈ RP d−1 | ρgR (L) = L, ∀gR ∈ G} = {span(u)}. By construction, P ∈ PA . We thus know that choosing a G-stable family will guarantee the recovery of the axis statistic by a unique minimiser. One such family is the von Mises–Fisher family. Each of its members is rotationally symmetric about a mean direction ν ∈ S d−1 , and indexed by a concentration parameter κ > 0 controlling how quickly the density decays as the angle from ν increases. Definition 5.2 (von Mises–Fisher family) A von Mises–Fisher (vMF) family Q ⊆ Pσ is a twoparameter family {Qν,κ | ν ∈ S d−1 , κ ∈ R>0 }, where for a normalising constant cd (κ) > 0, fQν,κ (x) = cd (κ) exp(κν ⊤ x), 8
for σ-a.e. x ∈ S d−1 .
Target with η < ηc
VI optimum with η < ηc
Target with η > ηc
ν=u
VI optimum with η > ηc
u
ν
Figure 2: Log-density contours of the target and variational posterior are shown using a Lambert azimuthal equal-area projection centred at the true direction u. Left: η < ηc and the reverse KL has a unique minimiser at ν = u, recovering the symmetry axis. Right: η > ηc and the minimisers form a latitude circle. None of them, including the one found during optimisation, recovers the statistic. Indeed, the von Mises–Fisher family is G-stable. For all (ν, κ) ∈ S d−1 × R>0 and all gR ∈ G, (gR )# Qν,κ = QRν,κ ∈ Q since by rotation invariance of σ, we have that for σ-a.e. x ∈ S d−1 f(gR )# Qν,κ (x) = cd (κ) exp κν ⊤ R−1 x = cd (κ) exp κ(Rν)⊤ x = fQRν,κ (x), (6) and since κ > 0, it follows that Q ⊆ PA . We thus obtain the following new recovery guarantee. Corollary 5.3 Let P ∈ PA be rotationally symmetric about u ∈ S d−1 , and let Q ⊆ PA be the von Mises–Fisher family. If Df (P∥Qν,κ ) has a unique minimiser (ν ⋆ , κ⋆ ) over S d−1 × R>0 , then A(Qν ⋆ ,κ⋆ ) = A(P) = span(u). Analogously to location and location–scale families, it is natural to consider a fixed-concentration von Mises–Fisher family, namely, the one-parameter subfamily of distributions that share the same concentration parameter κ0 > 0. Since the pushforward action leaves the concentration parameter unchanged as per (6), Corollary 5.3 equally applies to fixed-concentration vMF families. Remark 5.4 As with all corollaries of Theorem 3.1, uniqueness of the minimiser is necessary for the recovery to be guaranteed. In what follows, we highlight this with a concrete example. Fix λ ̸= 0 and u ∈ S d−1 and let (Pλ,η )η>0 ⊆ PA be a one-parameter family of target distributions such that for each η > 0, Pλ,η is rotationally symmetric about u with axial profile −1 ψλ,η (t) := Zλ,η exp λt − ηt2 , t ∈ [−1, 1], where Zλ,η > 0 is a normalising constant. Let Qκ0 be a fixed-concentration von Mises–Fisher family for some κ0 > 0 and consider the reverse KL objective DKL (Qν,κ0 ∥Pλ,η ). Since η > 0, Pλ,η is not a vMF distribution and the model is misspecified. Despite that, there exists a critical threshold ηc (d, λ, κ0 ) > 0 such that for all 0 < η ≤ ηc , the objective DKL (Qν,κ0 ∥Pλ,η ) has a unique minimiser ν ⋆ over S d−1 , and exact recovery holds with A(Qν ⋆ ,κ0 ) = A(Pλ,η ). By contrast, for every η > ηc , the minimiser set is non-singleton and no minimiser recovers the true axis statistic. This behaviour is visualised for d = 3 in Figure 2, with an explicit construction given in Appendix C.2.
6
Discussion
In this work, we have developed a general symmetry-based theory of statistic recovery in variational inference, and have used it both to recover existing results and to derive novel guarantees. With our framework providing a modular blueprint for establishing statistic recovery, we see significant opportunities to apply it across different domains, symmetry classes, and variational families. At the same time, we identify at least two important directions for future work. First, while our theory provides sufficient conditions for statistic recovery, we have shown that these conditions are not always necessary. Deriving necessary and sufficient conditions would therefore yield a complete characterisation of statistic recovery under symmetry in variational problems. Second, it would be valuable to extend the theory to settings where the target is only approximately symmetric. Such an extension could unlock new guarantees for a broad class of inference problems that exhibit substantial structure but do not exactly satisfy our assumptions — a direction we leave for future work. 9
Acknowledgments and Disclosure of Funding D. Marks acknowledges funding from a G-Research graduate scholarship. D. Paccagnan was partially supported by the EPSRC grant EP/Y001001/1, funded by the International Science Partnerships Fund (ISPF) and UKRI; and by an Imperial-MIT seed fund.
References Pierre Alquier and James Ridgway. Concentration of tempered posteriors and of their variational approximations. The Annals of Statistics, 48(3):1475–1497, June 2020. David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational Inference: A Review for Statisticians. Journal of the American Statistical Association, 112(518):859–877, April 2017. Kamélia Daudel, Joe Benton, Yuyang Shi, and Arnaud Doucet. Alpha-divergence Variational Inference Meets Importance Weighted Auto-Encoders: Methodology and Asymptotics. Journal of Machine Learning Research, 24(243):1–83, 2023. Adji Bousso Dieng, Dustin Tran, Rajesh Ranganath, John Paisley, and David Blei. Variational Inference via \textbackslash chi Upper Bound Minimization. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017. E. B. Dynkin. The maximal subgroups of the classical groups. volume 6, pages 245–378. American Mathematical Society, 1957. Ryan Giordano, Tamara Broderick, and Michael I. Jordan. Covariances, robustness and variational bayes. J. Mach. Learn. Res., 19(1):1981–2029, January 2018. Wei Han and Yun Yang. Statistical Inference in Mean-Field Variational Bayes, November 2019. Jose Hernandez-Lobato, Yingzhen Li, Mark Rowland, Thang Bui, Daniel Hernandez-Lobato, and Richard Turner. Black-Box Alpha Divergence Minimization. In Proceedings of The 33rd International Conference on Machine Learning, pages 1511–1520, June 2016. Geoffrey E. Hinton and Drew Van Camp. Keeping the neural networks simple by minimizing the description length of the weights. In Proceedings of the sixth annual conference on Computational learning theory - COLT ’93, pages 5–13, Santa Cruz, California, United States, 1993. ACM Press. Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of machine learning research, 2013. Michael I. Jordan, Zoubin Ghahramani, Tommi S. Jaakkola, and Lawrence K. Saul. An Introduction to Variational Methods for Graphical Models. Machine Learning, 37(2):183–233, November 1999. Anya Katsevich and Philippe Rigollet. On the approximation accuracy of Gaussian variational inference. The Annals of Statistics, 52(4):1384–1409, August 2024. Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David M. Blei. Automatic differentiation variational inference. J. Mach. Learn. Res., 18(1):430–474, January 2017. Yingzhen Li and Richard E Turner. Rényi Divergence Variational Inference. In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016. David J. C. MacKay. Bayesian Interpolation. Neural Computation, 4(3):415–447, May 1992. David J. C. MacKay. Information theory, inference, and learning algorithms. Cambridge University Press, Cambridge, UK ; New York, 2003. Kanti V. Mardia and Peter E. Jupp. Directional Statistics. John Wiley & Sons, September 2009. Charles Margossian and Lawrence K. Saul. Variational Inference in Location-Scale Families: Exact Recovery of the Mean and Correlation Matrix. In Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, pages 3466–3474, April 2025a. 10
Charles C. Margossian and Lawrence K. Saul. Generalized Guarantees for Variational Inference in the Presence of Even and Elliptical Symmetry, December 2025b. Dragoslav S. Mitrinović. Analytic Inequalities. Springer Berlin Heidelberg, Berlin, Heidelberg, 1970. Christian A. Naesseth, Fredrik Lindsten, and David Blei. Markovian score climbing: variational inference with KL(p||q). In Proceedings of the 34th International Conference on Neural Information Processing Systems, pages 15499–15510, Red Hook, NY, USA, December 2020. Curran Associates Inc. Constantin P. Niculescu and Lars-Erik Persson. Convex Functions and Their Applications. CMS Books in Mathematics. Springer International Publishing, Cham, 2018. Debdeep Pati, Anirban Bhattacharya, and Yun Yang. On Statistical Optimality of Variational Bayes. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, pages 1579–1588, March 2018. Carsten Peterson and James R Anderson. A mean field theory learning algorithm for neural networks. Complex Systems, 1:995–1019, 1987. Yury Polyanskiy and Yihong Wu. Information Theory: From Coding to Learning, January 2025. Rajesh Ranganath, Sean Gerrish, and David Blei. Black Box Variational Inference. In Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, pages 814–822, April 2014. Jean-Baptiste Regli and Ricardo Silva. Alpha-Beta Divergence For Variational Inference, May 2018. Christian P. Robert and George Casella. Monte Carlo Statistical Methods. Springer Texts in Statistics. Springer, New York, NY, 2004. Igal Sason and Sergio Verdu. f -Divergence Inequalities. IEEE Trans. Inf. Theor., 62(11):5973–6006, November 2016. Terence Tao. An introduction to measure theory. Number v. 126 in Graduate studies in mathematics. American Mathematical Society, Providence, R.I, 2011. Fuichi Uchida. Classification of real analytic SL(n, R) actions on n-sphere. Osaka Journal of Mathematics, 16(3):561–579, 1979. Martin J. Wainwright and Michael I. Jordan. Graphical Models, Exponential Families, and Variational Inference. Foundations and Trends in Machine Learning, 1(1-2):1–305, December 2008. Yixin Wang and David M. Blei. Frequentist Consistency of Variational Bayes. Journal of the American Statistical Association, 114(527):1147–1161, July 2019. Yun Yang, Debdeep Pati, and Anirban Bhattacharya. $\alpha $-variational inference with statistical guarantees. The Annals of Statistics, 48(2), April 2020. Fengshuo Zhang and Chao Gao. Convergence rates of variational posterior distributions. The Annals of Statistics, 48(4):2180–2207, August 2020.
11
A
Additional material for Section 3
A.1
Overview of f -divergences
In this appendix, we give a formal definition of f -divergences, which underpin our formulation of the variational inference problem throughout the paper. Let (X , B) be a measurable space, and denote by P the set of associated probability measures. Let P, Q ∈ P be two such probability measures, and let µ be any σ-finite measure on (X , B) such that P ≪ µ and Q ≪ µ, e.g. µ = P + Q. Write p = dP/dµ and q = dQ/dµ for the corresponding densities. Let f : R>0 → R be a convex function with f (1) = 0, and define the extensions: 1 f (0) := lim f (t), f ′ (∞) := lim tf . t↓0 t↓0 t The f -divergence generated by f , denoted Df : P × P → [0, ∞], is defined by Z p dµ + f ′ (∞) P({q = 0}), Df (P∥Q) := qf q {q>0} with the convention that 0 · ∞ = 0. This definition follows Polyanskiy and Wu (2025, Def. 7.1).6 A number of popular divergences can be recovered as special cases of f -divergences with different generators. These include the Kullback–Leibler (KL) divergence with f (t) = t log t, the χ2 divergence with f (t) = (t − 1)2 , the total variation distance with f (t) = 12 |t − 1|, and the squared √ Hellinger distance with f (t) = (1− t)2 . For any f -divergence, swapping the order of the arguments again yields an f -divergence, now generated by f˜(t) = tf (1/t). That is, Df (Q∥P) = Df˜(P∥Q). For example, since the forward KL divergence is generated by f (t) = t log t, the reverse KL divergence is generated by f˜(t) = tf (1/t) = − log t. For a detailed overview of f -divergences, see Sason and Verdu (2016). A.2
Proof of Theorem 3.1
In this appendix, we provide a proof of Theorem 3.1. The proof is divided into three parts. In the first part, we show that f -divergences are invariant under simultaneous pushforward by measurable automorphisms. In the second part, we use this fact to show that Conditions (1) and (2) alone suffice to guarantee that a unique variational minimiser inherits the G-invariance of the target. In the third part, we use the maps ρg on statistic values introduced in Section 3, whose well-definedness is guaranteed by Condition (3), and show that when all conditions hold, the statistic of both the target and a unique variational minimiser is constrained inside the fixed set ΓS (G). Setup. We formally recall the setup of Section 3. Let (X , B) be a measurable space, and denote by P the set of probability measures on (X , B). Let Aut(X , B) be the group of measurable automorphisms, acting on P by pushforward, and fix a subgroup G ≤ Aut(X , B). Throughout this section, we consider the variational inference problem (⋆) with f -divergences Q⋆ := argmin Df (P∥Q). Q∈Q
For the first two parts of the proof, which concern only the variational problem and the action of G on distributions, we work with a general P ∈ P and Q ⊆ P. For the third part, we additionally consider a statistic S : P ⇀ Y, with domain dom(S), and ultimately return to the assumption that P ∈ dom(S) and Q ⊆ dom(S) for the last step of the proof of Theorem 3.1. Whenever Conditions (1)–(3) are invoked below, we mean the three conditions stated in Section 3. We begin with the first part of the proof, and show that f -divergences are invariant under simultaneous pushforward by measurable automorphisms. Lemma A.1 (f -divergence invariance) For any P, Q ∈ P and g ∈ Aut(X , B), Df g# P g# Q = Df P Q . 6 Note that Polyanskiy and Wu (2025) work with the standing assumption that (X , B) is a standard Borel space. However,
their f -divergence definition extends verbatim to the case of an arbitrary measurable space, such as the one considered here.
12
Proof . Let µ be a σ-finite common dominating measure of P and Q, write p = dP/dµ and q = dQ/dµ for the corresponding densities, and define µ̃ = g# µ. It is easy to verify that µ̃ is σ-finite and that g# P, g# Q ≪ µ̃. Define p̃ := p ◦ g −1 and q̃ := q ◦ g −1 . Then for every A ∈ B, Z Z Z Z p̃ dµ̃ = 1A (p ◦ g −1 ) dµ̃ = (1A ◦ g) p dµ = p dµ = P(g −1 (A)) = (g# P)(A), X
A
g −1 (A)
X
meaning p̃ = d(g# P)/dµ̃, and similarly q̃ = d(g# Q)/dµ̃. We first reduce to the case f ≥ 0. For any a ∈ R, define ha (t) := f (t) + a(t − 1), t > 0. Then ha is convex, ha (1) = 0, and its extensions satisfy ha (0) = f (0) − a,
h′a (∞) = f ′ (∞) + a.
Moreover, Z
q f (p/q) + a(p/q − 1) dµ + f ′ (∞) + a P({q = 0})
Dha (P∥Q) = {q>0}
"Z = Df (P∥Q) + a
# (p − q) dµ + P({q = 0}) .
{q>0}
Now, by linearity of the integral Z (p − q) dµ + P({q = 0}) = P({q > 0}) − Q({q > 0}) + P({q = 0}), {q>0}
and since Q({q = 0}) = 0 we have Q({q > 0}) = 1. Hence Z (p − q) dµ + P({q = 0}) = 1 − 1 = 0, {q>0}
and therefore Dha (P∥Q) = Df (P∥Q). Now choose α ∈ ∂f (1). By convexity (Niculescu and Persson, 2018), f (t) ≥ f (1) + α(t − 1) = α(t − 1),
t > 0,
so h−α (t) = f (t) − α(t − 1) ≥ 0. Therefore, replacing f by the equivalent generator h−α if necessary, we may assume from the outset that f ≥ 0. With this reduction in hand, define q̃(x) f p̃(x)/q̃(x) , q̃(x) > 0, Φ(x) := 0, q̃(x) = 0. Since f ≥ 0, Φ is nonnegative and measurable, so a change of variables yields (Tao, 2011) Z Z Φ dµ̃ = Φ ◦ g dµ. X
X
Using p̃ ◦ g = p and q̃ ◦ g = q, we get Φ ◦ g = 1{q>0} q f (p/q). Therefore, Z Z q̃ f (p̃/q̃) dµ̃ = q f (p/q) dµ. {q̃>0}
{q>0}
For the singular term, since q̃ = q ◦ g −1 , we have g −1 ({q̃ = 0}) = {q = 0}, and hence (g# P)({q̃ = 0}) = P(g −1 ({q̃ = 0})) = P({q = 0}). Combining the two identities yields Z Df (g# P∥g# Q) =
q f (p/q) dµ + f ′ (∞) P({q = 0}) = Df (P∥Q).
{q>0}
□
This proves the claim.
We now proceed with the second part of the proof, where we isolate the role of Conditions (1) and (2), namely the G-invariance of the target and the G-stability of the variational family, to show that symmetries of the target propagate to unique variational minimisers. 13
Lemma A.2 (Symmetry propagation) Let G ≤ Aut(X , B). Assume P ∈ P is G-invariant, and let Q ⊆ P be G-stable. Then: (i) Q⋆ is G-stable, (ii) if Q⋆ = {Q⋆ } is a singleton, then Q⋆ is G-invariant. Proof . For any g ∈ G, the objective is invariant under simultaneous pushforward by Lemma A.1. Now if P is G-invariant, then g# P = P and hence Df (P∥g# Q) = Df (g# P ∥ g# Q) = Df (P ∥ Q),
(7)
for all Q ∈ Q. Thus the variational objective is constant along G-orbits. (i) If Q⋆ is empty, then it is trivially G-stable. Let Q⋆ ∈ Q⋆ and take g ∈ G. Since Q is G-stable, we have g# Q⋆ ∈ Q. Moreover, by (7), Df (P ∥ g# Q⋆ ) = Df (P ∥ Q⋆ ). Because Q⋆ is a minimiser over Q, the right-hand side equals inf Q∈Q Df (P∥Q), so g# Q⋆ also attains this infimum. Thus g# Q⋆ ∈ Q⋆ , proving that Q⋆ is G-stable. (ii) Assume Q⋆ = {Q⋆ }. By (i), for every g ∈ G we have g# Q⋆ ∈ Q⋆ . Since Q⋆ contains only Q⋆ , it follows that g# Q⋆ = Q⋆ for all g ∈ G, i.e. Q⋆ is G-invariant. □ We conclude the proof by finally turning to Condition (3). As shown in Section 3, this condition guarantees that for each g ∈ G, the map ρg : Im(S) → Im(S) is well defined. We now combine this observation with the preceding lemmas to obtain our statistic recovery guarantee. Theorem A.3 (Statistic recovery, restated) Let S : P ⇀ Y be a statistic, P ∈ dom(S) a target distribution, and Q ⊆ dom(S) a variational family. Assume the variational problem (⋆) has a unique minimiser, denoted by Q⋆ . If Conditions (1)–(3) are satisfied for some group of symmetries G, then S(P), S(Q⋆ ) ∈ ΓS (G). Proof . By Condition (3), for every g ∈ G, the map ρg introduced in Section 3 is well defined on Im(S). Since P is G-invariant by Condition (1), for every g ∈ G we have ρg (S(P)) = S(g# P) = S(P). Hence, by definition of ΓS (G), S(P) ∈ ΓS (G). By assumption, Q = {Q }. By Conditions (1) and (2), Lemma A.2(ii) implies that Q⋆ is G-invariant. Moreover, Q⋆ ∈ Q ⊆ dom(S), so S(Q⋆ ) is well defined. Therefore, for every g ∈ G, ⋆
⋆
ρg (S(Q⋆ )) = S(g# Q⋆ ) = S(Q⋆ ). Again by definition of ΓS (G), it follows that S(Q⋆ ) ∈ ΓS (G). Combining the two inclusions yields S(P), S(Q⋆ ) ∈ ΓS (G). □
B
Additional material for Section 4
B.1
Presentation differences relative to Margossian and Saul
In this appendix, we discuss the differences between our presentation of the recovery guarantees in location–scale families in Section 4 and that of Margossian and Saul (2025a,b), whose results we recover. The first difference is that we choose to work with probability measures rather than densities, 14
in order to follow the general blueprint of Section 3. Accordingly, our definitions of location–scale families, as well as of even and elliptical symmetry, are the measure-theoretic analogues of those presented in Margossian and Saul (2025a,b). This places the results in a more general setting, but at the same time forfeits some of the convenient properties that follow from the density-level formulation. For example, for covariance recovery, in addition to the necessary assumption that the target and the base distribution of the variational family have finite second moments, we explicitly assume that both have non-zero variance — a property that automatically follows from the existence of a density with respect to Lebesgue measure. Any similar minor differences in presentation are simply artefacts of our choice to work more generally with probability distributions rather than just probability densities. The other, and arguably more important, difference in presentation is that, whereas we state all our results under the assumption that the variational minimiser is unique, Margossian and Saul (2025a,b) take the extra step of imposing additional regularity conditions, such as differentiability and somewhere-strict log-concavity of the target, in order to guarantee uniqueness. Because our theory focuses on the fundamental mechanism underpinning symmetry-induced statistic recovery, we are less concerned with the specific assumptions needed to ensure uniqueness in concrete settings. Accordingly, we have chosen to state our results in the most general way possible, retaining uniqueness as a standing assumption. That said, any regularity conditions that guarantee uniqueness of the variational minimiser can be substituted for this assumption in all our corollaries. B.2
Correlation does not satisfy Condition (3)
Here, we demonstrate the problems that arise when trying to apply Theorem 3.1 directly to the correlation statistic. First, we show that the statistic does not satisfy Condition (3) because its domain is not G-stable. We then show that even after restricting the domain to distributions with finite second moment and positive-definite covariance to guarantee G-stability, the correlation statistic still fails to satisfy the condition. d Recall the setting of Section 4 and, in particular, the correlation statistic ρ : P2 ⇀ C+ , defined by
ρ(π) := D(π)−1/2 Σ(π)D(π)−1/2 ,
D(π) := diag(Σ(π)),
(8)
where dom(ρ) = P2◦ is the set of distributions with finite second moment and positive marginal d variances, and C+ is the set of positive semidefinite matrices with unit diagonal entries. Further, recall the relevant symmetry group: G := {gR | R ∈ O(d)},
gR (x) := m + AR (x − m),
AR := M 1/2 RM −1/2 .
Assume d ≥ 2 and choose v ∈ Rd with vi ̸= 0 for every i. Define 1 1 δv + δ−v . 2 2 Then π ∈ P2 since it is finitely supported, and moreover µ(π) = 0 and Σ(π) = vv ⊤ . Hence, for all i ∈ {1, . . . , d}, we have [Σ(π)]ii = vi2 > 0. π :=
So π ∈ P2◦ . Set z := M −1/2 v ̸= 0. Since d ≥ 2, the hyperplane (M 1/2 e1 )⊥ = {w ∈ Rd | ⟨w, M 1/2 e1 ⟩ = 0} has dimension d − 1 ≥ 1, and therefore intersects the sphere of radius ∥z∥. Choose w ∈ (M 1/2 e1 )⊥ such that ∥w∥ = ∥z∥. By transitivity of O(d) on the sphere, there exists R ∈ O(d) such that Rz = w. For this choice of R, 1/2 1/2 [AR v]1 = e⊤ RM −1/2 v = e⊤ w = ⟨M 1/2 e1 , w⟩ = 0. 1M 1M
Since covariance transforms under the affine map gR by congruence with AR , ⊤ ⊤ ⊤ Σ((gR )# π) = AR Σ(π)A⊤ R = AR vv AR = (AR v)(AR v) .
Its first diagonal entry is therefore [Σ((gR )# π)]11 = [AR v]21 = 0. 15
Hence (gR )# π ∈ / P2◦ . This proves that P2◦ is not G-stable, and thus ρ does not satisfy Condition (3). Since the problem already arises with the domain, one may try to suitably restrict it to ensure Gstability, with the hope that the remaining part of the condition will go through. We now show this not to be the case. Let P2□ be the set of distributions with finite second moment and positive-definite d covariance, and let C++ be the set of positive-definite matrices with unit diagonal entries. Clearly □ ◦ d d P2 ⊂ P2 and C++ ⊂ C+ , so we can define the correlation statistic in the same way as in (8). Let π ∈ P2□ and gR ∈ G. Since gR is affine, (gR )# π ∈ P2 . Moreover, by (4), Σ((gR )# π) = AR Σ(π)A⊤ R. d d Because M ∈ S++ and R ∈ O(d), AR is invertible. Therefore, since Σ(π) ∈ S++ , it follows that d AR Σ(π)A⊤ R ∈ S++ .
Thus (gR )# π ∈ P2□ , proving G-stability. However, the correlation statistic still does not satisfy Condition (3). It suffices to give a counterexample. We work in R2 and consider the spherical subcase m = 0, M = I2 , for which G = {gR | R ∈ O(2)}, Define
gR (x) = Rx.
1 π2 := N 0, 0
π1 := N (0, I2 ),
0 . 3
Clearly, both belong to P2□ , and ρ(π1 ) = I2 = ρ(π2 ), since both covariance matrices are diagonal with positive entries. Now let 1 1 −1 ∈ O(2). R := √ 2 1 1 Then Σ((gR )# π1 ) = I2 , so ρ((gR )# π1 ) = I2 . By contrast,
1 Σ((gR )# π2 ) = R 0
0 2 ⊤ R = 3 −1
−1 . 2
Therefore D((gR )# π2 ) = diag(Σ((gR )# π2 )) = 2I2 , and hence −1/2
ρ((gR )# π2 ) = (2I2 )
2 −1
−1 1 −1/2 (2I2 ) = 2 −1/2
−1/2 . 1
but
ρ((gR )# π1 ) ̸= ρ((gR )# π2 ).
Thus, there exists gR ∈ G such that ρ(π1 ) = ρ(π2 )
This violates the implication required by Condition (3).
C
Additional material for Section 5
C.1
Uniqueness of the axis of symmetry
In this appendix, we formally prove that if a non-uniform distribution is rotationally symmetric about two axis directions, then these determine the same line through the origin in Rd . Proposition C.1 (Uniqueness of the rotational axis) Fix d ≥ 3. Let π ∈ Pσ be a non-uniform distribution on S d−1 with density fπ , and suppose it is rotationally symmetric about both w1 , w2 ∈ S d−1 , i.e., there exist measurable functions ψ, φ : [−1, 1] → R≥0 such that fπ (x) = ψ(w1⊤ x) = φ(w2⊤ x)
for σ-a.e. x ∈ S d−1 .
Then span(w1 ) = span(w2 ). 16
Proof . Let T (π) := {R ∈ SO(d) | (gR )# π = π}, denote the subgroup of SO(d) under which π is invariant. T (π) is closed in SO(d): if Rn → R in SO(d) with (gRn )# π = π for all n, then for every h ∈ C(S d−1 ), the functions x 7→ h(Rn x) converge uniformly on S d−1 to x 7→ h(Rx) and hence Z Z Z Z h(x)dπ(x) = h(Rn x)dπ(x) −→ h(Rx)dπ(x) = h(x)d((gR )# π)(x) S d−1
S d−1
S d−1
S d−1
which guarantees that (gR )# π = π and therefore R ∈ T (π). Let T (π)0 ≤ T (π) denote the identity component of T (π), which is a connected and closed subgroup. As T (π) is closed in SO(d), T (π)0 is also closed in SO(d). Letting Gwi := {R ∈ SO(d) | Rwi = wi } for i = 1, 2, we know that since π is rotationally symmetric about both w1 , w2 ∈ S d−1 , we have Gw1 , Gw2 ⊆ T (π). Moreover, SO(d − 1) is a maximal connected closed subgroup of SO(d) (Dynkin, 1957; Uchida, 1979), and Gw1 , Gw2 are conjugate in SO(d) to the canonical embedding of SO(d − 1). Since conjugation is a group automorphism and a topological homeomorphism, it preserves maximal connected closed subgroups, and so Gw1 , Gw2 are themselves maximal connected closed subgroups of SO(d). As they are connected, and both contain the identity, it follows that Gw1 , Gw2 ≤ T (π)0 . Assume, towards a contradiction, that w1 ̸= ±w2 . Then there exists R ∈ Gw2 \ Gw1 . Indeed, otherwise every rotation fixing w2 would also fix w1 , so w1 would belong to the common fixed-point subspace of Gw2 , namely span(w2 ), which would force w1 = ±w2 . Since Gw2 ≤ T (π)0 , it follows that R ∈ T (π)0 \ Gw1 . Hence Gw1 is a proper subgroup of T (π)0 . But Gw1 is a maximal connected closed subgroup of SO(d) and so it must be that T (π)0 = SO(d). It follows that T (π) = SO(d). Therefore π is invariant under every rotation of S d−1 . The only rotation-invariant probability measure on S d−1 is the uniform measure, a contradiction. So span(w1 ) = span(w2 ). □ C.2
Explicit construction of the unique minimiser example
Here, we justify the claim made at the end of Section 5 for the reverse KL divergence objective. In particular, we show that there exists a threshold ηc (d, λ, κ0 ) > 0 such that the KL has a unique minimiser for 0 < η ≤ ηc (d, λ, κ0 ), whereas for η > ηc (d, λ, κ0 ) the minimiser set is non-singleton. Fix λ ̸= 0, u ∈ S d−1 , κ0 > 0, and η > 0. For ν ∈ S d−1 , let X ∼ Qν,κ0 , and write Lλ,η,κ0 (ν) for the reverse KL objective DKL (Qν,κ0 ∥Pλ,η ). Using log fQν,κ0 (x) = log cd (κ0 ) + κ0 ν ⊤ x and
log fPλ,η (x) = − log Zλ,η + λu⊤ x − η(u⊤ x)2 ,
we obtain fQν,κ0 (X) Lλ,η,κ0 (ν) = E log fPλ,η (X) = log cd (κ0 ) + log Zλ,η + κ0 E[ν ⊤ X] −λ E[u⊤ X] +η E[(u⊤ X)2 ] . | {z } | {z } | {z } (E1 )
(E2 )
(9)
(E3 )
We begin by reducing each expectation in (9) to a function of the single scalar c := u⊤ ν ∈ [−1, 1]. Set T := ν ⊤ X. With ρd (t) := (1 − t2 )(d−3)/2 , the density of T is (Mardia and Jupp, 2009, Eq. 9.3.12) fT (t) ∝ eκ0 t ρd (t),
t ∈ [−1, 1].
(10)
Now define
Ad,κ0 := E[T ], md,κ0 := E[T 2 ]. Then the mean and variance of X are given by (Mardia and Jupp, 2009, Eqs. 9.3.33–9.3.34) E[X] = Ad,κ0 ν,
Var(X) = Var(T ) νν ⊤ +
1 − md,κ0 (Id − νν ⊤ ). d−1
Therefore, the expectations (E1 ) and (E2 ) in (9) are E[ν ⊤ X] = Ad,κ0 , E[u⊤ X] = u⊤ E[X] = Ad,κ0 u⊤ ν = Ad,κ0 c. 17
(11)
Moreover, the second moment matrix is E[XX ⊤ ] = Var(X) + E[X]E[X]⊤ 1 − E[T 2 ] = Var(T ) + (E[T ])2 νν ⊤ + (Id − νν ⊤ ) d−1 1 − md,κ0 = md,κ0 νν ⊤ + (Id − νν ⊤ ), d−1 and hence, expectation (E3 ) in (9) reduces to E[(u⊤ X)2 ] = u⊤ E[XX ⊤ ]u 1 − md,κ0 1 − (u⊤ ν)2 d−1 d md,κ0 − 1 2 1 − md,κ0 + c . = d−1 d−1 | {z } = md,κ0 (u⊤ ν)2 +
(12)
Bd,κ0
Substituting (11) and (12) into (9) yields Lλ,η,κ0 (ν) = ηBd,κ0 c2 − λAd,κ0 c + Cλ,η,κ0 , where
1 − md,κ0 , d−1 is a constant independent of ν. We have thus reduced the optimisation problem to the scalar quadratic Cλ,η,κ0 := log cd (κ0 ) + log Zλ,η + κ0 Ad,κ0 + η
cλ,η,κ (c) := ηBd,κ c2 − λAd,κ c, L 0 0 0
c ∈ [−1, 1].
(13)
Since η > 0, once Bd,κ0 > 0 is established, the quadratic in (13) is strictly convex and therefore has a unique unconstrained minimiser. Because d ≥ 3, the denominator of Bd,κ0 is positive, so it suffices to show that d md,κ0 − 1 > 0, that is, md,κ0 > 1/d. To this end, let Y ∼ σ be uniform on S d−1 and define T0 := u⊤ Y . Then the density of T0 is (Mardia and Jupp, 2009, Eq. 9.3.1) fT0 (t) ∝ ρd (t),
t ∈ [−1, 1].
(14) 2
We use (14) to lower bound md,κ0 . Expanding the definition, observe that t ρd (t) is even while eκ0 t = cosh(κ0 t) + sinh(κ0 t), so the odd terms integrate to 0 on [−1, 1], and we obtain Z 1 Z 1 Z 1 2 t2 eκ0 t ρd (t) dt t cosh(κ0 t)ρd (t) dt t2 ρd (t) dt md,κ0 = Z−11 = Z0 1 > Z0 1 = E[T02 ], (15) κ0 t e ρd (t) dt cosh(κ0 t)ρd (t) dt ρd (t) dt −1
0
0
where the strict inequality follows from Chebyshev’s integral inequality (Mitrinović, 1970, Thm. 10, Ch. 2.5), since t 7→ t2 and t 7→ cosh(κ0 t) are strictly increasing on [0, 1]. To compute E[T02 ], note that since Y ∼ σ is uniform on S d−1 , its law is rotation invariant and therefore E[Y Y ⊤ ] = αId for some α ∈ R. Taking traces, tr(E[Y Y ⊤ ]) = E[tr(Y Y ⊤ )] = E[∥Y ∥2 ] = 1,
tr(αId ) = αd,
and hence α = 1/d. The second moment of T0 is thus 1 1 1 2 ⊤ ⊤ ⊤ Id u = ∥u∥2 = . E[T0 ] = u E[Y Y ]u = u d d d Combining this with (15) gives md,κ0 > 1/d, and therefore Bd,κ0 > 0. Consequently, the scalar objective in (13) has the unique unconstrained minimiser c⋆ =
λAd,κ0 . 2ηBd,κ0 18
Define
|λ|Ad,κ0 . 2Bd,κ0 Then, ηc (d, λ, κ0 ) > 0 since |λ|, Bd,κ0 > 0 and moreover Ad,κ0 > 0. Indeed, since ρd (t) is even, Z 1 Z 1 teκ0 t ρd (t) dt t eκ0 t − e−κ0 t ρd (t) dt Ad,κ0 = Z−11 = 0 > 0. Z 1 κ0 t κ0 t e ρd (t) dt e ρd (t) dt ηc (d, λ, κ0 ) :=
−1
−1 ⋆
If 0 < η ≤ ηc (d, λ, κ0 ), then |c | ≥ 1, and so the unique constrained minimiser of (13) over [−1, 1] is c = sgn(λ). Since u⊤ ν = sgn(λ) if and only if ν = sgn(λ)u, the original objective Lλ,η,κ0 (ν) has the unique minimiser ν ⋆ = sgn(λ)u. Consequently, exact statistic recovery holds with A(Qν ⋆ ,κ0 ) = span(u) = A(Pλ,η ). If η > ηc (d, λ, κ0 ), then c⋆ ∈ (−1, 1). Hence argmin Lλ,η,κ0 (ν) = {ν ∈ S d−1 | u⊤ ν = c⋆ }.
(16)
ν∈S d−1
Since d ≥ 3 and c⋆ ∈ (−1, 1), the set in (16) is a (d − 2)-sphere and hence is not a singleton. Moreover, no minimiser ν satisfies A(Qν,κ0 ) = A(Pλ,η ), because A(Qν,κ0 ) = span(ν), while span(ν) = span(u) would force ν = ±u and hence u⊤ ν ∈ {±1}, contradicting c⋆ ∈ (−1, 1). To visualise the axis recovery, we work with d = 3. For fixed λ > 0 and κ0 > 0, we can explicitly compute the constants A3,κ0 and B3,κ0 needed to evaluate ηc (3, λ, κ0 ). First observe that since ρ3 (t) = 1 for all t ∈ [−1, 1], the marginal density of T = ν T X in dimension 3 is proportional to eκ0 t . Writing Z 1 2 sinh κ0 Z(κ0 ) := eκ0 t dt = , κ0 −1 the first two moments of T can be obtained by differentiating the log-partition function: 2 d2 d d 2 log Z(κ0 ), m3,κ0 = E[T ] = 2 log Z(κ0 ) + log Z(κ0 ) . A3,κ0 = E[T ] = dκ0 dκ0 dκ0 Since log Z(κ0 ) = log 2 + log sinh κ0 − log κ0 , this yields A3,κ0 = coth κ0 −
1 , κ0
m3,κ0 = E[T 2 ] = 1 −
2 2 coth κ0 + 2 . κ0 κ0
Substituting m3,κ0 into the definition of B3,κ0 gives B3,κ0 = 1 −
3 3 coth κ0 + 2 κ0 κ0
and the critical threshold therefore becomes ηc (3, λ, κ0 ) =
λκ0 (κ0 coth κ0 − 1) . 2(κ20 − 3κ0 coth κ0 + 3)
In the example of Figure 2, we take λ = 1 and κ0 = 2.5, so A3,2.5 ≈ 0.6135, B3,2.5 ≈ 0.2637, ηc (3, 1, 2.5) ≈ 1.1632. In the left panel, we use η = 1 < ηc (3, 1, 2.5), and the unique minimiser ν ⋆ of Lλ,η,κ0 (ν) satisfies A(Qν ⋆ ,2.5 ) = A(P1,1 ) = span(u), while in the right panel, we use η = 2 > ηc (3, 1, 2.5), so the minimiser set is the latitude circle A3,2.5 2 ⊤ ν∈S :u ν= ≈ 0.5816 , 4B3,2.5 and A(Qν,2.5 ) ̸= A(P1,2 ) for all minimisers ν. 19