Theory Digest — Apr 2, 2026

Motivation: High-dimensional Bayesian clustering suffers from the curse of dimensionality where posterior distributions over partitions degenerate to trivial solutions (all observations in one cluster or all singletons), regardless of true structure. This fundamentally breaks mixture-model based clustering when $p \gg n$.

Mathematical setup: Given data matrix $X \in \mathbb{R}^{n \times p}$ with $n$ observations and $p$ dimensions. Let $z = (z_1, \ldots, z_n)$ denote a partition represented by cluster membership indicators where $z_i = h$ if observation $i$ belongs to cluster $h$. The space of all partitions is denoted $\mathcal{Z}$.

Traditional approach: posterior distribution on partitions

\[p(z | X) = \frac{p(z, X)}{p(X)}\]

VCI splits data vertically into $K$ shards: $X_i = (X_i^{(1)}, \ldots, X_i^{(K)})$ where $X_i^{(k)} \in \mathbb{R}^{p_k}$ and $\sum_{k=1}^K p_k = p$.

Assumptions:

1. Partition space $\mathcal{Z}$ is finite and equipped with ground metric $c: \mathcal{Z} \times \mathcal{Z} \to \mathbb{R}_+$

2. Shard posteriors $p_k(z

   X^{(k)})$ are well-defined for each shard $k$

3. Entropic regularization parameter $\epsilon_k > 0$ ensures numerical stability

   Toy example: For Old Faithful data with $n=272, p=2$ (eruption duration, waiting time), split into 2 one-dimensional shards: $X^{(1)}$ contains durations, $X^{(2)}$ contains waiting times. Each shard fits a truncated Dirichlet process mixture, producing shard posteriors $p_1(z

   X^{(1)})$ and $p_2(z

   X^{(2)})$.

   Formal objective: Construct consensus posterior as entropic regularized Wasserstein barycenter

   \[\bar{p}(z | X, \lambda) = \arg\min_{p \in \mathcal{P}(\mathcal{Z})} \sum_{k=1}^K \lambda_k W_{c,\epsilon_k}(p, p_k)\]

VCI constructs consensus posterior via Wasserstein barycenter of shard posteriors.

Method steps:

1. Split data vertically: $X_i = (X_i^{(1)}, \ldots, X_i^{(K)})$ for each observation $i$

2. Compute shard posteriors: $p_k(z

   X^{(k)})$ for each shard $k = 1, \ldots, K$

3. Define entropic Wasserstein distance between discrete measures:

   \[W_{c,\epsilon}(p_1, p_2) = \min_{\pi \in \Pi(p_1,p_2)} \sum_{i,j} c(z_i, z_j)\pi(z_i, z_j) - \epsilon H_\pi(z_1, z_2)\]

4. Construct barycenter weights $\lambda_k$ favoring meaningful partitions:

   \[\omega_k^P(p_k) = E\left[\frac{4(\tilde{H}-1)(n-\tilde{H})}{(n-1)^2} \Big| X^{(k)}\right] \cdot E[\exp(-aE(z)) | X^{(k)}] \cdot (1-4U)\]

   where $\tilde{H} = \exp(-H(z))$, $E(z) = -H(z)/\log(H(z))$, $U = \frac{2}{n(n-1)}\sum_{i<j} p_{ij}(1-p_{ij})$

5. Solve barycenter optimization:

   \[\min_{\alpha \in \Delta^{|\mathcal{Z}|}} \sum_{k=1}^K \lambda_k W_{c,\epsilon_k}(\alpha, \alpha^{(k)}; M)\]

   via iterative Bregman projections.

   Application to toy example: For Old Faithful with 2 shards, each produces ~1000 MCMC samples. Using VoI ground metric and $\epsilon_k = 0.05$, the barycenter optimization combines both shard posteriors. With $\lambda^P$ weights (penalizing trivial partitions), VCI achieves Wasserstein distance 0.2465 to ground truth versus 0.3030 and 0.3629 for individual shards.

Prior work:

1. Srivastava et al. (2015, 2018): Introduced Wasserstein consensus Monte Carlo for horizontal data splits (observations), using uniform barycenter weights
2. Chandra, Canale and Dunson (2023): Bayesian latent factor mixture model for high-dimensional clustering via dimensionality reduction

3. Scott et al. (2022): Consensus Monte Carlo framework for distributed Bayesian computation with horizontal splitting

   Delta: This paper introduces vertical consensus inference, splitting along feature dimensions rather than observations. Key additions:

