---
title: "A Random Matrix Theory Perspective on the Consistency of Diffusion Models"
authors:
- Binxu Wang
- Jacob A. Zavatone-Veth
- Cengiz Pehlevan
abstract: |
Diffusion models trained on different, non-overlapping subsets of a dataset often produce strikingly similar outputs when given the same noise seed. We trace this consistency to a simple linear effect: the shared Gaussian statistics across splits already predict much of the generated images. To formalize this, we develop a random matrix theory (RMT) framework that quantifies how finite datasets shape the expectation and variance of the learned denoiser and sampling map in the linear setting. For expectations, sampling variability acts as a renormalization of the noise level through a self-consistent relation $\sigma^2 \mapsto \kappa(\sigma^2)$, explaining why limited data overshrink low-variance directions and pull samples toward the dataset mean. For fluctuations, our variance formulas reveal three key factors behind cross-split disagreement: anisotropy across eigenmodes, inhomogeneity across inputs, and overall scaling with dataset size. Extending deterministic-equivalence tools to fractional matrix powers further allows us to analyze entire sampling trajectories. The theory sharply predicts the behavior of linear diffusion models, and we validate its predictions on UNet and DiT architectures in their non-memorization regime, identifying where and how samples deviate across training data splits. This provides a principled baseline for reproducibility in diffusion training, linking spectral properties of data to the stability of generative outputs.
arxiv: "2602.02908v2"
submitted: 2026-02-02
revised: 2026-07-04
venue: "ICML 2026 Oral; Outstanding Paper Honorable Mention"
---
# A Random Matrix Theory Perspective on the Consistency of Diffusion Models
**Binxu Wang, Jacob A. Zavatone-Veth, Cengiz Pehlevan**
Kempner Institute, Society of Fellows, John A. Paulson School of Engineering and Applied Sciences, and Center for Brain Science, Harvard University.
ArXiv v2, revised 2026-07-04. Accepted as an ICML 2026 oral and selected for an Outstanding Paper Honorable Mention.
# Abstract
Diffusion models trained on different, non-overlapping subsets of a dataset often produce strikingly similar outputs when given the same noise seed. We trace this consistency to a simple linear effect: the shared Gaussian statistics across splits already predict much of the generated images. To formalize this, we develop a random matrix theory (RMT) framework that quantifies how finite datasets shape the expectation and variance of the learned denoiser and sampling map in the linear setting. For expectations, sampling variability acts as a renormalization of the noise level through a self-consistent relation $\sigma^2 \mapsto \kappa(\sigma^2)$, explaining why limited data overshrink low-variance directions and pull samples toward the dataset mean. For fluctuations, our variance formulas reveal three key factors behind cross-split disagreement: anisotropy across eigenmodes, inhomogeneity across inputs, and overall scaling with dataset size. Extending deterministic-equivalence tools to fractional matrix powers further allows us to analyze entire sampling trajectories. The theory sharply predicts the behavior of linear diffusion models, and we validate its predictions on UNet and DiT architectures in their non-memorization regime, identifying where and how samples deviate across training data splits. This provides a principled baseline for reproducibility in diffusion training, linking spectral properties of data to the stability of generative outputs.
# Introduction
Diffusion models and their relatives such as flow matching have become the dominant generative modeling paradigm across diverse domains, including images, video, and proteins. By learning a time-dependent vector field, these models transform Gaussian noise into structured samples through an ordinary differential equation (ODE) or its stochastic variants [@song2021scorebased; @albergo2023stochasticinterpolants].
A distinctive feature of diffusion models is their striking *consistency across training runs* (Figure `\ref{fig:motivation_schematics}`{=latex}). When trained on the same distribution, even with disjoint datasets, different architectures, or repeated initializations, diffusion models often map the same noise seed to highly similar outputs under the deterministic probability flow [@kadkhodaie2024generalization; @zhang2024emergence]. This phenomenon contrasts with other generative modeling frameworks including GANs and VAEs, where the isotropic Gaussian latent space admits arbitrary rotations, leading to run-to-run variability in the mapping from latent codes to data [@Martinez2022ReproducibleRVAE].
#### Why does consistency matter?
Consistency across non-overlapping data splits suggests that diffusion models recover aspects of the underlying *data manifold* that are insensitive to the specific training set. This raises fundamental questions about how such models generalize beyond their training samples, to what extent they memorize idiosyncratic data, and whether their outputs reflect universal statistical regularities of the distribution. These issues connect to emerging theoretical and empirical debates on generalization, memorization, and creativity in diffusion models [@kamb2024creative; @niedoba2024towards; @kadkhodaie2024generalization; @chen2025interpolationEffectScoreSmoothing; @vastola2025generalization; @bonnaire2025diffusionRMTRegularize]; see also further discussion in App. `\ref{app:ext_related_works}`{=latex}.
#### Our approach.
We analyze this phenomenon through the lens of random matrix theory (RMT), beginning with the observation that the consistency effect can already be predicted by a linear Gaussian model (Fig. `\ref{fig:motivation_schematics}`{=latex}). Building on the linear denoiser framework, we develop a precise RMT analysis of how finite-sample variability in the empirical covariance affects both the expectation and fluctuation of denoisers and sampling maps (Fig. `\ref{fig:noise_renorm_eff}`{=latex}**A**). We then validate these theoretical predictions against deep diffusion models (CNNs and DiTs), showing that the same RMT principles still govern their inhomogeneity of consistency across data splits. Our **main contributions** are as follows:
- **Linear origin of consistency:** We show that shared Gaussian statistics i.e. linear denoiser already predict cross-split agreement.
- **Finite-sample RMT:** We show that randomness enters through a *renormalized* noise scale $\sigma^2\!\mapsto\!\kappa(\sigma^2)$, explaining overshrinkage of low-variance modes.
- **Variance law:** We derive a factorized form for fluctuation---explaining structures in cross-split deviation: anisotropy across eigenmodes, inhomogeneity across inputs, and global scaling with training set size.
- **Fractional-power equivalence:** We use integral formulae for fractional matrix powers to derive deterministic equivalents (DE) for full sampling trajectories.
- **Deep-net validation:** We qualitatively confirm overshrinkage, anisotropy, and inhomogeneity phenomenon in UNet and DiT models beyond the linear regime.
# Notation and Set up
#### Score-based Diffusion Models.
Let $p_0(\mathbf{x})$ be the target data distribution. For each noise scale $\sigma>0$, define the noised distribution as $p(\mathbf{x};\sigma)
= \bigl(p_0 * \mathcal{N}(0,\sigma^2\mathbf I)\bigr)(\mathbf{x})
= \int p_0(\mathbf{y})\,\mathcal{N}(\mathbf{x}\mid\mathbf{y},\sigma^2\mathbf I)\,d\mathbf{y}.$ The corresponding *score function* is $\nabla_{\mathbf x}\log p(\mathbf{x};\sigma),$ i.e. the gradient of the log--density. In the EDM formulation [@karras2022elucidatingDesignSp], the probability flow ODE (`\ref{eq:probflow_ode_edm}`{=latex}-ODE) reads, $$\label{eq:probflow_ode_edm}
\frac{d\mathbf{x}}{d\sigma}
= -\,\sigma\,\nabla_{\mathbf x}\log p(\mathbf{x};\sigma).\tag{PF}$$ This ODE transports samples from $p(\,\cdot\,;\sigma_2)$ to $p(\,\cdot\,;\sigma_1)$ when integrating $\sigma$ from $\sigma_2$ to $\sigma_1$. In particular, by starting from Gaussian noise $\mathcal N(0,\sigma_T^2 I)$ and integrating the PF-ODE from a sufficiently large $\sigma_T$ down to $\sigma=0$, one recovers clean samples from $p_0$. We adopt the EDM parametrization for its notational simplicity; other common diffusion formalisms are equivalent up to simple rescalings of time and space [@karras2022elucidatingDesignSp].
To estimate the score function of $p_0(\mathbf{x})$, we minimize the denoising score matching (DSM) objective [@vincent2011connection]. Parametrizing the denoiser $\mathbf{D}_\theta$ with a function approximator, then at noise level $\sigma$ the DSM objective reads `\small`{=latex} $$\begin{aligned}
\label{eq:dsm-obj}
\mathcal{L}_\sigma \;=\; \mathbb{E}_{\mathbf{x}_0\sim p_0,\; \mathbf{z} \sim \mathcal{N}(0,\mathbf{I})} \Bigl\|\mathbf{D}_\theta (\mathbf{x}_0 + \sigma \mathbf{z};\sigma)\;-\;\mathbf{x}_0\Bigr\|_2^2.\tag{DSM}
\end{aligned}$$ `\normalsize`{=latex} Using Tweedie's formula [@robbins1992empiricalTweedies], the score function can be obtained via the optimized \`denoiser' $\mathbf{s}_\theta(\mathbf{x},\sigma)=(\mathbf{D}_\theta(\mathbf{x},\sigma)-\mathbf{x})/\sigma^2$. In practice [@karras2022elucidatingDesignSp], diffusion models balance these scale-specific objectives with a weighting function $w(\sigma)$, yielding the overall training loss $\mathcal{L}=\int_\sigma d\sigma\,\, w(\sigma)\, \mathcal{L}_\sigma$.
#### Data distribution.
Consider a ground truth data distribution $p_0(\mathbf{x}),\ \mathbf{x}\in\mathbb{R}^d$, with population mean $\boldsymbol{\mu}$ and covariance $\boldsymbol{\Sigma}$. From this underlying distribution, we construct an empirical distribution $\{\mathbf{x}_i\}$ with $n$ samples, stacked as $X\in \mathbb{R}^{n\times d}$, then we denote the empirical mean $\hat{\boldsymbol{\mu}}$ and covariance $\hat{\boldsymbol{\Sigma}}$.
Here we are interested in the effect of the number of samples $n$, and different realizations of $X$ on the expectation (mean) and fluctuation (variance) of learned diffusion model. More specifically, we will study how randomness in the empirical covariance $\hat{\boldsymbol{\Sigma}}$ drives variability in the denoiser, relative to the population covariance $\boldsymbol{\Sigma}$.
#### Linear Denoiser.
A tractable setting for analytical study is the linear denoiser $$\begin{aligned}
\mathbf D(\mathbf{x};\sigma)=\mathbf{W}_\sigma\,\mathbf{x}+\mathbf b_\sigma,
\end{aligned}$$ which is an affine function of the noised state, independent across noise scales. As in linear regression, the training data enters the learned denoiser only through their first two moments [@Hastie2019Surprises]. More explicitly, minimizing `\ref{eq:dsm-obj}`{=latex} $\mathcal{L}_\sigma$ for the empirical dataset $p_0=\{\mathbf{x}_i\}$`\footnotemark`{=latex} yields the optimal empirical linear denoiser, depending on $\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\Sigma}}$. `\footnotetext{With $n$ samples, we average over infinite noise draws, so each sample is reused infinitely.}`{=latex} $$\begin{aligned}
\label{eq:empir_gauss_optim_denoiser}
\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^*(\mathbf{x};\sigma) \;=&\; \hat{\boldsymbol{\mu}} + (\hat{\boldsymbol{\Sigma}} + \sigma^2\mathbf{I})^{-1}\hat{\boldsymbol{\Sigma}}(\mathbf{x} - \hat{\boldsymbol{\mu}}).
\end{aligned}$$ For simplicity, we will later set $\hat{\boldsymbol{\mu}}=\boldsymbol{\mu}$ to isolate the effect of the empirical covariance $\hat{\boldsymbol{\Sigma}}$.
#### Sampling trajectory and sampling map.
Given an initial noise pattern $\mathbf{x}_{\sigma_{T}}\sim\mathcal N(0,\sigma_T^2I)$, the `\ref{eq:probflow_ode_edm}`{=latex} -ODE evolves it to a final sample $\mathbf{x}_0$. We refer to this mapping from $\mathbf{x}_{\sigma_{T}}$ to $\mathbf{x}_0$ as the *sampling map*; the phenomenon of consistency is precisely about the stability of this mapping across different realizations of training data. When the denoiser is linear and optimal at each noise scale, the PF-ODE can be solved in closed-form by projecting onto the eigenbasis of the data, yielding the analytic sampling trajectory [@wang_unreasonable_2024; @pierret2024diffusion4Gaussian]. $$\begin{aligned}
\label{eq:PFODE_sol_gauss}
&\mathbf{x}_{\hat{\boldsymbol{\Sigma}}}(\mathbf{x}_{\sigma_{T}},{\sigma}) \\
& =\hat{\boldsymbol{\mu}}+(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}(\mathbf{x}_{\sigma_{T}}-\hat{\boldsymbol{\mu}})\notag
\end{aligned}$$ Taking $\sigma\to0$ recovers the Wiener filter with Gaussian prior [@wiener1964extrapolation], which has been shown to be a strong predictor of the sampling map of the learned diffusion networks [@wang_unreasonable_2024; @li2024understanding; @Lukoianov2025Locality]. In the linear case, the mapping remains affine in the initial state, with the matrix $\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}$ emerging as the central object of analysis.
# Motivating Empirical Observation {#sec:motivation}
We begin with a simple experiment illustrating the consistency phenomenon. We train UNet-CNN [@song2019generative] and DiT [@Peebles_2023_ICCV] diffusion models under the EDM framework [@karras2022elucidatingDesignSp], each on two non-overlapping splits of FFHQ32 (30k images each; details in App. `\ref{app:DNN_method}`{=latex}). When sampling from the same noise seed with a deterministic solver, the outputs are visually similar across both splits and architectures (Fig. `\ref{fig:motivation_schematics}`{=latex}A). Quantification via pixel MSE confirms this effect: generated images are more similar across splits than to their nearest neighbors in the training set (Fig. `\ref{fig:motivation_schematics}`{=latex}B), ruling out memorization [@kadkhodaie2024generalization; @zhang2024emergence].
Strikingly, the linear Gaussian predictor (Wiener filter) [@wang_unreasonable_2024] already accounts for much of this behavior. Using the empirical mean and covariance $(\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\Sigma}})$ of each split in Eq. `\ref{eq:PFODE_sol_gauss}`{=latex}, the linear predictor yields nearly identical outputs across splits, also sharing visual similarities with CNN and DiT results (Fig. `\ref{fig:motivation_schematics}`{=latex}A,B). This suggests that consistency arises because different data splits share nearly identical Gaussian statistics, the only feature the linear diffusion can absorb [@Wang2025AnalySpec]. Pointwise, samples closer to the Gaussian solution are also more consistent across splits (Pearson $r=0.244, p=5\times 10^{-15}$), suggesting convergence toward the Gaussian predictor underlies consistency. More visual examples and quantitative comparisons for other datasets (CIFAR10, CIFAR100, FFHQ at 32 and 64 pixels, LSUN church and bedroom dataset at 32 and 64 pixels) are provided in Appendix `\ref{app:extended_visual_motivation}`{=latex} and [our website](https://animadversio.github.io/diffusion-consistency-rmt/).
In summary, (i) diffusion models trained on independent splits converge to nearly identical sampling maps, (ii) this property holds across architectures, and (iii) a simple Gaussian predictor already captures much of the effect. While linear diffusion is more consistent than deep networks---which can exploit higher-order statistics---it provides a necessary baseline: *if Gaussian statistics differ, deep models may not yield consistent samples*. To test this directly, we perform a counterfactual experiment in which the training data are intentionally partitioned to induce mismatches in mean and/or variance by stratifying samples along a chosen principal component. Diffusion models trained on these statistically misaligned splits produce markedly less consistent generations (Fig. `\ref{suppfig:counterfactual_visual}`{=latex}, `\ref{suppfig:counterfactual_quantitative}`{=latex}). These observations motivate our random matrix theory analysis of finite-sample effects.
# Theory of Diffusion Consistency Across Independent Data
The goal of the theoretical analysis is to calculate the expectation and covariance of various quantities in diffusion model under independent instantiation of dataset (Fig. `\ref{fig:noise_renorm_eff}`{=latex}**A**).
## Self consistency equation and renormalized noise scale
#### Deterministic equivalence of sample covariance.
Our central technical tool is deterministic equivalence [@Potters2020First; @Bun2015Rotational], which allows random matrices to be replaced by deterministic surrogates---an approximation that becomes exact in the large-dimensional limit. In particular, we rely on the deterministic equivalence relation for the empirical covariance matrix $\hat{\boldsymbol{\Sigma}}$ [@atanasov2024scalingRenormalization; @bach2024highDimRMT_doubleDescent], $$\label{eq:DE_cov_equiv}
\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+\lambda I)^{-1}\asymp \boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\lambda) I)^{-1} \tag{DE}$$ where $\kappa$ is the unique positive solution to the self-consistent equation [@silverstein1995strongConvergence; @Marcenko1967Distribution]. $$\begin{aligned}
\label{eq:silverstein_self_consistent}
\kappa(\lambda)-\lambda&=\gamma\kappa(\lambda)\int_{0}^{\infty}\frac{sd\nu(s)}{\kappa(\lambda)+s} \nonumber \\
&=\gamma\kappa(\lambda)\operatorname{tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\lambda)I)^{-1}],
\end{aligned}$$ where $\gamma=d/n$ is the aspect ratio, and $\nu$ is the (limiting) spectral measure of $\boldsymbol{\Sigma}$.[^1] Note we use $\operatorname{tr}$ to denote the *normalized trace*, such that $\operatorname{tr}[I]=1$, and $\operatorname{Tr}$ the unnormalized one. More elaborate two-point deterministic equivalences [@bach2024highDimRMT_doubleDescent; @Atanasov2024Risk; @Atanasov2025TwoPoint] are required to derive the variance results in the paper, which can be found in Appendix `\ref{app:deterministic_equivalence}`{=latex}.
#### Property of renormalized noise $\kappa(\sigma^2)$.
As Eq. `\ref{eq:DE_cov_equiv}`{=latex} suggests, with trace-like measurement, the stochastic effects of sample covariance $\hat{\boldsymbol{\Sigma}}$ can be absorbed into the scalar $\kappa(\lambda)$ leaving the population covariance $\boldsymbol{\Sigma}$ otherwise unchanged, similar to the renormalization of self-energy in field theory [@Azee2010quantumnutshell; @atanasov2024scalingRenormalization; @Hastie2019Surprises; @bach2024highDimRMT_doubleDescent]. In our context, $\lambda$ usually corresponds to noise variance $\sigma^2$, so we could understand $\kappa$ as the renormalized noise variance. To build intuition, we numerically evaluate this nonlinear mapping using the spectrum of natural images (FFHQ) (Fig. `\ref{fig:noise_renorm_eff}`{=latex}**B**, Method in `\ref{app:numerical_methods}`{=latex}). The renormalization effect $\kappa(\sigma^2)$ is most pronounced at low noise scales, and when the sample number is much fewer than the data dimension ($\gamma=d/n\gg1$).
#### Notation
We define the degrees-of-freedom functions $$\begin{aligned}
\begin{split}
\operatorname{df}_{1}(\lambda) &:=\operatorname{Tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\lambda I)^{-1}],
\\
\operatorname{df}_{2}(\lambda) &:=\operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda I)^{-2}],
\\
\operatorname{df}_{2}(\lambda,\lambda^\prime) &:=\operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda I)^{-1}(\boldsymbol{\Sigma}+\lambda^\prime I)^{-1}].
\end{split}
\end{aligned}$$ We have $0 \leq \operatorname{df}_{2}(\lambda) \leq \operatorname{df}_{1}(\lambda) \leq \operatorname{rank}(\boldsymbol{\Sigma}) \leq d$, and $\operatorname{df}_{2}(\lambda,\lambda') \leq \sqrt{\operatorname{df}_{2}(\lambda) \operatorname{df}_{2}(\lambda')}$. Moreover, on the solution of the Silverstein equation for $\lambda > 0$, we have the strict inequality $\operatorname{df}_{1}(\kappa) < n$.
## Expectation: Finite Data Renormalize Noise Scales
Next we apply these tools to compute the expectation and fluctuation of the denoiser under dataset realizations. The form of Eq. `\ref{eq:empir_gauss_optim_denoiser}`{=latex} naturally suggests the deterministic equivalence in Eq. `\ref{eq:DE_cov_equiv}`{=latex}, leading to the following result.
::: proposition
**Result 1** (Deterministic equivalent of the denoiser expectation). *`\label{prop:denoiser_exp_DE}`{=latex} Assuming $\hat{\boldsymbol{\mu}}=\boldsymbol{\mu}$, and given a fixed probe vector $\mathbf{v}\in \mathbb{R}^d$, then the optimal empirical linear denoiser has the following deterministic equivalent. (Derivation in App. `\ref{app:proof_denoiser_exp_de}`{=latex}). `\small`{=latex} $$\begin{aligned}
\label{eq:denoiser_onestep_mean_DE}
\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\Big] & \asymp\mathbf{v}^{\top}\mathbf{D}_{\boldsymbol{\Sigma}}^{*}(\mathbf{x};\sqrt{\kappa(\sigma^{2})})\\
&=\mathbf{v}^{\top}\Big[\boldsymbol{\mu}+\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}(\mathbf{x}-\boldsymbol{\mu})\Big]\notag.
\end{aligned}$$*
:::
#### Interpretation
In expectation, finite data act by renormalizing the noise scale, $\sigma^{2} \to \kappa(\sigma^{2})$, in the population denoiser. This is equivalent to adding an adaptive Ridge penalty to the `\ref{eq:dsm-obj}`{=latex} objective. Compared to the population solution $\mathbf{D}_{\boldsymbol{\Sigma}}^{*}$, the finite-sample denoiser shrinks low-variance directions more aggressively, treating them as noise and pulling outputs toward the dataset mean (Fig. `\ref{fig:noise_renorm_eff}`{=latex}**D**). Numerically, deviations are indeed most pronounced in the lower spectrum and at lower noise levels, where the renormalization effect is the strongest (Fig. `\ref{fig:noise_renorm_eff}`{=latex}**C**). Since smaller noise scale is associated with generation of high frequency details in image, this result suggests these detail eigenmodes take more samples to be learned correctly, which we will confirm in next section.
## Fluctuation: Anisotropic and Inhomogeneity of Denoiser Consistency
Next, we tackle the fluctuation due to dataset realizations, which addresses the consistency of diffusion models trained on independent data splits. We derive the following equivalence using two-point and one-point deterministic equivalence identities (Eq. `\ref{eq:DE_mod_resolv_quad_two_z}`{=latex},`\ref{eq:DE_mod_resolv_lin_one_z}`{=latex}, @bach2024highDimRMT_doubleDescent).
::: proposition
**Result 2** (Deterministic equivalent of the denoiser variance). *`\label{prop:denoiser_var_DE}`{=latex} Assuming $\hat{\boldsymbol{\mu}}=\boldsymbol{\mu}$, across dataset realizations of size $n$, the covariance of the optimal empirical linear denoiser is denoted as $\mathcal{S}_{D}$; at point $\mathbf{x}$, along direction $\mathbf{v}$, the variance is given by $\mathbf{v}^\top \mathcal{S}_{D}(\mathbf{x}) \mathbf{v}$, which admits the following deterministic equivalence. $$\begin{aligned}
\label{eq:denoiser_var_deter_equiv}
&\mathbf{v}^{\top}\mathcal{S}_{D}(\mathbf{x})\mathbf{v} :=\operatorname{Var}_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)] \\
&\asymp {\frac{\kappa(\sigma^{2})^{2}}{n-\operatorname{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}}
\underbrace{\Diamond(\mathbf{v},\kappa(\sigma^{2}),\boldsymbol{\Sigma})}_{\text{anisotropy}}\
\underbrace{\Diamond(\mathbf{x}-\boldsymbol{\mu},\kappa(\sigma^{2}),\boldsymbol{\Sigma})}_{\text{inhomogeneity}} \notag
\end{aligned}$$ where $\Diamond(\mathbf{u},\kappa,\boldsymbol{\Sigma}):=\mathbf{u}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{u}$. Derivation in App. `\ref{app:proof_denoiser_var_de}`{=latex}.[^2]*
:::
#### Interpretation.
The variance of denoiser across dataset realizations factorizes into three interpretable components: a dependence on probe direction (*anisotropy*), a dependence on noised sample location (*inhomogeneity*), and an overall scale with $n$ and $\sigma$ (*global scaling*). Note, given the relation between score and denoiser, the covariance of score is $\sigma^{-4}\mathcal{S}_{D}(\mathbf{x})$, *i.e.* all results translate by scaling.
#### Anisotropy in probe direction.
The anisotropy of consistency is governed by $\Diamond(\mathbf{v},\kappa,\boldsymbol{\Sigma})$. When the probe $\mathbf{v}$ aligns with a principal component (PC) $\mathbf u_k$ of $\boldsymbol{\Sigma}$ with eigenvalue $\lambda_k$, this reduces to $\chi(\lambda_k,\kappa):=\lambda_k/(\lambda_k+\kappa)^2$. The function $\chi(\lambda,\kappa)$ is bell-shaped in $\lambda$, uniquely maximized at $\lambda=\kappa$ with value $1/(4\kappa)$. Thus, for each noise scale, the directions of greatest uncertainty are precisely those whose variances match the renormalized noise $\kappa(\sigma^2)$ (Fig. `\ref{fig:denoiser_var_anisotropy_scaling}`{=latex}**B**). This effect is evident visually. For linear denoisers trained on non-overlapping splits of human face dataset (FFHQ), their differences follow the spectral structure of natural images [@ruderman1994statNatImage]: at high noise the deviations appear as low-frequency facial envelopes, while at low noise they shift to high-frequency specular patterns (Fig. `\ref{fig:denoiser_var_anisotropy_scaling}`{=latex}**B**). Quantitatively, the MSE between two denoisers along each PC matches the variance prediction of Eq. `\ref{eq:denoiser_var_deter_equiv}`{=latex}, with the expected factor of two from independent sampling (Lemma `\ref{lemma:dev_double_var}`{=latex})`\footnotemark`{=latex}. `\footnotetext{Note that the empirical analyses in our paper, as well as in prior works, examine consistency across diffusion models trained on \textit{non-overlapping} data splits. This differs slightly from \textit{independently sampled} subsets when the full dataset is finite. Nevertheless, the agreement between theory and linear denoiser suggests this correction is negligible.}`{=latex}
#### Inhomogeneity in input location.
The inhomogeneity of denoiser variance across input space is governed by $\Diamond(\mathbf{x}-\boldsymbol{\mu},\kappa,\boldsymbol{\Sigma})$. While structurally similar to the anisotropy factor, here $\mathbf{x}-\boldsymbol{\mu}$ is drawn from the noised data distribution rather than a unit probe. Approximating $\mathbf{x}-\boldsymbol{\mu}$ as lying on the ellipsoidal shell of $\mathcal{N}(0,\boldsymbol{\Sigma}+\sigma^2I)$, its displacement along eigenvector $\mathbf u_k$ has typical radius $\sqrt{\sigma^2+\lambda_k}$. Substituting gives $\Diamond(\sqrt{\sigma^2+\lambda_k}\,\mathbf u_k,\kappa,\boldsymbol{\Sigma})
= (\sigma^2+\lambda_k)\,\chi(\lambda_k,\kappa)$. Unlike the pure anisotropy factor, this expression grows monotonically with $\lambda_k$. Thus, denoiser variability is amplified for inputs displaced along high-variance modes, yielding larger uncertainty for such locations (Fig. `\ref{fig:denoiser_var_anisotropy_scaling}`{=latex}**C**), which agree quantitatively with numerical results. Based on this factor, denoiser consistency can be predicted for noisy images point by point (e.g. Pearson $r=0.94$ across noised images, at $\sigma^2=1$, $n=1000$, Fig. `\ref{suppfig:inhomogeneity_theory_vs_linear_denoiser}`{=latex}).
#### Global scaling with sample size.
Finally, marginalizing over all directions and noised samples yields a closed-form expression for the overall denoiser variance (Eq. `\ref{eq:denoiser_var_overall_scaling}`{=latex}, Fig. `\ref{fig:denoiser_var_anisotropy_scaling}`{=latex}**D**). At large $n$ limit, denoiser variance scale inversely with sample number $n^{-1}$, reminiscent of classic statistical laws; while at smaller $n$, the renormalization effects modify the scaling.
#### Summary.
In sum, the variance structure reveals three key effects. *Anisotropy*: uncertainty is maximized along eigenmodes whose variance $\lambda_k$ is comparable to the renormalized noise $\kappa(\sigma^2)$. *Inhomogeneity*: noised points displaced along high-variance directions experience larger uncertainty. *Scaling*: the overall variance shrinks with dataset size $n$, recovering the population model in the large-sample limit. Together, these predictions yield a detailed spatial and spectral map of where denoisers trained on different data splits are most likely to disagree.
# Consistency of Diffusion Samples for Linear Denoisers
Beyond the consistency of single-step denoiser output or score, we are interested in the final diffusion sample from the same initial noise seed $\mathbf{x}_{\sigma_T}$. For linear denoisers, sampling map from initial noise to generated sample is captured by Wiener filter (Eq. `\ref{eq:PFODE_sol_gauss}`{=latex}, $\sigma=0$). However, unlike one-step denoiser, this mapping involves fractional power of covariances $\boldsymbol{\Sigma}^{1/2}(\boldsymbol{\Sigma}+\sigma^2 I)^{-1/2}$, for which the deterministic equivalence is not readily available. Here, we leveraged the integral representation of fractional power (@Balakrishnan1960FractionalPower's formula) and deterministic equivalence, and arrived at a few novel equivalence of these matrices (Result `\ref{prop_app:deter_equiv_cov_mat_half}`{=latex}, `\ref{prop_app:deter_equiv_cov_half_resolv}`{=latex}, Derivation in App. `\ref{app:Balakrishnan_identity}`{=latex}). Using these developments, we can calculate the expectation and fluctuation of sampling map and, consequently, of samples generated from a fixed seed.
## Expectation of diffusion sample: over-shrinkage to the mean
We note that when the initial noise scale $\sigma_T$ is large, the sampling map admits the approximation $$\begin{aligned}
\mathbf{x}_{\hat{\boldsymbol{\Sigma}}}(\mathbf{x}_{\sigma_{T}},0)&=\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}(\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu})\notag \\
&\approx\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}},
\end{aligned}$$ where we denote the shift and normalized noise $\bar{\mathbf{x}}:=\frac{\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu}}{\sigma_{T}}$. At the $\sigma_T\to\infty$ limit, this approximation becomes exact, and $\bar{\mathbf{x}}\sim\mathcal N(0,I)$. For clarity, we present results under this infinite-$\sigma_T$ approximation; the expressions accounting for finite $\sigma_T$ effects are provided in App. `\ref{app:proof_samp_mapping_exp_full}`{=latex}.
::: proposition
**Result 3** (Deterministic equivalence for expectation of diffusion sampling map). *`\label{prop:gen_sample_exp_DE}`{=latex} The sample generated from initial state $\mathbf{x}_{\sigma_{T}}$ has the following deterministic equivalence. Derivation in App. `\ref{app:proof_samp_mapping_exp_approx}`{=latex}. $$\begin{aligned}
\label{eq:gen_sample_exp_DE}
\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}[\mathbf{x}_{\hat{\boldsymbol{\Sigma}}}(\mathbf{x}_{\sigma_{T}},0)] &\approx\boldsymbol{\mu}+\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\big[\hat{\boldsymbol{\Sigma}}^{1/2}\big]\frac{\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu}}{\sigma_{T}}\\
&\asymp\boldsymbol{\mu}+\frac{2}{\pi}\int_{0}^{\infty}\boldsymbol{\Sigma}\Big(\boldsymbol{\Sigma}+\kappa(u^{2})I\Big)^{-1}\bar{\mathbf{x}}du.\notag
\end{aligned}$$*
:::
#### Interpretation
This expression mirrors the deterministic equivalence of denoisers (Eq. `\ref{eq:denoiser_onestep_mean_DE}`{=latex}), but with an integration over effective noise scales. Comparing to the population sampling map, where $\kappa(u^2)$ reduces to $u^2$, the finite data case integrates over a stronger shrink factor $\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}$ (since $\kappa(u^2)>u^2$), especially on the lower eigenmodes. This effect is confirmed with numerics of empirical covariance (Fig. `\ref{fig:diff_sample_RMT_expectation_variance}`{=latex}**A**). This leads to a systematic overshrinkage toward the dataset mean along these modes, reducing the generated variance along lower-variance directions [^3].
## Variance of diffusion sample: Anisotropy and inhomogeneity
::: proposition
**Result 4** (Deterministic equivalence for variance of diffusion sampling map). *`\label{prop:gen_sample_var_DE}`{=latex} Due to dataset realization, the variance of generated sample starting from initial state $\mathbf{x}_{\sigma_T}$, along vector $\mathbf{v}$ admits the following deterministic equivalence, `\small`{=latex} $$\begin{aligned}
\label{eq:gen_sample_var_DE}
&\operatorname{Var}_{\hat{\boldsymbol{\Sigma}}}\!\big[\mathbf v^{\top}\mathbf{x}_{\hat{\boldsymbol{\Sigma}}}(\mathbf{x}_{\sigma_T},0)\big]
= \operatorname{Var}_{\hat{\boldsymbol{\Sigma}}}\!\big[\mathbf v^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf x}\big] \\
&\asymp \frac{4}{\pi^{2}}\!\int_{0}^{\infty}\!\!\!\int_{0}^{\infty}
\frac{\kappa\,\kappa'}{n-\operatorname{df}_{2}(\kappa,\kappa')}\,
\underbrace{\Psi(\mathbf v;\kappa,\kappa',\boldsymbol{\Sigma})}_{\text{anisotropy}}
\underbrace{\Psi(\bar{\mathbf x};\kappa,\kappa',\boldsymbol{\Sigma})}_{\text{inhomogeneity}}\notag
\,du\,dv
\end{aligned}$$ `\normalsize`{=latex} where $\Psi(\mathbf a;\kappa,\kappa',\boldsymbol{\Sigma})
:= \mathbf a^{\top}\,\boldsymbol{\Sigma}\,(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa' I)^{-1}\,\mathbf a$ , and $\kappa:=\kappa(u^{2}),\kappa^{\prime}:=\kappa(v^{2})$ are variables to be integrated over. Derivation in App. `\ref{app:proof_samp_mapping_var_approx}`{=latex}.*
:::
#### Interpretation
The variance of sampling map Eq. `\ref{eq:gen_sample_var_DE}`{=latex} simplifies to a double integral of the denoiser-variance (Eq. `\ref{eq:denoiser_var_deter_equiv}`{=latex}). The integrand factorizes into a direction-dependent term (*anisotropy*), an initial noise-dependent term (*inhomogeneity*), and a scaling term. Note the anisotropy and inhomogeneity factors rely on the same $\Psi(.;\kappa,\kappa',\boldsymbol{\Sigma})$ function, showing that dependency on $\mathbf{v}$ and $\bar\mathbf{x}$ has the same spectral structure. As in Result `\ref{prop:denoiser_var_DE}`{=latex}, this result is the generic factorized contribution to the variance, which dominates for high-dimensional Gaussian $\bar{\mathbf{x}}$ when the data resolvent has high effective rank.
We resort to numerical simulation to provide more intuition. We note that integrals in Eqs. `\ref{eq:gen_sample_exp_DE}`{=latex},`\ref{eq:gen_sample_var_DE}`{=latex} are nontrivial to evaluate; we describe our numerical scheme in App. `\ref{app:numerical_methods}`{=latex}. Using this procedure, the theoretical predictions align closely with direct computations of linear diffusion (Fig. `\ref{fig:diff_sample_RMT_expectation_variance}`{=latex}). **Inhomogeneity.** Spatially, when initial noise $\bar\mathbf{x}$ deviates more along the top eigenspace of $\boldsymbol{\Sigma}$, there will be larger uncertainty (Fig. `\ref{fig:diff_sample_RMT_expectation_variance}`{=latex}**C**), this enables us to predict the sample difference point by point. **Anisotropy.** Directionally, the dependency on $\mathbf{v}$ has the same structure, in absolute term, the deviation is larger at higher eigenspace (Fig. `\ref{fig:diff_sample_RMT_expectation_variance}`{=latex}**B**). Note that when scaling up the dataset size, the variance in the top eigenspace decay immediately from small sample size; while the deviation in lower eigenspace will stay put and start decaying only later at larger dataset size (Fig. `\ref{fig:diff_sample_RMT_expectation_variance}`{=latex}**D**). This shows that the fine detail of the samples needs a larger dataset size to be consistency across training.
# Validating Predictions on Deep Networks
Finally, given that linear diffusion behavior is well captured by our random matrix theory (RMT), we test the applicability of its prediction to practical deep diffusion networks.
#### Setup.
We trained UNet- and DiT-based denoisers under the EDM framework on FFHQ64, FFHQ32, AFHQ32 [@choi2020starganv2], LSUN church and bedroom at 32 and 64 pixels [@yu2015lsun], CIFAR10, and CIFAR100 (UNet on all; DiT on a subset). For each dataset we trained on two non-overlapping splits at sizes $n \in \{300,\,1000,\,3000,\,10^4,\,3\!\cdot\!10^4\}$ (10 runs total per architecture). Sampling was performed with the same random seed using the Heun solver [@karras2022elucidatingDesignSp]. We train for 50,000 steps with Adam optimizer, further details are provided in App. `\ref{app:DNN_method}`{=latex}.
#### Expectation: from memorization to renormalization.
We observe a clear two-phase behavior as dataset size increases. *Memorization phase* ($n\leq1000$): models largely reproduce training samples (Fig. `\ref{fig:DNN_validation_results}`{=latex}**A,B**), and samples are much closer to the nearest neighbor in their training split than the control split, consistent with prior observations. This regime is outside the scope of linear theory, since linear score models cannot memorize individual points [@Wang2025AnalySpec]. The same transition occurs across architectures, though the dataset size $n$ at which it occurs depends on model capacity and image resolution (Fig. `\ref{suppfig:DNN_val_NN_AFHQ}`{=latex},`\ref{suppfig:DNN_val_NN_FFHQ32}`{=latex},`\ref{suppfig:DNN_val_NN_CIFAR}`{=latex}, see also App. `\ref{app:nearest_neighbor_across_dataset_size}`{=latex}). *Renormalization phase* ($n \geq 3000$): the samples have comparable distance to the neighbor in the training split and control split, showing generalization. Further, as $n$ grows, samples increasingly resemble and approach the linear predictors (Fig. `\ref{suppfig:DNN_approach_Gaussian}`{=latex}) [@li2024understanding]. In this regime, the overshrinkage predicted by Result `\ref{prop:gen_sample_exp_DE}`{=latex} becomes visible: generated face samples resemble the average face [@langlois1994avgface], with smoother textures and background (Fig. `\ref{fig:DNN_validation_results}`{=latex}**A**, $n=3000$). Quantitatively, we observe reduced variance along low- and mid-spectrum eigenmodes of the generated samples (Fig. `\ref{fig:DNN_validation_results}`{=latex}**D**, see also App. `\ref{app:spectral_anisotropy_var_consistency}`{=latex}). This bias decreases as dataset size increases, and vanishes when learned and population spectra coincide at $n\sim30000$.
#### Fluctuations: inhomogeneity of consistency.
Within the renormalization phase, RMT further predicts which initial noise and direction exhibit the largest discrepancies across data splits, due to their alignment with data covariance (Result `\ref{prop:gen_sample_var_DE}`{=latex}). Spectrally, measuring the cross-split deviation along population eigenbases, we can see the characteristic anisotropy profile. Further the decrease of MSE majorly occurs in top eigenspace, while the middle or lower eigenspace remains unchanged or becomes less consistent when sample size increases (Fig. `\ref{fig:DNN_validation_results}`{=latex}**E**). This is consistent with the prediction of the theory that lower eigenmodes needs more training samples to be consistent (Fig. `\ref{fig:diff_sample_RMT_expectation_variance}`{=latex}**B**, see also App. `\ref{app:spectral_anisotropy_var_consistency}`{=latex}). Spatially, the inhomogeneity effect is borne out: RMT predictions correlate with observed cross-split deviations for each initial noise noise; e.g., UNets trained on FFHQ64 with $n=30000$ achieve a Spearman correlation of $0.33$ ($p=2.5\times10^{-26}$) over 1000 seeds (Fig. `\ref{fig:DNN_validation_results}`{=latex}**F**). Remarkably, the prediction requires only the population covariance and dataset size, with no knowledge of split identities or network architecture. The absolute deviation magnitudes, however, are much larger in deep networks than predicted by linear theory, reflecting nonlinear source of idiosyncrasy across splits. As controls, correlations collapse in the memorization regime (Fig. `\ref{suppfig:DNN_val_inhomogeneity_FFHQ32}`{=latex},`\ref{suppfig:DNN_val_inhomogeneity_FFHQ64}`{=latex}, see App. `\ref{app:anisotropy_init_noise}`{=latex}) and disappear when mismatched noise seeds are used.
#### Summary.
Across architectures and datasets, the predictions of our linear RMT framework extend to deep diffusion models: limited data induce overshrinkage toward the mean, and the variance structure across splits exhibits inhomogeneity and anisotropy predicted by theory.
# Discussion
Our analysis shows that much of the consistency in diffusion models across training data is already captured by Gaussian statistics, which inspires an in-depth analysis of consistency in linear diffusion. Random matrix theory sharpens this picture by showing that finite data act through a renormalized noise scale $\sigma^2 \mapsto \kappa(\sigma^2)$, and that fluctuations across splits factor into anisotropy over eigenmodes, inhomogeneity across inputs, and a global scaling with dataset size. These results extend deterministic-equivalence tools to fractional matrix powers, allowing closed-form predictions for both denoisers and sampling trajectories, and align qualitatively with deep networks in terms of where deviation accentuates, even if nonlinear effects amplify the magnitudes of deviation.
#### Why are diffusion models special?
Previous work suggests that consistency and reproducibility are distinctive features of diffusion and flow-based models with deterministic samplers [@zhang2023emergenceRepro]. Why should this be the case? In diffusion models, the function approximator learns a score vector field $\nabla \log p(\mathbf{x};\sigma)$, which is uniquely determined by the data distribution[^4]. The lowest-order approximation to this score is a linear vector field determined by the Gaussian statistics of the data, namely its mean and covariance. Our random matrix analysis shows that these Gaussian statistics, and observables derived from them (e.g. denoiser and sampling map), are highly stable across independent data splits. Moreover, at high noise levels, the linear Gaussian score approximation is often accurate; since the high-variance, low-frequency aspects of samples are determined primarily at these high-noise scales, even nonlinear score networks can inherit substantial stability from this shared Gaussian structure. This provides a mechanism by which independently trained diffusion models can share a common mapping from initial noise to generated samples.
This mechanism is less directly applicable to other generative modeling paradigms. In autoregressive models, randomness typically enters through sequential sampling from conditional token distributions at nonzero temperature, so repeated generation depends on a chain of stochastic discrete choices rather than a deterministic flow from a fixed initial noise. In VAEs and GANs, samples are generated by a learned map from a usually lower-dimensional latent space, $x=f(z)$ with $z\sim\mathcal{N}(0,I)$. This latent space is not identifiable: composing $f$ with an orthogonal transformation of $z$ preserves the latent and generated distributions, but changes which latent vector corresponds to which generated sample. Thus, two independently trained models may learn rotated latent coordinate systems, so the same latent vector need not produce similar outputs. One could in principle align the latent spaces of VAEs or GANs by fixing this orthogonal freedom. For example, latent vectors can be expressed in a meaningful local basis, such as the right singular vector basis of the averaged generator Jacobian $\partial f/\partial z$. Previous work showed that, in this basis, axes of different GANs trained on human face datasets can have relatively consistent semantics, such as face orientation or background changes [@wang2021aGANGeom Fig. 10]. A full analysis of consistency in VAEs or GANs is left for future investigation.
#### Limitations
At the same time, our framework has limitations. Linear surrogates underestimate variability in expressive models and do not capture architecture-specific inductive biases. Extending the theory to random-feature models would better explain the transition from memorization to renormalization [@bonnaire2025diffusionRMTRegularize; @george2025diffusionRandFeatRMT], and help quantify how model capacity shifts the required dataset size for consistency. Another promising direction is to study the *anisotropy of the initial noise space* and its alignment with the data manifold. The seemingly unstructured noise space is already anchored by the data covariance before generation. Manipulating noise components along dataset principal-component subspaces leads to interpretable changes in the sampled outputs, such as suppressing background structure or reducing specific facial variations (Fig. `\ref{suppfig:anisotropic_structure_initial_noise_top50PC}`{=latex}, `\ref{suppfig:anisotropic_structure_initial_noise_top200PC}`{=latex}). Our analysis shows that different directions in the noise space contribute unequally to cross-split deviations (Result `\ref{prop:gen_sample_var_DE}`{=latex}); consequently, the alignment between the initial noise and the covariance eigenframe predicts the degree of generation consistency across dataset splits. Such anisotropic structures might explain why certain "magic" random seeds may consistently yield better generations [@xu2025goodseed]. This echoes anisotropic effects observed in GANs' latent space, where noise vectors aligned strongly with top eigenspaces of Jacobian can lead to degraded generations [@wang2021aGANGeom]. Such connections suggest that spectral geometry of the input space deserves closer attention as a unifying factor across generative models.
# Impact Statement {#impact-statement .unnumbered}
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
## Acknowledgments {#acknowledgments .unnumbered}
B.W. was supported by the Kempner Fellowship from the Kempner Institute at Harvard. J.Z.-V. was supported by a Junior Fellowship from the Harvard Society of Fellows. C.P. is supported by an NSF CAREER Award (IIS-2239780), DARPA grants DIAL-FP-038 and AIQ-HR00112520041, the Simons Collaboration on the Physics of Learning and Neural Computation, and the William F. Milton Fund from Harvard University. We thank David Alvarez-Melis for his insightful suggestions on the counterfactual moment manipulation experiments. This work has been made possible in part by a gift from the Chan Zuckerberg Initiative Foundation to establish the Kempner Institute for the Study of Natural and Artificial Intelligence, and by the generous computing resources provided by the Kempner Institute and Harvard Research Computing.
# Extended Related Work {#app:ext_related_works}
#### Consistency and Reproducibility in Diffusion
As a motivating observation, [@kadkhodaie2024generalization] found that diffusion models trained on non-overlapping splits of training data could produce visually highly similar images. The seminal paper studying this effect is [@zhang2024emergence]; there, the authors found that models trained on the same dataset have a consistent mapping from noise to sample across architectures (transformer vs. UNet), objectives, training runs, samplers, and noising kernels, provided that an ODE deterministic sampler is used. In their Appendix B, they also discuss the lack of reproducibility in VAEs and GANs. The consistency studied in our paper is more closely related to reproducibility in the generalization regime.
#### Hidden Linear Score Structure in Diffusion Models
Recent work has shown that, for much of diffusion time (*i.e.*, across a broad range of signal-to-noise ratios), the learned neural score is closely approximated by the linear score of a Gaussian fit to the data, which is usually the best linear approximation [@Binxu_2023_hidden; @li2024understanding]. Crucially, this Gaussian linear score admits a closed-form solution to the probability‐flow ODE, which can be exploited to accelerate sampling and improve its quality [@WangVastola_unreasonable]. Moreover, this same linear structure has been linked to the generalization--memorization transition in diffusion models [@li2024understanding]. In sum, across many noise levels, the Gaussian linear approximation captures many salient aspects of the learned score. Here, we leverage it to explain the observed consistency across splits and as a tractable setup for random matrix theory analysis.
#### Memorization, Generalization and Creativity in Diffusion
The question of when diffusion models can generate genuinely novel samples matters both scientifically and for mitigating data leakage. From the score-matching perspective, if the learned score exactly matches that of the empirical data distribution, then the reverse process reproduces that empirical distribution, and thus does not create new samples beyond the training set [@kamb2024creative; @li2024goodscoredoeslead; @wang_unreasonable_2024]. Yet high-quality diffusion models routinely generate images that are not identical copies of images from the training set. @kamb2024creative take an important step toward reconciling this: when the score network is a simple CNN, its inductive biases (locality and translation equivariance) favor patch wise composition, enabling global samples that are novel while remaining locally consistent "mosaics." Similarly, @Wang2025AnalySpec observed that score networks with different architectural constraints learn different approximations of the dataset and therefore generalize differently: e.g., linear networks learn the Gaussian approximation, and circular convolutional networks learn the stationary Gaussian process approximation. @finn2025originscreativity provided evidence that adding a final self-attention layer promotes global consistency across distant regions, organizing locally plausible features into coherent layouts that move beyond purely patch-level mosaics. This result is consistent with preliminary observations by @kamb2024creative regarding cases in which their purely convolutional models fail to generate coherent images, while models including attention succeed. Related theoretical work further probes why well-trained diffusion models can generalize despite apparent memorization pressures [@bonnaire2025diffusionRMTRegularize; @vastola2025generalization; @chen2025interpolationEffectScoreSmoothing]. These results suggest that departures from exact empirical-score fitting---mediated by inductive biases (both architectural and training dynamics) can explain how diffusion models avoid pure memorization while maintaining visual plausibility [@ambrogioni2023searchdispersedmemoriesgenerative].
# Extended Results and Figures
## Extended visual examples for the motivating observation {#app:extended_visual_motivation}
![**Motivating observation and the linear theory for FFHQ64 dataset**. Similar format to Fig. `\ref{fig:motivation_schematics}`{=latex}, but for FFHQ64 dataset. **A.** Examples of generated samples from the same noise seed, for UNet, DiT, and a linear denoiser on split 1 and split 2 of the data, each with 30k non-overlapping samples. The closest 4 samples in the corresponding training set are shown above and below the generated sample. One can appreciate the visual similarity of samples generated from models trained on separate splits, even with different neural architectures, and also with the linear denoiser on each split. Admittedly, the generated outcomes of linear denoisers at 64-pixel resolution look worse, especially around edges, showing signatures of non-Gaussian statistics, as [@wang_unreasonable_2024] has pointed out. **B.** Quantification of **A**, paired image distances (MSE) averaging from 512 initial noises. ](figs/Figure_consistency_motivation-02.png){#suppfig:motivation_schematics width="0.99\\linewidth"}
{#suppfig:motivation_schematics_CIFAR10 width="0.99\\linewidth"}
{#suppfig:motivation_schematics_CIFAR100 width="0.99\\linewidth"}
{#suppfig:motivation_schematics_LSUNbedroom32 width="0.99\\linewidth"}
{#suppfig:motivation_schematics_LSUNchurch32 width="0.99\\linewidth"}
{#suppfig:motivation_schematics_LSUNbedroom64 width="0.99\\linewidth"}
{#suppfig:motivation_schematics_LSUNchurch64 width="0.99\\linewidth"}
## Counterfactual moment manipulation: when the first two moments do not match {#app:counterfactual_mean_cov_manip}
{#suppfig:counterfactual_visual width="\\linewidth"}
{#suppfig:counterfactual_quantitative width="\\linewidth"}
## Additional validation with linear denoisers
{#suppfig:inhomogeneity_theory_vs_linear_denoiser width="0.65\\linewidth"}
{#suppfig:diff_sample_RMT_expectation_variance width="0.8\\linewidth"}
## Additional validation with deep networks
### Nearest-neighbor distances across dataset sizes {#app:nearest_neighbor_across_dataset_size}
{#suppfig:DNN_val_NN_AFHQ width="0.8\\linewidth"}
{#suppfig:DNN_val_NN_FFHQ32 width="0.8\\linewidth"}
DNN validation experiments (CIFAR10 and CIFAR100), nearest neighbor in training and control set
DNN validation experiments (LSUN bedroom 32 and 64), nearest neighbor in training and control set
DNN validation experiments (LSUN church 32 and 64), nearest neighbor in training and control set
### Cross-split consistency across dataset sizes {#app:consistency_scaling_dataset_size}
{#suppfig:DNN_val_scaling width="0.95\\linewidth"}
### DNN samples approach the linear predictor at larger dataset sizes {#app:DNN_approach_Gaussian}
{#suppfig:DNN_approach_Gaussian width="0.95\\linewidth"}
### Spectral structure of generated variance and cross-split disagreement {#app:spectral_anisotropy_var_consistency}
In this section, we provide further validation of the theoretically predicted spectral structure on additional datasets. We examine the variance of the generated samples along each population PC (blue traces), which ideally should match the population covariance exactly. This reveals effects similar to those in Fig. `\ref{fig:DNN_validation_results}`{=latex}**D**: insufficient variance and over-shrinkage along mid-range eigenvalues when the dataset size lies between the memorization and full-generalization regimes. We also examine the anisotropy of the consistency effect (red traces), which mirrors the behavior in Fig. `\ref{fig:DNN_validation_results}`{=latex}**E**: the top eigenspaces exhibit larger cross-split deviations. However, as the dataset size increases, the top eigenspaces become consistent faster than the lower ones.
{#suppfig:DNN_val_var_MSE_AFHQ32 width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_FFHQ32 width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_CIFAR width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_FFHQ64 width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_CIFAR100 width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_LSUNbedroom32 width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_LSUNbedroom64 width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_LSUNchurch32 width="0.95\\linewidth"}
{#suppfig:DNN_val_var_MSE_LSUNchurch64 width="0.95\\linewidth"}
### Spatial inhomogeneity of cross-split disagreement {#app:spatial_inhomogeneity_consistency_rmt}
Here, we present the per-seed cross split disagreement of DNN samples with corresponding RMT prediction, with similar format as Fig. `\ref{fig:DNN_validation_results}`{=latex}**F**. In general, we observe that, the RMT theory exhibits stronger and significant prediction power at larger dataset size, non-memorizing regime; in small data, memorization regime the prediction power is weak or non-significant. This corroborates prior finding that the linear score approximation is valid only in the non-memorization regime [@li2024understanding]. This holds across datasets and architectures.
DNN validation experiments, RMT predicting spatial inhomogeneity across seeds (FFHQ32)
DNN validation experiments, RMT predicting spatial inhomogeneity across seeds (FFHQ64)
{#suppfig:DNN_val_inhomogeneity_CIFAR100 width="0.99\\linewidth"}
{#suppfig:DNN_val_inhomogeneity_LSUNbedroom32 width="0.99\\linewidth"}
{#suppfig:DNN_val_inhomogeneity_LSUNbedroom64 width="0.99\\linewidth"}
{#suppfig:DNN_val_inhomogeneity_LSUNchurch32 width="0.99\\linewidth"}
{#suppfig:DNN_val_inhomogeneity_LSUNchurch64 width="0.99\\linewidth"}
### Anisotropy of the initial noise space {#app:anisotropy_init_noise}
In this section, we provide additional evidence supporting the point that the initial noise space already has well defined axes structure, i.e. the covariance eigenbasis, when using a deterministic sampler. We show that by amplifying or shrinking the random noise within certain eigenspaces with top or bottom eigenvalues. One see that, manipulating noise projection in top eigenspace has semantically interesting effect (e.g. removing backgrounds, Fig. `\ref{suppfig:anisotropic_structure_initial_noise_top50PC}`{=latex}), which can inspire future training-free adaptations of diffusion models.
{#suppfig:anisotropic_structure_initial_noise_top50PC width="0.85\\linewidth"}
{#suppfig:anisotropic_structure_initial_noise_top200PC width="0.85\\linewidth"}
{#suppfig:anisotropic_structure_initial_noise_bottom200PC width="0.85\\linewidth"}
### Consistency and proximity to linear predictor across training time {#app:longer_training_runs}
We next examine how the consistency effect evolves over longer training times. Here, we train FFHQ64 UNet models for 250k steps instead of 50k steps. Across training, the samples generated by the DNN first move closer to the samples predicted by linear theory, but later gradually drift away from them (Fig. `\ref{suppfig:longtrain_mse2linpred_s1model_own_vs_heldout_panels_250k}`{=latex}). Notably, the MSE distances to the linear predictors constructed from the training split and from the held-out split are nearly identical, highlighting the robustness of the underlying Gaussian statistics. This agrees with the previous observation that, during training, the learned score and denoiser first approach their linear counterparts before gradually deviating from them (@wang_unreasonable_2024, Fig. 12). In the memorization regime, the DNN samples begin to deviate from the linear predictor earlier and more strongly, consistent with smaller dataset sizes leading to earlier memorization. By contrast, for larger dataset sizes, the generated samples remain relatively close to the linear predictor throughout $10^5$ training steps, suggesting that the linear approximation remains valid for a long training window in the non-memorizing regime.
Examining cross-split consistency over training time (Fig. `\ref{suppfig:longtrain_consistency_panels_250k}`{=latex}), we find a similar trend to the MSE-to-linear-predictor curves: when models drift away from the linear predictor, they also become less consistent with models trained on the other split. This supports the interpretation that cross-split consistency arises from learning shared Gaussian statistics.
Finally, we examine FID during training (Fig. `\ref{suppfig:longtrain_FID_s1model_vs_s1_s2_full_panels_250k}`{=latex}), which measures the distributional distance between generated samples and the dataset in a deep feature space. In general, FID improves over training, i.e. decreases. In the memorization regime, however, FID decreases more strongly when evaluated against the training split than against the held-out split. In some cases, it also exhibits a \`\`double-descent"-like behavior: FID initially worsens when memorization begins, but later improves again as the model memorizes the training samples more accurately.
{#suppfig:longtrain_mse2linpred_s1model_own_vs_heldout_panels_250k width="0.99\\linewidth"}
{#suppfig:longtrain_consistency_panels_250k width="0.99\\linewidth"}
{#suppfig:longtrain_FID_s1model_vs_s1_s2_full_panels_250k width="0.99\\linewidth"}
# Proofs and Derivations
## Deterministic equivalence relations {#app:deterministic_equivalence}
Here we collect the one-point and two-point deterministic equivalent relationships from [@atanasov2024scalingRenormalization; @Atanasov2025TwoPoint; @bach2024highDimRMT_doubleDescent], rewritten in a unified notation.
#### Setup
Using similar notation as @bach2024highDimRMT_doubleDescent, we consider data matrix $X\in\mathbb{R}^{n\times d}$, where each row is an i.i.d. sample $\mathbf{x}_i$. The population covariance of these samples is denoted as $\boldsymbol{\Sigma}$. The key object of analysis is their empirical covariance $$\hat{\boldsymbol{\Sigma}}=\frac{1}{n}X^{\top}X$$, which we will assume lies in a Wishart-like universality class [@bach2024highDimRMT_doubleDescent; @silverstein1995strongConvergence]. We use $\gamma=d/n$ to denote the aspect ratio of $X$, i.e. the dimensionality of data over the number of data points. We use $\mathrm{Tr}$ to denote the normal trace of a matrix, and $\operatorname{tr}$ to denote the normalized trace, i.e. $\operatorname{tr}[M]=\frac 1d\mathrm{Tr}[M]$ for a $d\times d$ matrix $M$.
#### Self-consistency equation for the renormalized noise scale
The spectral properties of a matrix are determined by the Stieltjes transform. We consider the Stieltjes transform of the kernel matrix $\frac 1n XX^{\top}$, defined as $\hat \varphi(z):=\operatorname{Tr}[(XX^{\top}-nzI)^{-1}]$. In the large-matrix limit, the limiting variable $\varphi(z)$ satisfies the following self-consistent equation, $$\frac{1}{\varphi(z)}+z=\gamma\int_{0}^{\infty}\frac{sd\nu(s)}{1+s\varphi(z)}$$ where $\nu(s)$ is the limiting spectral measure of the population covariance $\boldsymbol{\Sigma}$. This follows from the leave-one-out arguments in the Appendix of @bach2024highDimRMT_doubleDescent, as well as @bai2010spectralbook [@ledoit2011eigenvectors].
This can be translated to the self-consistent equation of the renormalized ridge variable $\kappa(z):=\frac{1}{\varphi(-z)}$, which is used throughout the paper, $$\begin{aligned}
\frac{1}{\varphi(-z)}-z & =\gamma\int_{0}^{\infty}\frac{sd\nu(s)}{1+s\varphi(-z)}\\
\kappa(z)-z & =\gamma\int_{0}^{\infty}\frac{sd\nu(s)}{1+s\frac{1}{\kappa(z)}}\\
\kappa(z)-z & =\gamma\kappa(z)\int_{0}^{\infty}\frac{sd\nu(s)}{\kappa(z)+s}\\
z & =\kappa(z)\Big[1-\gamma\int_{0}^{\infty}\frac{sd\nu(s)}{\kappa(z)+s}\Big]
\end{aligned}$$
Practically, when solving such equations, given a finite size population covariance matrix, the integral over the spectral measure can be represented as normalized trace, leading to the Silverstein equation (Eq. `\ref{eq:silverstein_self_consistent}`{=latex}). $$\kappa(\lambda)-\lambda=\gamma\kappa(\lambda)\operatorname{tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\lambda)I)^{-1}]\tag{Silverstein}$$ The self-consistent solution to this equation obeys $\kappa(\lambda) \geq \lambda$ for all $\lambda \geq 0$, as can be seen from the fact that the right-hand-side is non-negative.
#### Degrees-of-freedom functions
We define the degrees-of-freedom functions with unnormalized trace, following the convention in @bach2024highDimRMT_doubleDescent, unlike @Atanasov2025TwoPoint, which used normalized trace. $$\begin{aligned}
\mathrm{df}_{1}(\lambda) & :=\operatorname{Tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\lambda I)^{-1}]\\
\mathrm{df}_{2}(\lambda) & :=\operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda I)^{-2}].
\end{aligned}$$ Note that $$\begin{aligned}
\mathrm{df}_{2}(\lambda)-\mathrm{df}_{1}(\lambda) & =\operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda I)^{-2}]-\operatorname{Tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\lambda I)^{-1}]\\
& =\operatorname{Tr}[\big(\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\lambda I)^{-1}-I\big)\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\lambda I)^{-1}]\\
& = - \lambda\operatorname{Tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\lambda I)^{-2}]\\
& \leq 0.
\end{aligned}$$ Note that both $\mathrm{df}_{2}(\lambda),\mathrm{df}_{1}(\lambda)$ are smaller than the number of non-zero eigenvalues of $\boldsymbol{\Sigma}$, i.e. $\mathrm{rank}(\boldsymbol{\Sigma})$. Thus, we have the chain of inequalities $$0 \leq \mathrm{df}_{2}(\lambda) \leq \mathrm{df}_{1}(\lambda) \leq \operatorname{rank}(\boldsymbol{\Sigma}) \leq d.$$ We will also use the mixed degrees-of-freedom function $$\begin{aligned}
\operatorname{df}_{2}(\lambda,\lambda') := \operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda I)^{-1}(\boldsymbol{\Sigma}+\lambda^\prime I)^{-1}] .
\end{aligned}$$ Clearly, $\operatorname{df}_{2}(\lambda,\lambda') \geq 0$. By the Cauchy-Schwarz inequality, $$\begin{aligned}
\operatorname{df}_{2}(\lambda,\lambda') = \operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda I)^{-1}(\boldsymbol{\Sigma}+\lambda^\prime I)^{-1}] \leq \sqrt{\operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda I)^{-2}] \operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\lambda^\prime I)^{-2}]} = \sqrt{\operatorname{df}_{2}(\lambda) \operatorname{df}_{2}(\lambda') }.
\end{aligned}$$
Now, by re-writing the Silverstein equation, we have that $$\begin{aligned}
\operatorname{df}_{1}(\kappa) = \left(1 - \frac{\lambda}{\kappa} \right) n < n,
\end{aligned}$$ where the inequality is strict if both $\lambda$ and $\kappa$ are strictly positive. This leads to an additional upper bound $$\begin{aligned}
\operatorname{df}_{2}(\kappa,\kappa') < n.
\end{aligned}$$
#### Basic deterministic equivalences
Following Proposition 1 of @bach2024highDimRMT_doubleDescent, we use the shorthand $\kappa(z):=1/\varphi(-z)$ to express the deterministic equivalences in the more convenient forms below. In what follows, $A$ and $B$ are test matrices of bounded spectral norm. For the resolvent, we have: $$\label{eq:DE_resolv_lin_one_z}
\operatorname{Tr}\bigl[A\,(\hat{\boldsymbol{\Sigma}}+\lambda I)^{-1}\bigr]\asymp\frac{\kappa(\lambda )}{\lambda }\operatorname{Tr}\bigl[A\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-1}\bigr]$$ and $$\begin{aligned}
\label{eq:DE_resolv_quad_one_z}
\begin{split}
&\operatorname{Tr}\bigl[A(\hat{\boldsymbol{\Sigma}}+\lambda I)^{-1}B(\hat{\boldsymbol{\Sigma}}+\lambda I)^{-1}\bigr]
\\ & \quad\asymp\frac{\kappa(\lambda )^{2}}{\lambda ^{2}\,}\operatorname{Tr}\bigl[A\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-1}\,B\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-1}\bigr]\\
& \qquad+\frac{\kappa(\lambda )^{2}}{\lambda ^{2}}\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(\lambda )\bigr)}\operatorname{Tr}\bigl[A\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-2}\boldsymbol{\Sigma}\bigr]\operatorname{Tr}\bigl[B\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-2}\boldsymbol{\Sigma}\bigr].
\end{split}
\end{aligned}$$ Equivalently, $$\label{eq:DE_mod_resolv_lin_one_z}
\operatorname{Tr}\bigl[A\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+\lambda I)^{-1}\bigr]\asymp\operatorname{Tr}\bigl[A\,\boldsymbol{\Sigma}\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-1}\bigr]$$ and $$\begin{aligned}
\label{eq:DE_mod_resolv_quad_one_z}
\begin{split}
&\operatorname{Tr}\bigl[A\,\hat{\boldsymbol{\Sigma}}\,(\hat{\boldsymbol{\Sigma}}+\lambda I)^{-1}\,B\,\hat{\boldsymbol{\Sigma}}\,(\hat{\boldsymbol{\Sigma}}+\lambda I)^{-1}\bigr] \\&\quad \asymp\operatorname{Tr}\bigl[A\,\boldsymbol{\Sigma}\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-1}\,B\,\boldsymbol{\Sigma}\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-1}\bigr]\\
& \qquad+\frac{\kappa^{2}(\lambda )}{n-\mathrm{df}_{2}\bigl(\kappa(\lambda )\bigr)}\operatorname{Tr}\bigl[A\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-2}\boldsymbol{\Sigma}\bigr]\operatorname{Tr}\bigl[B\,\bigl(\boldsymbol{\Sigma}+\kappa(\lambda )I\bigr)^{-2}\boldsymbol{\Sigma}\bigr]
\end{split}
\end{aligned}$$ where $\kappa(\lambda)$ can be solved from the self-consistent equation (`\ref{eq:silverstein_self_consistent}`{=latex}) above. Note that given the unnormalized trace, the equivalence $\asymp$ shall be understood through convergence of ratio.
#### Two-point equivalence for resolvents with different arguments
This can be further generalized to equivalence with two \"Ridge\" variables $\lambda,\lambda^\prime$, $$\begin{aligned}
\label{eq:DE_mod_resolv_quad_two_z}
\operatorname{Tr}\Big[A\hat{\boldsymbol{\Sigma}}(\lambda+\hat{\boldsymbol{\Sigma}})^{-1}B\hat{\boldsymbol{\Sigma}}(\lambda^{\prime}+\hat{\boldsymbol{\Sigma}})^{-1}\Big]\asymp\
&\operatorname{Tr}\big[AT_{\boldsymbol{\Sigma}}BT_{\boldsymbol{\Sigma}}^{\prime}\big]+ \\
&\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\operatorname{Tr}\big[AG_{\boldsymbol{\Sigma}}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}^{\prime}\big]\operatorname{Tr}\big[G_{\boldsymbol{\Sigma}}^{\prime}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}B\big] \notag
\end{aligned}$$ where $T_{\boldsymbol{\Sigma}}:=\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa)^{-1},T_{\boldsymbol{\Sigma}}^{\prime}:=\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa^{\prime})^{-1}$, $G_{\boldsymbol{\Sigma}}:=(\boldsymbol{\Sigma}+\kappa)^{-1},G_{\boldsymbol{\Sigma}}^{\prime}:=(\boldsymbol{\Sigma}+\kappa^{\prime})^{-1}$. and $\mathrm{df}_{2}(\kappa,\kappa'):=\operatorname{Tr}[\boldsymbol{\Sigma}^{2}G_{\boldsymbol{\Sigma}}G_{\boldsymbol{\Sigma}}^{\prime}]$. When $\kappa=\kappa^\prime$ it recovers Eq. `\ref{eq:DE_mod_resolv_quad_one_z}`{=latex}.
As a brief note for derivation, this follows from the deterministic equivalence for free product of matrices $A*B$ presented in Appendix A of @Atanasov2025TwoPoint. Let $A=\boldsymbol{\Sigma}$ be population covariance, $B=\frac{1}{n}ZZ^{T}$ be whitened data covariance, then $A*B=\hat{\boldsymbol{\Sigma}}$. Thus, $$\hat{\boldsymbol{\Sigma}}(\lambda+\hat{\boldsymbol{\Sigma}})^{-1}M\hat{\boldsymbol{\Sigma}}(\lambda^{\prime}+\hat{\boldsymbol{\Sigma}})^{-1}\asymp T_{\boldsymbol{\Sigma}}MT_{\boldsymbol{\Sigma}}^{\prime}+\kappa\kappa'G_{\boldsymbol{\Sigma}}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}^{\prime}\frac{\operatorname{Tr}\left[G_{\boldsymbol{\Sigma}}^{\prime}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}M\right]}{n-\mathrm{df}_{2}(\kappa,\kappa')}$$ Note that $q$ in their notation corresponds to our $\gamma$ and that their definition of $\mathrm{df}$ used normalized trace. This two-point equivalence was derived in @Atanasov2025TwoPoint [@Atanasov2024Risk] using a diagrammatic moment-method argument; it could also be derived by extending the leave-one-out arguments used by @bach2024highDimRMT_doubleDescent to prove the deterministic equivalents with a single $\lambda$ listed above.
#### Assumptions for the results
The deterministic-equivalence calculations in the paper are stated for the standard general-covariance sample-covariance model $\mathbf{x}_i=\boldsymbol{\mu}+\boldsymbol{\Sigma}^{1/2}\mathbf{z}_i$, with $\mathbf{z}_i$ i.i.d. standardized Gaussian or in a universality class for which the one- and two-point deterministic equivalents hold. Probe vectors and input displacements are deterministic or conditionally independent of $\hat{\boldsymbol{\Sigma}}$. Expectation formulas use anisotropic one-point deterministic equivalents; variance formulas use the factorized part of the second-order deterministic equivalent. When $\mathbf{v}$ and the input displacement $\bar\mathbf{x}$ have non-negligible weighted alignment, additional exchange terms may appear. Moreover, we will assume that all deterministic equivalents hold uniformly enough in the argument of the resolvent (namely, $\lambda$) so that we can use them to obtain the limiting behavior of integrals.
## Derivation of the denoiser expectation equivalence (Result `\ref{prop:denoiser_exp_DE}`{=latex}) {#app:proof_denoiser_exp_de}
::: proposition
**Result 5** (Main result, deterministic equivalence of the expectation of score and denoiser). *The optimal linear score and denoiser using empirical covariance has the following deterministic equivalence. $$\begin{aligned}
\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\Big] & \asymp\mathbf{v}^{\top}\mathbf{x}+\mathbf{v}^{\top}\kappa(\sigma^{2})(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}(\boldsymbol{\mu}-\mathbf{x})\\
& =\mathbf{v}^{\top}\boldsymbol{\mu}+\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}(\mathbf{x}-\boldsymbol{\mu})\\
\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\Big] & \asymp\frac{\kappa(\sigma^{2})}{\sigma^{2}}\mathbf{v}^{\top}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}(\boldsymbol{\mu}-\mathbf{x})
\end{aligned}$$*
:::
```{=latex}
\begin{IEEEproof}
Per assumption, assume the sample mean $\hat{\boldsymbol{\mu}}=\boldsymbol{\mu}$, consider
only the effect of empirical covariance $\hat{\boldsymbol{\Sigma}}$,
\[
\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)=\mathbf{x}+\sigma^{2}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})
\]
Using the deterministic equivalence Eq. \ref{eq:DE_resolv_lin_one_z},\ref{eq:DE_mod_resolv_lin_one_z}, in the sense that the trace with any independent matrix converge in ratio at limit.
\begin{align*}
(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1} & \asymp\frac{\kappa(\sigma^{2})}{\sigma^{2}}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}\\
\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1} & \asymp\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}
\end{align*}
Then, given a fixed measurement vector $\mathbf{v}$, and a noised
input $\mathbf{x}$, the projection of score onto a vector can be
framed as trace. The equivalence reads,
\begin{align*}
\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\Big] & =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})\Big]\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\mathrm{Tr}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})\mathbf{v}^{\top}\Big]\\
& \asymp\frac{\kappa(\sigma^{2})}{\sigma^{2}}\mathrm{Tr}\Big[(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}(\boldsymbol{\mu}-\mathbf{x})\mathbf{v}^{\top}\Big]\\
& =\frac{\kappa(\sigma^{2})}{\sigma^{2}}\mathbf{v}^{\top}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}(\boldsymbol{\mu}-\mathbf{x})
\end{align*}
Similarly, use the other equivalence, the denoiser projection has
equivalence,
\begin{align*}
\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\Big] & =\mathbf{v}^{\top}\boldsymbol{\mu}+\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\mathbf{x}-\boldsymbol{\mu})\Big]\\
& \asymp\mathbf{v}^{\top}\boldsymbol{\mu}+\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}(\mathbf{x}-\boldsymbol{\mu})\\
& =\mathbf{v}^{\top}\mathbf{D}_{\boldsymbol{\Sigma}}^{*}(\mathbf{x};\kappa^{1/2})
\end{align*}
Thus, in the expectation sense, the effect of empirical data covariance
(finite data) on the denoiser, is equivalent to renormalizing and
increasing the effective noise scale $\sigma^{2}\to\kappa(\sigma^{2})$,
similar to adding an adaptive Ridge regularization.
\end{IEEEproof}
```
#### Interpretation
Measuring the deviation of the empirical covariance denoiser from the population covariance denoiser, at the same noise scale, $$\begin{aligned}
& \mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\big(\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)-\mathbf{D}_{\boldsymbol{\Sigma}}^{*}(\mathbf{x};\sigma)\Big)\Big]\\
\asymp & \mathbf{v}^{\top}\Big[\kappa(\sigma^{2}) (\boldsymbol{\Sigma}+\kappa(\sigma^{2})I)^{-1}-\sigma^{2}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1}\Big](\boldsymbol{\mu}-\mathbf{x})
\end{aligned}$$ Using push through identity $A^{-1}-B^{-1}=A^{-1}(B-A)B^{-1}$, $$\begin{aligned}
& \kappa(\boldsymbol{\Sigma}+\kappa I)^{-1}-\sigma^{2}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1}\\
= & \kappa\sigma^{2}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1}\Big(\frac{1}{\sigma^{2}}(\boldsymbol{\Sigma}+\sigma^{2}I)-\frac{1}{\kappa}(\boldsymbol{\Sigma}+\kappa I)\Big)\\
= & (\kappa(\sigma^{2})-\sigma^{2})\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1}
\end{aligned}$$ We can represent the deviation as a resolvent product. This makes it clear that the deviation is proportional to the effect of renormalization $(\kappa(\sigma^{2})-\sigma^{2})$. $$\begin{aligned}
& \mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\big(\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)-\mathbf{D}_{\boldsymbol{\Sigma}}^{*}(\mathbf{x};\sigma)\Big)\Big]\\
\asymp & (\kappa(\sigma^{2})-\sigma^{2})\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})
\end{aligned}$$ Setting the measurement vector along population eigenvector $\mathbf{u}_k$, with eigenvalue $\lambda_{k}$, then the deviation reads $$\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{u}_k^{\top}\big(\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)-\mathbf{D}_{\boldsymbol{\Sigma}}^{*}(\mathbf{x};\sigma)\Big)\Big] \asymp \frac{\lambda_{k}(\kappa-\sigma^{2})}{(\lambda_{k}+\sigma^{2})(\lambda_{k}+\kappa)}\mathbf{u}_k^{\top}(\boldsymbol{\mu}-\mathbf{x})$$ Viewed as a function of $\lambda_k$, the coefficient $$f(\lambda_k)=\frac{\lambda_k(\kappa-\sigma^2)}{(\lambda_k+\sigma^2)(\lambda_k+\kappa)}$$ is unimodal, peaking at the geometric mean of the two scales, $\lambda_k^\star=\sqrt{\sigma^2\kappa}$, where it attains $$f(\lambda_k^\star)=\frac{\kappa-\sigma^2}{\big(\sigma+\sqrt{\kappa}\big)^2}=\frac{\sqrt{\kappa}-\sigma}{\sqrt{\kappa}+\sigma}.$$ Since $\kappa(\sigma^2)>\sigma^2$ (renormalization inflates the noise scale), this peak sits somewhat above $\sigma^2$ but still in the low-to-moderate part of the spectrum, and the weight decays to zero on both sides:
- For top eigenmodes ($\lambda_k\gg\kappa$), $f(\lambda_k)\approx(\kappa-\sigma^2)/\lambda_k\to 0$: large-variance directions are essentially unaffected, because the resolvents $(\boldsymbol{\Sigma}+\kappa I)^{-1}$ and $(\boldsymbol{\Sigma}+\sigma^2 I)^{-1}$ both act like $\lambda_k^{-1}$ and the additive shifts $\kappa,\sigma^2$ are negligible relative to $\lambda_k$.
- For very small eigenmodes ($\lambda_k\to 0$), $f(\lambda_k)\to 0$ as well, suppressed by the leading $\lambda_k$ factor.
- In between, modes with $\lambda_k$ comparable to $\sqrt{\sigma^2\kappa}$---the lower, signal-to-noise-marginal part of the spectrum---receive the largest correction, of magnitude $(\sqrt{\kappa}-\sigma)/(\sqrt{\kappa}+\sigma)$.
Thus the deviation is the most salient when $\kappa$ deviates from $\sigma^2$; and along the spectrum, it peaked at the eigenmodes with variance around $\sigma^2$ and $\kappa$.
## Derivation of the denoiser fluctuation equivalence (Result `\ref{prop:denoiser_var_DE}`{=latex}) {#app:proof_denoiser_var_de}
::: proposition
**Result 6** (Main result, deterministic equivalence of the denoiser variance). *Assuming $\hat{\boldsymbol{\mu}}=\boldsymbol{\mu}$, across dataset realizations of size $n$, the variance of the optimal empirical linear denoiser at point $\mathbf{x}$ in direction $\mathbf{v}$, given by $\mathbf{v}^\top \mathcal{S}_{D}(\mathbf{x}) \mathbf{v}$, admits the following deterministic equivalence. $$\begin{aligned}
\label{eq:denoiser_var_deter_equiv_full}
\mathbf{v}^{\top}\mathcal{S}_{D}(\mathbf{x})\mathbf{v} &=\operatorname{Var}_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\mathbf{D}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)] \\
&\asymp {\frac{\kappa(\sigma^{2})^{2}}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}}%_{\hexagon(n,\sigma^{2},\kappa)}\
\underbrace{\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{v}\Big)}_{\text{anisotropy: }\Diamond(\mathbf{v},\kappa,\boldsymbol{\Sigma})}\
\underbrace{\Big((\mathbf{x}-\boldsymbol{\mu})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}(\mathbf{x}-\boldsymbol{\mu})\Big)}_{\text{inhomogeneity: }\Diamond(\mathbf{x}-\boldsymbol{\mu},\kappa,\boldsymbol{\Sigma})} \notag
\end{aligned}$$*
:::
```{=latex}
\begin{IEEEproof}
Next, we examine the covariance of denoiser due to dataset realization,
the score variance reads,
\begin{align*}
\mathcal{S}_{s}:=Cov_{\hat{\boldsymbol{\Sigma}}}[\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)] & =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)^{\top}-\Big(\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\Big)\Big(\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)\Big)^{\top}\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big]-\\
&\qquad\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})\Big]\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[(\boldsymbol{\mu}-\mathbf{x})^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big]\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big]-\\
&\qquad \mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big](\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big]
\end{align*}
Note that the variance of denoiser and that of score has the simple
scaling relationship, so we just need to study the score.
\[
\mathcal{S}_{\mathbf{D}}=\sigma^{4}\mathcal{S}_{s}
\]
We are interested in the variance of score vector along a fixed probe
vector $\mathbf{v}$,
\begin{align*}
\mathbf{v}^{\top}\mathcal{S}_{s}\mathbf{v} & =Var_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)]\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\mathbf{v}\Big]-\Big(\mathbf{v}^{\top}\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big](\boldsymbol{\mu}-\mathbf{x})\Big)^{2}\\
& =\underbrace{\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\operatorname{Tr}\Big[\mathbf{v}\mathbf{v}^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big]}_{\text{2nd moment}}-\Big(\underbrace{\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\operatorname{Tr}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})\mathbf{v}^{\top}\Big]}_{\text{1st moment}}\Big)^{2}
\end{align*}
The two terms can be tackled by one-point and two-point equivalence
Eq. \ref{eq:DE_resolv_quad_one_z},\ref{eq:DE_resolv_lin_one_z}.
Abbreviating $A:=\mathbf{v}\mathbf{v}^{\top}$,
$B:=(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}$, $z:=\sigma^{2}$.
\[
\operatorname{Tr}\bigl[A\,(\hat{\boldsymbol{\Sigma}}+zI)^{-1}\bigr]\sim\frac{\kappa(z)}{z}\operatorname{Tr}\bigl[A\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}\bigr]
\]
\begin{align*}
\operatorname{Tr}\bigl[A(\hat{\boldsymbol{\Sigma}}+zI)^{-1}B(\hat{\boldsymbol{\Sigma}}+zI)^{-1}\bigr] & \sim\frac{\kappa(z)^{2}}{z^{2}\,}\operatorname{Tr}\bigl[A\,\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}\,B\,\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}\bigr]\\
& +\frac{\kappa(z)^{2}}{z^{2}}\operatorname{Tr}\bigl[A\,\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}\bigr]\operatorname{Tr}\bigl[B\,\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}\bigr]\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(z)\bigr)}
\end{align*}
The 2nd moment term is equivalent to,
\begin{align*}
&\operatorname{Tr}\bigl[A(\hat{\boldsymbol{\Sigma}}+zI)^{-1}B(\hat{\boldsymbol{\Sigma}}+zI)^{-1}\bigr] \\ &\nonumber \sim\frac{\kappa(z)^{2}}{z^{2}\,}\operatorname{Tr}\bigl[\mathbf{v}\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}\bigr]\\
& \quad +\frac{\kappa(z)^{2}}{z^{2}}\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(z)\bigr)}\operatorname{Tr}\bigl[\mathbf{v}\mathbf{v}^{\top}\,\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}\bigr]\operatorname{Tr}\bigl[(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}\,\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}\bigr]\\
& =\frac{\kappa(z)^{2}}{z^{2}\,}\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}(\boldsymbol{\mu}-\mathbf{x})\Big)^{2}\\
&\quad +\frac{\kappa(z)^{2}}{z^{2}}\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(z)\bigr)}\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{v}\Big)\ \Big((\boldsymbol{\mu}-\mathbf{x})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}(\boldsymbol{\mu}-\mathbf{x})\Big)
\end{align*}
The first moment term is equivalent to,
\begin{align*}
\operatorname{Tr}\Big[(\hat{\boldsymbol{\Sigma}}+zI)^{-1}(\boldsymbol{\mu}-\mathbf{x})\mathbf{v}^{\top}\Big] & \sim\frac{\kappa(z)}{z}\operatorname{Tr}\bigl[\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}(\boldsymbol{\mu}-\mathbf{x})\mathbf{v}^{\top}\bigr]\\
& =\frac{\kappa(z)}{z}\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}(\boldsymbol{\mu}-\mathbf{x})
\end{align*}
Thus, combining the two terms, we obtain the variance of score at
noised datapoint $\mathbf{x}$, along direction $\mathbf{v}$,
\begin{align*}
\mathbf{v}^{\top}\mathcal{S}_{s}(\mathbf{x})\mathbf{v} & =Var_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\mathbf{s}_{\hat{\boldsymbol{\Sigma}}}^{*}(\mathbf{x};\sigma)]\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\operatorname{Tr}\Big[\mathbf{v}\mathbf{v}^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})(\boldsymbol{\mu}-\mathbf{x})^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}\Big] \nonumber\\&\quad -\Big(\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\operatorname{Tr}\Big[(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1}(\boldsymbol{\mu}-\mathbf{x})\mathbf{v}^{\top}\Big]\Big)^{2}\\
& \sim\frac{\kappa(z)^{2}}{z^{2}\,}\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}(\boldsymbol{\mu}-\mathbf{x})\Big)^{2}\\
& +\frac{\kappa(z)^{2}}{z^{2}}\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{v}\Big)\ \Big((\boldsymbol{\mu}-\mathbf{x})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}(\boldsymbol{\mu}-\mathbf{x})\Big)\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(z)\bigr)}\\
& -\Big(\frac{\kappa(z)}{z}\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-1}(\boldsymbol{\mu}-\mathbf{x})\Big)^{2}\\
& =\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(z)\bigr)}\frac{\kappa(z)^{2}}{z^{2}}\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{v}\Big)\ \Big((\boldsymbol{\mu}-\mathbf{x})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(z)I\bigr)^{-2}\boldsymbol{\Sigma}(\boldsymbol{\mu}-\mathbf{x})\Big)\\
(z\mapsto\sigma^{2}) & =\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}\frac{\kappa(\sigma^{2})^{2}}{\sigma^{4}}\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{v}\Big)\ \Big((\boldsymbol{\mu}-\mathbf{x})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}(\boldsymbol{\mu}-\mathbf{x})\Big)
\end{align*}
Per simple scaling, the variance of denoisers reads,
\begin{align*}
\mathbf{v}^{\top}\mathcal{S}_{D}(\mathbf{x})\mathbf{v} & =\sigma^{4}\mathbf{v}^{\top}\mathcal{S}_{s}(\mathbf{x})\mathbf{v}\\
& \sim\frac{\kappa(\sigma^{2})^{2}}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}\ \underbrace{\Big(\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{v}\Big)}_{\Diamond(\mathbf{v},\kappa,\boldsymbol{\Sigma})}\!\underbrace{\Big((\boldsymbol{\mu}-\mathbf{x})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}(\boldsymbol{\mu}-\mathbf{x})\Big)}_{\Diamond(\boldsymbol{\mu}-\mathbf{x},\kappa,\boldsymbol{\Sigma})}
\end{align*}
\end{IEEEproof}
```
::: remark
*Remark 7*. `\label{rmk:generic_vector}`{=latex} The equivalence above keeps the factorized, generic contribution $(\mathbf{v}^\top \mathbf{K} \mathbf{v})(\mathbf{y}^\top \mathbf{K} \mathbf{y})$ to the variance, where $\mathbf{K}=\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-2}$ and $\mathbf{y}=\mathbf{x}-\boldsymbol{\mu}$. However, a fully anisotropic real-Wishart second-order formula can also contain an exchange/alignment contribution proportional to $(\mathbf{v}^\top \mathbf{K} \mathbf{y})^2$. This term is negligible for independent isotropic or high-effective-rank generic inputs, and for the eigenband averaged or fully marginalized settings used in the figures, but can be comparable for deliberately aligned probes [@bach2024highDimRMT_doubleDescent; @Atanasov2024Risk; @Atanasov2025TwoPoint]. For $\mathbf{y}\sim N(0,\tau^2 I)$ independent of $v$, the ratio of the expected exchange contribution to the factorized contribution is $$\begin{aligned}
\frac{\mathbb{E}_{\mathbf{y}} (\mathbf{v}^\top \mathbf{K}\mathbf{y})^2}{\mathbb{E}_{\mathbf{y}}[(\mathbf{v}^\top \mathbf{K} \mathbf{v})(\mathbf{y}^{\top}\mathbf{K}\mathbf{y})]}
&= \frac{\mathbf{v}^\top \mathbf{K}^2 \mathbf{v}}{(\mathbf{v}^\top \mathbf{K} \mathbf{v})\operatorname{Tr}(\mathbf{K})}
\le
\frac{\|\mathbf{K}\|_{\rm op}}{\operatorname{Tr}\mathbf{K}},
\end{aligned}$$ so the exchange term is suppressed when the effective rank of $\mathbf{K}$ is large.
:::
### Interpretation and derivations
#### Dependence on the probe direction $\mathbf{v}$
This dependency on $\mathbf{v}$ describes the *anisotropy of uncertainty*, or the variance of the score/denoiser prediction along different directions. $$\begin{aligned}
\Diamond(\mathbf{v},\kappa,\boldsymbol{\Sigma}) & :=\mathbf{v}^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{v}\\
& =\mathbf{v}^{\top}U\frac{\Lambda}{(\Lambda+\kappa(\sigma^{2}))^{2}}U^{\top}\mathbf{v}
\end{aligned}$$ By assumption, the probe vector $\mathbf{v}$ is a unit vector. Then this dependency is determined by the diagonal matrix $\frac{\Lambda}{(\Lambda+\kappa(\sigma^{2}))^{2}}=\text{diag}\left(\frac{\lambda_{k}}{(\lambda_{k}+\kappa(\sigma^{2}))^{2}}\right)$ and the alignment between $\mathbf{v}$ and the eigenbasis $U$.
Consider when the probing vector is aligned exactly with the $k$-th eigenvector $\mathbf{u}_k$, this term reads $$\begin{aligned}
\Diamond(\mathbf{u}_k,\kappa,\boldsymbol{\Sigma}) & =\mathbf{u}_k^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{u}_k\\
& =\frac{\lambda_{k}}{(\lambda_{k}+\kappa(\sigma^{2}))^{2}}\\
& =:\chi(\lambda_{k},\kappa(\sigma^{2}))
\end{aligned}$$
We can discuss the different regimes of $\chi(\lambda_{k},\kappa)$ depending on $\lambda_{k}$ and $\kappa(\sigma^{2})$
- High noise regime $\lambda_{k}\ll\kappa$: $\chi(\lambda_{k},\kappa)\approx\frac{\lambda_{k}}{\kappa^{2}}$, so $\chi(\lambda_{k},\kappa)$ increases with $\lambda_{k}$. Higher variance directions have larger uncertainties.
- Low noise regime $\lambda_{k}\gg\kappa$: $\chi(\lambda_{k},\kappa)\approx\frac{1}{\lambda_{k}}$, so $\chi(\lambda_{k},\kappa)$ will decrease with $\lambda_{k}$. Lower variance directions have larger uncertainties.
- Regarding $\kappa$, $\chi(\lambda_{k},\kappa)$ decreases monotonically with $\kappa$, i.e. the higher the noise scale, the smaller the variance.
- Regarding $\lambda_{k}$, $\chi(\lambda_{k},\kappa)$ has one unique maximum, where $\arg\max_{\lambda}\chi(\lambda,\kappa)=\kappa$, and $\max_{\lambda}\chi(\lambda,\kappa)=\frac{1}{4\kappa}$. Thus, it is a bell shaped function of $\lambda_{k}$. (Proof below.)
- This shows that at different noise level or $\kappa(\sigma^{2})$, there is always some direction with variance comparable to $\kappa(\sigma^{2})$ which will have the largest variance!
- Further, the largest variance will be inversely proportional to $\kappa(\sigma^{2})$, i.e. generally larger at lower noise.
This result shows that the anisotropy of the uncertainty depends on the renormalized noise scale $\kappa(\sigma^{2})$, and the maximal uncertainty is focused around the PC dimensions with variance similar to $\kappa(\sigma^{2})$.
`\noindent`{=latex}`\fbox{%\begin{minipage}[t]{1\columnwidth \fboxsep \fboxrule}%
\begin{minipage}[t]{\columnwidth-2\fboxsep-2\fboxrule}
\textbf{Proof of unique maximum of $\chi(\lambda,\kappa)$}
\begin{proof}
Given
\[
\chi(\lambda,\kappa)=\frac{\lambda}{(\lambda+\kappa)^{2}}
\]
Then
\begin{align*}
\frac{d\chi(\lambda,\kappa)}{d\lambda} & =\frac{(\lambda+\kappa)^{2}-2(\lambda+\kappa)\lambda}{(\lambda+\kappa)^{4}}\\
& =\frac{\kappa-\lambda}{(\lambda+\kappa)^{3}}
\end{align*}
Setting the gradient to zero yields the unique stationary point, $\kappa=\lambda$.
Given $\kappa,\lambda>0$, we have the unique maximum w.r.t. $\lambda$.
\begin{align*}
\arg\max_{\lambda}\chi(\lambda,\kappa) & =\kappa\\
\max_{\lambda}\chi(\lambda,\kappa) & =\frac{1}{4\kappa}
\end{align*}
\end{proof}
\end{minipage}}`{=latex}
#### Dependence on the probe point $\mathbf{x}$.
The dependency on probe point $\mathbf{x}$ tells us about the spatial inhomogeneity of the uncertainty. $$\begin{aligned}
\Diamond(\mathbf{x}-\boldsymbol{\mu},\kappa,\boldsymbol{\Sigma}) & =(\mathbf{x}-\boldsymbol{\mu})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}(\mathbf{x}-\boldsymbol{\mu})\\
& =(\mathbf{x}-\boldsymbol{\mu})^{\top}U\frac{\Lambda}{(\Lambda+\kappa(\sigma^{2}))^{2}}U^{\top}(\mathbf{x}-\boldsymbol{\mu})\\
& =\sum_{k}\frac{\lambda_{k}}{(\lambda_{k}+\kappa(\sigma^{2}))^{2}}\Big(\mathbf{u}_k^{\top}(\mathbf{x}-\boldsymbol{\mu})\Big)^{2}\\
& =\sum_{k}\chi(\lambda_{k},\kappa(\sigma^{2}))\ \Big(\mathbf{u}_k^{\top}(\mathbf{x}-\boldsymbol{\mu})\Big)^{2}
\end{aligned}$$
This is similar to the dependency above, except that now the argument $\mathbf{x}-\boldsymbol{\mu}$ is no longer unit normed, but any probing direction in the sample space.
Note, generally the noised sample $\mathbf{x}$ from a certain realization of dataset is distributed like $\mathcal{N}(\boldsymbol{\mu},\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)$ (under Gaussian data assumption), so $$\mathbf{v}^{\top}(\mathbf{x}-\boldsymbol{\mu})\sim\mathcal{N}(0,\mathbf{v}^{\top}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)\mathbf{v})$$
Consider a probe point on the hyperelliptical shell defined by $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma}+\sigma^{2}I)$, then if the point falls on the line $\mathbf{x}=\boldsymbol{\mu}+c\mathbf{u}_k$. $\|\mathbf{x}-\boldsymbol{\mu}\|^{2}=c^{2}\approx\sigma^{2}+\lambda_{k}$
Then $$\begin{aligned}
\Diamond(\mathbf{x}-\boldsymbol{\mu},\kappa,\boldsymbol{\Sigma}) & =\Diamond(c\mathbf{u}_k,\kappa,\boldsymbol{\Sigma})\\
& =c^{2}\frac{\lambda_{k}}{(\lambda_{k}+\kappa(\sigma^{2}))^{2}}\\
& \approx\frac{(\lambda_{k}+\sigma^{2})\lambda_{k}}{(\lambda_{k}+\kappa(\sigma^{2}))^{2}}\\
& =(\sigma^{2}+\lambda_{k})\ \chi(\lambda_{k},\kappa(\sigma^{2}))\\
& =\xi(\lambda_{k},\sigma^{2})
\end{aligned}$$
- **High noise regime**, $\kappa>\sigma^{2}\gg\lambda$, then $\xi(\lambda,\sigma^{2})\approx\frac{\sigma^{2}\lambda}{\kappa^{2}(\sigma^{2})}<\frac{\lambda}{\kappa(\sigma^{2})}\ll1$, which scale linearly with PC variance $\lambda$, higher the PC, larger the variance.
- **Low noise regime**, $\lambda\gg\kappa>\sigma^{2}$, then $\xi(\lambda,\sigma^{2})\approx1$. Then all points on the ellipsoid have large variance.
- At any fixed $\sigma^{2}$, this function increases monotonically with $\lambda_{k}$.
- The score or denoiser variance is larger when the probing point $\boldsymbol{\mu}+c\mathbf{u}_k$ is deviating along those higher variance directions $\mathbf{u}_k$.
- When the probing point is deviating along low variance directions, the variance is lower.
`\noindent`{=latex}`\fbox{\begin{minipage}[t]{\columnwidth-2\fboxsep-2\fboxrule}
\paragraph{Derivation of properties of $\xi(\lambda,\sigma^{2})$
\[
\xi(\lambda,\sigma^{2})=\frac{(\sigma^{2}+\lambda)\lambda}{(\lambda+\kappa(\sigma^{2}))^{2}}
\]
}
Derivative
\begin{align*}
\frac{d\xi(\lambda,\sigma^{2})}{d\lambda} & =\frac{(\sigma^{2}+2\lambda)(\lambda+\kappa(\sigma^{2}))^{2}-2(\lambda+\kappa(\sigma^{2}))(\sigma^{2}+\lambda)\lambda}{(\lambda+\kappa(\sigma^{2}))^{4}}\\
& =\frac{(\sigma^{2}+2\lambda)(\lambda+\kappa(\sigma^{2}))-2(\sigma^{2}+\lambda)\lambda}{(\lambda+\kappa(\sigma^{2}))^{3}}\\
& =\frac{(\sigma^{2}+2\lambda)\kappa(\sigma^{2})-\lambda\sigma^{2}}{(\lambda+\kappa(\sigma^{2}))^{3}}\\
& =\frac{\sigma^{2}\kappa(\sigma^{2})+(2\kappa(\sigma^{2})-\sigma^{2})\lambda}{(\lambda+\kappa(\sigma^{2}))^{3}}
\end{align*}
Note that through the self-consistent equation $\kappa(\sigma^{2})-\sigma^{2}>0$,
thus $\frac{d\xi(\lambda,\sigma^{2})}{d\lambda}>0,\forall\lambda$.
The function is monotonically increasing for $\lambda$.
Given that $\kappa(\sigma^{2})>\sigma^{2}>0$, we have bounds
\[
\xi(\lambda,\sigma^{2})=\frac{(\sigma^{2}+\lambda)\lambda}{(\lambda+\kappa(\sigma^{2}))^{2}}<\frac{\lambda}{\lambda+\kappa(\sigma^{2})}<1
\]
\end{minipage}}`{=latex}
#### Overall scaling with sample size
Finally, we marginalize over space and direction, obtaining an overall quantification of consistency of denoiser, and study its scaling property.
First, *marginalizing (summing) all directions*, we have $$\begin{aligned}
\sum_{k}\Diamond(\mathbf{u}_k,\kappa,\boldsymbol{\Sigma}) & =\sum_{k}\mathbf{u}_k^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{u}_k\\
& =\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\Big]
\end{aligned}$$
This can be further abbreviated as following, $$\begin{aligned}
\sum_{k}\Diamond(\mathbf{u}_k,\kappa,\boldsymbol{\Sigma}) & =\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\!I\Big]\\
& =\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\!\frac{1}{\kappa(\sigma^{2})}(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I-\boldsymbol{\Sigma})\Big]\\
& =\frac{1}{\kappa(\sigma^{2})}\Big(\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-1}\boldsymbol{\Sigma}\Big]-\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}^{2}\Big]\Big)\\
& =\frac{\mathrm{df}_{1}(\kappa)-\mathrm{df}_{2}(\kappa)}{\kappa}
\end{aligned}$$
Next, *marginalize (averaging) over space*. Here we consider the noised distribution starting from the true target distribution $\mathbf{x}\sim p(\mathbf{x};\sigma)=p_{0}(\mathbf{x})*\mathcal{N}(0,\sigma^{2}I)$. For us, the only thing that matters is the second moment, so for an arbitrary distribution we have, $$\mathbb{E}_{\mathbf{x}}[(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{\top}]=\boldsymbol{\Sigma}+\sigma^{2}I$$
Thus, $$\begin{aligned}
\mathbb{E}_{\mathbf{x}}\Diamond(\mathbf{x}-\boldsymbol{\mu},\kappa,\boldsymbol{\Sigma}) & =\mathbb{E}_{\mathbf{x}}\Big[(\mathbf{x}-\boldsymbol{\mu})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}(\mathbf{x}-\boldsymbol{\mu})\Big]\\
& =\mathrm{Tr}\Big[(\boldsymbol{\Sigma}+\sigma^{2}I)\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\Big]
\end{aligned}$$
This can also be abbreviated using degree of freedom, $$\begin{aligned}
\mathbb{E}_{\mathbf{x}}\Diamond(\mathbf{x}-\boldsymbol{\mu},\kappa,\boldsymbol{\Sigma}) & =\mathrm{Tr}\Big[(\boldsymbol{\Sigma}+\sigma^{2}I)\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\Big]\\
& =\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}^{2}\Big]+\sigma^{2}\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\Big]\\
& =\mathrm{df}_{2}(\kappa)+\frac{\sigma^{2}}{\kappa}(\mathrm{df}_{1}(\kappa)-\mathrm{df}_{2}(\kappa))\\
& =\frac{\sigma^{2}}{\kappa}\mathrm{df}_{1}(\kappa)+(1-\frac{\sigma^{2}}{\kappa})\mathrm{df}_{2}(\kappa)
\end{aligned}$$
Thus, we have $$\begin{aligned}
\mathbb{E}_{\mathbf{x}}\sum_{k}\mathbf{u}_k^{\top}\mathcal{S}_{D}(\mathbf{x})\mathbf{u}_k & \asymp\frac{\kappa(\sigma^{2})^{2}}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}\ \sum_{k}\underbrace{\Big(\mathbf{u}_k^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\mathbf{u}_k\Big)}_{\Diamond(\mathbf{v},\kappa,\boldsymbol{\Sigma})}\!\mathbb{E}_{\mathbf{x}}\underbrace{\Big((\boldsymbol{\mu}-\mathbf{x})^{\top}\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}(\boldsymbol{\mu}-\mathbf{x})\Big)}_{\Diamond(\boldsymbol{\mu}-\mathbf{x},\kappa,\boldsymbol{\Sigma})}\\
& =\frac{\kappa(\sigma^{2})^{2}}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}\mathrm{Tr}\Big[\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\Big]\mathrm{Tr}\Big[(\boldsymbol{\Sigma}+\sigma^{2}I)\bigl(\boldsymbol{\Sigma}+\kappa(\sigma^{2})I\bigr)^{-2}\boldsymbol{\Sigma}\Big]\\
& =\frac{\kappa(\sigma^{2})^{2}}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}\times\frac{\mathrm{df}_{1}(\kappa)-\mathrm{df}_{2}(\kappa)}{\kappa}\times\Big(\frac{\sigma^{2}}{\kappa}\mathrm{df}_{1}(\kappa)+(1-\frac{\sigma^{2}}{\kappa})\mathrm{df}_{2}(\kappa)\Big)\\
& =\frac{\big(\mathrm{df}_{1}(\kappa)-\mathrm{df}_{2}(\kappa)\big)\times\big(\sigma^{2}\mathrm{df}_{1}(\kappa)+(\kappa-\sigma^{2})\mathrm{df}_{2}(\kappa)\big)}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}\\
& =:\Delta(n,\sigma^{2},\Lambda)
\end{aligned}$$
$$\label{eq:denoiser_var_overall_scaling}
\Delta(n,\sigma^2,\Lambda)
\;=\;
\frac{(\mathrm{df}_1(\kappa)-\mathrm{df}_2(\kappa))\;\big(\sigma^2\mathrm{df}_1(\kappa)+(\kappa-\sigma^2)\mathrm{df}_2(\kappa)\big)}{n-\mathrm{df}_2(\kappa)}$$
Now, marginalized over space and direction, this is only a function of the population spectrum, sample number and noise scale. Note $n$ is the sample number, so it makes sense when $n$ goes to infinity, then $\hat{\boldsymbol{\Sigma}}\to\boldsymbol{\Sigma}$ and $\kappa\to\sigma^{2}$, the variance reduces to zero.
Basically the higher the $\kappa$, the smaller the $\mathrm{df}_{2}(\kappa)$, so $n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)$ will be larger, which scales down $\frac{1}{n-\mathrm{df}_{2}\bigl(\kappa(\sigma^{2})\bigr)}\frac{\kappa(\sigma^{2})^{2}}{\sigma^{4}}$.
When we compare our theory with empirical measurements of denoiser or sample deviations between two non-overlapping splits, we use the following lemma to predict the expected MSE deviation from the variance.
::: lemma
**Lemma 8** (Expected MSE between two i.i.d. samples doubles the variance). *`\label{lemma:dev_double_var}`{=latex} Let $X,Y$ be i.i.d. random variables with variance $S=\operatorname{Var}(X)$. Then their mean squared error (MSE) is double the variance. $$\mathbb{E}\!\left[(X-Y)^2\right] \;=\; 2S.$$*
:::
::: proof
*Proof.* Expanding and using independence, $$\mathbb{E}\!\left[(X-Y)^2\right]
= \mathbb{E}[X^2] + \mathbb{E}[Y^2] - 2\,\mathbb{E}[XY]
= 2\,\mathbb{E}[X^2] - 2\,\mathbb{E}[X]\mathbb{E}[Y].$$ Since $\operatorname{Var}(X)=\mathbb{E}[X^2]-\big(\mathbb{E}[X]\big)^2=S$, this simplifies to $2S$. ◻
:::
## Integral representation of matrix fractional powers {#app:Balakrishnan_identity}
::: lemma
**Lemma 9** (Scalar beta integral identity [@Balakrishnan1960FractionalPower]). *`\label{lem:Scalar_beta_identity}`{=latex}We have the integral identity $$\int_{0}^{\infty}\frac{t^{-\alpha}dt}{\lambda+t}=\frac{\pi}{\sin(\pi\alpha)}\lambda^{-\alpha},\quad\alpha\in(0,1)$$*
:::
```{=latex}
\begin{IEEEproof}
Recall the definition of Beta function,
\begin{align*}
B(p,q) & =\int_{0}^{1}u^{p-1}(1-u)^{q-1}du =\frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)}
\end{align*}
We can turn it into beta function via change of variable $u=\frac{t}{\lambda+t}$,
then $t\in[0,\infty)$ maps to $u\in[0,1)$.
\begin{align*}
t & =\frac{u\lambda}{1-u}\\
dt & =\frac{\lambda}{(1-u)^{2}}du
\end{align*}
\begin{align*}
\int_{0}^{\infty}\frac{t^{-\alpha}dt}{\lambda+t} & =\int_{0}^{\infty}(\frac{t}{t+\lambda})t^{-1-\alpha}dt\\
& =\int_{0}^{1}u(\frac{u\lambda}{1-u})^{-1-\alpha}\frac{\lambda}{(1-u)^{2}}du\\
& =\lambda^{-\alpha}\int_{0}^{1}u^{-\alpha}(1-u)^{\alpha-1}du\\
& =\lambda^{-\alpha}B(1-\alpha,\alpha)
\end{align*}
and using Euler's reflection formula, we have
\[
B(1-\alpha,\alpha)=\frac{\pi}{\sin(\pi\alpha)}
\]
Thus,
\[
\int_{0}^{\infty}\frac{t^{-\alpha}dt}{\lambda+t}=\frac{\pi}{\sin(\pi\alpha)}\lambda^{-\alpha}
\]
\end{IEEEproof}
```
::: cor
**Corollary 1** (Integral formula for power one half). *In the special case of $\alpha=1/2$ $$\pi\lambda^{-1/2}=\int_{0}^{\infty}\frac{t^{-1/2}dt}{\lambda+t}=2\int_{0}^{\infty}\frac{ds}{\lambda+s^{2}}$$*
:::
```{=latex}
\begin{IEEEproof}
Use simple change of variable $t\to s^{2}$,
\begin{align*}
\pi\lambda^{-1/2} & =\int_{0}^{\infty}\frac{t^{-1/2}dt}{\lambda+t}
=\int_{0}^{\infty}\frac{s^{-1}ds^{2}}{\lambda+s^{2}}
=\int_{0}^{\infty}\frac{2ds}{\lambda+s^{2}}
\end{align*}
\end{IEEEproof}
```
::: cor
**Corollary 2** (Integral representation of fractional matrix power). *`\label{cor:int_repr_mat_frac_power}`{=latex} The matrix version of such identity, for self-adjoint, positive semi definite matrix $A\succeq0$, $$\int_{0}^{\infty}(A+tI)^{-1}t^{-\alpha}dt=\frac{\pi}{\sin(\pi\alpha)}A^{-\alpha},\quad\alpha\in(0,1)$$*
*Similarly, for $z\geq0,z\in\mathbb{R}$, $$\int_{0}^{\infty}(A+(z+t)I)^{-1}t^{-\alpha}dt=\frac{\pi}{\sin(\pi\alpha)}(A+zI)^{-\alpha},\quad\alpha\in(0,1)$$*
:::
::: cor
**Corollary 3** (Integral representation of matrix one half). *`\label{cor:int_repr_mat_one_half}`{=latex} The matrix version of such identity, for self-adjoint, positive semi definite matrix $A\succeq0$, $$A^{-1/2}=\frac{1}{\pi}\int_{0}^{\infty}(A+tI)^{-1}t^{-1/2}dt=\frac{2}{\pi}\int_{0}^{\infty}(A+s^{2}I)^{-1}ds$$*
:::
::: lemma
**Lemma 10** (Resolvent Identity). *`\label{lem:resolvant_identity}`{=latex}When $u\neq s$, we have the identity $$(A+sI)^{-1}(A+uI)^{-1}=\frac{1}{s-u}\Big((A+uI)^{-1}-(A+sI)^{-1}\Big)$$ $$\begin{aligned}
A(A+sI)^{-1}(A+uI)^{-1} & =\frac{1}{s-u}\Big(A(A+uI)^{-1}-A(A+sI)^{-1}\Big)\\
& =\frac{s(A+sI)^{-1}-u(A+uI)^{-1}}{s-u}
\end{aligned}$$*
:::
```{=latex}
\begin{IEEEproof}
Note that
\begin{align*}
& \Big((A+sI)-(A+uI)\Big)(A+sI)^{-1}(A+uI)^{-1}\\
= & (A+uI)^{-1}-(A+sI)^{-1}\\
= & (s-u)(A+sI)^{-1}(A+uI)^{-1}
\end{align*}
Thus,
\begin{align*}
(A+sI)^{-1}(A+uI)^{-1} & =\frac{1}{(s-u)}\Big((A+uI)^{-1}-(A+sI)^{-1}\Big)
\end{align*}
as corollary
\begin{align*}
A(A+sI)^{-1}(A+uI)^{-1} & =\frac{1}{s-u}\Big(A(A+uI)^{-1}-A(A+sI)^{-1}\Big)\\
& =\frac{1}{s-u}\Big(I-u(A+uI)^{-1}-I+s(A+sI)^{-1}\Big)\\
& =\frac{s(A+sI)^{-1}-u(A+uI)^{-1}}{s-u}
\end{align*}
Note that this formula has no real pole, and it behaves well when
denominator vanishes, and the RHS becomes a derivative.
\begin{align*}
\lim_{s\to u}\frac{s(A+sI)^{-1}-u(A+uI)^{-1}}{s-u} & =\frac{d}{ds}s(A+sI)^{-1}\\
& =(A+sI)^{-1}-s(A+sI)^{-2}\\
& =A(A+sI)^{-2}
\end{align*}
\begin{align*}
\lim_{s\to u}\frac{1}{(s-u)}\Big((A+uI)^{-1}-(A+sI)^{-1}\Big) & =-\frac{d}{du}(A+uI)^{-1}\\
& =(A+uI)^{-2}
\end{align*}
\end{IEEEproof}
```
## Derivation of the sampling-map expectation equivalence: infinite-$\sigma_T$ approximation (Result `\ref{prop:gen_sample_exp_DE}`{=latex}) {#app:proof_samp_mapping_exp_approx}
Using empirical covariance and mean to realize the sampling, we have $$\mathbf{x}(\mathbf{x}_{\sigma_{T}},\sigma_{0})=\hat{\boldsymbol{\mu}}+(\hat{\boldsymbol{\Sigma}}+\sigma_{0}^{2}I)^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}(\mathbf{x}_{\sigma_{T}}-\hat{\boldsymbol{\mu}})$$
For the final sampling outcome $\sigma_{0}\to0$, this reads $$\mathbf{x}(\mathbf{x}_{\sigma_{T}},0)=\hat{\boldsymbol{\mu}}+\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}(\mathbf{x}_{\sigma_{T}}-\hat{\boldsymbol{\mu}})$$
As before, assume the sample mean equals the population one, then the finite sample effect comes from the matrix $\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}$ $$\mathbf{x}(\mathbf{x}_{\sigma_{T}},0)=\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}(\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu})$$
Note that for sampling, under the EDM convention, the initial noise is sampled with variance $\sigma_{T}^{2}I$, $\mathbf{x}_{\sigma_{T}}\sim\mathcal{N}(0,\sigma_{T}^{2}I)$, notably for practical diffusion models, initial noise variances are large, $\sigma_{T}^{2}\sim6000$. Thus we can define a normalized initial noise $\bar{\mathbf{x}}=(\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu})/\sigma_{T}$.
As a large initial noise limit, given that $\boldsymbol{\Sigma}$ has finite spectral norm, $$\lim_{\sigma\to\infty}\sigma\boldsymbol{\Sigma}^{1/2}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1/2}=\boldsymbol{\Sigma}^{1/2}$$ and when $\sigma_{T}\to\infty$ the normalized initial noise is sampled from standard Gaussian, $\bar{\mathbf{x}}\sim\mathcal{N}(0,I)$.
Equivalently, we can consider expansion as orders of $1/\sigma$, $$\begin{aligned}
\sigma\boldsymbol{\Sigma}^{1/2}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1/2} & =\boldsymbol{\Sigma}^{1/2}(I+\frac{1}{\sigma^{2}}\boldsymbol{\Sigma})^{-1/2}\\
& \approx\boldsymbol{\Sigma}^{1/2}(I-\frac{1}{2}\frac{1}{\sigma^{2}}\boldsymbol{\Sigma}+...)\\
& \approx\boldsymbol{\Sigma}^{1/2}-\frac{1}{2}\frac{1}{\sigma^{2}}\boldsymbol{\Sigma}^{3/2}+...
\end{aligned}$$
If we keep the zeroth-order term, then we get the approximation $$\sigma\boldsymbol{\Sigma}^{1/2}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1/2}\approx\boldsymbol{\Sigma}^{1/2}$$
Consider approximation, $$\begin{aligned}
\mathbf{x}(\mathbf{x}_{\sigma_{T}},0) & =\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}(\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu})\\
& \approx\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}(\frac{\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu}}{\sigma_{T}})\\
& =\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}
\end{aligned}$$ then we can study the finite-sample effect on the sampling mapping via the matrix $\hat{\boldsymbol{\Sigma}}^{1/2}$.
::: proposition
**Result 11**. *`\label{prop_app:deter_equiv_cov_mat_half}`{=latex} Deterministic equivalence of empirical covariance matrix square root $$\begin{aligned}
\hat{\boldsymbol{\Sigma}}^{1/2}= & \frac{2}{\pi}\int_{0}^{\infty}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+u^{2}I)^{-1}du \label{eq:matrix_half_int_repr}\\
\asymp & \frac{2}{\pi}\int_{0}^{\infty}\boldsymbol{\Sigma}\Big(\boldsymbol{\Sigma}+\kappa(u^{2})I\Big)^{-1}du\label{eq:emp_cov_half_equivalence-app}
\end{aligned}$$*
:::
```{=latex}
\begin{IEEEproof}
Combining Lemma \ref{cor:int_repr_mat_one_half} with deterministic
equivalence of one point (Eq.~\ref{eq:DE_cov_equiv}).
\end{IEEEproof}
```
This result can be compared to population covariance half, when renormalization effect vanish $\kappa(u^{2})\to u^{2}$. $$\boldsymbol{\Sigma}^{1/2}=\frac{2}{\pi}\int_{0}^{\infty}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+u^{2}I)^{-1}du$$
Since $\kappa(u^{2})>u^{2}$ point by point in the integral, the sample version leads to larger shrinkage. $$\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\mathbf{v}<\mathbf{v}^{\top}\boldsymbol{\Sigma}^{1/2}\mathbf{v}$$
Concretely, if we measure along spectral modes $\mathbf{u}_k$ of population covariance, $$\begin{aligned}
\mathbf{u}_k^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\mathbf{u}_k & \asymp\frac{2}{\pi}\int_{0}^{\infty}\mathbf{u}_k^{\top}\boldsymbol{\Sigma}\Big(\boldsymbol{\Sigma}+\kappa(u^{2})I\Big)^{-1}\mathbf{u}_k\,du\\
& =\frac{2}{\pi}\int_{0}^{\infty}\frac{\lambda_{k}}{\lambda_{k}+\kappa(u^{2})}du\\
& <\frac{2}{\pi}\int_{0}^{\infty}\frac{\lambda_{k}}{\lambda_{k}+u^{2}}du\\
& =\lambda_{k}^{1/2}
\end{aligned}$$
## Derivation of the sampling-map expectation equivalence: finite $\sigma_T$ {#app:proof_samp_mapping_exp_full}
Next, we consider the finite $\sigma_{T}$ case, which involves two matrix square roots and their equivalence. To prove this, we proceed in two steps: 1) use an integral identity to represent matrices of the form $A^{1/2}(A+zI)^{-1/2}$, and 2) apply one-point deterministic equivalence.
::: lemma
**Lemma 12**. *`\label{prop:int_repr_half_resolv}`{=latex}Integral representation, for self-adjoint, positive semidefinite matrix $A\succeq0$, $$\begin{aligned}
A^{1/2}(A+zI)^{-1/2} & =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}A(A+u^{2}I)^{-1}(A+(z+v^{2})I)^{-1}dudv\\
& =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{A(A+(z+u^{2})I)^{-1}-A(A+v^{2}I)^{-1}}{v^{2}-u^{2}-z}dudv
\end{aligned}$$*
:::
```{=latex}
\begin{IEEEproof}
We can study matrix of this form,
\[
A^{1/2}(A+zI)^{-1/2}
\]
using the integral representation above twice, we have
\begin{align*}
& A^{1/2}(A+zI)^{-1/2}\\
= & AA^{-1/2}(A+zI)^{-1/2}\\
= & \frac{1}{\pi^{2}}A\int_{0}^{\infty}(A+sI)^{-1}s^{-1/2}ds\int_{0}^{\infty}(A+(z+t)I)^{-1}t^{-1/2}dt\\
= & \frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}A(A+sI)^{-1}(A+(z+t)I)^{-1}t^{-1/2}s^{-1/2}dsdt\\
= & \frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}A(A+u^{2}I)^{-1}(A+(z+v^{2})I)^{-1}dudv
\end{align*}
To deal with this \emph{product of resolvents}, we can turn it into
a \emph{difference of resolvents} via Lemma \ref{lem:resolvant_identity},
\[
(A+sI)^{-1}(A+tI)^{-1}=\frac{1}{(s-t)}\Big((A+tI)^{-1}-(A+sI)^{-1}\Big)
\]
Now using the identity, we have
\begin{align*}
& A^{1/2}(A+zI)^{-1/2}\\
= & \frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}A(A+sI)^{-1}(A+(z+t)I)^{-1}t^{-1/2}s^{-1/2}dsdt\\
= & \frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{A(A+(z+t)I)^{-1}-A(A+sI)^{-1}}{s-z-t}t^{-1/2}s^{-1/2}dsdt
\end{align*}
Putting it together,
\begin{align*}
A^{1/2}(A+zI)^{-1/2} & =\frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{A(A+(z+t)I)^{-1}-A(A+sI)^{-1}}{s-t-z}t^{-1/2}s^{-1/2}dsdt\\
& =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{A(A+(z+u^{2})I)^{-1}-A(A+v^{2}I)^{-1}}{v^{2}-u^{2}-z}dudv
\end{align*}
\end{IEEEproof}
```
Next we are ready to use the one-point deterministic equivalence.
::: proposition
**Result 13**. *`\label{prop_app:deter_equiv_cov_half_resolv}`{=latex} For sample covariance matrix $\hat{\boldsymbol{\Sigma}}$, the following expression has deterministic equivalent to the double integral of population covariance, $$\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1/2}\asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\kappa(\sigma^{2}+u^{2})-\kappa(v^{2})}{(\sigma^{2}+u^{2})-v^{2}}\boldsymbol{\Sigma}\big(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+u^{2})I\big)^{-1}\big(\boldsymbol{\Sigma}+\kappa(v^{2})I\big)^{-1}dudv$$*
:::
```{=latex}
\begin{IEEEproof}
Using Result \ref{prop:int_repr_half_resolv}, set $A\to\hat{\boldsymbol{\Sigma}}$.
We can then apply deterministic equivalence for
resolvents:
\begin{align*}
\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1/2} & =\hat{\boldsymbol{\Sigma}}\hat{\boldsymbol{\Sigma}}^{-1/2}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1/2}\\
& =\frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+sI)^{-1}(\hat{\boldsymbol{\Sigma}}+(\sigma^{2}+t)I)^{-1}t^{-1/2}s^{-1/2}dsdt\\
& =\frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+(\sigma^{2}+t)I)^{-1}-\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+sI)^{-1}}{s-t-\sigma^{2}}t^{-1/2}s^{-1/2}dsdt\\
& \asymp\frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+t)I)^{-1}-\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(s)I)^{-1}}{s-t-\sigma^{2}}t^{-1/2}s^{-1/2}dsdt
\end{align*}
Note there is no pole in this double integral, i.e. when $s=t+\sigma^{2}$,
$\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+t)I)^{-1}=\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(s)I)^{-1}$,
thus both numerator and denominator vanish, and the limit is well
defined as a derivative!
\begin{align*}
RHS= & \frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+t)I)^{-1}-\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(s)I)^{-1}}{s-t-\sigma^{2}}t^{-1/2}s^{-1/2}dsdt\\
= & \frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{(\kappa(s)-\kappa(\sigma^{2}+t))\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+t)I)^{-1}(\boldsymbol{\Sigma}+\kappa(s)I)^{-1}}{s-t-\sigma^{2}}t^{-1/2}s^{-1/2}dsdt\\
= & \frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\kappa(s)-\kappa(\sigma^{2}+t)}{s-(\sigma^{2}+t)}\boldsymbol{\Sigma}\big(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+t)I\big)^{-1}\big(\boldsymbol{\Sigma}+\kappa(s)I\big)^{-1}t^{-1/2}s^{-1/2}dsdt
\end{align*}
This formulation shows that there are no real poles.
We can remove the singularity at $0$ via the change of variables $t\to u^{2},s\to v^{2}$:
\begin{align*}
RHS & =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\kappa(\sigma^{2}+u^{2})-\kappa(v^{2})}{(\sigma^{2}+u^{2})-v^{2}}\!\boldsymbol{\Sigma}\!\big(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+u^{2})I\big)^{-1}\big(\boldsymbol{\Sigma}+\kappa(v^{2})I\big)^{-1}\!dudv
\end{align*}
Thus, we obtain the desired equivalence,
\[
\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma^{2}I)^{-1/2}\asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{\kappa(\sigma^{2}+u^{2})-\kappa(v^{2})}{(\sigma^{2}+u^{2})-v^{2}}\!\boldsymbol{\Sigma}\!\big(\boldsymbol{\Sigma}+\kappa(\sigma^{2}+u^{2})I\big)^{-1}\big(\boldsymbol{\Sigma}+\kappa(v^{2})I\big)^{-1}\!dudv
\]
Note that the coefficient $\frac{\kappa(\sigma^{2}+u^{2})-\kappa(v^{2})}{(\sigma^{2}+u^{2})-v^{2}}$
behaves well when $(\sigma^{2}+u^{2})-v^{2}\to0$, i.e. it becomes
a derivative of $\kappa$ (Lemma \ref{lem:resolvant_identity}). Thus,
there is no singularity in the integrand.
\end{IEEEproof}
```
#### Interpretation
We can compare it to the sampling mapping with the population covariance, i.e. infinite data limit. Using Prop `\ref{prop:int_repr_half_resolv}`{=latex}, setting $A\to\boldsymbol{\Sigma}$, the double integral representation of the denoiser mapping reads, $$\begin{aligned}
\boldsymbol{\Sigma}^{1/2}(\boldsymbol{\Sigma}+\sigma^{2}I)^{-1/2} & =\frac{1}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\boldsymbol{\Sigma}\Big(\boldsymbol{\Sigma}+(\sigma^{2}+t)I\Big)^{-1}(\boldsymbol{\Sigma}+sI)^{-1}t^{-1/2}s^{-1/2}\!dsdt\\
& =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\boldsymbol{\Sigma}\Big(\boldsymbol{\Sigma}+(\sigma^{2}+u^{2})I\Big)^{-1}(\boldsymbol{\Sigma}+v^{2}I)^{-1}\!dudv
\end{aligned}$$
Indeed, since $\kappa(\sigma^{2}+u^{2})>(\sigma^{2}+u^{2})$ and $\kappa(v^{2})>v^{2}$, this creates a larger shrinkage, especially at small eigen dimensions.
## Derivation of the sampling-map fluctuation equivalence: infinite-$\sigma_T$ approximation (Result `\ref{prop:gen_sample_var_DE}`{=latex}) {#app:proof_samp_mapping_var_approx}
Now let us consider the variance of the generated outcome with the infinite $\sigma_{T}$ approximation, ignoring estimation error in $\boldsymbol{\mu}$, $$\begin{aligned}
\mathbf{x}_{\sigma_{0}} & =\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}(\hat{\boldsymbol{\Sigma}}+\sigma_{T}^{2}I)^{-1/2}(\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu})\\
& \approx\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}(\frac{\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu}}{\sigma_{T}})\\
& =\boldsymbol{\mu}+\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}
\end{aligned}$$ Thus, the variance of the generated output comes from the sample estimate of the covariance. Let $\bar{\mathbf{x}}:=\frac{\mathbf{x}_{\sigma_{T}}-\boldsymbol{\mu}}{\sigma_{T}}$, i.e. the normalized deviation from the center.
::: proposition
**Result 14** (Main result, variance of generated sample under empirical data covariance.). *$$\begin{aligned}
Var_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}] & \asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\mathbf{v}\big]\!\\
& \quad\quad\times\big[\bar{\mathbf{x}}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\bar{\mathbf{x}}\big]\!\Big\}\!dudv
\end{aligned}$$ where $\kappa:=\kappa(u^{2}),\kappa^{\prime}:=\kappa(v^{2})$ are the variables to be integrated over.*
:::
```{=latex}
\begin{IEEEproof}
Representing the variance by moments,
\begin{align*}
& \quad Var_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}]\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}[(\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}})^{2}]-\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}]^{2}\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}\bar{\mathbf{x}}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\mathbf{v}]-\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}]\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}[\mathbf{\bar{\mathbf{x}}}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\mathbf{v}]\tag*{\textit{using Eq.~\ref{eq:matrix_half_int_repr}}}\\
& =\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big\{\mathbf{v}^{\top}\Big[\frac{2}{\pi}\int_{0}^{\infty}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+u^{2}I)^{-1}du\Big]\bar{\mathbf{x}}\bar{\mathbf{x}}^{\top}\Big[\frac{2}{\pi}\int_{0}^{\infty}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+v^{2}I)^{-1}dv\Big]\mathbf{v}\Big\}\\
& \quad-\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big\{\mathbf{v}^{\top}\Big[\frac{2}{\pi}\int_{0}^{\infty}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+u^{2}I)^{-1}du\Big]\bar{\mathbf{x}}\Big\}\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big\{\bar{\mathbf{x}}^{\top}\Big[\frac{2}{\pi}\int_{0}^{\infty}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+v^{2}I)^{-1}dv\Big]\mathbf{v}\Big\}
\end{align*}
Using the integral representation and exchanging the integral with
expectation,
\begin{align*}
RHS & =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big\{\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+u^{2}I)^{-1}\bar{\mathbf{x}}\bar{\mathbf{x}}^{\top}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+v^{2}I)^{-1}\mathbf{v}\Big\}\\
& \quad-\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+u^{2}I)^{-1}\bar{\mathbf{x}}\Big]\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big[\bar{\mathbf{x}}^{\top}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+v^{2}I)^{-1}\mathbf{v}\Big]\ \Big\}\!dudv\tag*{\textit{ integral representation of matrix half}}\\
& \asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\mathbb{E}_{\hat{\boldsymbol{\Sigma}}}\Big\{\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+u^{2}I)^{-1}\bar{\mathbf{x}}\bar{\mathbf{x}}^{\top}\hat{\boldsymbol{\Sigma}}(\hat{\boldsymbol{\Sigma}}+v^{2}I)^{-1}\mathbf{v}\Big\}\\
& \quad-\Big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(u^{2})I)^{-1}\bar{\mathbf{x}}\Big]\ \Big[\bar{\mathbf{x}}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(v^{2})I)^{-1}\mathbf{v}\Big]\ \Big\}\!dudv\tag*{\textit{ using one point equivalence}}\\
& \asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\operatorname{Tr}\big[\mathbf{v}\mathbf{v}^{\top}T_{\boldsymbol{\Sigma}}\bar{\mathbf{x}}\bar{\mathbf{x}}^{\top}T_{\boldsymbol{\Sigma}}^{\prime}\big]+\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\operatorname{Tr}\big[\mathbf{v}\mathbf{v}^{\top}G_{\boldsymbol{\Sigma}}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}^{\prime}\big]\operatorname{Tr}\big[\bar{\mathbf{x}}\bar{\mathbf{x}}^{\top}G_{\boldsymbol{\Sigma}}^{\prime}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}\big]\\
& \quad-\Big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(u^{2})I)^{-1}\bar{\mathbf{x}}\Big]\ \Big[\bar{\mathbf{x}}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(v^{2})I)^{-1}\mathbf{v}\Big]\ \Big\}\!dudv\tag*{\textit{ using two point equivalence}}\\
& =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\operatorname{Tr}\big[\mathbf{v}\mathbf{v}^{\top}G_{\boldsymbol{\Sigma}}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}^{\prime}\big]\operatorname{Tr}\big[\bar{\mathbf{x}}\bar{\mathbf{x}}^{\top}G_{\boldsymbol{\Sigma}}^{\prime}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}\big]\!\Big\}\!dudv\tag*{\textit{ first trace cancels out.}}\\
& =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\big[\mathbf{v}^{\top}G_{\boldsymbol{\Sigma}}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}^{\prime}\mathbf{v}\big]\!\big[\bar{\mathbf{x}}^{\top}G_{\boldsymbol{\Sigma}}^{\prime}\boldsymbol{\Sigma} G_{\boldsymbol{\Sigma}}\bar{\mathbf{x}}\big]\!\Big\}\!dudv\\
& =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(v^{2})I)^{-1}(\boldsymbol{\Sigma}+\kappa(u^{2})I)^{-1}\mathbf{v}\big]\!\\
& \quad\quad\times\big[\bar{\mathbf{x}}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa(v^{2})I)^{-1}(\boldsymbol{\Sigma}+\kappa(u^{2})I)^{-1}\bar{\mathbf{x}}\big]\!\Big\}\!dudv\\
& =\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\mathbf{v}\big]\!\\
& \quad\quad\times\big[\bar{\mathbf{x}}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\bar{\mathbf{x}}\big]\!\Big\}\!dudv
\end{align*}
Thus, we arrive at our result
\begin{align*}
Var[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}] & \asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\mathbf{v}\big]\!\\
& \quad\quad\times\big[\bar{\mathbf{x}}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\bar{\mathbf{x}}\big]\!\Big\}\!dudv
\end{align*}
\end{IEEEproof}
```
As in Remark `\ref{rmk:generic_vector}`{=latex}, this retains the generic, factorized contribution and discards the possible exchange/alignment term.
### Interpretation and derivations
#### Anisotropy: effect of the probe vector
If we marginalize over $\bar{\mathbf{x}}$, assuming $\bar{\mathbf{x}}\sim\mathcal{N}(0,I)$ from white noise, and consider only the effect of probe direction $\mathbf{v}$, $$\begin{aligned}
\mathbb{E}_{\bar{\mathbf{x}}}Var[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}] & \asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa'}{n-\mathrm{df}_{2}(\kappa,\kappa')}\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\mathbf{v}\big]\!\\
& \quad\quad\times\mathbb{E}_{\bar{\mathbf{x}}}\big[\bar{\mathbf{x}}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\bar{\mathbf{x}}\big]\!\Big\}\!dudv\\
& \asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa\kappa^{\prime}\operatorname{Tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}]}{n-\mathrm{df}_{2}(\kappa,\kappa^{\prime})}\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\mathbf{v}\big]\!\Big\}\!dudv\\
\end{aligned}$$
The mixed resolvent term can be simplified, $$\begin{aligned}
\operatorname{Tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}] & =\frac{1}{\kappa}\operatorname{Tr}[(\boldsymbol{\Sigma}+\kappa I-\boldsymbol{\Sigma})\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}]\\
& =\frac{1}{\kappa}\operatorname{Tr}[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}]-\frac{1}{\kappa}\operatorname{Tr}[\boldsymbol{\Sigma}^{2}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}]\\
& =\frac{1}{\kappa}\mathrm{df}_{1}(\kappa^{\prime})-\frac{1}{\kappa}\mathrm{df}_{2}(\kappa,\kappa^{\prime})\\
& =\frac{1}{\kappa^{\prime}}\mathrm{df}_{1}(\kappa)-\frac{1}{\kappa^{\prime}}\mathrm{df}_{2}(\kappa,\kappa^{\prime})
\end{aligned}$$
Using this identity $$\begin{aligned}
\mathbb{E}_{\bar{\mathbf{x}}}Var[\mathbf{v}^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}] & \asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa^{\prime}\big(\mathrm{df}_{1}(\kappa^{\prime})-\mathrm{df}_{2}(\kappa,\kappa^{\prime})\big)}{n-\mathrm{df}_{2}(\kappa,\kappa^{\prime})}\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa^{\prime}I)^{-1}\mathbf{v}\big]\!\Big\}\!dudv
\end{aligned}$$ Setting the direction $\mathbf{v}$ to the eigenvector $\mathbf{u}_k$, with corresponding eigenvalue $\lambda_{k}$, the variance along the eigen direction reads, $$\mathbb{E}_{\bar{\mathbf{x}}}Var[\mathbf{u}_k^{\top}\hat{\boldsymbol{\Sigma}}^{1/2}\bar{\mathbf{x}}]\asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa^{\prime}\big(\mathrm{df}_{1}(\kappa^{\prime})-\mathrm{df}_{2}(\kappa,\kappa^{\prime})\big)}{n-\mathrm{df}_{2}(\kappa,\kappa^{\prime})}\ \frac{\lambda_{k}}{(\lambda_{k}+\kappa)(\lambda_{k}+\kappa^{\prime})}\!\Big\}\!dudv\$$
#### Inhomogeneity: effect of initial noise
Since the variance is symmetric in $\bar{\mathbf{x}}$ and $\mathbf{v}$, we can also marginalize over $\mathbf{v}$ while keeping the $\bar{\mathbf{x}}$ dependency. Note that we assume $\mathbf{v}$ is unit norm, so summation over $\mathbf{u}_k$ eigenvectors (instead of expectation) is equivalent to trace. $$\sum_{k}Var[\mathbf{u}_{k}^{\top}\hat{\Sigma}^{1/2}\bar{\mathbf{x}}] =\operatorname{Tr} \operatorname{Var}[\hat{\Sigma}^{1/2}\bar{\mathbf{x}}]\asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\kappa^{\prime}\big(\mathrm{df}_{1}(\kappa^{\prime})-\mathrm{df}_{2}(\kappa,\kappa^{\prime})\big)}{n-\mathrm{df}_{2}(\kappa,\kappa^{\prime})}\big[\bar{\mathbf{x}}^{\top}\Sigma(\Sigma+\kappa I)^{-1}(\Sigma+\kappa^{\prime}I)^{-1}\bar{\mathbf{x}}\big]\!\Big\}\!dudv$$
#### Scaling: effect of sample size
Finally, marginalizing over both factors, we have the overall scaling. $$\mathbb{E}_{\bar{\mathbf{x}}}\sum_{k}\operatorname{Var}[\mathbf{u}_{k}^{\top}\hat{\Sigma}^{1/2}\bar{\mathbf{x}}]=\mathbb{E}_{\bar{\mathbf{x}}}\operatorname{Tr} \operatorname{Var}[\hat{\Sigma}^{1/2}\bar{\mathbf{x}}]\asymp\frac{4}{\pi^{2}}\int_{0}^{\infty}\int_{0}^{\infty}\Big\{\frac{\big(\mathrm{df}_{1}(\kappa^{\prime})-\mathrm{df}_{2}(\kappa,\kappa^{\prime})\big)\!\big(\mathrm{df}_{1}(\kappa)-\mathrm{df}_{2}(\kappa,\kappa^{\prime})\big)}{n-\mathrm{df}_{2}(\kappa,\kappa^{\prime})}\Big\}\!dudv$$
# Experimental Details
## Numerical methods {#app:numerical_methods}
#### Numerical evaluation of renormalized Ridge $\kappa(z)$.
We computed $\kappa(z)$ as the solution to the self-consistent Silverstein equation $$\kappa(z) - z \;=\; \gamma \sum_{k=1}^p w_k \,\frac{\kappa(z)\,\lambda_k}{\kappa(z) + \lambda_k},$$ where $\{\lambda_k\}$ are the eigenvalues of $\boldsymbol{\Sigma}$ and $\{w_k\}$ are their normalized weights. For scalar $z$, we solved this nonlinear equation using Newton's method with analytical derivative. Using implicit derivative with respect to $\kappa$, we have $$1-\frac {dz}{d\kappa}
=\gamma \sum_{k=1}^p w_k \,\frac {d}{d\kappa}\Big[\frac{\kappa(z)\,\lambda_k}{\kappa(z) + \lambda_k}\Big]$$ $$\big(\kappa'(z)\big)^{-1} \;=\; 1 \;-\; \gamma \sum_{k=1}^p w_k \,\frac{\lambda_k^2}{\big(\kappa(z)+\lambda_k\big)^2},$$ This falls back to a robust root-finder for purely real inputs. For a sequence of $z$ values along a path, we used an \`\`analytic continuation" procedure in which the solution at the previous $z$ served as the initial guess for the next, ensuring branch continuity and numerical stability, particularly for small $z$. Further, we generally start the path from $z$ with high norm and solve with continuation back to small $z$. A caching mechanism stored previously computed $(z,\,\kappa)$ pairs, with nearest-neighbor retrieval for initial guesses, further accelerating repeated evaluations. This approach yields accurate and smooth $\kappa(z)$ profiles suitable for downstream quadrature-based integration.
#### Numerical evaluation of the integral over deterministic equivalence
The analytical results in Eqs. `\ref{eq:gen_sample_exp_DE}`{=latex},`\ref{eq:gen_sample_var_DE}`{=latex}, which involve integrals to infinity, are nontrivial to evaluate. To avoid truncation error, we used the following scheme, which maps the integration variable onto a finite domain.
We approximated the double integral $$\frac{4}{\pi^{2}}\int_{0}^{\infty}\!\!\int_{0}^{\infty}
\frac{\kappa\,\kappa'\,\operatorname{Tr}\!\big[\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa' I)^{-1}\big]}
{n-\operatorname{df}_{2}(\kappa,\kappa')}
\big[\mathbf{v}^{\top}\boldsymbol{\Sigma}(\boldsymbol{\Sigma}+\kappa I)^{-1}(\boldsymbol{\Sigma}+\kappa' I)^{-1}\mathbf{v}\big]\,du\,dv$$ using a Gauss--Legendre quadrature scheme combined with the tangent mapping $u = \tan\theta$ to transform the semi-infinite domain $[0,\infty)$ to a finite interval $[0,\pi/2)$.
We first generated $n_{\mathrm{nodes}}$ Gauss--Legendre nodes $\theta_i$ and weights $w_i$ on $[0,\pi/2]$, then applied the transformation $u = \tan\theta$ with Jacobian $J(\theta) = 1/\cos^2\theta$ to obtain quadrature points on $[0,\infty)$. This was performed independently for $u$ and $v$, and their 2D tensor product provided the integration grid.
The $\kappa$ values were computed at each $u^2$ and $v^2$ using a numerically stable, vectorized evaluation of the spectral mapping function $\kappa(z)$ derived from the eigenspectrum of $\boldsymbol{\Sigma}$. The integrand was then assembled by evaluating the trace term, the scalar bilinear form $\mathbf{v}^{\top}(\cdot)\mathbf{v}$, and the denominator $n - \operatorname{df}_2(\kappa,\kappa')$ on the full 2D grid. Quadrature weights and Jacobians were applied multiplicatively, and the sum over all grid points yielded the numerical approximation to the integral.
A similar quadrature is used for the single-integral equivalence Eq. `\ref{eq:gen_sample_exp_DE}`{=latex}, where we integrate over a 1D grid.
This approach yields high accuracy while avoiding explicit truncation of the infinite domain, as the nonlinear mapping concentrates quadrature nodes where the integrand varies most rapidly.
## Linear denoiser experiments {#app:linear_denoiser_method}
To cross-validate against our theory and numerical scheme, we performed extensive validation using empirical linear denoisers.
We computed the empirical covariance of a dataset and then used the following functions implementing the linear one-step denoiser and the full sampling map (Wiener filter).
```{=latex}
\begin{RoundedListing}
def dnoised_X(x, Xmean, sample_cov, sigma2,):
# single step denoiser
return x + sigma2 * (Xmean - x ) @ torch.inverse(sample_cov + torch.eye(sample_cov.shape[0], device=x.device) * sigma2)
def wiener_gen_X(x, Xmean, wiener_matrix, sigmaT,):
if x.dim() == 1:
# Single vector case
return Xmean + wiener_matrix @ (x * sigmaT - Xmean)
else:
# Batched vector case - x should be shape (batch_size, ndim)
return Xmean[None,:] + (x * sigmaT - Xmean[None,:]) @ wiener_matrix.T
def build_wiener_matrix(eigvals, eigvecs, sigmaT=80.0, sigma0=0.0, EPS=1E-16, clip=True):
if clip:
eigvals = torch.clamp(eigvals, min=EPS)
scaling = ((eigvals + sigma0**2) / (eigvals + sigmaT**2)).sqrt()
return eigvecs @ torch.diag(scaling) @ eigvecs.T
\end{RoundedListing}
```
We set $\sigma_0=0$ in theory. In practice, it is usually set to a small positive number, e.g., 0.002. We tested this in a few cases and report the results in the appendix. Generally, a positive $\sigma_0$ acts as a floor for generated variance, thereby reducing the overshrinking effect.
We found that when the dataset size is insufficient, e.g., when $\hat{\boldsymbol{\Sigma}}$ is rank deficient, the eigendecomposition can be unstable and sometimes produces negative eigenvalues, which affects the matrix square-root operation in the Wiener matrix. Even if we clip them, numerical artifacts often remain in small eigenspaces. As a remedy, using higher precision `float64` numbers yields results that closely match the theory.
## Deep neural network experiments {#app:DNN_method}
We used the following preconditioning scheme, inspired by [@karras2022elucidatingDesignSp], for all architectures in our comparison.
```{=latex}
\begin{RoundedListing}
class EDMPrecondWrapper(nn.Module):
def __init__(self, model, sigma_data=0.5, sigma_min=0.002, sigma_max=80, rho=7.0):
super().__init__()
self.model = model
self.sigma_data = sigma_data
self.sigma_min = sigma_min
self.sigma_max = sigma_max
self.rho = rho
def forward(self, X, sigma, cond=None, ):
sigma[sigma == 0] = self.sigma_min
## edm preconditioning for input and output
## https://github.com/NVlabs/edm/blob/main/training/networks.py#L632
# unsqueeze sigma to have same dimension as X (which may have 2-4 dim)
sigma_vec = sigma.view([-1, ] + [1, ] * (X.ndim - 1))
c_skip = self.sigma_data ** 2 / (sigma_vec ** 2 + self.sigma_data ** 2)
c_out = sigma_vec * self.sigma_data / (sigma_vec ** 2 + self.sigma_data ** 2).sqrt()
c_in = 1 / (self.sigma_data ** 2 + sigma_vec ** 2).sqrt()
c_noise = sigma.log() / 4
model_out = self.model(c_in * X, c_noise, cond=cond)
return c_skip * X + c_out * model_out
\end{RoundedListing}
```
#### EDM Loss Function
We employ the loss function $\mathcal{L}_\text{EDM}$ introduced in the Elucidated Diffusion Model (EDM) paper [@karras2022elucidatingDesignSp], which is one specific weighting scheme for training diffusion models with $x_0$ prediction with variance-exploding noise schedule.
For each data point $\mathbf{x}\in \mathbb R^d$, the loss is computed as follows. The noise level for each data point is sampled from a log-normal distribution with hyperparameters $P_{\text{mean}}$ and $P_{\text{std}}$ (e.g., $P_{\text{mean}} = -1.2$ and $P_{\text{std}} = 1.2$). Specifically, the noise level $\sigma$ is sampled via $$\sigma = \exp\big(P_{\text{mean}} + P_{\text{std}}\,\epsilon\big),\quad \epsilon \sim \mathcal{N}(0,1).$$ The weighting function per noise scale is defined as: $$w(\sigma) = \frac{\sigma^2 + \sigma_{\text{data}}^2}{\left(\sigma\,\sigma_{\text{data}}\right)^2},$$ with hyperparameter $\sigma_{\text{data}}$ (e.g., $\sigma_{\text{data}} = 0.5$). The noisy input $\mathbf{y}$ is created by the following, $$\mathbf{y} = \mathbf{x} + \sigma \mathbf{n},\quad \mathbf{n} \sim \mathcal{N}\big(\mathbf{0},\, \mathbf{I}_d\big),$$ Let $\mathbf{D}_\theta(\mathbf{y}, \sigma, \text{labels})$ denote the output of the denoising network when given the noisy input $\mathbf{y}$, the noise level $\sigma$, and optional conditioning labels. The EDM loss per data point can be computed as: $$\mathcal{L}(\mathbf{x}) = w(\sigma) \, \| \mathbf{D}_\theta(\mathbf{x}+\sigma \mathbf{n}, \sigma, \text{labels}) - \mathbf{x} \|^2.$$ Taking expectation over the data points and noise scales, the overall loss reads $$\begin{aligned}
\mathcal{L}_{EDM}=\mathbb{E}_{\mathbf{x}\sim p_{data}}\mathbb{E}_{\mathbf{n}\sim \mathcal{N}(0,\mathbf{I}_d)}\mathbb{E}_{\sigma} \Big[w(\sigma) \, \| \mathbf{D}_\theta(\mathbf{x}+\sigma \mathbf{n}, \sigma, \text{labels}) - \mathbf{x} \|^2\Big]
\end{aligned}$$
#### Hyperparameter Settings: DiT
All experiments use DiT backbones with consistent architectural and optimization settings unless otherwise specified. Key hyperparameters:
- **Model architecture:** patch size 2 or 4 (used once for FFHQ64, discarded for worse performance), hidden size 384, depth 6 layers, 6 attention heads, MLP ratio 4.
- **Datasets:** FFHQ-32, AFHQ-32, CIFAR-32, and FFHQ-64; subsampled at varying sizes (300, 1k, 3k, 10k, 30k) with two non-overlapping splits per size.
- **Training objective:** Denoising Score Matching (DSM) under EDM parametrization.
- **Training schedule:** 50000 steps with batch size 256, Adam optimizer with learning rate $1\times 10^{-4}$.
- **Evaluation:** fixed-noise seed, sampling with 35 steps with Heun sampler; evaluation sample size 1000, batch size 512.
#### Hyperparameter Settings: UNet
All CNN-UNet experiments follow consistent architectural and optimization settings unless noted. Key hyperparameters:
- **Model architecture:** UNet with base channels $128$ ; channel multipliers $\{1,2,2,2\}$; self-attention at resolution $8$.
- **Datasets:** FFHQ-32/64, AFHQ-32, CIFAR10-32, CIFAR100-32, LSUN-church-32/64, LSUN-bedroom-32/64; subsampled at varying sizes (300, 1k, 3k, 10k, 30k) with two non-overlapping splits per size.
- **Training objective:** Denoising Score Matching (DSM) under EDM parametrization.
- **Training schedule:** 50000 steps, batch size $256$, Adam with learning rate $1\times10^{-4}$.
- **Evaluation:** fixed-noise seed, sampling with 35 steps with Heun sampler; evaluation sample size 1000, batch size 512.
#### Computation Cost
All experiments were conducted on NVIDIA A100 or H100 GPUs. Training DiT and CNN models on $32\times 32$ resolution datasets typically required $5$--$8$ hours to complete. In contrast, DiT models trained on FFHQ64 were substantially more expensive, taking approximately $24$ hours per run.
# Usage of LLMs
We used LLMs in three ways. First, as a research assistant, to look up tools related to deterministic equivalence and to point us toward integral identities for fractional matrix powers, which we then verified and derived independently. Second, as a coding agent to help us generate plotting and analysis code for our results. Third, as a writing aid, for polishing technical text and providing feedback on clarity and presentation of the whole paper.
[^1]: We write $A_n \asymp B_n$ for *deterministic equivalence*: for any sequence of deterministic matrices $C_n$ with uniformly bounded spectral norm, $\operatorname{tr}\!\left[C_n (A_n - B_n)\right] \to 0 \quad \text{as } d,n\to\infty,\ d/n\to\gamma.$. For the case of rank-one observables formed of sequences of bounded-norm vectors $\mathbf{u}_{n}$ and $\mathbf{v}_{n}$, the corresponding equivalence statement is that $\mathbf{u}_{n}^{\top} (A_{n}-B_{n}) \mathbf{v}_{n} \to 0$. Equivalences of scalar expressions are denoted similarly with $\asymp$. We state the detailed technical assumptions of our results in Appendix `\ref{app:deterministic_equivalence}`{=latex}.
[^2]: *We note that this result assumes $\mathbf{v}$ and $\mathbf{x}-\boldsymbol{\mu}$ are generic. When they are strongly aligned, additional exchange terms proportional to $[\mathbf{v}^{\top}\boldsymbol{\Sigma} (\boldsymbol{\Sigma}+\kappa I)^{-2} (\mathbf{x}-\boldsymbol{\mu})]^{2}$ can contribute at leading order; see Appendix `\ref{app:proof_denoiser_var_de}`{=latex} for details.*
[^3]: Though the sample covariance $\hat{\boldsymbol{\Sigma}}$ is an unbiased estimator of the population covariance $\boldsymbol{\Sigma}$, taking the square root introduces finite sample bias: $\boldsymbol{\Sigma} = \mathbb{E}[\hat{\boldsymbol{\Sigma}}]=\mathbb{E}[\hat{\boldsymbol{\Sigma}}^{1/2}\hat{\boldsymbol{\Sigma}}^{1/2}]\neq(\mathbb{E}[\hat{\boldsymbol{\Sigma}}^{1/2}])^2$.
[^4]: For finite data, the underlying population distribution is unknown, so there remains ambiguity in how the model generalizes beyond the empirical distribution; nevertheless, the score of the empirical distribution itself is uniquely defined.