4. vertical data sharding for high-dimensional clustering
5. structured barycenter weights favoring non-trivial partitions

6. theoretical justification as variational approximation to hierarchical generalized Bayes model.

   Theory-specific assessment:

   • Main result (Proposition 2.1) shows VCI corresponds to upper bound of NELBO in variational inference framework - this connection is somewhat predictable given existing variational-optimal transport links
   • Proof technique is routine assembly of entropy inequalities and Wasserstein barycenter properties
   • No lower bounds established; optimality gaps in Proposition 2.1 depend on restrictive “star structure” assumptions

   • The hierarchical model interpretation (equation 5) with generalized Bayes prior $p(z^{(k)}

     z) \propto \exp(-\zeta_k c(z^{(k)}, z))$ provides post-hoc justification but wasn’t the original motivation

   Verdict: INCREMENTAL — solid extension of existing Wasserstein consensus methods from horizontal to vertical splitting, with reasonable theoretical backing but no fundamental breakthroughs in understanding high-dimensional clustering.

Main proof strategy for Proposition 2.1 uses entropy decomposition and convex optimization bounds.

Key steps:

1. Decompose joint entropy using chain rule:

   \[H_q(z, z^{(1):(K)}) = H_q(z) + H_q(z^{(1):(K)} | z)\] \[= H_q(z) + \frac{1}{K}\sum_{k=1}^K [H_q(z^{(k)} | z) + H_q(z^{(-k)} | z, z^{(k)})]\]

2. Apply non-negativity of conditional entropy:

   \[H_q(z^{(-k)} | z, z^{(k)}) \geq 0\]

   to obtain lower bound:

   \[H_q(z, z^{(1):(K)}) \geq \frac{1}{K}\sum_{k=1}^K H_q(z, z^{(k)})\]

3. Define upper bound for NELBO:

   \[\mathcal{L}(q) = E_{(z,z^{(1):(K)}) \sim q}[-\log p(z, z^{(1):(K)})] - \frac{1}{K}\sum_{k=1}^K H_q(z, z^{(k)}) \geq L(q)\]

4. Key inequality: Restrict to star couplings $\Pi^*(q_0, q_1, \ldots, q_K)$ where:

   \[q(z, z^{(1)}, \ldots, z^{(K)}) = q_0(z)\prod_{k=1}^K q(z^{(k)} | z)\]

5. Key equation: Express log-likelihood using hierarchical model:

   \[\log p(z, z^{(1):(K)}, X) = \sum_{k=1}^K [-\zeta_k c(z^{(k)}, z) + \log p_k(X^{(k)} | z^{(k)})] - C\]

6. Final bound: Combine steps to show:

   \[\min_{q \in \Pi(q_0,q_1,\ldots,q_K)} L(q) \leq \sum_{k=1}^K \zeta_k W_{c,\frac{1}{K\zeta_k}}(q_0, q_k) - \sum_{k=1}^K E_{z^{(k)} \sim q_k}[\log p_k(X^{(k)} | z^{(k)})] + C\]

   Technical insight: The proof exploits that entropic Wasserstein distance emerges naturally when upper-bounding KL divergence terms in the NELBO via Jensen’s inequality and optimal transport duality.

Datasets and setup:

1. Scenario 1: Old Faithful Geyser (n=272, p=2) - validation in low dimensions
2. Scenario 2: Old Faithful + 18 noise dimensions (total p=20) - robustness to irrelevant features

3. Scenario 3: Single-cell data (n=800, p=25,348 genes, K=100 shards) - high-dimensional performance

   Baselines: Shard posteriors, uniform mixture of shards, entropy-weighted mixture, full-data posterior (when feasible)

   Key numbers:

   • Scenario 1: VCI achieves 0.2465 VoI distance to ground truth vs 0.3030-0.3629 for individual shards
   • Scenario 2: VCI with structured weights λ^P achieves 0.0031 distance vs 0.8847-3.4301 for noisy shards
   • Scenario 3: VCI posteriors achieve 1.5496-1.7756 expected VoI to pseudo-labels vs 2.3026 for full posterior

   Computational details: Uses iterative Bregman projections with ε_k = 0.05, 100-1000 MCMC samples per shard, Dirichlet process mixtures for continuous data, Poisson-Gamma mixtures for count data.

Limitations:

1. TECHNICAL: Requires pre-specification of number of shards K and shard assignments - likely removable via adaptive splitting schemes

2. TECHNICAL: Entropic regularization parameter ε_k needs manual tuning for numerical stability - standard in optimal transport literature

3. RESTRICTIVE: Theoretical guarantees in Proposition 2.1 depend on “star coupling” structure where z^{(-k)} is deterministic given (z, z^{(k)}) - significantly narrows applicability

4. TECHNICAL: Barycenter computation complexity O(

   Z

   ^2) limits scalability to moderate numbers of partition samples - addressable via free support methods

5. RESTRICTIVE: Ground truth evaluation relies on pseudo-labels in high-dimensional experiments - fundamental challenge in unsupervised learning validation

6. NATURAL: Method assumes partition structure exists within feature subsets - reasonable for many applications but may fail when clustering requires full feature interaction

   Open problems:

7. Develop principled methods for optimal vertical shard construction beyond manual splitting or prior knowledge
8. Establish finite-sample convergence rates for VCI approximation quality relative to full-data posterior under specific data-generating assumptions

Motivation: High-dimensional Bayesian inverse problems require expressive priors that can capture complex data distributions while maintaining interpretable posterior summaries. Deep generative models provide flexibility but typically produce entangled latent representations that are difficult to interpret in terms of physical parameters.

Mathematical setup: Let $x \in \mathbb{R}^p$ denote an unknown high-dimensional field and $y \in \mathbb{R}^m$ denote observations. The forward model is:

\[y = F(x) + \epsilon, \quad \epsilon \sim p_\epsilon\]

The structured hierarchical prior uses a disentangled generator $G_\theta$ with partitioned latent space $z = (z_{aux}, z_{rec})$ where $z_{aux} \in \mathbb{R}^d$ aligns with auxiliary variables $u \in \mathbb{R}^d$:

\[u \sim p(u)\] \[z_{aux} | u \sim N(u, \tau^2 I_d)\] \[z_{rec} \sim N(0, I_{d_Z-d})\] \[x = G_\theta(z_{aux}, z_{rec})\]

Assumptions:

1. The generator $G_\theta$ is differentiable and trained via auxiliary-variable guidance
2. Auxiliary variables $u$ represent interpretable physical parameters
3. Noise model $p_\epsilon$ has known covariance $\Sigma_\epsilon$

4. Prior $p(u)$ encodes domain knowledge about interpretable parameters

   Toy example: For $d = 2$ auxiliary variables representing mean and variance of a 1D conductivity field, and $d_Z - d = 4$ residual dimensions, the generator maps $(z_{aux}, z_{rec}) \in \mathbb{R}^6$ to $x \in \mathbb{R}^{100}$ representing log-conductivity values on a discretized domain.

   Formal objective: The posterior inference problem is to characterize:

   \[p(u, z_{rec} | y) \propto p_\epsilon(y - F(G_\theta(u, z_{rec}))) p(u) N(z_{rec}; 0, I)\]

The method consists of offline training of a disentangled generator followed by online latent-space inference.

Offline Training: Use Aux-VAE or DL-CFM with correlation-based penalties. For Aux-VAE, minimize:

\[L_{AuxVAE} = -L_{VAE} + \lambda_1 \sum_{j=1}^d [\text{Align}(u_j, \mu_{\phi,aux,j}) + \text{Decorr}(u_j, \mu_{\phi,aux,-j})] + \lambda_2 \text{Decorr}(u, \mu_{\phi,rec})\]

where the penalties use polynomial correlation features:

\[R_0^{(K)}(v,w) := \frac{1}{K^2 m_v m_w} \sum_{k,k'=1}^K \sum_{i,j} |\text{Corr}(v^{(k)}, w^{(k')})_{ij}|\]

Online Inference: In the hard-constrained limit $\tau \to 0$, perform inference on $(u, z_{rec})$ with log-posterior:

\[\log p(u, z_{rec} | y) = -\frac{1}{2} \|y - F(G_\theta(u, z_{rec}))\|^2_{\Sigma_\epsilon^{-1}} + \log p(u) - \frac{1}{2}\|z_{rec}\|_2^2\]

Use gradient-based MAP estimation or HMC sampling in the low-dimensional latent space.

Toy example application: For the 1D conductivity example, the method would learn to map auxiliary variables $(u_1, u_2)$ representing field mean/variance to the guided latent block, while residual variations are captured by $z_{rec}$. During inference, posterior samples of $(u_1, u_2)$ directly quantify uncertainty in interpretable physical parameters.

Prior work:

1. “Deep generative priors for inverse problems” (Asim et al. 2020) - used VAE/GAN priors but without interpretable structure
2. “Aux-VAE: Auxiliary Variables for Variational Autoencoders” (Larsen et al. 2020) - introduced auxiliary-guided disentanglement but not for inverse problems

3. “Physics-informed neural networks” (Raissi et al. 2019) - incorporated physics constraints but without structured Bayesian priors

   Delta: This paper formulates auxiliary-guided disentanglement as a hierarchical prior for Bayesian inverse problems. The key additions are:

4. treating auxiliary variables as random with their own prior rather than fixed supervision
5. theoretical analysis of when latent disentanglement translates to posterior decoupling via linearization

6. latent-space inference procedures.

   Theory-specific assessment:

   • Main results (Lemma 4.2, Proposition 4.5) are predictable applications of linearization and block-diagonal structure analysis
   • Proof techniques are routine: Taylor expansions, QR decompositions, standard Gaussian posterior calculations
   • No fundamental lower bounds are established; tightness of linearization approximations is not quantified
   • The connection between encoder-side correlation penalties and generator-side tangent space orthogonality remains loose

   The theoretical analysis provides useful local characterizations but relies on linearization assumptions that may not hold globally. The correlation-based disentanglement penalties are heuristic proxies that don’t guarantee independence.

   Verdict: INCREMENTAL — solid combination of existing disentanglement techniques with inverse problem formulation, but theoretical insights are straightforward applications of standard tools.

The main theoretical results use standard linearization and Gaussian analysis:

Lemma 4.2 (Linearized covariance): Uses Taylor expansion around reference point:

\[x_{lin} := G_\theta(m_u, 0) + J_{aux}(z_{aux} - m_u) + J_{rec} z_{rec}\]

where $J_{aux} = \frac{\partial G_\theta}{\partial z_{aux}}

_{(m_u,0)}$ and $J_{rec} = \frac{\partial G_\theta}{\partial z_{rec}}

_{(m_u,0)}$.

Key insight: Independence in latent prior yields additive covariance structure:

\[\text{Cov}(x_{lin}) = J_{aux}(\Sigma_u + \tau^2 I) J_{aux}^T + J_{rec} J_{rec}^T\]

Proposition 4.5 (Posterior decoupling): Under linearization $y \approx H(u_0, z_0) + A_u(u-u_0) + A_{rec}(z_{rec} - z_0) + \epsilon$, the key decoupling condition is:

\[A_u^T \Sigma_\epsilon^{-1} A_{rec} = 0\]

When satisfied, the Gaussian posterior precision matrix becomes block-diagonal, implying independence between $u$ and $z_{rec}$.

Corollary 4.4 (Block structure): Uses QR decomposition and submultiplicative norm bounds to show that cross-covariance between projected coordinates satisfies:

\[\|\text{Cov}(x_{aux}, x_{rec})\|_2 \leq \|Q_{aux}^T Q_{rec}\|_2 (\|R_{aux}(\Sigma_u + \tau^2 I)R_{aux}^T\|_2 + \|R_{rec}R_{rec}^T\|_2)\]

The proofs are elementary applications of matrix algebra and Gaussian calculations. No non-standard techniques or sharp concentration bounds are employed.

Datasets: Two elliptic PDE problems on 28×28 grids (676 DOFs):

1. Conductivity identification: Unknown log-conductivity field, GP-generated training data (N=12,000)

2. Source identification: Unknown log-source field, Gaussian bump structure (N=5,000)

   Baselines:

   • Plain VAE (no auxiliary structure)
   • GP-Fixed (oracle hyperparameters)
   • GP-Hier (learnable GP hyperparameters)

   Key numbers:

   • Conductivity (well-specified): AuxVAE achieves RMSE=0.358, coverage=0.969, ESS=1046 vs GP-Fixed RMSE=0.340, coverage=0.978, ESS=3000
   • Source (misspecified): AuxVAE achieves RMSE=1.455, coverage=0.944 vs GP-Fixed RMSE=4.484, coverage=0.198

   Architecture: Conv1D encoder, MLP decoder, latent dims d_aux=4, d_rec=16-48. HMC with 3000 samples, dual-averaging adaptation.

   Results show AuxVAE matches oracle GP under correct specification and substantially outperforms GP baselines under misspecification, while recovering interpretable physical parameters.

Limitations:

1. Linearization assumptions for theoretical analysis - TECHNICAL (needed for tractable analysis but likely removable with more sophisticated techniques)

2. Correlation-based disentanglement penalties don’t guarantee independence - TECHNICAL (practical proxy but theoretically incomplete)

3. Generator tangent space orthogonality not directly controlled by training objective - TECHNICAL (could be enforced with additional penalties)

4. Limited to moderate-dimensional latent spaces and specific PDE problems - RESTRICTIVE (scalability to high-dimensional latents unclear)

5. Requires paired (x,u) training data with known auxiliary variables - RESTRICTIVE (limits applicability when interpretable factors are unknown)

6. No convergence guarantees for disentanglement training or posterior sampling - TECHNICAL (standard limitation of VAE training and MCMC)

   Open problems:

7. Develop non-linear analysis beyond first-order Taylor approximations to characterize when posterior decoupling occurs without linearization assumptions

8. Design training objectives that directly control generator-side tangent space properties rather than relying on encoder-side correlation penalties

1. Motivation: Estimating smooth functionals f(θ) of a mean parameter θ in high-dimensional or infinite-dimensional spaces is challenging because naive plug-in estimators suffer from bias that can dominate the target estimation rate. This problem matters for precision matrix estimation, covariance operator analysis, and nonparametric functional estimation where optimal rates depend on both functional smoothness and space complexity.

2. Mathematical setup: Let B be a separable Banach space with norm

   ·

   . Given distribution P on B with mean θ = E[W], consider functional f: Θ → ℝ on open set Θ ⊆ B. For m-smooth functional with m = s + ρ (s = ⌈m⌉ - 1, ρ ∈ (0,1]), f admits expansion:

   \[f(x + h) = f(x) + \sum_{k=1}^{s} \frac{D^k f(x)[h,...,h]}{k!} + o(||h||^{s+ρ})\]

   Assumptions:

   1. Distribution P has finite weak variance ν = sup_{

      u

      *≤1} E[

      ⟨W-θ,u⟩

      ²] and strong variance V = E[

      W-θ

      ²]

   2. Preliminary estimators θ̂_n satisfy probability bounds with rate r_n
   3. Functional derivatives D^k f(θ) are bounded in operator norm

3. Toy example: When B = ℝ², θ = (θ₁,θ₂), and f(θ) = θ₁²θ₂, the bias of plug-in estimator f(θ̂) involves second-order terms like E[θ̂₁²θ̂₂] - θ₁²θ₂ ≈ 2θ₁θ₂Var(θ̂₁) + θ₁²Var(θ̂₂) which can be O(1/n) rather than o(1/√n).

4. Formal objective: Estimate

   \[f(θ) = E_P[f(W)]\]

   achieving √n-rate and asymptotic normality under minimal regularity conditions.

The method uses cross-fitted debiasing with sample splitting and higher-order corrections.

Algorithm:

1. Split sample W₁,…,W₂ₙ into S₁ = {W₁,…,Wₙ} and S₂ = {Wₙ₊₁,…,W₂ₙ}
2. Compute pilot estimator θ̂_S₂ from S₂

3. Form U-statistic corrections from S₁:

   \[T_k[\bar{U}^{(k)}(b)] = \binom{n}{k}^{-1} \sum_{1≤j₁<...<j_k≤n} T_k[W_{j₁}-b,...,W_{j_k}-b]\]

4. One-sided estimator:

   \[\hat{f}_s(S₁,S₂) = f(\hat{θ}_{S₂}) + \sum_{k=1}^{s} \frac{1}{k!} D^k f(\hat{θ}_{S₂})[\bar{U}^{(k)}(\hat{θ}_{S₂})]\]

5. Final cross-fitted estimator:

   \[\hat{f}_s = \frac{1}{2}[\hat{f}_s(S₁,S₂) + \hat{f}_s(S₂,S₁)]\]

   Toy example application: For f(θ) = θ₁²θ₂ with θ̂ = (θ̂₁,θ̂₂), the correction terms involve:

   • First-order: D¹f(θ̂)[Ū⁽¹⁾(θ̂)] = 2θ̂₁θ̂₂(W̄₁-θ̂₁) + θ̂₁²(W̄₂-θ̂₂)
   • Second-order: D²f(θ̂)[Ū⁽²⁾(θ̂)] involving products of centered observations

   The cross-fitting ensures the U-statistics remain degenerate conditional on the independent pilot.

Prior work:

1. Zhou et al. (2021): Developed higher-order U-statistic corrections in Hilbert spaces without sample splitting, achieving rate 1/√n + (√d/n)^m
2. Koltchinskii & Li (2026): Bootstrap bias correction in Banach spaces using independent blocks, requiring 1 + s(s+1)/2 base estimators

3. Koltchinskii series (2018-2023): Iterative bootstrap bias correction achieving minimax rates but with computational burden

   Delta: This paper introduces single sample-split cross-fitted debiasing that:

   • Uses only 2 blocks vs. O(s²) blocks in prior work
   • Extends to infinitely differentiable (Gevrey) functionals
   • Achieves dimension regime d log²(en) = o(n) vs. previous d = o(n^(m-1)/m)
   • Provides computational relaxations for matrix functionals

   Theory assessment:

   • Main theorem provides moment bounds E^(1/2)[

     f̂_s - f(θ)

     ²] ≤ √(ν/n) + (√d/n)^m which matches minimax lower bounds

   • Proof technique combining cross-fitting with degenerate U-statistics is novel
   • Berry-Essén bounds with explicit constants are non-trivial
   • Extension to Gevrey classes with truncation level s ≍ log(n) is genuinely new

   The bounds are tight against known lower bounds in Gaussian shift models. The Gevrey extension required developing new truncation analysis.

   Verdict: SIGNIFICANT — The cross-fitted construction cleverly preserves degeneracy while using minimal sample splitting, and the Gevrey extension substantially broadens the functional classes amenable to efficient estimation.

Main proof strategy combines three components:

1. Cross-fitting analysis: Key insight is that sample splitting preserves conditional degeneracy. For independent pilot θ̂_S₂, the corrected terms remain degenerate U-statistics:

   \[E[T_k[W_{j₁}-\hat{θ}_{S₂},...,W_{j_k}-\hat{θ}_{S₂}] | \hat{θ}_{S₂}] = 0\]

2. Moment bounds via Rosenthal inequality: For U-statistic of order k, the variance scaling uses:

   \[\text{Var}[T_k[\bar{U}^{(k)}(θ)]] = \binom{n}{k}^{-1} \text{Var}[T_k[W_1-θ,...,W_k-θ]]\]

   Combined with Rosenthal inequality for sums of independent random variables.

3. Remainder control: The key remainder bound uses functional smoothness:

   \[|R_s| ≤ \sum_{k=0}^{s} \frac{1}{k!} |J^{s-k}\Delta^{(s)}[(\theta-\tilde{θ})^{\otimes(s-k)}, \bar{U}^{(k)}(θ)]|\]

   where J^{s-k} is Riemann-Liouville operator and Δ^{(s)} captures derivative differences.

4. Berry-Essén bounds: Use Stein’s method for multilinear forms, establishing:

   \[\sup_t |P(\sqrt{n}(\hat{f}_s - f(θ))/σ_f ≤ t) - Φ(t)| ≤ C[\text{moment terms} + \text{bias terms}]\]

5. Gevrey analysis: For infinitely differentiable functionals, balance stochastic error O(√(ν/n)) against remainder with truncation s = ⌊log(en)⌋:

   \[\|D^s f\|_{\text{Lip}^1} ≤ \|f\|_{G^α(U),R} R^{s+1} \{(s+1)!\}^α\]

   The choice s ≍ log(n) optimally balances these terms under Gevrey regularity.

Purely theoretical with applications to:

1. Precision matrix estimation: η₁ᵀΣ⁻¹η₂ under d log²(en) = o(n) without sparsity
2. Linear regression projections: c ᵀ Σ⁻¹Γ estimation

3. Covariance operator functionals: tr(Σ), ‖Σ‖²_HS in Hilbert spaces

   Empirical validation would require:

   • Finite-sample performance vs. plug-in and bootstrap methods
   • Computational timing for polynomial-time relaxations
   • Dimension regimes where d ≈ n/log²(n) in practice
   • Robustness to heavy-tailed distributions under moment assumptions

Limitations:

1. Sample splitting reduces effective sample size by factor 2 - NATURAL (standard tradeoff in cross-fitting)

2. Fourth moment conditions for Berry-Essén bounds - NATURAL (optimal for CLT-type results)

3. Positive definiteness assumption for matrix applications - NATURAL (inherent to precision matrix problems)

4. Computational complexity O(n^s) for exact U-statistics - TECHNICAL (addressed by polynomial relaxations)

5. Gevrey regularity assumption for infinitely smooth case - RESTRICTIVE (substantial class but not all smooth functionals)

6. Pilot estimator must satisfy probability tail bounds - TECHNICAL (achievable by robust estimators)

   Open problems:

7. Extend dimension regime beyond d log²(en) = o(n) for specific functional classes
8. Develop adaptive methods that automatically choose truncation level s without knowing functional smoothness