Computational Bottlenecks for Denoising Diffusions

## Abstract

Abstract Denoising diffusions sample from a probability distribution μ 𝜇 \mu italic_μ in ℝ d superscript ℝ 𝑑 \mathbb{R}^{d} blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT by constructing a stochastic process ( 𝒙 ^ t : t ≥ 0 ) : subscript ^ 𝒙 𝑡 𝑡 0 (\hat{\boldsymbol{x}}_{t}:t\geq 0) ( over^ start_ARG bold_italic_x end_ARG start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT : italic_t ≥ 0 ) in ℝ d superscript ℝ 𝑑 {\mathbb{R}}^{d} blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT such that 𝒙 ^ 0 subscript ^ 𝒙 0 \hat{\boldsymbol{x}}_{0} over^ start_ARG bold_italic_x end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is easy to sample, but the distribution
of 𝒙 ^ T subscript ^ 𝒙 𝑇 \hat{\boldsymbol{x}}_{T} over^ start_ARG bold_italic_x end_ARG start_POSTSUBSCRIPT italic_T end_POSTSUBSCRIPT at large T 𝑇 T italic_T approximates μ 𝜇 \mu italic_μ . The drift 𝒎 : ℝ d × ℝ → ℝ d : 𝒎 → superscript ℝ 𝑑 ℝ superscript ℝ 𝑑 \boldsymbol{m}:{\mathbb{R}}^{d}\times{\mathbb{R}}\to{\mathbb{R}}^{d} bold_italic_m : blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT × blackboard_R → blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT of this diffusion process is learned my minimizing a
score-matching objective. Is every probability distribution μ 𝜇 \mu italic_μ , for which sampling is tractable, also amenable to sampling via diffusions?
We provide evidence to the contrary by studying a probability distribution μ 𝜇 \mu italic_μ for which sampling is easy, but the drift of the diffusion process is intractable —under a popular conjecture on information-computation gaps in statistical estimation.
We show that there exist drifts that are superpolynomially close to the optimum value
(among polynomial time drifts) and yet yield samples with distribution that is very far from the target one.

## 1 Introduction

### 1.1 Background

Diffusion sampling (DS) (Song and Ermon, 2019; Ho et al., 2020)
has emerged as a central paradigm in generative artificial intelligence (AI).
Given a target distribution $\mu$ on ${\mathbb{R}}^{d}$, we want to sample $\boldsymbol{x}\sim\mu$.
Diffusions achieve this goal by generating trajectories of a stochastic process
$(\hat{\boldsymbol{x}}_{t}:0\leq t\leq T)$ whose state $\hat{\boldsymbol{x}}_{T}$ at large $T$ is approximately distributed according to $\mu$.
This suggests a natural question:

- Q:
Are there distributions $\mu$
for which sampling via diffusions fails even if sampling from $\mu$ is easy?

In order to explain how DS might fail, it is useful to recall the setup and introduce some notations(^1^11We follow the formulation of Montanari (2023), which does not require time reversal.).
The basic approach implements an approximation of the following
stochastic differential equation (SDE), with initialization $\boldsymbol{y}_{0}=\boldsymbol{0}$

$$ $\displaystyle{\rm d}{\boldsymbol{y}_{t}}$ $\displaystyle=\boldsymbol{m}(\boldsymbol{y}_{t};t){\rm d}{t}+{\rm d}{ \boldsymbol{B}_{t}},$ (1) $\displaystyle\boldsymbol{m}(\boldsymbol{y},t)$ $\displaystyle:=\mathbb{E}\{\boldsymbol{x}|t\boldsymbol{x}+\sqrt{t}\boldsymbol{ g}=\boldsymbol{y}\}\,,$ (2) $$

where $(\boldsymbol{B}_{t})_{t\geq 0}$ is Brownian motion and in the second line
$\boldsymbol{x}\sim\mu$ is independent of $\boldsymbol{g}\sim{\mathsf{N}}(\boldsymbol{0},\boldsymbol{I}_{d})$.

It is not hard to show that, if $\boldsymbol{y}_{t}$ is generated according the the above SDE,
then there exists $\boldsymbol{x}\sim\mu$ and an independent standard Brownian motion $(\boldsymbol{W}_{t})_{t\geq 0}$ (different from $(\boldsymbol{B}_{t})_{t\geq 0}$)
such that

$$ $\displaystyle\boldsymbol{y}_{t}=t\,\boldsymbol{x}+\boldsymbol{W}_{t}\,.$ (3) $$

Therefore, running the diffusion ([1](https://arxiv.org/html/2503.08028v2#S1.E1)) until some large time $T$,
and returning $\boldsymbol{y}_{T}/T$ or $\boldsymbol{m}(\boldsymbol{y}_{T},T)$ yields a sample approximately distributed according to $\mu$.

A standard Euler discretization of the SDE
([1](https://arxiv.org/html/2503.08028v2#S1.E1)) is given by
(for $(\hat{\boldsymbol{z}}_{t})_{t\in{\mathbb{N}}\Delta}\sim_{iid}{\mathsf{N}}(
\boldsymbol{0},\boldsymbol{I}_{d})$):

$$ $\displaystyle\hat{\boldsymbol{y}}_{t+\Delta}$ $\displaystyle=\hat{\boldsymbol{y}}_{t}+\hat{\boldsymbol{m}}(\hat{\boldsymbol{y }}_{t},t)\Delta+\sqrt{\Delta}\,\hat{\boldsymbol{z}}_{t}\,,$ (4) $$

with $\Delta$ a small stepsize.
The algorithm outputs $\hat{\boldsymbol{x}}_{T}=\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{T},T)$ for some suitably large $T$.
Here $\hat{\boldsymbol{m}}(\,\cdot\,,t)$ is the estimated drift.

Diffusions reduce the problem of sampling to the problem of
approximating the conditional expectation ([2](https://arxiv.org/html/2503.08028v2#S1.E2)). The mapping
$\boldsymbol{y}\mapsto\boldsymbol{m}(\boldsymbol{y},t)$ is the Bayes-optimal estimator
of $\boldsymbol{x}$ in Gaussian noise:

$$ $\displaystyle\boldsymbol{m}(\,\cdot\,,t)=\underset{{\boldsymbol{\varphi}}:{ \mathbb{R}}^{n}\to{\mathbb{R}}^{n}}{\operatorname{argmin}}\;\mathbb{E}\big{\{} \|{\boldsymbol{\varphi}}(\boldsymbol{y}_{t})-\boldsymbol{x}\|^{2}\big{\}}\,.$ (5) $$

In words, we are given a Gaussian observation $\boldsymbol{y}_{t}\sim{\mathsf{N}}(t\boldsymbol{x},t{\boldsymbol{I}}_{d})$
(for a single $t$) and want to estimate $\boldsymbol{x}$ as to minimize
mean square error (MSE). This is also known as the ‘score-matching objective’.
Minimization is over all measurable functions ${\boldsymbol{\varphi}}$.

In practice the minimization in Eq. ([5](https://arxiv.org/html/2503.08028v2#S1.E5)) has to be modified
for two reasons:

*First,* in general we do not know the distribution of $\boldsymbol{x}$ over which
the expectation in Eq. ([5](https://arxiv.org/html/2503.08028v2#S1.E5)) is taken, and instead
we are given a finite sample $(\boldsymbol{x}_{i})_{i\leq N}\sim_{iid}\mu$. We thus replace
the MSE by its empirical version:

$$ minimize $\displaystyle\;\;\;\frac{1}{N}\sum_{i=1}^{N}\big{\|}{\boldsymbol{\varphi}}( \boldsymbol{y}_{i,t})-\boldsymbol{x}_{i}\big{\|}^{2}\,,\;\;\;\;\mbox{subj. to} \;\;{\boldsymbol{\varphi}}\in\mathscrsfs{N}\,,$ (6) $$

where $\boldsymbol{y}_{i,t}=t\boldsymbol{x}_{i}+\sqrt{t}\boldsymbol{g}_{i}$ for i.i.d. standard Gaussian vectors $\boldsymbol{g}_{i}$.
Note that, for statistical reasons, the minimization over all functions must be replaced by by minimization
over a suitable function class $\mathscrsfs{N}$ (e.g. neural nets). Changing the population problem
([5](https://arxiv.org/html/2503.08028v2#S1.E5)) into the empirical risk minimization problem of Eq. ([6](https://arxiv.org/html/2503.08028v2#S1.E6))
results into a statistical error that (ideally) vanishes as $N\to\infty$.

*Second,* the evaluation
$\boldsymbol{y}\mapsto\hat{\boldsymbol{m}}(\boldsymbol{y},t)$ (a.k.a. ‘inference’ in deep learning jargon)
must admit a polytime algorithm,
in order to
be able to implement the generative process ([4](https://arxiv.org/html/2503.08028v2#S1.E4)).
For this reason, $\mathscrsfs{N}$ must be a set of functions that can be evaluated in polynomial time.
This is a purely computational constraint, and is present even if we have access to
the population distribution (i.e., for $N=\infty$).

It is natural to suspect that the second constraint might lead to incorrect
DS if near Bayes-optimal denoising
cannot be achieved in polynomial time. The recent literature on ‘information-computation gaps’
provides a wealth of examples of statistical estimation problems for which
evidence exists of a gap between Bayes optimal error and the optimal error achievable in polynomial time
(Brennan et al., 2018; Bandeira et al., 2022; Celentano and Montanari, 2022; Schramm and Wein, 2022).
Several of these problems take the form of estimating $\boldsymbol{x}\sim\mu$ from Gaussian observations
$\boldsymbol{y}_{t}=t\boldsymbol{x}+\sqrt{t}\,\boldsymbol{g}$.

Koehler and Vuong (2024) already pointed out informally that
denoising problems presenting an information-computation gap can result into a failure of DS.
As a concrete example, they suggested the sparse spike Wigner model (see below).
While this informal remark is natural, making it mathematically precise is far from obvious.
In fact –strictly speaking– the remark is false. Roughly speaking, if sampling
from $\mu$ is easy, then the drift $\boldsymbol{m}(\boldsymbol{y},t)$ can be constructed to return (for all $t\geq t_{0}$)
a fixed random sample $\boldsymbol{x}\sim\mu$. Then the diffusion will sample correctly.
However such $\boldsymbol{m}$ will be very far from an optimal denoiser. (See below
for formal version of this argument.)

We also note that several earlier paper pointed at gaps between
general sampling and DS or sequential sampling (Montanari et al., 2007; Ricci-Tersenghi and Semerjian, 2009; Ghio et al., 2024; Huang et al., 2024).
None of these papers presented a formal claim either.

The present work fills this gap by providing mathematically rigorous
versions of this claim in the setting of the sparse spiked Wigner model.
We prove that there exists denoisers that are approximate optimizers (or sometimes optimizers)
of the score matching objective ([6](https://arxiv.org/html/2503.08028v2#S1.E6))
among polynomial time algorithms (up to an exponentially small error)
and yet lead to completely incorrect sampling.

We emphasize that this failure of DS is computational of nature and
purely related to the requirement to approximate the Bayes optimal denoiser $\boldsymbol{m}(\boldsymbol{y},t)$
by a polytime computable function.

### 1.2 Summary of results

The construction of the distribution $\mu$
is fairly simple. Let $B_{n,k}:=\{\boldsymbol{u}\in\{0,\pm 1/\sqrt{k}\}^{n}|\|\boldsymbol{u}\|_{0}=k\}$ be the set of $0/\pm(1/\sqrt{k})$ unit vectors with $k$ nonzero entries ($\|\boldsymbol{u}\|_{0}$ denotes the number of nonzeros in $\boldsymbol{u}$). We define the target distribution $\mu=\mu_{n,k}$ to be the distribution of

$$ $\displaystyle\boldsymbol{x}=\boldsymbol{u}\boldsymbol{u}^{\sf T},\;\;\;\;\;\; \;\boldsymbol{u}\sim\operatorname{Unif}(B_{n,k}).$ (7) $$

Note that $\boldsymbol{x}\in{\mathbb{R}}^{n\times n}$ is a matrix that we identify with a vector in ${\mathbb{R}}^{d}$ for $d=n^{2}$.

Sampling from the target distribution $\mu$ is trivial. Just sample a vector
with entries in $\{0,1/\sqrt{k},-1/\sqrt{k}\}$ and exactly $k$ non-zero entries, and let $\boldsymbol{x}=\boldsymbol{u}\boldsymbol{u}^{\sf T}$. However, DS requires to approximate the denoiser

$$ $\displaystyle\boldsymbol{m}(\boldsymbol{y},t):=\mathbb{E}\big{\{}\boldsymbol{u }\boldsymbol{u}^{\sf T}\Big{|}\boldsymbol{y}_{t}=\boldsymbol{y}\big{\}},\;\;\; \;\boldsymbol{y}_{t}=t\boldsymbol{u}\boldsymbol{u}^{\sf T}+\sqrt{t}\boldsymbol {g},$ (8) $$

where $\boldsymbol{g}\sim{\mathsf{N}}({\boldsymbol{0}},{\boldsymbol{I}}_{n\times n})$. In other words,
it requires to estimate a sparse rank-one matrix from a single
observation in Gaussian noise, sometimes referred to ‘sparse spiked Wigner,’
and
closely related to the hidden clique problem or sparse PCA.
Substantial rigorous evidence supports the claim that —for certain
scalings of $k$, $t$ with $n$— polynomial-time algorithms cannot
approach the Bayes-optimal error (Butucea et al., 2015; Ma and Wu, 2015; Cai et al., 2017; Brennan et al., 2018; Schramm and Wein, 2022).

We will prove three sets of results:

1. Near optimizers of score-matching can sample incorrectly. We prove that there
exists $\hat{\boldsymbol{m}}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}\to{\mathbb{R}}
^{n\times n}$
such that:

- M1.
$\hat{\boldsymbol{m}}(\,\cdot\,)$ can be evaluated in polynomial time.
- M2.
The estimation error achieved by $\hat{\boldsymbol{m}}$
(namely, $\mathbb{E}\{\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\}$) is close to the conjectured optimal estimation error achieved by polynomial-time algorithms. Hence $\hat{\boldsymbol{m}}(\,\cdot\,,t)$ will be
a near minimizer of the score-matching objective ([5](https://arxiv.org/html/2503.08028v2#S1.E5)).
- M3.
Samples $\hat{\boldsymbol{x}}_{T}$ generated by the discretized diffusions ([4](https://arxiv.org/html/2503.08028v2#S1.E4)) with
drift approximation $\hat{\boldsymbol{m}}(\,\cdot\,,t)$
at some large time $T$ have distribution that is very far from the target
$\mu$.

2. Reduction from estimation to
DS. We prove that if accurate DS is possible in polynomial time,
then near Bayes optimal estimation of $\boldsymbol{x}$ from $\boldsymbol{y}_{t}=t\boldsymbol{x}+\sqrt{t}\boldsymbol{g}$
must also be possible in polynomial time for all $t$. The contrapositive of this statement implies
that if an information-computation gap exists, then diffusion sampling is impossible in polynomial time.

3. *All* Lipschitz score-matching optimizers sample incorrectly.
We prove that any denoiser that optimizes the score matching among polytime algorithms for $t\leq t_{1}$, and is $C/t$-Lipschitz for $t>t_{1}$ (for any constant $C$ a suitable $t_{1}$), samples incorrectly (under a modified version of $\mu_{n,k}$.)

The rest of the paper is organized as follows.
In Section [2](https://arxiv.org/html/2503.08028v2#S2) we motivate our setting and assumptions, and discuss some limitations of our results.
In Section [3](https://arxiv.org/html/2503.08028v2#S3) we
state formally our results
(for technical reason we state two separate results depending
on the growth of $k$ with $n$.)
Section [4](https://arxiv.org/html/2503.08028v2#S4) presents the general reduction from estimation to
diffusion sampling. Section [5](https://arxiv.org/html/2503.08028v2#S5) proves that all Lipschitz score matching optimizers fail.
A summary of notationss and proofs are presented in the appendices.

Throughout the paper, $C,C_{i},c_{i},\dots$ denote constants, whose value can change from line to line.

## 2 Discussion

##### Setting.

Our results indicate that
a standard application of denoising diffusions methodology
will fail to sample from $\mu_{n,k}$ despite the sampling problem in itself is easy.
On the other hand, the latent structure of the distribution
$\mu_{n,k}$ can be exploited to construct a better algorithm.
Namely, one can use diffusions to sample $\boldsymbol{u}\sim\operatorname{Unif}(B_{n,k})$ (this is easy since
the posterior expectation $\boldsymbol{m}(\boldsymbol{y},t)=\mathbb{E}[\boldsymbol{u}|t\boldsymbol{u}+
\sqrt{t}\boldsymbol{g}=\boldsymbol{y}]$ is easy to compute)
and then generate $\boldsymbol{x}=\boldsymbol{u}\boldsymbol{u}^{{\sf T}}$. However, identifying such latent structures from data
is (beyond the present example) a hard problem, both from a statistical and
from a computational viewpoint.

##### Limitations.

We prove that there exists drifts $\hat{\boldsymbol{m}}(\,\cdot\,,t)$ that lead
to poor generative sampling, despite being nearly optimal (among poly-time algorithms)
in terms of the score matching objective ([5](https://arxiv.org/html/2503.08028v2#S1.E5)).
In particular, these bad drifts will be near optimal solutions of the empirical risk minimization problem ([6](https://arxiv.org/html/2503.08028v2#S1.E6)), as long as $\mathscrsfs{N}$ pnly contains polytime methods and
is rich enough to approximate
them. We further exclude the existence of Lipschitz drifts $\hat{\boldsymbol{m}}(\,\cdot\,,t)$
that also satisfy conditions M1 and M2 but yield good generative sampling.

In principle there could still be non-Lipschitz polytime drifts that are
near score matching optimizers and sample well.
However if such drift existed, our results suggest that minimizing the score matching objective is not the right approach to find them (since the difference in value with bad drifts will be
superpolynomially small).

##### Correct samplers violating M2 .

As mentioned abovre,
of we drop condition M2, i.e. we accept drifts that are bad
from the point of view of the score matching, then it is possible to construct drifts that
can be evaluated in polynomial time and yield good sampling. However,
these are constructions that completely ignore the score-matching prescription ([2](https://arxiv.org/html/2503.08028v2#S1.E2)). This is stated formally bnelow and proved in Appendix [E](https://arxiv.org/html/2503.08028v2#A5).

###### Proposition 2.1 .

Suppose that a discretized SDE per Equation ([4](https://arxiv.org/html/2503.08028v2#S1.E4)) is generated, with step size $\Delta>0$ and noise stream $\hat{\boldsymbol{z}}_{t}\overset{i.i.d.}{\sim}{\mathsf{N}}(0,\boldsymbol{I}_{n
\times n})$. Then for every $n,k$, there exists a function $\hat{\boldsymbol{m}}(\boldsymbol{y},t)=\hat{\boldsymbol{m}}(\boldsymbol{y},t;
\hat{\boldsymbol{z}}_{1})$ parametrized by $\hat{\boldsymbol{z}}_{1}$ (with no additional randomness) such that:
$(i)$ $\mathbb{E}[\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}]=
2(1-o(1))$ uniformly for every $t\geq 0$ (sub-optimal score-matching error);
$(ii)$ $W_{1}(\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{\ell\Delta},\ell\Delta),
\boldsymbol{x})=0$ for all $\ell\geq 0$ ($\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{t},t)$ is an approximate sample of $\boldsymbol{x}$ at every point);
$(iii)$ $\lim_{\ell\to\infty}W_{1}(\hat{\boldsymbol{y}}_{\ell\Delta}/(\ell\Delta),
\boldsymbol{x})=0$ ($\hat{\boldsymbol{y}}_{t}/t$ is an approximate sample of $\boldsymbol{x}$ at large time).

Again, the drift constructed in this proposition
is artificial and has very poor value of the score-matching objective.

###### Proposition 2.1 .

##### Hardness assumption.

Condition M2 refers to the *conjectured* optimal estimation error
achieved by polynomial-time algorithms because lower bounds on this error do necessarily rely either
on complexity-theoretic assumptions or on optimality or certain restricted classes of algorithms.

##### Further related work.

A number of groups proved positive results on diffusion sampling.
Alaoui et al. (2022); Chen et al. (2023b); Montanari and Wu (2023); Lee et al. (2023); Benton et al. (2023) provide reductions
from diffusion sampling to score estimation.
Chen et al. (2023a); Shah et al. (2023); Mei and Wu (2025); Li et al. (2024)
give end-to-end guarantees for classes of distributions $\mu$.

The computational bottleneck that we study here has been
observed before in the context of certain Gibbs measures and Bayes posterior distributions Ghio et al. (2024); Alaoui et al. (2023); Huang et al. (2024),
and random constraint satisfaction problems Montanari et al. (2007); Ricci-Tersenghi and Semerjian (2009) (the latter papers use however
sequential sampling rather than diffusion sampling).

Our work provides a simpler example and a rigorous formalization of
the latter line of work. (Notice that despite superficial similarities,
our construction is different from the examples of
Ghio et al. (2024); Alaoui et al. (2023); Huang et al. (2024), and allows us to build existing hardness results.)

## 3 Main results

We will state two separate results for the probability distribution $\mu=\mu_{n,k}$
described in the introduction,
depending on the scaling of $k$ with $n$:
in Section [3.1](https://arxiv.org/html/2503.08028v2#S3.SS1) we assume $\sqrt{n}\ll k\ll n$;
while in Section [3.3](https://arxiv.org/html/2503.08028v2#S3.SS3) we assume $k\ll\sqrt{n}$.
In fact, the nature of the problem
changes at the threshold $k\asymp\sqrt{n}$.

In what follows $(\boldsymbol{x},\boldsymbol{y}_{t})$ will always be distributed according to the ideal
diffusion process of Eq. ([3](https://arxiv.org/html/2503.08028v2#S1.E3)), which also satisfies Eq. ([2](https://arxiv.org/html/2503.08028v2#S1.E2)),
and $\hat{\boldsymbol{y}}_{t}$ will denote the process generated with the implemented procedure
([4](https://arxiv.org/html/2503.08028v2#S1.E4)). A crucial role will be played by the
following algorithmic threshold

$$ $\displaystyle t_{\mbox{\tiny\rm alg}}(n,k):=\begin{cases}k^{2}\log\left(\dfrac {n}{k^{2}}\right)&\text{ if $k\ll\sqrt{n}$}\\ \dfrac{n}{2}&\text{ if $\sqrt{n}\ll k\ll n$}\end{cases}$ (9) $$

It is expected that no polynomial-time algorithm can estimate $\boldsymbol{x}$ significantly better than the null estimator $\hat{\boldsymbol{m}}_{\mbox{\tiny\rm null}}=\boldsymbol{0}$ for $\underline{k}_{n}\leq k\ll n$, with $\underline{k}_{n}\ll n$ a suitable sequence
(see Conjecture [3.2](https://arxiv.org/html/2503.08028v2#S3.Thmclaim2)).

Since $\|\boldsymbol{x}\|_{F}=1$ for $\boldsymbol{x}\sim\mu$, the Bayes denoiser
$\boldsymbol{m}(\boldsymbol{y},t)=\boldsymbol{m}(\boldsymbol{y})$
does not depend on $t$.

### 3.1 Moderately sparse regime: n ≪ k ≪ n much-less-than 𝑛 𝑘 much-less-than 𝑛 \sqrt{n}\ll k\ll n square-root start_ARG italic_n end_ARG ≪ italic_k ≪ italic_n

Given an arbitrary polynomial time computable drift $\hat{\boldsymbol{m}}_{0}$,
we construct a different poly-time drift $\hat{\boldsymbol{m}}$, with nearly equal score matching objective and yet incorrect sampling. In the next subsection we state a related result that establishes incorrect sampling from an explicit denoiser with near-optimal denoising
properties.

###### Assumption 1 (Small norm below threshold) .

Let $\boldsymbol{y}_{t}=t\boldsymbol{x}+\boldsymbol{B}_{t}$, and assume $\sqrt{n}\ll k\ll n$.
Then, there exists (known) constants $\gamma,\varepsilon\in(0,1)$, independent of $n$
such that, for every $D>0$,

$$ $\displaystyle\int_{0}^{(1-\gamma)t_{\mbox{\tiny\rm alg}}}\mathbb{P}\big{(}\| \hat{\boldsymbol{m}}_{0}(\boldsymbol{y}_{t},t)\|\geq 1-\varepsilon\big{)}\,{ \rm d}t=O(n^{-D})\,.$ $$

We expect this condition to be satisfied by any reasonable polynomial time algorithm,
under the well-accepted
Conjecture [3.2](https://arxiv.org/html/2503.08028v2#S3.Thmclaim2) below on information-computation gaps.
Indeed, the conjecture states that
$\inf_{0\leq t\leq(1-\gamma)t_{\mbox{\tiny\rm alg}}}\mathbb{E}\left[\|\hat{
\boldsymbol{m}}_{0}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\right]\geq 1-o_
{n}(1)$. A simple calculation shows that
any polytime estimator matching this error must also satisfy
$\mathbb{E}\{\|\hat{\boldsymbol{m}}_{0}(\boldsymbol{y}_{t},t)\|^{2}\}=o_{n}(1)$ (because otherwise $c_{t}\hat{\boldsymbol{m}}_{0}(\,\cdot\,)$ would violate the lower bound for a
suitable deterministic $c_{t}$).

###### Theorem 1 .

Assume $\sqrt{n}\ll k\ll n$, and note that in this case $t_{\mbox{\tiny\rm alg}}(n,k):=n/2$
per Eq. ([9](https://arxiv.org/html/2503.08028v2#S3.E9)).
Let $\hat{\boldsymbol{m}}_{0}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}\to{\mathbb
{R}}^{n\times n}$ be an arbitrary poly-time algorithm such that $\sup_{\boldsymbol{y},t}\|\hat{\boldsymbol{m}}_{0}(\boldsymbol{y},t)\|_{F}\leq 1$ and Assumption
[1](https://arxiv.org/html/2503.08028v2#Thmassumption1) holds. Then there exists an estimator
$\hat{\boldsymbol{m}}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}\to{\mathbb{R}}
^{n\times n}$ such that:

- M1.
$\hat{\boldsymbol{m}}(\,\cdot\,)$ can be evaluated in polynomial time.
- M2.
If $\boldsymbol{y}_{t}=t\boldsymbol{x}+\boldsymbol{B}_{t}$ is the true diffusion (equivalently given by Eq. ([1](https://arxiv.org/html/2503.08028v2#S1.E1))), then, for every $D>0$,
$\displaystyle\int_{0}^{\infty}\mathbb{E}[\|\hat{\boldsymbol{m}}(\boldsymbol{y}
_{t},t)-\hat{\boldsymbol{m}}_{0}(\boldsymbol{y}_{t},t)\|^{2}]\,{\rm d}t=O(n^{-
D})\,.$
- M3.
There exists $\delta=o_{n}(1)$ such that, for
any step size $\Delta=\Delta_{n}>0$, we have
$\displaystyle\inf_{t\in\mathbb{N}\cdot\Delta,t\geq(1+\delta)t_{\mbox{\tiny\rm
alg
}}}W_{1}(\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{t},t),\boldsymbol{x})\geq 1
-o_{n}(1)\,.$
(10)

In order to connect this theorem with the discussion in the introduction,
we recall two facts from the literature on submatrix estimation:
$(i)$ The ideal posterior mean estimator $\boldsymbol{m}(\boldsymbol{y}_{t})$ achieves small mean square error
in a large interval below $t_{\mbox{\tiny\rm alg}}$ (Proposition [3.1](https://arxiv.org/html/2503.08028v2#S3.Thmclaim1)); $(ii)$ No polynomial-time estimator is expected to
perform substantially better than the null estimator below $t_{\mbox{\tiny\rm alg}}$ (Conjecture [3.2](https://arxiv.org/html/2503.08028v2#S3.Thmclaim2)).

Regarding $(i)$, we state a characterization of the Bayes optimal error.
The proof is essentially the same (indeed simpler) as the main result in Butucea et al. (2015),
which considers the case of asymmetric matrices. An alternative proof
(for $k\leq n^{a}$, $a<5/6$) also follows from Barbier et al. (2020).

###### Proposition 3.1 (Modification of Butucea et al. ( 2015 ) ) .

Let $\boldsymbol{m}(\,\cdot\,)$ be the posterior mean estimator
in Eq. ([8](https://arxiv.org/html/2503.08028v2#S1.E8)). Assume $1\ll k\ll n$,
and define $t_{\mbox{\tiny\rm Bayes}}(n,k):=2k\log(n/k)$.
Then, for any $\delta>0$, we have

$$ $\displaystyle\inf_{t\leq(1-\delta)t_{\mbox{\tiny\rm Bayes}}}\mathbb{E}\big{\{} \|\boldsymbol{m}(\boldsymbol{y}_{t})-\boldsymbol{x}\|^{2}\}=1-o_{n}(1)\,,\;\; \;\;\;\;\sup_{t\geq(1+\delta)t_{\mbox{\tiny\rm Bayes}}}\mathbb{E}\big{\{}\| \boldsymbol{m}(\boldsymbol{y}_{t})-\boldsymbol{x}\|^{2}\}=o_{n}(1)\,.$ $$

Several groups have proved hardness results for estimating $\boldsymbol{u}$ in the model
considered here, for $t$ sufficiently smaller than the spectral threshold $n/2$.
In particular Brennan et al. (2018) establishes hardness conditional on
hardness on the planted clique problem for $t\leq n^{1-\varepsilon}$, $\varepsilon>0$,
Hopkins et al. (2017) proved failure of sum-of-squares relaxations also for
$t\leq n^{1-\varepsilon}$, Schramm and Wein (2022) proved
hardness for low-degree polynomials for the related hidden submatrix problem for $t\leq cn$, $c>0$ a sufficiently small constant. Finally,
Kunisky et al. (2019) prove hardness for low-degree polynomials algorithms when $t\leq(1-\varepsilon)n/2$, albeit
their results require $k$ to be of order $n$.

###### Conjecture 3.2 .

There exists $\underline{k}_{n}\ll n$ such that the following holds for any
$k=k_{n}$, with $\underline{k}_{n}\leq k_{n}\ll n$. Let $\{\hat{\boldsymbol{m}}_{n}\}_{n\geq 1}$,
$\hat{\boldsymbol{m}}_{n}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}\to{\mathbb
{R}}^{n\times n}$ be any sequence of polynomial time algorithms indexed by $n$ (namely $\hat{\boldsymbol{m}}_{n}$ can be evaluated in time bounded by $n^{C}$ for some constant $C$). Then for any $\delta>0$, we have

$$ $\displaystyle\inf_{t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}}\mathbb{E}\big{\{}\| \hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\}$ $\displaystyle\geq 1-o_{n}(1)\,.$ (11) $$

Summarizing, in the moderately sparse regime,
for $2k\log(n/k)\ll t\ll n$ the Bayes optimal estimator can
estimate the signal $\boldsymbol{x}$ accurately, but we expect that no polynomial-time algorithm can achieve the same.

It makes sense to assume that $\|\hat{\boldsymbol{m}}_{0}(\cdot,\cdot)\|_{F}\leq 1$. Since $\boldsymbol{x}\sim\mu_{n,k}$
and $\text{supp}(\mu_{n,k})\subset\{\boldsymbol{a}:\|\boldsymbol{a}\|_{F}\leq 1\}$ and the latter is a convex set, projecting any estimator
$\hat{\boldsymbol{m}}_{0}$ onto this set yields a smaller MSE.

###### Assumption 1 (Small norm below threshold) .

###### Theorem 1 .

###### Proposition 3.1 (Modification of Butucea et al. ( 2015 ) ) .

###### Conjecture 3.2 .

### 3.2 Moderately sparse regime: A concrete example

Figure: Algorithm 1 Submatrix Estimation Algorithm (moderately sparse regime)

Since Theorem [1](https://arxiv.org/html/2503.08028v2#Thmtheorem1) is somewhat abstract, we complement it
with an explicit example. Namely, we prove that a modification of a standard spectral estimator,
while achieving near optimal estimation error (among polytime algorithms) fails to generate
samples from the correct distribition.

###### Theorem 2 .

Assume $\sqrt{n}\ll k\ll n$, and note that in this case $t_{\mbox{\tiny\rm alg}}(n,k):=n/2$
per Eq. ([9](https://arxiv.org/html/2503.08028v2#S3.E9)).
Then the estimator $\hat{\boldsymbol{m}}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}\to{\mathbb{R}}
^{n\times n}$
defined in Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1) satisfies the following:

- M1.
$\hat{\boldsymbol{m}}(\,\cdot\,)$ can be evaluated in polynomial time.
- M2.
For any $\delta>0$, there exists $c=c(\delta)$, $C=C(\delta)$ such that
$\displaystyle\phantom{A}\inf_{t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}}\mathbb{E
}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\big
{\}}=1-o_{n}(1)\,,\;\;\;\;\sup_{t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}}\mathbb
{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}
\big{\}}\leq C\,e^{-cn/k}\,.$
- M3.
For any $\Delta>0$,
$\displaystyle\inf_{t\in{\mathbb{N}}\cdot\Delta}W_{1}(\hat{\boldsymbol{m}}(\hat
{\boldsymbol{y}}_{t},t),\boldsymbol{x})=1-o_{n}(1)\,.$

As mentioned, Algorithm
[1](https://arxiv.org/html/2503.08028v2#alg1) uses a spectral approach.
We compute the leading
eigenvector of (the symmetrized version of) $\boldsymbol{y}_{t}$, call it $\boldsymbol{v}_{t}\in{\mathbb{R}}^{n}$.
We then estimate the support $S$
of the latent rank-one matrix using the largest entries of $\boldsymbol{v}_{t}$.

###### Theorem 2 .

### 3.3 Very sparse regime: k ≪ n much-less-than 𝑘 𝑛 k\ll\sqrt{n} italic_k ≪ square-root start_ARG italic_n end_ARG

Figure: Algorithm 2 Submatrix Estimation Algorithm (very sparse regime)

###### Theorem 3 .

Assume $(\log n)^{5/2}\lesssim k\ll\sqrt{n}$, and note that in this case $t_{\mbox{\tiny\rm alg}}(n,k)=k^{2}\log(n/k^{2})$,
per Eq. ([9](https://arxiv.org/html/2503.08028v2#S3.E9)).
Then the randomized estimeator $\hat{\boldsymbol{m}}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}\to{\mathbb{R}}
^{n\times n}$
of Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2) satisfies the following:

- M1.
$\hat{\boldsymbol{m}}(\,\cdot\,)$ can be evaluated in polynomial time.
- M2.
For any $\delta>0$, $D>0$, there exists $C(\delta,D)>0$ such that:
$\displaystyle\phantom{A}\inf_{t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}}\mathbb{E
}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\big
{\}}=1-o_{n}(1)\,,\;\;\;\sup_{t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}}\mathbb{E
}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\big
{\}}\leq C(\delta,D)n^{-D}\,.$
- M3.
For any $\Delta>0$,
$\displaystyle\inf_{t\in{\mathbb{N}}\cdot\Delta}W_{1}(\hat{\boldsymbol{m}}(\hat
{\boldsymbol{y}}_{t},t),\boldsymbol{x})=1-o_{n}(1)\,.$

The pseudocode for the estimator $\hat{\boldsymbol{m}}(\,\cdot\,)$ that
is constructed in the above is given as Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2).
This is based on a standard approach in the literature Deshpande and Montanari (2016); Cai et al. (2017),
with some modifications
to allow for its analysis in the diffusion setting.
The main steps are as follows:
$(1)$ Perform Gaussian data splitting of $\boldsymbol{y}_{t}$ into $\boldsymbol{y}_{t,+}$, $\boldsymbol{y}_{t,-}$, see Line 3 of
lgorithm [2](https://arxiv.org/html/2503.08028v2#alg2), with most of the information preserved in $\boldsymbol{y}_{t,+}$.
$(2)$ Use entrywise soft thresholding $\eta_{s}(x)=(|x|-s)_{+}\operatorname{sign}(x)$ to reduce the noise in
the symmetrized version of $\boldsymbol{y}_{t,+}$.
$(3)$ Compute a first estimate of the latent vector $\boldsymbol{1}_{S}$ by the principal eigenvector
of the above matrix.
$(4)$ Refine this estimate using the remaining information $\boldsymbol{y}_{t,-}$.

We point out that Proposition [3.1](https://arxiv.org/html/2503.08028v2#S3.Thmclaim1) remains true in the regime
$\sqrt{n}\ll k\ll n$, and hence we observe a gap between $t_{\mbox{\tiny\rm alg}}(n,k)$ and $t_{\mbox{\tiny\rm Bayes}}(n,k)$
in this regime as well.

###### Theorem 3 .

## 4 Reduction of estimation to diffusion-based sampling

In order to clarify our conclusions,
we prove a general reduction: if diffusion sampling can be performed in polynomial time with sufficient accuracy, then
we can sample in polynomial time from the posterior distribution of
$\boldsymbol{x}\sim\mu$ given noisy information $\boldsymbol{y}=\boldsymbol{x}+\sigma\boldsymbol{g}$ (where $\boldsymbol{g}$ is standard Gaussian noise independent of $\boldsymbol{x}$). By repeating this sampling, we can also perform estimation.

Here we consider a general sampling problem for $\boldsymbol{x}\sim\mu$ a probability distribution in ${\mathbb{R}}^{d}$.
For the sake of simplicity, we will assume $\mu$ to be supported on the unit sphere ${\mathbb{S}}^{d-1}=\{\boldsymbol{x}:\|\boldsymbol{x}\|=1\}$.
We denote by ${\rm P}_{\boldsymbol{y}}^{T}$ the law of $(\boldsymbol{y}_{t})_{0\leq t\leq T}$ where $\boldsymbol{y}_{t}$ given by Eq. ([1](https://arxiv.org/html/2503.08028v2#S1.E1))
and by ${\rm P}_{\hat{\boldsymbol{y}}}^{T,\Delta}$ the law of $(\hat{\boldsymbol{y}}_{t})_{0\leq t\leq T}$, which is the discretized diffusion trajectory defined in ([4](https://arxiv.org/html/2503.08028v2#S1.E4)) (interpolated linearly outside ${\mathbb{N}}\cdot\Delta$).

It is further useful to define $\overline{\rm P}_{\hat{\boldsymbol{y}}}^{T,\Delta}$ to be the law of the the solution following SDE

$$ $\displaystyle{\rm d}\hat{\boldsymbol{y}}_{t}=\hat{\boldsymbol{m}}(\hat{ \boldsymbol{y}}_{\lfloor t\rfloor_{\Delta}},\lfloor t\rfloor_{\Delta})\,{\rm d }t+{\rm d}\boldsymbol{B}_{t}\,,$ (12) $$

where $\lfloor t\rfloor_{\Delta}:=\max\{s\in{\mathbb{N}}\cdot\Delta:s\leq t\}$. This definition interpolates $(\hat{\boldsymbol{y}}_{t}:t\in{\mathbb{N}}\cdot\Delta)$
defined in ([4](https://arxiv.org/html/2503.08028v2#S1.E4)).

###### Theorem 4 .

Assume that $\hat{\boldsymbol{m}}(\,\cdot\,,\,\cdot\,)$ has complexity $\chi$
and that for any $T\leq\theta d$,
$D_{\mbox{\tiny\rm KL}}(\overline{\rm P}_{\hat{\boldsymbol{y}}}^{T,\Delta}\|{
\rm P}_{\boldsymbol{y}}^{T})\leq\varepsilon$
(where $\overline{\rm P}_{\hat{\boldsymbol{y}}}^{T,\Delta}$ is the continuous time process obtained by Brownian-linear interpolation of Eq. ([4](https://arxiv.org/html/2503.08028v2#S1.E4))).

Then for any $\sigma>0$
there exists an algorithm with complexity $O(\chi\cdot T/\Delta)$, that takes as input
$\boldsymbol{y}=\boldsymbol{x}+\sigma\boldsymbol{g}$, $(\boldsymbol{x},\boldsymbol{g})\sim\mu\otimes{\mathsf{N}}(0,{\boldsymbol{I}})$, and outputs
$\hat{\boldsymbol{x}}$, such that

$$ $\displaystyle\mathbb{E}\|{\rm P}_{\boldsymbol{x}|\boldsymbol{y}}-{\rm P}_{\hat {\boldsymbol{x}}|\boldsymbol{y}}\|_{\mbox{\tiny\rm TV}}\leq\sqrt{2\varepsilon} +\varepsilon_{0}(\theta)=:\overline{\varepsilon}\,,$ (13) $$

where $\varepsilon_{0}(\theta):=\mathbb{E}\|{\rm P}_{\boldsymbol{x}|\boldsymbol{y}}-{
\mathsf{N}}(\boldsymbol{0},(\theta d)^{-1}\boldsymbol{I}_{d})*{\rm P}_{
\boldsymbol{x}|\boldsymbol{y}}\|_{\mbox{\tiny\rm TV}}$
is the expected TV distance between ${\rm P}_{\boldsymbol{x}|\boldsymbol{y}}$ and the convolution of ${\rm P}_{\boldsymbol{x}|\boldsymbol{y}}$
with a Gaussian with variance $1/(\theta d)$.
As a consequence, there exists a randomized algorithm $\hat{\boldsymbol{m}}_{+}$ with complexity $(N\chi\cdot T/\Delta)$
that approximates the posterior expectation:

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}_{+}(\boldsymbol{y})- \mathbb{E}(\boldsymbol{x}|\boldsymbol{y})\|^{2}\big{\}}\leq 2\overline{ \varepsilon}+2N^{-1}\,.$ (14) $$

The next statement makes a weaker assumption on the accuracy of the diffusion sampler (transportation instead of KL
distance), but in exchnage assumes the approximate drift $\hat{\boldsymbol{m}}$ to be Lipschitz.
We note that ${\rm Lip}(\boldsymbol{m}(\,\cdot\,,t))=\sup_{\boldsymbol{y}}\|{\rm Cov}(
\boldsymbol{x}|\boldsymbol{y}_{t}=\boldsymbol{y})\|_{{\mbox{\rm\tiny op}}}$,
and the latter is of $O(1/d)$ (for instance) if the coordinates of $\boldsymbol{x}$ are weakly dependent under
the posterior.

###### Theorem 5 .

Assume that $\hat{\boldsymbol{m}}(\,\cdot\,,\,\cdot\,)$ has computational complexity $\chi$ and satisfies the following: $(a)$ For every $t\geq 1/\sigma^{2}$, $\boldsymbol{y}\mapsto\hat{\boldsymbol{m}}(\boldsymbol{y},t)$ is $L/d$-Lipschitz.
$(b)$ There is a stepsize $\Delta$ such that $W_{1}({\rm P}_{\hat{\boldsymbol{y}}}^{T,\Delta},{\rm P}_{\boldsymbol{y}}^{T})\leq\varepsilon$
for any $T\leq\theta d$.

Then for any $\sigma>0$
there exists an algorithm with complexity $O(\chi\cdot T/\Delta)$, that takes as input
$\boldsymbol{y}=\boldsymbol{x}+\sigma\boldsymbol{g}$, $(\boldsymbol{x},\boldsymbol{g})\sim\mu\otimes{\mathsf{N}}(0,{\boldsymbol{I}})$, and outputs
$\hat{\boldsymbol{x}}$, such that

$$ $\displaystyle\mathbb{E}_{\boldsymbol{y}}W_{1}({\rm P}_{\boldsymbol{x}| \boldsymbol{y}},{\rm P}_{\hat{\boldsymbol{x}}|\boldsymbol{y}})\leq 2e^{\theta L }\varepsilon+\frac{1}{\sqrt{\theta}}=:\overline{\varepsilon}\,.$ (15) $$

As a consequence, Eq. ([14](https://arxiv.org/html/2503.08028v2#S4.E14)) holds also in this case with the new definition of $\overline{\varepsilon}$.

###### Theorem 4 .

###### Theorem 5 .

## 5 All Lipschitz polytime algorithms fail

In Section [3](https://arxiv.org/html/2503.08028v2#S3) (Theorems [1](https://arxiv.org/html/2503.08028v2#Thmtheorem1), [2](https://arxiv.org/html/2503.08028v2#Thmtheorem2)
and [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3))
we proved that there exist near-optimizers of the
score matching objective that perform poorly. This means that, unless we optimize
this objective to very high precision (and the sample size is large enough to approximate
the population objective with similarly high precision), DS will fail.

However, we did not rule out the possibility that the optimal polytime
drift $\hat{\boldsymbol{m}}$ will perform well.

We next show that this is not the case, under an additional assumption,
namely that the drift $\hat{\boldsymbol{m}}(\,\cdot\,;t)$ is Lipschitz continuous for $t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}(n,k)$, where $t_{\mbox{\tiny\rm alg}}(n,k)$ is defined as in Eq. ([9](https://arxiv.org/html/2503.08028v2#S3.E9)).
We do so under a modification of the target distribution.
Letting $\mu^{0}_{n,k}$ be the centered version of $\mu_{n,k}$ we define

$$ $\displaystyle\overline{\mu}_{n,k}=\frac{1}{2}\,\delta_{\boldsymbol{0}}+\frac{1 }{2}\,\mu^{0}_{n,k}\,.$ (16) $$

In words, with probability $1/2$ we let $\boldsymbol{x}=\boldsymbol{0}$ and with probability $1/2$ we
draw $\boldsymbol{x}=\tilde{\boldsymbol{x}}-\mathbb{E}[\tilde{\boldsymbol{x}}]$, $\tilde{\boldsymbol{x}}\sim\mu_{n,k}$, a sparse rank-one matrix, as in previous sections.

###### Theorem 6 .

Assume $\underline{k}_{n}$ exists as in Conjecture [3.2](https://arxiv.org/html/2503.08028v2#S3.Thmclaim2), and let $k=k_{n}$
be such that $\underline{k}_{n}\vee\sqrt{n}\leq k_{n}\ll n$ (moderately sparse regime).
Let $\hat{\boldsymbol{m}}_{n}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}_{\geq 0}
\to{\mathbb{R}}^{n\times n}$ be
a polynomial time denoiser for such that the following holds for some $\delta=\delta_{n}\leq n^{-c_{*}}$ with $c_{*}>0$:

- 1.
$\hat{\boldsymbol{m}}_{n}$ is nearly optimal, namely (for $\boldsymbol{y}_{t}=t\boldsymbol{x}+\boldsymbol{W}_{t}$, $\boldsymbol{W}_{t}$ standard Brownian motion)
$\displaystyle\sup_{t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}}\Big{|}\mathbb{E}
\big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}
\big{\}}-\frac{1}{2}\Big{|}=o_{n}(n^{-1})\,,$
(17)
$\displaystyle\sup_{t\leq(1+\delta)t_{\mbox{\tiny\rm alg}}}\mathbb{E}\big{\{}\|
\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\big{\}}=o_
{n}(1)\,,$
(18)
and further, for any $t\geq 0$ the MSE of $\hat{\boldsymbol{m}}_{n}$ smaller or equal than the MSE
of $c(\lambda_{1}(\boldsymbol{y}_{t}))\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)$ for any function $c(\,\cdot\,)$
of the maximum eigenvalue of $(\boldsymbol{y}_{t}+\boldsymbol{y}_{T}^{{\sf T}})/\sqrt{2}$, and than the MSE of $\boldsymbol{P}_{{\sf B}}\hat{\boldsymbol{m}}_{n}$, for $\boldsymbol{P}_{{\sf B}}$
the projection onto the unit ball.
- 2.
$\hat{\boldsymbol{m}}_{n}(\,\cdot\,,t):{\mathbb{R}}^{n\times n}\to{\mathbb{R}}^
{n\times n}$ is $C/t$-Lipschitz for some constant $C$
and all $t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}$.

Further assume that for any polynomial time computable test
$\phi:{\mathbb{R}}^{n\times n}\to[0,1]$, $\sup_{t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}}\big{|}\mathbb{E}\{\phi(
\boldsymbol{y}_{t})\}-\mathbb{E}\{\phi(\boldsymbol{W}_{t})\}|\leq o_{n}(1)$.

Then, there exists a constant $\gamma_{0}>0$ such that and any $\Delta>0$

$$ $\displaystyle\inf_{t\in{\mathbb{N}}\cdot\Delta\cap[0,n^{1+\gamma_{0}}]}W_{1}( \hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{t},t),\boldsymbol{x})\geq\frac{1}{2 }-o_{n}(1)\,.$ (19) $$

We note that the error of polytime denoisers in Eq. ([17](https://arxiv.org/html/2503.08028v2#S5.E17))
(and correspondingly, the sampling error of Eq. ([19](https://arxiv.org/html/2503.08028v2#S5.E19)))
is $1/2$ instead of $1$ because the trivial denoiser (which always returns $\boldsymbol{0}$)
achieves error $1/2$.

Theorem [6](https://arxiv.org/html/2503.08028v2#Thmtheorem6) does not rule out the possiblility that
there exists a near-optimizer of score matching that violates the Lipschitz condition
and samples well. However, in for $t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}$ accurate estimation is possible with Lipschitz algorithms, and indeed many natural methods follow in this class (e.g. neural nets with bounded number of layers and suitable operator norm bounds on the weigths.)

###### Theorem 6 .

## 6 Numerical illustration

Figure: Figure 1: Generating sparse rank-one matrices $\boldsymbol{x}\sim\tilde{\mu}_{n,k}$ using denoising diffusions, for $n=350$, $k=20$. Left: MSE of various denoisers (vertical line corresponds to the algorithmic threshold $t_{\mbox{\tiny\rm alg}}$.) Right: Frobenius norms of generated samples.
Refer to caption: x1.png

The theory developed in the previous section yields a concrete prediction of the failure mode of DS when applied to the distribution $\tilde{\mu}_{n,k}=(1/2)\delta_{\boldsymbol{0}}+(1/2)\mu_{n,k}$ (with $\mu_{n,k}$ the law of $\boldsymbol{x}=\boldsymbol{u}\boldsymbol{u}^{{\sf T}}$, $\boldsymbol{u}\sim\operatorname{Unif}(B_{n,k})$). Namely
(for large $n$, and $\sqrt{n}\ll k\ll n$):

1. Given sufficient model complexity and training samples,
we expect the learnt denoiser $\hat{\boldsymbol{m}}_{n}(t,\,\cdot\,)$ to achieve MSE
close to $1/2$ for $t<(1-\delta)t_{\mbox{\tiny\rm alg}}$, and close to 0 0 for $t>(1+\delta)t_{\mbox{\tiny\rm alg}}$.

2. We expect DS based on such a denoiser to generate samples concentrated around $\boldsymbol{0}$.

We tested these predictions in a numerical experiment. We considered three polytime denoisers:

- $(a)$
The spectral-plus-projection denoiser of Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1);
- $(b)$
A modification of the latter whereby the eigenvector calculation is replaced by $25$ iterations of power method;
- $(c)$
A learned graph neural network (GNN)
(Scarselli et al., 2008; Kipf and Welling, 2016).

We carry out experiments with denoiser $(b)$ because
$\ell$ iterations of power method can be approximated by an $\ell$-layers
GNN. Hence, method $(b)$ provides a baseline for GNN denoisers.

Our GNN architecture uses node embeddings that are generated by $3$ iterations
of the power method and $10$ message passing layers.
Each message passing layer comprises of the ‘message’ and ‘node-update’
multi-layer perceptrons (MLPs), both of which are $2$-layer neural networks with LeakyReLU nonlinearity. We simply use the complete graph with self-loops for node embedding updates. We find that ‘seeding’ the node embeddings with iterations of
power method is crucial for effective training.

During training of the denoiser, we sample time points $t$ as follows: choose a time threshold $t_{\star}$, and sample so that times $t>t_{\star}$ are picked with total probability $0.95$ (and times $t\leq t_{\star}$ are picked with total probability $0.05$). Within each interval $(0,t_{\star}]$ and $(t_{\star},T)$, times are chosen at random. This allows the neural network to initially prioritize learning in a low-noise regime. Several fine-tuning steps are taken, for which $t_{\star}$ is gradually decreased to refine the network on lower SNR.

Empirically, training directly with $10$ layers is difficult, due to its depth. We find that training initially with $7$ layers, then subsequently introducing the later layers results in more stable training.

We train such a network using $N=10000$ samples $\boldsymbol{x}_{i}$ from the distribution $\overline{\mu}_{n,k}$,
and evaluate their MSE on $N_{\text{test}}=15000$ samples.

Figure [1](https://arxiv.org/html/2503.08028v2#S6.F1), left frame, reports the MSE achieved by the three denoisers
$(a)$, $(b)$, $(c)$ as a function of $t/n$, for $n=350$, $k=20$. As GNNs are permutation-equivariant, we are training on $\approx 1\%$ of all possible outcomes, for $n=350$ and $k=20$.
We observe that the GNN denoiser outperforms both the spectral algorithm and its approximation via power iteration. However, none of the three approaches can overcome the barrier at $t_{\mbox{\tiny\rm alg}}=n/2$, while they perform reasonably well above that threshold. This confirms the prediction at point 1 above.

On the right, we plot the histogram of Frobenius norms of samples generated
with the GNN denoiser. These values are close to 0 0, which confirms
the prediction at point 2 above. By using $\|\,\cdot\,\|_{F}$ as a $1$-Lipschitz test function, we obtain that the Wasserstein distance between diffusion samples and the target distribution is at least $0.48$ (the asymptotic prediction from theory is $0.50$).

## Acknowledgements

This work was supported by the NSF through award DMS-2031883, the Simons Foundation through
Award 814639 for the Collaboration on the Theoretical Foundations of Deep Learning,
and the ONR grant N00014-18-1-2729.

## References

- Alaoui et al. [2022]
Ahmed El Alaoui, Andrea Montanari, and Mark Sellke.
Sampling from the sherrington-kirkpatrick gibbs measure via
algorithmic stochastic localization.
In *2022 IEEE 63rd Annual Symposium on Foundations of Computer
Science (FOCS)*, pages 323–334. IEEE, 2022.
- Alaoui et al. [2023]
Ahmed El Alaoui, Andrea Montanari, and Mark Sellke.
Sampling from mean-field gibbs measures via diffusion processes.
*arXiv:2310.08912*, 2023.
- Baik et al. [2005]
Jinho Baik, Gérard Ben Arous, and Sandrine Péché.
Phase transition of the largest eigenvalue for nonnull complex sample
covariance matrices.
*Ann. Probab.*, 33(1):1643–1697, 2005.
- Bandeira and van Handel [2016]
Afonso S. Bandeira and Ramon van Handel.
Sharp nonasymptotic bounds on the norm of random matrices with
independent entries.
*The Annals of Probability*, 44(4):2479 –
2506, 2016.
doi: 10.1214/15-AOP1025.
URL [https://doi.org/10.1214/15-AOP1025](https://doi.org/10.1214/15-AOP1025).
- Bandeira et al. [2022]
Afonso S Bandeira, Ahmed El Alaoui, Samuel Hopkins, Tselil Schramm, Alexander S
Wein, and Ilias Zadik.
The franz-parisi criterion and computational trade-offs in high
dimensional statistics.
*Advances in Neural Information Processing Systems*,
35:33831–33844, 2022.
- Barbier et al. [2020]
Jean Barbier, Nicolas Macris, and Cynthia Rush.
All-or-nothing statistical and computational phase transitions in
sparse spiked matrix estimation.
*Advances in Neural Information Processing Systems*,
33:14915–14926, 2020.
- Benton et al. [2023]
Joe Benton, Valentin De Bortoli, Arnaud Doucet, and George Deligiannidis.
Linear convergence bounds for diffusion models via stochastic
localization.
*arXiv:2308.03686*, 2023.
- Brennan et al. [2018]
Matthew Brennan, Guy Bresler, and Wasim Huleihel.
Reducibility and computational lower bounds for problems with planted
sparse structure.
In *Conference On Learning Theory*, pages 48–166. PMLR, 2018.
- Butucea et al. [2015]
Cristina Butucea, Yuri I Ingster, and Irina A Suslina.
Sharp variable selection of a sparse submatrix in a high-dimensional
noisy matrix.
*ESAIM: Probability and Statistics*, 19:115–134,
2015.
- Cai et al. [2017]
T Tony Cai, Tengyuan Liang, and Alexander Rakhlin.
Computational and statistical boundaries for submatrix localization
in a large noisy matrix.
*The Annals of Statistics*, 45(4):1403–1430, 2017.
- Celentano and Montanari [2022]
Michael Celentano and Andrea Montanari.
Fundamental barriers to high-dimensional regression with convex
penalties.
*The Annals of Statistics*, 50(1):170–196,
2022.
- Chen et al. [2023a]
Hongrui Chen, Holden Lee, and Jianfeng Lu.
Improved analysis of score-based generative modeling: User-friendly
bounds under minimal smoothness assumptions.
In *International Conference on Machine Learning*, pages
4735–4763. PMLR, 2023a.
- Chen et al. [2023b]
Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R Zhang.
Sampling is as easy as learning the score: theory for diffusion
models with minimal data assumptions.
In *International Conference on Learning Representations*,
2023b.
- Deshpande and Montanari [2016]
Yash Deshpande and Andrea Montanari.
Sparse pca via covariance thresholding.
*Journal of Machine Learning Research*, 17(141):1–41, 2016.
- Ghio et al. [2024]
Davide Ghio, Yatin Dandi, Florent Krzakala, and Lenka Zdeborová.
Sampling with flows, diffusion, and autoregressive neural networks
from a spin-glass perspective.
*Proceedings of the National Academy of Sciences*, 121(27):e2311810121, 2024.
- Ho et al. [2020]
Jonathan Ho, Ajay Jain, and Pieter Abbeel.
Denoising diffusion probabilistic models.
*Advances in Neural Information Processing Systems*,
33:6840–6851, 2020.
- Hopkins et al. [2017]
Samuel B Hopkins, Pravesh K Kothari, Aaron Potechin, Prasad Raghavendra, Tselil
Schramm, and David Steurer.
The power of sum-of-squares for detecting hidden structures.
In *2017 IEEE 58th Annual Symposium on Foundations of Computer
Science (FOCS)*, pages 720–731. IEEE, 2017.
- Huang et al. [2024]
Brice Huang, Andrea Montanari, and Huy Tuan Pham.
Sampling from spherical spin glasses in total variation via
algorithmic stochastic localization.
*arXiv:2404.15651*, 2024.
- Kipf and Welling [2016]
Thomas N Kipf and Max Welling.
Semi-supervised classification with graph convolutional networks.
*arXiv preprint arXiv:1609.02907*, 2016.
- Koehler and Vuong [2024]
Frederic Koehler and Thuy-Duong Vuong.
Sampling multimodal distributions with the vanilla score: Benefits of
data-based initialization.
In *The Twelfth International Conference on Learning
Representations*, 2024.
- Kunisky et al. [2019]
Dmitriy Kunisky, Alexander S Wein, and Afonso S Bandeira.
Notes on computational hardness of hypothesis testing: Predictions
using the low-degree likelihood ratio.
In *ISAAC Congress (International Society for Analysis, its
Applications and Computation)*, pages 1–50. Springer, 2019.
- Lee et al. [2023]
Holden Lee, Jianfeng Lu, and Yixin Tan.
Convergence of score-based generative modeling for general data
distributions.
In *International Conference on Algorithmic Learning Theory*,
pages 946–985. PMLR, 2023.
- Li et al. [2024]
Gen Li, Zhihan Huang, and Yuting Wei.
Towards a mathematical theory for consistency training in diffusion
models.
*arXiv:2402.07802*, 2024.
- Ma and Wu [2015]
Zongming Ma and Yihong Wu.
Computational barriers in minimax submatrix detection.
*The Annals of Statistics*, pages 1089–1116, 2015.
- Mei and Wu [2025]
Song Mei and Yuchen Wu.
Deep networks as denoising algorithms: Sample-efficient learning of
diffusion models in high-dimensional graphical models.
*IEEE Transactions on Information Theory*, 2025.
- Montanari [2023]
Andrea Montanari.
Sampling, diffusions, and stochastic localization.
*arXiv*, 2023.
- Montanari and Wu [2023]
Andrea Montanari and Yuchen Wu.
Posterior Sampling in High Dimension via Diffusion Processes.
*arXiv:2304.11449*, 2023.
- Montanari et al. [2007]
Andrea Montanari, Federico Ricci-Tersenghi, and Guilhem Semerjian.
Solving constraint satisfaction problems through belief
propagation-guided decimation.
*arXiv:0709.1667*, 2007.
- Nguyen et al. [2017]
Hoi Nguyen, Terence Tao, and Van Vu.
Random matrices: tail bounds for gaps between eigenvalues.
*Probability Theory and Related Fields*, 167, 04 2017.
doi: 10.1007/s00440-016-0693-5.
- Peng [2012]
Minyu Peng.
Eigenvalues of deformed random matrices.
*arXiv:1205.0572*, 2012.
- Ricci-Tersenghi and Semerjian [2009]
Federico Ricci-Tersenghi and Guilhem Semerjian.
On the cavity method for decimated random constraint satisfaction
problems and the analysis of belief propagation guided decimation algorithms.
*Journal of Statistical Mechanics: Theory and Experiment*,
2009(09):P09001, 2009.
- Scarselli et al. [2008]
Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele
Monfardini.
The graph neural network model.
*IEEE transactions on neural networks*, 20(1):61–80, 2008.
- Schramm and Wein [2022]
Tselil Schramm and Alexander S Wein.
Computational barriers to estimation from low-degree polynomials.
*The Annals of Statistics*, 50(3):1833–1858, 2022.
- Shah et al. [2023]
Kulin Shah, Sitan Chen, and Adam Klivans.
Learning mixtures of gaussians using the ddpm objective.
*Advances in Neural Information Processing Systems*,
36:19636–19649, 2023.
- Song and Ermon [2019]
Yang Song and Stefano Ermon.
Generative modeling by estimating gradients of the data distribution.
*Advances in neural information processing systems*, 32, 2019.

## Appendix A Notations

Throughout the paper it will be understood that we are considering sequences of problems indexed by $n$, where $\boldsymbol{x}\in{\mathbb{R}}^{n\times n}$ and the sparsity index $k=k_{n}$ diverges as well. We write $f(n)\ll g(n)$ or $f(n)=o(g(n))$ if $f(n)/g(n)\to 0$ and $f(n)\lesssim g(n)$ or
$f(n)=O(g(n))$ if $f(n)/g(n)\leq C$ for a constant $C$. Finally $f(n)=\Theta(g(n))$ or $f(n)\asymp g(n)$
if $1/C\leq f(n)/g(n)\leq C$.

Given two measures $\mu,\nu$ on a metric space $(M,{\sf dist})$, their Wasserstein-$1$ distance is

$$ $W_{1}(\mu,\nu):=\inf_{\gamma\in\mathcal{C}(\mu,\nu)}\int_{M\times M}{\sf dist} (\boldsymbol{x},\boldsymbol{y})\,\gamma({\rm d}\boldsymbol{x},{\rm d} \boldsymbol{y})\,,$ (20) $$

where the infimum is over couplings of $\mu$, $\nu$. We will apply this definitions in two cases:
$(i)$ $M={\mathbb{R}}^{d}$, with ${\sf dist}(\boldsymbol{x},\boldsymbol{y})=\|\boldsymbol{x}-\boldsymbol{y}\|_{2}$; $(i)$ $M=C([0,T],{\mathbb{R}}^{d})$, with
${\sf dist}(\boldsymbol{x},\boldsymbol{y})=\sup_{t\in[0,T]}\|\boldsymbol{x}(t)-
\boldsymbol{y}(t)\|_{2}$.
Given random variables $\boldsymbol{X}$, $\boldsymbol{Y}$ with laws ${\rm P}_{\boldsymbol{X}}$, ${\rm P}_{\boldsymbol{Y}}$,
we overload notations $W_{1}(\boldsymbol{X},\boldsymbol{Y}):=W_{1}({\rm P}_{\boldsymbol{X}},{\rm P}_{
\boldsymbol{Y}})$.

We write $\boldsymbol{W}\sim{\sf GOE}(n)$ if $\boldsymbol{W}=\boldsymbol{W}^{{\sf T}}$ is a random symmetric matrix with $(W_{ij})_{i\leq j\leq n}$
independent entries $W_{ii}\sim{\mathsf{N}}(0,2)$, and $W_{ij}\sim{\mathsf{N}}(0,1)$ for $i<j$.
We say that $(\boldsymbol{W}_{t}:t\geq 0)$is a ${\sf GOE}(n)$ process if $\boldsymbol{W}_{t}\in{\mathbb{R}}^{n\times n}$ is a symmetric matrix
with entries above and on the diagonal $(W_{t}(i,j):i<j\leq n;W_{t}(i,i)/\sqrt{2}:i\leq n;t\geq 0)$ forming
a collection of $n(n+1)/2$ independent Brownian motions.

We use $C,C_{i},c_{i},\dots$ to denote absolute constants, whose value can change from line to line.

## Appendix B Proof of Theorem 2

### B.1 Properties of the estimator 𝒎 ^ ⁢ ( ⋅ ) ^ 𝒎 ⋅ \hat{\boldsymbol{m}}(\,\cdot\,) over^ start_ARG bold_italic_m end_ARG ( ⋅ )

###### Proposition B.1 .

Assume $\sqrt{n}\ll k\ll n$,
and note that in this case $t_{\mbox{\tiny\rm alg}}(n,k)=n/2$. Let $\hat{\boldsymbol{m}}(\,\cdot\,)$ be the estimator of
Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1) with input parameter $\varepsilon$.
For every $\delta>0$, there exists $\varepsilon>0$ such that

$$ $\displaystyle\sup_{t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}}\mathbb{E}\left[\| \hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^{2}\right]\leq C\, e^{-n\varepsilon^{2}/64k}\,.$ (21) $$

The proof of this proposition is standard, and will be presented in Appendix [H](https://arxiv.org/html/2503.08028v2#A8). We note that the rate in Equation ([21](https://arxiv.org/html/2503.08028v2#A2.E21)) gets slower the closer $k$ is to $n$; it is super-polynomial if $n\gg k\log n$.

By definition, when $\sqrt{n}\ll k\ll n$ and $t<t_{\mbox{\tiny\rm alg}}(n,k)$, the Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1) returns $\hat{\boldsymbol{m}}(\boldsymbol{y},t)=\boldsymbol{0}$, so we automatically have the following result.

###### Proposition B.2 .

For any fixed $\delta>0$, and $t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}$, we have $\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|=1$.

###### Proposition B.1 .

###### Proposition B.2 .

### B.2 Auxiliary lemmas

The following lemmas are needed for the analysis of the generated diffusion. Their proofs are deferred to
Appendices [J](https://arxiv.org/html/2503.08028v2#A10), [K](https://arxiv.org/html/2503.08028v2#A11), [L](https://arxiv.org/html/2503.08028v2#A12), [M](https://arxiv.org/html/2503.08028v2#A13).

###### Lemma B.3 .

Let $\boldsymbol{W}_{t}$ be a ${\sf GOE}$ process. Then for each time $t_{0}\geq 0$,

$$ $\mathbb{P}\left(\max_{0\leq t\leq t_{0}}\|\boldsymbol{W}_{t}\|_{{\mbox{\rm \tiny op}}}\geq 16\sqrt{t_{0}n}\right)\leq 2\exp\left(-32n\right).$ $$

###### Lemma B.4 .

Let $\boldsymbol{W}_{t}$ be a ${\sf GOE}$ process, and let $\boldsymbol{v}_{t}$ be any eigenvector of $\boldsymbol{W}_{t}$ for every $t\geq 0$. Define the set

$$ $A(\boldsymbol{v}_{t};C)=\left\{i:1\leq i\leq n,|v_{ti}|\geq\dfrac{C\sqrt{\log( n/k)}}{\sqrt{n}}\right\}$ $$

Then for any $C>4$, we have

$$ $\mathbb{P}\left(|A(\boldsymbol{v}_{t};C)|\geq\max\{\sqrt{k},k^{2}/n\}\right)=O \left(\exp(-(1/3)n^{1/4})\right)$ $$

As a consequence, using this eigenvector, $\hat{\boldsymbol{m}}$ will evaluate to $\boldsymbol{0}$ per line 8 of Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1).

###### Lemma B.5 .

Let $\boldsymbol{W}_{t}$ be a ${\sf GOE}$ process, and for each $t$, let $\boldsymbol{v}_{t}$ be a top eigenvector of $\boldsymbol{W}_{t}$. Then for any times $t_{0}\leq t_{1}$, with probability at least $1-2\exp(-32n)$,

$$ $\sup_{t_{0}\leq t\leq t_{1}}|\langle\boldsymbol{v}_{t},\boldsymbol{W}_{t_{0}} \boldsymbol{v}_{t}\rangle-\lambda_{1}(\boldsymbol{W}_{t_{0}})|\leq 32\sqrt{n(t _{1}-t_{0})}.$ $$

###### Lemma B.6 (Concentration for deformed GOE model) .

Consider the model $\boldsymbol{Y}=\theta\boldsymbol{v}\boldsymbol{v}^{{\sf T}}+\boldsymbol{W}$ for $\boldsymbol{W}\sim{\sf GOE}(n)/\sqrt{n}$ and $\theta>1$ a constant, $\boldsymbol{v}$ a unit vector. Let $\boldsymbol{v}_{1}(\boldsymbol{Y})$ be the top eigenvector of $\boldsymbol{Y}$. Define $(x^{\star},u^{\star})=(\theta+1/\theta,1-1/\theta^{2})$. For any closed set $F$ such that $d((x^{\star},u^{\star}),F)>0$, there exists a constant $c>0$ such that

$$ $\mathbb{P}\left((\lambda_{1}(\boldsymbol{Y}),\langle\boldsymbol{v}_{1}( \boldsymbol{Y}),\boldsymbol{v}\rangle^{2})\in F\right)\leq\exp(-cn)$ $$

for all $n$ large enough.

We only use Lemma [B.6](https://arxiv.org/html/2503.08028v2#A2.Thmclaim6) for the alignment $\langle\boldsymbol{v}_{1}(\boldsymbol{Y}),\boldsymbol{v}\rangle^{2}$.

###### Lemma B.3 .

###### Lemma B.4 .

###### Lemma B.5 .

###### Lemma B.6 (Concentration for deformed GOE model) .

### B.3 Analysis of the diffusion process: Proof of Theorem 2

We will prove Theorem [2](https://arxiv.org/html/2503.08028v2#Thmtheorem2) for $1\gg\varepsilon\geq C\sqrt{\log(n/k)/(n/k)}$ for some sufficiently large constant $C$.

Suppose that we generate the following diffusion, with $(\boldsymbol{z}_{t})_{t\geq 0}$ a standard $n^{2}$-dimensional Brownian motion, and $\hat{\boldsymbol{y}}_{0}=\boldsymbol{0}$:

$$ $\hat{\boldsymbol{y}}_{\ell\Delta}=\hat{\boldsymbol{y}}_{(\ell-1)\Delta}+\Delta \cdot\hat{\boldsymbol{m}}\left(\hat{\boldsymbol{y}}_{(\ell-1)\Delta},(\ell-1) \Delta\right)+\left(\boldsymbol{z}_{\ell\Delta}-\boldsymbol{z}_{(\ell-1)\Delta }\right).$ $$

We will prove that the generated diffusion never passes the termination conditions (c.f. Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1), lines 3, 5, 8).

#### B.3.1 Analysis up to an intermediate time

Define $t_{\mbox{\tiny\rm between}}=n^{2}$. Following the same strategy with Section [D.2](https://arxiv.org/html/2503.08028v2#A4.SS2), we will first show that $\hat{\boldsymbol{m}}=\boldsymbol{0}$ up to $t_{\mbox{\tiny\rm between}}$ with high probability by analyzing only the noise process (in short, if $\hat{\boldsymbol{m}}=\boldsymbol{0}$ always, our generated diffusion coincides with the noise process).
Our strategy is of the same nature as that of Section [D.2](https://arxiv.org/html/2503.08028v2#A4.SS2). Indeed, we will attempt to prove that $\hat{\boldsymbol{m}}=\boldsymbol{0}$ simultaneously for all $t$, with high probability. In this phase ($0\leq t\leq t_{\mbox{\tiny\rm between}}$), we will show that $|\hat{S}|<k/2$ (c.f. definition in Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1), lines 7, 8) for $0\leq t\leq t_{\mbox{\tiny\rm between}}$, with high probability ($\boldsymbol{v}_{t}$ is the top eigenvector of $\boldsymbol{A}_{t}$, c.f. Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1)). Note that line 5 of Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1) is not relevant in this phase. We first show this for a sequence of time points $\{t_{\ell}\}_{\ell\geq 1}$, then control the in-between fluctuations. We can set $t_{1}$ to be any value in $[0,n/2)$, as the algorithm returns 0 0 if $t<t_{\mbox{\tiny\rm alg}}=n/2$ anyway. We denote the ${\sf GOE}$ process

$$ $\boldsymbol{W}_{t}=\dfrac{\boldsymbol{B}_{t}+\boldsymbol{B}_{t}^{{\sf T}}}{2}= \sqrt{t}\boldsymbol{A}_{t}.$ $$

It is clear that the eigenvectors of $\boldsymbol{W}_{t}$ and $\boldsymbol{A}_{t}$ coincide.

We choose the following time points:

$$ $t_{\ell}=\dfrac{n}{2}-1+\dfrac{\ell}{n^{4}}.$ $$

To exceed $t_{\mbox{\tiny\rm between}}=n^{2}$, we will need $n^{6}$ values of $\ell$. By union bound from Lemma [B.4](https://arxiv.org/html/2503.08028v2#A2.Thmclaim4) (recall also the definition of the set $A(\boldsymbol{v};C)$ from this Lemma),

$$ $\displaystyle\mathbb{P}\left(\exists 1\leq\ell\leq n^{6}:|A(\boldsymbol{v}_{t_ {\ell}};C)|\geq\max\{\sqrt{k},k^{2}/n\}\right)\leq O\left(\exp(-(1/3)n^{1/4}+6 \log n)\right)$ (22) $$

Next, we will control the in-between fluctuations; specifically, we would like to show that $\max_{t_{\ell}\leq t\leq t_{\ell+1}}|A(\boldsymbol{v}_{t};C)|\leq C_{0}\max\{
\sqrt{k},k^{2}/n\}$ simultaneously for many values of $\ell$ (with high probability), for some constant $C_{0}>0$. Our approach is as follows.

- (i)
Let $\boldsymbol{v}_{t}$ be a top eigenvector of $\boldsymbol{W}_{t}$. If $t$ is close to $t_{\ell}$, then $\boldsymbol{v}_{t}$ is an approximate solution to the equation (in $\boldsymbol{v}$):
$\boldsymbol{v}^{{\sf T}}\boldsymbol{W}_{t_{\ell}}\boldsymbol{v}=\lambda_{1}(
\boldsymbol{W}_{t_{\ell}})$
- (ii)
$\boldsymbol{v}_{t}$ can be written in the coordinate system of the orthonormal eigenvectors $\boldsymbol{U}_{t_{\ell}}=[\boldsymbol{u}_{1}|\cdots|\boldsymbol{u}_{n}]$ of $\boldsymbol{W}_{t_{\ell}}$, corresponding to decreasing eigenvalues $\lambda_{1}(\boldsymbol{W}_{t_{\ell}})\geq\cdots\geq\lambda_{n}(\boldsymbol{W}
_{t_{\ell}})$. Namely, $\boldsymbol{v}_{t}=\boldsymbol{U}_{t_{\ell}}\boldsymbol{U}_{t_{\ell}}^{{\sf T}
}\boldsymbol{v}_{t}=\boldsymbol{U}_{t_{\ell}}\boldsymbol{w}$ with $\|\boldsymbol{w}\|=1$.
- (iii)
Let $m$ be a (sufficiently large) constant integer with $1\leq m\leq n$. The first $m$ components of $\boldsymbol{w}$ take up $1-o(1)$ in $L_{2}$-norm by (i) and (ii) with overwhelming probability, from which we can simply use triangle inequality to upper bound $|A(\boldsymbol{v}_{t};C)|$ according to Lemma [B.4](https://arxiv.org/html/2503.08028v2#A2.Thmclaim4) for $\boldsymbol{u}_{1},\cdots,\boldsymbol{u}_{m}$, which incurs only a constant factor of error probability, by union bound.

Define
$p_{n}=P\left(|\lambda_{1}(\boldsymbol{W}_{1})-\lambda_{7}(\boldsymbol{W}_{1})|
\leq n^{-C^{\prime}-1/2}\right)$
for any $C^{\prime}>0$ (here we take $m=7$). We use the following result (we have accounted for the scaling).

###### Lemma B.7 (Corollary 2.5, Nguyen et al. [ 2017 ] ) .

Let $\boldsymbol{W}_{1}\sim(1/\sqrt{2}){\sf GOE}(n)$. For any fixed $l\geq 1,C^{\prime}>0$, there exists a constant $c_{0}=c_{0}(l,C^{\prime})$ such that

$$ $\mathbb{P}\left(\lambda_{1}(\boldsymbol{W}_{1})-\lambda_{1+l}(\boldsymbol{W}_{ 1})\leq\dfrac{1}{2}n^{-C^{\prime}-1/2}\right)\leq c_{0}n^{-C^{\prime}\cdot \frac{l^{2}+2l}{3}}.$ $$

We materialize our approach above. We can write, with $\boldsymbol{v}_{t}=\boldsymbol{U}_{t_{\ell}}\boldsymbol{w}$:

$$ $\boldsymbol{v}_{t}^{{\sf T}}\boldsymbol{W}_{t_{\ell}}\boldsymbol{v}_{t}= \boldsymbol{w}^{{\sf T}}\boldsymbol{D}_{t_{\ell}}\boldsymbol{w}=\sum_{i=1}^{n} (D_{t_{\ell}})_{ii}w_{i}^{2}.$ $$

We then obtain that

$$ $\boldsymbol{v}_{t}^{{\sf T}}\boldsymbol{W}_{t_{\ell}}\boldsymbol{v}_{t}- \lambda_{1}(\boldsymbol{W}_{t_{\ell}})\leq\sum_{i=8}^{n}(\lambda_{i}( \boldsymbol{W}_{t_{\ell}})-\lambda_{1}(\boldsymbol{W}_{t_{\ell}}))w_{i}^{2}<- \dfrac{1}{2}\sqrt{t_{\ell}}n^{-3/2}\sum_{i=8}^{n}w_{i}^{2}$ $$

with probability at least $1-p_{n}\geq 1-c_{0}n^{-8}$, from Lemma [B.7](https://arxiv.org/html/2503.08028v2#A2.Thmclaim7) and $\boldsymbol{W}_{t_{\ell}}\sim\sqrt{t_{\ell}}\boldsymbol{W}_{1}$. Now from Lemma [B.5](https://arxiv.org/html/2503.08028v2#A2.Thmclaim5), we know that with probability at least $1-2\exp(-32n)$,

$$ $\boldsymbol{v}_{t}^{{\sf T}}\boldsymbol{W}_{t_{\ell}}\boldsymbol{v}_{t}- \lambda_{1}(\boldsymbol{W}_{t_{\ell}})\geq-32\sqrt{n(t_{\ell+1}-t_{\ell})}.$ $$

With probability at least $1-c_{0}n^{-8}-2\exp(-32n)$, both of these statements are true, uniformly over $t_{\ell}\leq t\leq t_{\ell+1}$, leading to

$$ $64\sqrt{\dfrac{t_{\ell+1}-t_{\ell}}{t_{\ell}}}\geq n^{-3/2}\sum_{i=8}^{n}w_{i} ^{2}.$ $$

A simple bit of algebra shows that

$$ $\sqrt{\dfrac{t_{\ell+1}-t_{\ell}}{t_{\ell}}}\leq 2n^{-5/2}\Rightarrow\sum_{i=8 }^{n}w_{i}^{2}\leq 128n^{-1}.$ $$

Consider the first $7$ eigenvectors $\{\boldsymbol{u}_{t_{\ell},i}\}_{i=1}^{7}$ of $\boldsymbol{W}_{t_{\ell}}$. Let

$$ $A=\bigcup_{i=1}^{7}A(\boldsymbol{u}_{t_{\ell},i};C).$ $$

From Lemma [B.4](https://arxiv.org/html/2503.08028v2#A2.Thmclaim4) and a union bound, that $|A|\leq 7\max\{\sqrt{k},k^{2}/n\}$ with probability at least $1-O(\exp(-(1/3)n^{1/4}))$. For every $j\in A^{c}$, we have

$$ $|v_{tj}|\leq\sum_{i=1}^{n}|w_{i}|\cdot|u_{t_{\ell},i,j}|<\sum_{i=1}^{7}|w_{i}| \cdot\dfrac{C\sqrt{\log(n/k)}}{\sqrt{n}}+\sum_{i=8}^{n}|w_{i}|\leq\dfrac{C^{ \prime}\sqrt{\log(n/k)}}{\sqrt{n}}+\dfrac{128}{\sqrt{n}}<C^{\prime\prime}\sqrt {\dfrac{\log(n/k)}{n}}$ $$

for a large enough constant $C^{\prime\prime}>0$. This means that with probability at least $1-O(n^{-8})$,

$$ $\sup_{t_{\ell}\leq t\leq t_{\ell+1}}|A(\boldsymbol{v}_{t};C^{\prime\prime})| \leq 7\max\{\sqrt{k},k^{2}/n\}<k/2$ $$

From Equation ([22](https://arxiv.org/html/2503.08028v2#A2.E22)) and a union bound over $\ell\geq 1$, we know that with high probability,

$$ $\sup_{n/2-1\leq t\leq n^{2}}|A(\boldsymbol{v}_{t};C)|\leq 7\max\{\sqrt{k},k^{2 }/n\}$ $$

for some absolute constant $C>0$, meaning that $\hat{\boldsymbol{m}}=0$ up to $t_{\mbox{\tiny\rm between}}$, as long as

$$ $\varepsilon>\dfrac{C\sqrt{\log(n/k)}}{\sqrt{n/k}}$ $$

###### Lemma B.7 (Corollary 2.5, Nguyen et al. [ 2017 ] ) .

#### B.3.2 Analysis to the infinite horizon

We will prove that simultaneously for all $t\geq n^{2}$, Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1) always terminates at line 5, or that $\lambda_{1}(\boldsymbol{W}_{t})\leq t/2$. Similar to Subsection [B.3.1](https://arxiv.org/html/2503.08028v2#A2.SS3.SSS1), we choose the following sequence of time points for all $\ell\geq 1$:

$$ $t_{\ell}^{(2)}=n^{2}+\ell-1$ $$

By standard Gaussian concentration and the Bai-Yin theorem, we get, for instance, the following tail bound (constants are loose) for all $x\geq 0$:

$$ $\mathbb{P}\left(\lambda_{1}\left(\dfrac{\boldsymbol{W}_{t}}{\sqrt{t}}\right) \geq 4\sqrt{n}+x\right)\leq 2\exp\left(-\dfrac{x^{2}}{2}\right)$ $$

Set $x=\dfrac{\sqrt{t}}{8}$. Since $t\geq n^{2}$, we have $x\geq n/8$, and so $x+4\sqrt{n}\leq 2x$ for $n$ large enough. Consequently,

$$ $\mathbb{P}\left(\lambda_{1}\left(\dfrac{\boldsymbol{W}_{t}}{\sqrt{t}}\right) \geq\dfrac{\sqrt{t}}{4}\right)\leq 2\exp\left(-c_{1}t\right)$ $$

for some universal constant $c_{1}>0$. A union bound for the chosen points gives:

$$ $\displaystyle\mathbb{P}\left(\exists\ell\geq 1:\lambda_{1}\left(\boldsymbol{W} _{t_{\ell}^{(2)}}\right)\geq t_{\ell}^{(2)}/4\right)\leq 2\sum_{\ell=1}^{ \infty}\exp\left(-c_{1}t_{\ell}^{(2)}\right)=2\exp(-c_{1}n^{2})\sum_{\ell=1}^{ \infty}\exp(-c_{1}(\ell-1))\lesssim\exp(-c_{1}n^{2})$ (23) $$

Next we control the in-between fluctuations. From a simple modification of Lemma [B.3](https://arxiv.org/html/2503.08028v2#A2.Thmclaim3), we have

$$ $\mathbb{P}\left(\sup_{t_{\ell}^{(2)}\leq t\leq t_{\ell+1}^{(2)}}\left|\lambda_ {1}\left(\boldsymbol{W}_{t_{\ell}^{(2)}}\right)-\lambda_{1}\left(\boldsymbol{W }_{t}\right)\right|\geq 16\sqrt{(t_{\ell+1}^{(2)}-t_{\ell}^{(2)})\cdot t_{\ell }^{(2)}n}\right)\leq 2\exp\left(-32nt_{\ell}^{(2)}\right)$ $$

so that by union bound

$$ $\displaystyle\mathbb{P}\left(\exists\ell\geq 1:\sup_{t_{\ell}^{(2)}\leq t\leq t _{\ell+1}^{(2)}}\left|\lambda_{1}\left(\boldsymbol{W}_{t_{\ell}^{(2)}}\right)- \lambda_{1}\left(\boldsymbol{W}_{t}\right)\right|\geq 16\sqrt{(t_{\ell+1}^{(2) }-t_{\ell}^{(2)})\cdot t_{\ell}^{(2)}n}\right)\lesssim\exp(-32n^{3})$ (24) $$

Consider the intersection of events described in Equations ([23](https://arxiv.org/html/2503.08028v2#A2.E23)) and ([24](https://arxiv.org/html/2503.08028v2#A2.E24)):

$$ $A=\left\{\exists\ell\geq 1:\lambda_{1}\left(\boldsymbol{W}_{t_{\ell}^{(2)}} \right)\geq t_{\ell}^{(2)}/4\right\}\cup\left\{\exists\ell\geq 1:\sup_{t_{\ell }^{(2)}\leq t\leq t_{\ell+1}^{(2)}}\left|\lambda_{1}\left(\boldsymbol{W}_{t_{ \ell}^{(2)}}\right)-\lambda_{1}\left(\boldsymbol{W}_{t}\right)\right|\geq 16 \sqrt{(t_{\ell+1}^{(2)}-t_{\ell}^{(2)})\cdot t_{\ell}^{(2)}n}\right\}$ $$

For each $t\geq n^{2}$, let $t_{\ell}$ be largest such that $t_{\ell}\leq t<t_{\ell+1}$. On $A$, we have

$$ $\lambda_{1}(\boldsymbol{W}_{t})\leq\dfrac{t_{\ell}^{(2)}}{4}+16\sqrt{(t_{\ell+ 1}^{(2)}-t_{\ell}^{(2)})\cdot t_{\ell}^{(2)}n}=\dfrac{t_{\ell}^{(2)}}{4}+16 \sqrt{t_{\ell}^{(2)}n}\leq\dfrac{t_{\ell}^{(2)}}{2}\leq\dfrac{t}{2}$ $$

for $n$ large enough, since $t_{\ell}^{(2)}\geq n^{2}\gg n$. Hence the algorithm always returns $\boldsymbol{0}$ with high probability, and we are done.

## Appendix C Proof of Theorem 1

### C.1 Auxiliary lemmas

We will use the following lemmas, whose proofs are deferred to Appendix
[G](https://arxiv.org/html/2503.08028v2#A7).

###### Lemma C.1 .

Let $\boldsymbol{W}\sim{\sf GOE}(n,1/2)$, and $C>\sqrt{2}$ some positive constant. Then we have

$$ $\displaystyle\mathbb{P}\left(\max_{\boldsymbol{v}\in\Omega_{n,k}}|\langle \boldsymbol{v},\boldsymbol{W}\boldsymbol{v}\rangle|\geq C\sqrt{\log{\binom{n}{ k}}}\right)\leq 2\binom{n}{k}^{-C^{2}/2+2}\,.$ $$

We state the following non-asymptotic result from Peng [2012].

###### Lemma C.2 (Theorem 3.1, Peng [ 2012 ] .) .

Let $\boldsymbol{y}=\theta\boldsymbol{u}\boldsymbol{u}^{{\sf T}}+{\sf GOE}(n,1/n)$, and denote by $\lambda_{1}(\boldsymbol{y})$ the top eigenvalue of $\boldsymbol{y}$. Letting $\xi(\theta):=\theta+\theta^{-1}$, the following holds for
every $x\geq 0$ and $\theta>1$:

$$ $\displaystyle\mathbb{P}\left(\lambda_{1}(\boldsymbol{y})\leq\xi(\theta)-x- \frac{2}{n}\right)\leq\exp\left(-\dfrac{(n-1)(\theta-1)^{4}}{16\theta^{2}} \right)+8\exp\left(-\dfrac{1}{4}\cdot\dfrac{(n-1)(\theta-1)^{5}x^{2}}{(\theta+ 1)^{3}}\right)\,.$ $$

In Appendix [G](https://arxiv.org/html/2503.08028v2#A7), we will use the last lemma to prove the bound below.

###### Lemma C.3 (Alignment bound) .

Let $\boldsymbol{y}=\theta\boldsymbol{u}\boldsymbol{u}^{{\sf T}}+{\sf GOE}(n,1/n)$, where $\theta=\sqrt{1+\delta}$. Let $\boldsymbol{v}_{1}$ be the top eigenvector of $\boldsymbol{y}$. Then,
there exist constants $C,c>0$ such that the following holds for any $\delta=\delta_{n}$
with $n^{-c}\ll\delta\ll 1$ for some $c>0$ small enough:

$$ $\displaystyle\mathbb{P}\left(|\langle\boldsymbol{v}_{1},\boldsymbol{u}\rangle| \leq c\delta^{2}\right)\leq C\,e^{-cn^{1/3}}\,.$ (25) $$

We also use the following lemma, which is implied in the proof of Lemma [B.4](https://arxiv.org/html/2503.08028v2#A2.Thmclaim4).

###### Lemma C.4 .

Let $\boldsymbol{g}\sim{\mathsf{N}}(\boldsymbol{0},\boldsymbol{I}_{n})$ and define the set

$$ $\displaystyle{\mathcal{L}}(\boldsymbol{g};C)=\left\{i:1\leq i\leq n,|g_{i}| \geq C\sqrt{\log(n/k)}\right\}\,.$ $$

Assume $\sqrt{n}\ll k\ll n$.
Then, for any $C$ large enough, there exists $C_{*}$ such that

$$ $\displaystyle\mathbb{P}\left(|{\mathcal{L}}(\boldsymbol{g};C)|\geq\Big{(}\sqrt {k}\vee\frac{k^{2}}{n}\Big{)}\right)\leq C_{*}e^{-n^{1/4}}\,.$ $$

###### Lemma C.1 .

###### Lemma C.2 (Theorem 3.1, Peng [ 2012 ] .) .

###### Lemma C.3 (Alignment bound) .

###### Lemma C.4 .

### C.2 Proof of Theorem 1

Let $\delta=o_{n}(1)$ be a parameter to be chosen later
and recall that $t_{\mbox{\tiny\rm alg}}=n/2$. We define $\hat{\boldsymbol{m}}$ as follows:

$$ $\displaystyle\hat{\boldsymbol{m}}(\boldsymbol{y},t)=\begin{cases}\hat{ \boldsymbol{m}}_{0}(\boldsymbol{y},t)\boldsymbol{1}_{\|\hat{\boldsymbol{m}}_{0 }(\boldsymbol{y},t)\|\leq 1-\varepsilon}&\text{ if $t\leq(1-\gamma)t_{\mbox{ \tiny\rm alg}}$,}\\ \hat{\boldsymbol{m}}_{0}(\boldsymbol{y},t)&\text{ if $(1-\gamma)t_{\mbox{\tiny \rm alg}}<t<(1+\delta)t_{\mbox{\tiny\rm alg}}$,}\\ \hat{\boldsymbol{m}}_{0}(\boldsymbol{y},t)\phi_{1}(\boldsymbol{y},t)&\text{ if $(1+\delta)t_{\mbox{\tiny\rm alg}}\leq t<n^{4}$,}\\ \hat{\boldsymbol{m}}_{0}(\boldsymbol{y},t)\phi_{2}(\boldsymbol{y},t)&\text{ if $t\geq n^{4}$,}\end{cases}$ $$

where $\phi_{1},\phi_{2}:{\mathbb{R}}^{n\times n}\times{\mathbb{R}}_{\geq 0}\to\{0,1\}$ are
defined below and $\gamma$ is given in Assumption [1](https://arxiv.org/html/2503.08028v2#Thmassumption1).
It will be clear from the constructions below that $\phi_{1},\phi_{2}$
can be evaluated in polynomial time.

In order to prove claim M2, we need to bound $J_{n}(0,\infty)$, whereby, for $t_{a}\leq t_{b}$,

$$ $\displaystyle J_{n}(t_{a},t_{b}):=\int_{t_{a}}^{t_{b}}\mathbb{E}\left[\|\hat{ \boldsymbol{m}}(\boldsymbol{y}_{t},t)-\hat{\boldsymbol{m}}_{0}(\boldsymbol{y}_ {t},t)\|^{2}\right]{\rm d}t$ $$

Setting $t_{1}=(1+\delta)t_{\mbox{\tiny\rm alg}}$ and $t_{2}=n^{4}$, we write

$$ $\displaystyle J_{n}(0,\infty)=J_{n}(0,t_{1})+J_{n}(t_{1},t_{2})+J_{n}(t_{2}, \infty)\,,$ (26) $$

and will bound each of the three terms separately.

##### Bounding J n ⁢ ( 0 , t 1 ) subscript 𝐽 𝑛 0 subscript 𝑡 1 J_{n}(0,t_{1}) italic_J start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( 0 , italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ) .

By Assumption [1](https://arxiv.org/html/2503.08028v2#Thmassumption1), we know that

$$ $\displaystyle J_{n}(0,t_{1})\leq\int_{0}^{(1-\gamma)t_{0}}\mathbb{P}(\|\hat{ \boldsymbol{m}}_{0}(\boldsymbol{y}_{t},t)\|>1-\varepsilon)\,{\rm d}t=O(n^{-D})\,,$ (27) $$

for every $D>0$.

##### Bounding J n ⁢ ( t 1 , t 2 ) subscript 𝐽 𝑛 subscript 𝑡 1 subscript 𝑡 2 J_{n}(t_{1},t_{2}) italic_J start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) .

For a matrix $\boldsymbol{y}\in{\mathbb{R}}^{n\times n}$ and time point $t$,
we define $\phi_{1}(\boldsymbol{y},t)$
according to the following procedure:

- 1.
Compute sthe symmetrized matrix $\boldsymbol{A}=(1/2)(\boldsymbol{y}+\boldsymbol{y}^{{\sf T}})$;
- 2.
Compute its top eigenvector $\boldsymbol{v}$ (choose at random if this is not unique).
- 3.
Let $S\subseteq[n]$ be the $k$ positions in $\boldsymbol{v}$ with the largest magnitude, and
define $\hat{\boldsymbol{v}}\in{\mathbb{R}}^{n}$ with $\hat{v}_{i}=(1/\sqrt{k})\operatorname{sign}(v_{i})\,\boldsymbol{1}_{i\in S}$;
- 4.
Compute the test statistic $s:=\langle\hat{\boldsymbol{v}},\boldsymbol{A}\hat{\boldsymbol{v}}\rangle$, and return $1$ if $s\geq\beta t$; 0 0 otherwise.

Here $\beta\in(0,1)$ is a fixed constant to be chosen later.

We will show that $\phi_{1}(\boldsymbol{y}_{t},t)=1$ with overwhelming probability for the true model
$\boldsymbol{y}_{t}=t\boldsymbol{x}+\boldsymbol{B}_{t}$. Define

$$ $\boldsymbol{A}_{t}=\dfrac{\boldsymbol{y}_{t}+\boldsymbol{y}_{t}^{{\sf T}}}{2}= t\boldsymbol{u}\boldsymbol{u}^{{\sf T}}+\boldsymbol{W}_{t}$ $$

where we recall that $\boldsymbol{u}\sim\operatorname{Unif}(\Omega_{n,k})$
and $\boldsymbol{W}_{t}$ is a ${\sf GOE}(n)$ process.
Let $\boldsymbol{v}_{t},\hat{\boldsymbol{v}}_{t}$ be the top eigenvector and thresholded vector of $\boldsymbol{A}_{t}$, respectively.

We have

$$ $\displaystyle s_{t}:=\langle\hat{\boldsymbol{v}}_{t},\boldsymbol{A}_{t}\hat{ \boldsymbol{v}}_{t}\rangle=t\cdot\langle\boldsymbol{u},\hat{\boldsymbol{v}}_{t }\rangle^{2}+\langle\hat{\boldsymbol{v}}_{t},\boldsymbol{W}_{t}\hat{ \boldsymbol{v}}_{t}\rangle\,.$ (28) $$

Using Lemma [C.1](https://arxiv.org/html/2503.08028v2#A3.Thmclaim1), we know that $|\langle\hat{\boldsymbol{v}}_{t},\boldsymbol{W}\hat{\boldsymbol{v}}_{t}\rangle
|\leq 4\sqrt{k\log(n/k)}\cdot\sqrt{t}$ with probability at least $1-\binom{n}{k}^{-6}$, say. Further

$$ $\boldsymbol{v}_{t}=\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle\boldsymbol{ u}+\sqrt{1-\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle^{2}}\,\boldsymbol{w }\,,$ $$

where $\boldsymbol{w}$ is a uniformly random unit vector orthogonal to $\boldsymbol{u}$, independent of $\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle$.
Alternatively, there exists $\boldsymbol{g}\sim{\mathsf{N}}(\boldsymbol{0},\boldsymbol{I}_{n})$, , independent of $\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle$
so that:

$$ $\boldsymbol{w}=\dfrac{(\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{{\sf T }})\boldsymbol{g}}{\|(\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{{\sf T} })\boldsymbol{g}\|}\,.$ $$

For every $1\leq i\leq n$, we have

$$ $|((\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{{\sf T}})\boldsymbol{g})_{ i}|\leq|g_{i}|+\dfrac{|\langle\boldsymbol{u},\boldsymbol{g}\rangle|}{\sqrt{k}}\,.$ $$

Since by assumption $k\log(n/k)\ll n$, with probability at least
$1-C_{1}\exp(-k/2)$, $|\langle\boldsymbol{u},\boldsymbol{g}\rangle|\leq\sqrt{k\log(n/k)}$ and $\|(\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{{\sf T}})\boldsymbol{g}\|
\geq\sqrt{n}/2$, so that

$$ $\displaystyle i\in{\mathcal{L}}(\boldsymbol{g};C)\;\;\Rightarrow\;\;|w_{i}| \leq(2C+2)\cdot\sqrt{\frac{1}{n}\log(n/k)}\,.$ $$

Therefore, by using Lemma [C.4](https://arxiv.org/html/2503.08028v2#A3.Thmclaim4)
and $k\log(n/k)\gg n^{1/2}$, we get that with probability at
least $1-C_{*}\exp(-n^{1/4})$, $|{\mathcal{A}}(\boldsymbol{w};C)|\leq\max\{\sqrt{k},k^{2}/n\}$ for some constant $C>0$, where

$$ $\displaystyle{\mathcal{A}}(\boldsymbol{w};C):=\left\{i:1\leq i\leq n,|w_{i}| \geq C\sqrt{\frac{1}{n}\log(n/k)}\right\}\,.$ $$

By Lemma [C.3](https://arxiv.org/html/2503.08028v2#A3.Thmclaim3), we know that with probability at
least $1-C_{*}\exp(-cn^{1/3})$ for some $C_{*},c>0$, $|\langle\boldsymbol{v}_{t_{0}(1+\delta)},\boldsymbol{u}\rangle|\geq c\delta^{2}$. Since $|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|$ is stochastically increasing with $t$ (Fact [H.2](https://arxiv.org/html/2503.08028v2#A8.Thmclaim2)), we actually have
that for any $t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}$, with the same probability
$|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|\geq c\delta^{2}$.

On the event $\mathcal{E}_{\delta}:=\{|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|\geq c
\delta^{2}\}$,
further suppose without loss of generality that
$\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle>0$. As a consequence of the result above, if $i\in\text{supp}(\boldsymbol{u})$ and $i\notin{\mathcal{A}}(\boldsymbol{w};C)$, then

$$ $\displaystyle u_{i}>0,i\notin{\mathcal{A}}(\boldsymbol{w};C)$ $\displaystyle\;\;\Rightarrow\;\;v_{t,i}\geq\frac{c\delta^{2}}{\sqrt{k}}-C\sqrt {\frac{1}{n}\log(n/k)}\,,$ $\displaystyle u_{i}<0,i\notin{\mathcal{A}}(\boldsymbol{w};C)$ $\displaystyle\;\;\Rightarrow\;\;v_{t,i}\leq-\frac{c\delta^{2}}{\sqrt{k}}+C \sqrt{\frac{1}{n}\log(n/k)}\,.$ $$

Similarly, if $i\notin\text{supp}(\boldsymbol{u})$, then

$$ $\displaystyle u_{i}=0,\;\;i\notin{\mathcal{A}}(\boldsymbol{w};C)\;\; \Rightarrow\;\;|v_{t,i}|\leq C\sqrt{\frac{1}{n}\log(n/k)}\,.$ $$

We next choose $\delta=\delta_{n}$
such that $\sqrt{(k/n)\log(n/k)}\ll\delta_{n}\ll 1$.
Hence, we conclude that

$$ $\displaystyle i\in{\rm supp}(\boldsymbol{u})\setminus{\mathcal{A}}(\boldsymbol {w};C)\;\;\Rightarrow\;\;\hat{v}_{t_{i}}=\operatorname{sign}(\langle \boldsymbol{v}_{t},\boldsymbol{u}\rangle)\cdot u_{i}\,,$ (29) $$

whence

$$ $\displaystyle|\langle\boldsymbol{u},\hat{\boldsymbol{v}}_{t}\rangle|\geq 1- \frac{2}{k}|{\mathcal{A}}(\boldsymbol{w};C)|=1-o_{n}(1)\,,$ (30) $$

with probability at least $1-C_{*}\exp(-n^{1/4}/3)$. On this event, by Eq. ([28](https://arxiv.org/html/2503.08028v2#A3.E28)) we obtain that

$$ $s_{t}\geq t\cdot(1-o_{n}(1))^{2}-4\cdot\sqrt{k\log(n/k)}\cdot\sqrt{t}$ $$

For $t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}=(1+\delta)n/2$, we this implies that, for any fixed $\beta\in(0,1)$,
with probability at least $1-C_{*}\exp(-n^{1/4}/3)$ (possibly adjusting the constant $C_{*}$).

Recalling that $\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)=\hat{\boldsymbol{m}}_{0}(
\boldsymbol{y}_{t},t)\phi_{1}(\boldsymbol{y}_{t},t)$ for $t\in[t_{1},t_{2}]$, and $\|\hat{\boldsymbol{m}}_{0}(\boldsymbol{y},t)\|\leq 1$,
we have, for $t_{1}=(1+\delta)t_{\mbox{\tiny\rm alg}}$, $t_{2}=n^{4}$ as defined above,

$$ $\displaystyle J_{n}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}\mathbb{E}[\|\hat{ \boldsymbol{m}}(\boldsymbol{y}_{t},t)-\hat{\boldsymbol{m}}_{0}(\boldsymbol{y}_ {t},t)\|^{2}]{\rm d}t\leq\int_{t_{1}}^{t_{2}}\mathbb{P}\left(\phi_{1}( \boldsymbol{y}_{t},t)=0\right)\leq C_{*}n^{4}\,e^{-cn^{1/3}}\,.$ (31) $$

##### Bounding J n ⁢ ( t 2 , ∞ ) subscript 𝐽 𝑛 subscript 𝑡 2 J_{n}(t_{2},\infty) italic_J start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ( italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , ∞ ) .

When $t\geq t_{2}$, we use a simple eigenvalue test.
For a matrix $\boldsymbol{y}$, and time point $t$, we define
$\phi_{2}(\boldsymbol{y},t)$ according to the following procedure:

- 1.
Compute symmetrized matrix $\boldsymbol{A}=(1/2)(\boldsymbol{y}+\boldsymbol{y}^{{\sf T}})$.
- 2.
Compute top eigenvalue $\lambda_{1}(\boldsymbol{A})$.
- 3.
Return $1$ if $\lambda_{1}(\boldsymbol{A})\geq t/2$, and 0 0 otherwise.

Under the true model $\boldsymbol{y}_{t}=t\boldsymbol{x}-\boldsymbol{B}_{t}$, we have:

$$ $\lambda_{1}(\boldsymbol{A}_{t})\geq t-\|\boldsymbol{W}_{t}\|_{{\mbox{\rm\tiny op }}}\geq t-t^{2/3}$ $$

with error probability given by

$$ $\mathbb{P}\left(\|\boldsymbol{W}_{t}\|_{{\mbox{\rm\tiny op}}}\geq t^{2/3} \right)\leq C\exp\left(-ct^{1/3}\right)\,,$ $$

for constants $C,c>0$. Thus, we get that

$$ $\displaystyle J_{n}(t_{2},\infty)$ $\displaystyle\leq 2\int_{t_{2}}^{\infty}\mathbb{P}(\|\boldsymbol{W}_{t}\|_{{ \mbox{\rm\tiny op}}}\geq t^{2/3})\,{\rm d}t\leq C_{*}\int_{t_{2}}^{\infty}e^{- ct^{1/3}}{\rm d}t$ $\displaystyle\leq C_{**}e^{-ct_{2}^{1/3}}\leq C_{**}e^{-cn^{4/3}}\,.$ (32) $$

Claim M2 of Theorem [1](https://arxiv.org/html/2503.08028v2#Thmtheorem1) follows from Eqs. ([27](https://arxiv.org/html/2503.08028v2#A3.E27)), ([31](https://arxiv.org/html/2503.08028v2#A3.E31)), ([32](https://arxiv.org/html/2503.08028v2#A3.E32)).

We finally consider claim M3.
Equation ([4](https://arxiv.org/html/2503.08028v2#S1.E4)) yields, for every $\ell\geq 1$:

$$ $\displaystyle\hat{\boldsymbol{y}}_{\ell\Delta}$ $\displaystyle=\Delta\sum_{i=1}^{\ell}\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}} _{(i-1)\Delta},(i-1)\Delta)+\sqrt{\Delta}\sum_{i=1}^{\ell}\boldsymbol{g}_{i \Delta}\,.$ (33) $$

We define

$$ $\displaystyle\bar{\boldsymbol{m}}_{t_{1}}$ $\displaystyle:=\Delta\sum_{i=1}^{\ell_{1}}\hat{\boldsymbol{m}}(\hat{ \boldsymbol{y}}_{(i-1)\Delta},(i-1)\Delta)\,,\,,\;\;\;\;\ell_{1}:=\lfloor t_{ \mbox{\tiny\rm alg}}(1+\delta)/\Delta\rfloor\,.$ (34) $$

We further define the following auxiliary process for $t\geq\ell_{1}\Delta=\lfloor t_{1}/\Delta\rfloor\Delta$:

$$ $\displaystyle\tilde{\boldsymbol{y}}_{t}=\bar{\boldsymbol{m}}_{t_{1}}+ \boldsymbol{B}_{t}\,,\;\;\;\;$ (35) $$

where $(\boldsymbol{B}_{t}:\,t\geq 0)$ is a Brownian motion such that $\boldsymbol{B}_{j\Delta}=\sqrt{\Delta}\sum_{i=1}^{k}\boldsymbol{g}_{i\Delta}$.
In particular, $\tilde{\boldsymbol{y}}_{\ell_{1}\Delta}=\hat{\boldsymbol{y}}_{\ell_{1}\Delta}$.

From triangle inequality, using Assumption [1](https://arxiv.org/html/2503.08028v2#Thmassumption1), the condition $\|\hat{\boldsymbol{m}}_{0}(\boldsymbol{y},t)\|_{F}\leq 1$, and that by construction $\|\hat{\boldsymbol{m}}(\boldsymbol{y},t)\|_{F}\leq 1-\varepsilon$ for $t\leq(1-\gamma)t_{\mbox{\tiny\rm alg}}$, we get

$$ $\displaystyle\|\bar{\boldsymbol{m}}\|_{F}\leq(1-\gamma)t_{\mbox{\tiny\rm alg}} (1-\varepsilon)+(\gamma+\delta)t_{\mbox{\tiny\rm alg}}\,.$ $$

We claim that (with high probability) $\phi_{1}(\tilde{\boldsymbol{y}}_{t},t)=0$ simultaneously for all $t\in[t_{1},t_{2}]$.
As a consequence, recalling the definition of $\hat{\boldsymbol{m}}$, we obtain that $\tilde{\boldsymbol{y}}_{t}=\hat{\boldsymbol{y}}_{t}$ for all $t\in[t_{1},t_{2}]$.
In order to prove the claim, define, for $\ell\geq 1$:

$$ $\displaystyle\mathcal{T}_{n}:=\big{\{}t^{+}_{\ell}=t_{1}+\dfrac{\ell-1}{n}: \ell\in{\mathbb{N}}\big{\}}\cap[t_{1},t_{2}]\,.$ $$

For every $t\in[t_{1},t_{2}]$, consider the thresholded vector $\hat{\boldsymbol{v}}_{t}$ in the definition of thest $\phi_{1}$.
We know that

$$ $\displaystyle\langle\hat{\boldsymbol{v}}_{t},\tilde{\boldsymbol{y}}_{t}\hat{ \boldsymbol{v}}_{t}\rangle=\langle\hat{\boldsymbol{v}}_{t},\bar{\boldsymbol{m} }\hat{\boldsymbol{v}}_{t}\rangle+\langle\hat{\boldsymbol{v}}_{t},\boldsymbol{B }_{t}\hat{\boldsymbol{v}}_{t}\rangle\leq(1-\gamma)t_{\mbox{\tiny\rm alg}}(1- \varepsilon)+(\gamma+\delta)t_{\mbox{\tiny\rm alg}}+\max_{\boldsymbol{v}\in \Omega_{n,k}}\langle\boldsymbol{v},\boldsymbol{B}_{t}\boldsymbol{v}\rangle\,.$ (36) $$

Using Lemma [C.1](https://arxiv.org/html/2503.08028v2#A3.Thmclaim1), we get that

$$ $\mathbb{P}\left(\max_{\boldsymbol{v}\in\Omega_{n,k}}|\langle\boldsymbol{v}, \boldsymbol{B}_{t}\boldsymbol{v}\rangle|\geq C\sqrt{\log{\binom{n}{k}}}\sqrt{t }\right)\leq 2\binom{n}{k}^{-C^{2}/2+2}\,.$ $$

Note that $|\mathcal{T}_{n}|\leq n^{5}$, whence

$$ $\mathbb{P}\left(\exists t\in\mathcal{T}_{n}:\max_{\boldsymbol{v}\in\Omega_{n,k }}|\langle\boldsymbol{v},\boldsymbol{W}_{t_{\ell}}\boldsymbol{v}\rangle|\geq C \sqrt{\log{\binom{n}{k}}}\sqrt{t_{\ell}}\right)\leq 2\binom{n}{k}^{-C^{2}/2+2} n^{5}\,.$ $$

For $t\in[t_{\ell},t_{\ell+1}]$, we have

$$ $\max_{t_{\ell}\leq t\leq t_{\ell+1}}\left\{\max_{\boldsymbol{v}\in\Omega_{n,k} }|\langle\boldsymbol{v},\boldsymbol{W}_{t}\boldsymbol{v}\rangle|-\max_{ \boldsymbol{v}\in\Omega_{n,k}}|\langle\boldsymbol{v},\boldsymbol{W}_{t_{\ell}} \boldsymbol{v}\rangle|\right\}\leq\max_{t_{\ell}\leq t\leq t_{\ell+1}}\| \boldsymbol{W}_{t}-\boldsymbol{W}_{t_{\ell}}\|_{{\mbox{\rm\tiny op}}}\,.$ $$

Using Lemma [B.3](https://arxiv.org/html/2503.08028v2#A2.Thmclaim3), we get that $\max_{t_{\ell}\leq t\leq t_{\ell+1}}\|\boldsymbol{W}_{t}-\boldsymbol{W}_{t_{
\ell}}\|_{{\mbox{\rm\tiny op}}}\leq 16\sqrt{(t_{\ell+1}-t_{\ell})n}=16$ with probability at least $1-2\exp(-32n)$.
Taking a union bound over $\mathcal{T}_{n}$, we get that

$$ $\mathbb{P}\left(\exists t\in[t_{1},t_{2}]:\max_{\boldsymbol{v}\in\Omega_{n,k}} |\langle\boldsymbol{v},\boldsymbol{W}_{t}\boldsymbol{v}\rangle|\geq C\sqrt{ \log\binom{n}{k}}\sqrt{t}\right)\leq 2\binom{n}{k}^{-C^{2}/2+2}n^{5}+2n^{5}\,e ^{-32n}\,,$ $$

for possibly a different constant $C>0$.

Using Eq. ([36](https://arxiv.org/html/2503.08028v2#A3.E36)) and the last estimate, we obtain that

$$ $\displaystyle\mathbb{P}\left(\exists t\in[t_{1},t_{2}]:\langle\hat{\boldsymbol {v}}_{t},\tilde{\boldsymbol{y}}_{t}\hat{\boldsymbol{v}}_{t}\rangle\geq b_{n}t_ {\mbox{\tiny\rm alg}}+C\sqrt{\log\binom{n}{k}}\sqrt{t}\right)\leq 2\binom{n}{k }^{-C^{2}/2+2}n^{5}+2n^{5}\exp(-32n)\,,$ (37) $$

where $b_{n}:=(1-\gamma)(1-\varepsilon)+(\gamma+\delta_{n})$.
Since $\delta_{n}=o_{n}(1)$, we have $b_{n}\to(1-\gamma)(1-\varepsilon)+\gamma<1$. Hence, we get that for $t\geq t_{1}$, and all $n$ large enough

$$ $b_{n}t_{\mbox{\tiny\rm alg}}+C\sqrt{\log\binom{n}{k}}\sqrt{t}<b_{*}t$ $$

for some constant $b_{*}\in(0,1)$ large enough. Therefore Eq. ([37](https://arxiv.org/html/2503.08028v2#A3.E37))
implies that $\phi_{1}(\tilde{\boldsymbol{y}}_{t},t)=0$ simultaneously for all $t\in[t_{1},t_{2}]$, with high probability.
We conclude that, with high probability $\hat{\boldsymbol{y}}_{t}=\tilde{\boldsymbol{y}}_{t}$ throughout $t\in[t_{1},t_{2}]$.

Finally, we extend the analysis to $t\in[t_{2},\infty)$ by proving that, with high probability,
$\phi_{2}(\tilde{\boldsymbol{y}}_{t},t)=0$ and hence $\hat{\boldsymbol{y}}_{t}=\tilde{\boldsymbol{y}}_{t}$ for all $t\in[t_{2},\infty)$.
We use (for $\boldsymbol{A}_{t}=(\tilde{\boldsymbol{y}}_{t}+\tilde{\boldsymbol{y}}_{t}^{{
\sf T}})/2$, $\boldsymbol{W}_{t}=(\boldsymbol{B}_{t}+\boldsymbol{B}_{t}^{{\sf T}})/2$)

$$ $\displaystyle\lambda_{1}(\boldsymbol{A}_{t})\leq(1-\gamma)t_{\mbox{\tiny\rm alg }}(1-\varepsilon)+(\delta+\gamma)t_{\mbox{\tiny\rm alg}}+\lambda_{1}( \boldsymbol{W}_{t})$ (38) $$

Following exactly the argument as in the proof of Theorem [2](https://arxiv.org/html/2503.08028v2#Thmtheorem2) (in particular, Subsection [B.3.2](https://arxiv.org/html/2503.08028v2#A2.SS3.SSS2)),
we get that $\lambda_{1}(\boldsymbol{W}_{t})\leq t/3$ simultaneously for all $t\geq t_{2}$, with high probability. On this event,
$\lambda_{1}(\boldsymbol{A}_{t})<t/2$ with high probability
(because $t/6\gg(1-\gamma)t_{\mbox{\tiny\rm alg}}(1-\varepsilon)+(\gamma+\delta)t_{\mbox
{\tiny\rm alg}}$). Hence, with high probability $\phi_{2}(\tilde{\boldsymbol{y}}_{t},t)=0$ and hence $\hat{\boldsymbol{y}}_{t}=\tilde{\boldsymbol{y}}_{t}$ for all $t\in[t_{2},\infty)$.

We thus proved that, with high probability, $\hat{\boldsymbol{y}}_{t}=\tilde{\boldsymbol{y}}_{t}$ for all $t\geq t_{1}=(1+\delta)t_{\mbox{\tiny\rm alg}}$, whence
$\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{t},t)=0$ as well (because $\phi_{1}(\tilde{\boldsymbol{y}}_{t},t)=0$ for $t\in[t_{1},t_{2}]$ and $\phi_{2}(\tilde{\boldsymbol{y}}_{t},t)=0$ for $t\in[t_{2},\infty)$).
Claim M3 thus follows.

## Appendix D Proof of Theorem 3

### D.1 Properties of the estimator 𝒎 ^ ⁢ ( ⋅ ) ^ 𝒎 ⋅ \hat{\boldsymbol{m}}(\,\cdot\,) over^ start_ARG bold_italic_m end_ARG ( ⋅ )

###### Proposition D.1 .

Assume $(\log n)^{2}\ll k\ll\sqrt{n}$ and let $\hat{\boldsymbol{m}}(\,\cdot\,)$ be the estimator defined in Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2). Recall that in this case, $t_{\mbox{\tiny\rm alg}}=k^{2}\log(n/k^{2})$. Then for any $\delta>0$ there exists $\varepsilon>0$ such that, letting
$s=\sqrt{(1+\varepsilon)\log(n/k^{2})}$, we have

$$ $\displaystyle\sup_{t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}}\mathbb{P}(\hat{ \boldsymbol{m}}(\boldsymbol{y}_{t},t)\neq\boldsymbol{x})=O(n^{-D})$ (39) $$

for any fixed $D>0$.

The proof of this proposition is a modification of the one in Cai et al. [2017], and will be presented in Appendix [F](https://arxiv.org/html/2503.08028v2#A6). Note that Proposition [D.1](https://arxiv.org/html/2503.08028v2#A4.Thmclaim1) directly implies the first inequality in
Condition M2v of Theorem [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3),
as $\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|\leq 2$.

By definition, when $t<t_{\mbox{\tiny\rm alg}}$, the algorithm will return $\hat{\boldsymbol{m}}=\boldsymbol{0}$, so we automatically have the following, which implies the second inequality
in Condition M2v.

###### Proposition D.2 .

For any fixed $\delta>0$, and $t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}$, we have $\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|=1$.

In the proof of Theorem [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3),
we will also make use of the following estimates, whose
proof is deferred to the Appendix [I](https://arxiv.org/html/2503.08028v2#A9).

###### Lemma D.3 .

Let $(\boldsymbol{w}_{t}:t\geq 0)$ be a process defined as

$$ $\boldsymbol{w}_{t}=\dfrac{1}{2}\left\{(\boldsymbol{B}_{t}+\sqrt{\varepsilon} \boldsymbol{g}_{t})+(\boldsymbol{B}_{t}+\sqrt{\varepsilon}\boldsymbol{g}_{t})^ {{\sf T}}\right\}$ $$

where $\boldsymbol{B},\boldsymbol{g}$ are independent $n^{2}$-dimensional Brownian motions, and $0\leq\varepsilon<1$. Then, for any $0\leq t_{0}\leq t_{1}$, and $t\geq 0$, $s\geq 1$,

$$ $\displaystyle\mathbb{P}\Big{(}\max_{t_{0}\leq t\leq t_{1}}\|\boldsymbol{w}_{t} -\boldsymbol{w}_{t_{0}}\|_{F}\geq 4\sqrt{(t_{1}-t_{0})s}\cdot n\Big{)}$ $\displaystyle\leq 2\,e^{-n^{2}s/4}\,,$ (40) $\displaystyle\mathbb{P}\big{(}\|\boldsymbol{w}_{t}\|_{F}\geq 4\sqrt{ts}\cdot n \big{)}$ $\displaystyle\leq 2\,e^{-n^{2}s/4}\,.$ (41) $$

###### Proposition D.1 .

###### Proposition D.2 .

###### Lemma D.3 .

### D.2 Analysis of the diffusion process: Proof of Theorem 3

We are left to prove that Condition M3v of
Theorem [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3) holds.

For that purpose, we make the following choices about
Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2):

- (C1)
We select the constants in the algorithm to be
$\varepsilon_{n}=o_{n}(1)$ and $s_{n}=\sqrt{(1+\varepsilon_{n})\log(n/k^{2})}$. We will use the shorthands $s=s_{n}$ and $\varepsilon=\varepsilon_{n}$, unless there is ambiguity.
- (C2)
The process $(\boldsymbol{g}_{t})_{t\geq 0}$ used in Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2) follows a $n^{2}$-dimensional Brownian motion.

Note that Propositions [D.1](https://arxiv.org/html/2503.08028v2#A4.Thmclaim1), [D.2](https://arxiv.org/html/2503.08028v2#A4.Thmclaim2)
hold under these choices, and in particular $\boldsymbol{g}_{t}\sim{\mathsf{N}}({\boldsymbol{0}},t{\boldsymbol{I}}_{n\times
n})$ at all times. Also the sequence
of random vectors $\boldsymbol{g}_{\ell\Delta}$, $\ell\in{\mathbb{N}}$
can be generated via $\boldsymbol{g}_{\ell\Delta}=\boldsymbol{g}_{(\ell-1)\Delta}+\sqrt{\Delta}\hat{
\boldsymbol{g}}_{\ell}$, for some i.i.d. standard normal vectors $\{\hat{\boldsymbol{g}}_{\ell}\}_{\ell\geq 0}$.

Letting $(\boldsymbol{z}_{t})_{t\geq 0}$ a standard Brownian motion in ${\mathbb{R}}^{n\times n}$, and $\hat{\boldsymbol{y}}_{0}=\boldsymbol{0}$
we can rewrite the approximate diffusion ([4](https://arxiv.org/html/2503.08028v2#S1.E4)) as follows (for $t\in{\mathbb{N}}\cdot\Delta$)

$$ $\displaystyle\hat{\boldsymbol{y}}_{t+\Delta}=\hat{\boldsymbol{y}}_{t}+\Delta \cdot\hat{\boldsymbol{m}}\left(\hat{\boldsymbol{y}}_{t},t\right)+\left( \boldsymbol{z}_{t+\Delta}-\boldsymbol{z}_{t}\right)\,.$ (42) $$

We further define

$$ $\displaystyle\boldsymbol{w}_{t}=\dfrac{1}{2}\big{\{}(\boldsymbol{z}_{t}+\sqrt{ \varepsilon}\boldsymbol{g}_{t})+(\boldsymbol{z}_{t}+\sqrt{\varepsilon} \boldsymbol{g}_{t})^{{\sf T}}\big{\}}\,.$ (43) $$

It is easy to see that $(c(\varepsilon)\boldsymbol{w}_{t}:t\geq 0)$ is a ${\sf GOE}(n)$ process for $c(\varepsilon):=((1+\varepsilon)/2)^{-1/2}$.
The key technical estimate in the proof of Theorem [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3), Condition M3v
is stated in the next proposition.

###### Proposition D.4 .

Let $(\boldsymbol{w}_{t}:t\geq 0)$ be defined as per Eq. ([43](https://arxiv.org/html/2503.08028v2#A4.E43)), and assume
$\varepsilon_{n}=o_{n}(1)$ and $s_{n}=\sqrt{(1+\varepsilon_{n})\log(n/k^{2})}$. Further, assume
$k\geq C(\log n)^{5/2}$ for some sufficiently large absolute constant $C>0$. Then

$$ $\displaystyle\lim_{n\to\infty}\mathbb{P}\left\{\|\eta_{s}(\boldsymbol{w}_{t}/ \sqrt{t})\|_{{\mbox{\rm\tiny op}}}\leq k+\sqrt{t}/s\;\forall t\geq 1\right\}=1\,.$ (44) $$

Before proving this proposition, let us show that it implies Condition M3v of
Theorem [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3). Indeed we claim that, with high probability,
for all $\ell\in{\mathbb{N}}$, $\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{\ell\Delta},\ell\Delta)={
\boldsymbol{0}}$ and $\hat{\boldsymbol{y}}_{\ell\Delta}=\boldsymbol{z}_{\ell{\Delta}}$.
This is proven by induction over $\ell$. Indeed, if it holds up to a certain $\ell-1\in{\mathbb{N}}$,
then we have $\hat{\boldsymbol{y}}_{\ell\Delta}=\boldsymbol{z}_{\ell\Delta}$ by Eq. ([42](https://arxiv.org/html/2503.08028v2#A4.E42))
whence it follows that $\boldsymbol{A}_{t,+}=\boldsymbol{w}_{t}/\sqrt{t}$, for $t=\ell\Delta$ (c.f. Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2), line 4)
and therefore $\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{t},t)={\boldsymbol{0}}$ by Proposition [D.4](https://arxiv.org/html/2503.08028v2#A4.Thmclaim4)
(because the condition in Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2), line 6, is never passed).

We therefore have

$$ $\displaystyle\inf_{\ell\geq 0}W_{1}(\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_ {\ell\Delta},\ell\Delta),\boldsymbol{x})\geq\inf_{\ell\geq 0}\mathbb{P}(\hat{ \boldsymbol{m}}(\hat{\boldsymbol{y}}_{\ell\Delta},\ell\Delta)={\boldsymbol{0}} )=1-o_{n}(1)\,.$ (45) $$

This concludes the proof of Theorem [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3). We next turn to proving Proposition [D.4](https://arxiv.org/html/2503.08028v2#A4.Thmclaim4).

###### Proof of Proposition D.4 .

We follow a strategy analogous to the proof of the Law of Iterated Logarithm. We choose a sparse sequence of time points $\{t_{\ell}\}_{\ell=1}^{\infty}$, and $(i)$ establish the statement jointly for these time points, and $(ii)$ control deviations in between. In particular, we consider

$$ $t_{\ell}=\left(1+\dfrac{\ell-1}{n^{3}}\right)^{2}$ $$

for all $\ell\geq 1$.
We first show that simultaneously for all $\ell\geq 1$, we have $\max_{i,j}\left|(\boldsymbol{w}_{t_{\ell}})_{ij}/\sqrt{t_{\ell}}\right|\leq 8t
_{\ell}^{1/4}\sqrt{\log n}$. We have, by sub-gaussianity of $(\boldsymbol{w}_{t_{\ell}})_{ij}$ and a union bound (here we account also for the case where $i=j$, in which there is an inflated variance), along with $\varepsilon=o_{n}(1)$:
using the bound $2xy\geq x+y$ when $x,y\geq 1$ for $x=t_{\ell}^{1/2},y=\log n$. Taking a union bound once again over $\ell$, we have

$$ $\mathbb{P}\left(\exists\ell\geq 1,1\leq i,j\leq n:|(\boldsymbol{w}_{t_{\ell}}) _{ij}|\geq 8t_{\ell}^{3/4}\sqrt{\log n}\right)\leq n^{-6}\cdot\sum_{\ell=0}^{ \infty}\exp\left(-8\cdot\left(1+\dfrac{\ell}{n^{3}}\right)\right)$ $$

We have, as the summands form a decreasing function of $\ell$ integer:

$$ $\displaystyle\sum_{\ell=0}^{\infty}\exp\left(-8\cdot\left(1+\dfrac{\ell}{n^{3} }\right)\right)\leq$ $\displaystyle C+\int_{0}^{\infty}\exp\left(-\dfrac{8x}{n^{3}}\right){\rm d}x \leq Cn^{3}\,.$ (46) $$

We thus obtain that

$$ $\displaystyle\mathbb{P}\left(\exists\ell\geq 1,1\leq i,j\leq n:|(\boldsymbol{w }_{t_{\ell}})_{ij}|\geq 8t_{\ell}^{3/4}\sqrt{\log n}\right)=O(n^{-3})\,.$ (47) $$

The point of this calculation is that simultaneously for all $\ell\geq 1$, we can truncate the entries of $\eta_{s}(\boldsymbol{w}_{t_{\ell}}/\sqrt{t_{\ell}})$ by $8t_{\ell}^{1/4}\sqrt{\log n}$ without worry.

Namely, for each $\ell\geq 1$, $\vartheta_{\ell}=8t_{\ell}^{1/4}\sqrt{\log n}$
we define $\tilde{\boldsymbol{w}}_{t_{\ell}}\in{\mathbb{R}}^{n\times n}$
by

$$ $\displaystyle(\tilde{\boldsymbol{w}}_{t_{\ell}})_{ij}:=\begin{cases}\eta_{s}( \boldsymbol{w}_{t_{\ell}}/\sqrt{t_{\ell}})&\mbox{ if }\;\;|\eta_{s}( \boldsymbol{w}_{t_{\ell}}/\sqrt{t_{\ell}})|\leq\vartheta_{\ell}\,,\\ \vartheta_{\ell}&\mbox{ if }\;\;\eta_{s}(\boldsymbol{w}_{t_{\ell}}/\sqrt{t_{ \ell}})>\vartheta_{\ell}\,,\\ -\vartheta_{\ell}&\mbox{ if }\eta_{s}(\boldsymbol{w}_{t_{\ell}}/\sqrt{t_{\ell} })<-\vartheta_{\ell}\,.\end{cases}$ (48) $$

By Eq. ([47](https://arxiv.org/html/2503.08028v2#A4.E47)), we have

$$ $\displaystyle\mathbb{P}\left(\exists\ell\geq 1\;:\eta_{s}(\boldsymbol{w}_{t_{ \ell}}/\sqrt{t_{\ell}})\neq\tilde{\boldsymbol{w}}_{t_{\ell}}\right)=O(n^{-3})\,.$ (49) $$

We have from Bandeira and van Handel [2016], for every $x\geq 0$:

$$ $\displaystyle\mathbb{P}\left(\left\|\tilde{\boldsymbol{w}}_{t_{\ell}}\right\|_ {\text{op}}\geq 4\sigma+x\right)\leq n\exp\left(-\dfrac{cx^{2}}{\sigma_{\star} ^{2}}\right)\,,$ (50) $$

for some absolute constant $c>0$, where

$$ $\displaystyle\sigma^{2}$ $\displaystyle:=\max_{i\leq n}\sum_{j=1}^{n}\mathbb{E}\left[(\tilde{\boldsymbol {w}}_{t_{\ell}})_{ij}^{2}\right]\leq\sum_{j=1}^{n}\mathbb{E}\left[\eta_{s} \left(\dfrac{\boldsymbol{w}_{t_{\ell}}}{\sqrt{t_{\ell}}}\right)_{ij}^{2}\right ]\,,$ (51) $\displaystyle\sigma_{\star}$ $\displaystyle:=\max_{i,j\leq n}|(\tilde{\boldsymbol{w}}_{t_{\ell}})_{ij}|\leq 8 t_{\ell}^{1/4}\sqrt{\log n}\,.$ (52) $$

It can be seen from an immediate Gaussian calculation that, for $i\neq j$ and $Z\sim{\mathsf{N}}(0,1)$:

$$ $\displaystyle\mathbb{E}\left[\eta_{s}\left(\dfrac{\boldsymbol{w}_{t_{\ell}}}{ \sqrt{t_{\ell}}}\right)_{ij}^{2}\right]=$ $\displaystyle\int_{0}^{\infty}4z\cdot\mathbb{P}\left(\sqrt{\frac{1+\varepsilon }{2}}Z\geq z+s\right){\rm d}z$ $\displaystyle\overset{(a)}{\leq}$ $\displaystyle\int_{0}^{\infty}4z\cdot\dfrac{1}{z+s}\cdot\exp\left(-\dfrac{(z+s )^{2}}{1+\varepsilon}\right){\rm d}z$ $\displaystyle\overset{(b)}{\ll}$ $\displaystyle\dfrac{1}{s}\exp\left(-\dfrac{s^{2}}{1+\varepsilon}\right)\leq \exp\left(-\dfrac{s^{2}}{1+\varepsilon}\right)$ $$

Here in $(a)$ we employ the Mill’s ratio bound, and $(b)$ follows from $z+s\geq s$ and $s\to\infty$.

Proceeding analogously for the diagonal entries of $\eta_{s}(\boldsymbol{w}_{t_{\ell}}/\sqrt{t_{\ell}})$, we obtain that $\sigma\ll\sqrt{n}\exp(-s^{2}/(2(1+\varepsilon)))=k$ by definition of $s$.

We set $x=k/3+\sqrt{t_{\ell}}/(3s)$. Since $x\gg\sigma$, we have $4\sigma+x\leq(3/2)x$ if $n,k$ are sufficiently large.
Using Eq. ([50](https://arxiv.org/html/2503.08028v2#A4.E50)) we obtain that, for some universal constants $c,c^{\prime},c^{\prime\prime}>0$:

$$ $\mathbb{P}\left(\|\tilde{\boldsymbol{w}}_{t_{\ell}}\|_{\text{op}}\geq\dfrac{k} {2}+\dfrac{\sqrt{t_{\ell}}}{2s}\right)\leq n\exp\left(-\dfrac{c\left(\dfrac{k} {3}+\dfrac{\sqrt{t_{\ell}}}{3s}\right)^{2}}{64t_{\ell}^{1/2}\log n}\right) \overset{(a)}{\leq}n\exp\left(-\dfrac{c^{\prime}k}{s\log n}-\dfrac{c^{\prime \prime}t_{\ell}^{1/2}}{s^{2}\log n}\right)\,.$ $$

In step $(a)$, we simply expand the squared term in the numerator and drop the quadratic term in $k$. Now, taking a union bound over $\ell\geq 1$, we get that (similar to Eq. ([46](https://arxiv.org/html/2503.08028v2#A4.E46)))

$$ $\displaystyle\mathbb{P}\left(\exists\ell\geq 1:\|\tilde{\boldsymbol{w}}_{t_{ \ell}}\|_{\text{op}}\geq\dfrac{k}{2}+\dfrac{\sqrt{t_{\ell}}}{2s}\right)\leq$ $\displaystyle\;n\exp\left(-\dfrac{c^{\prime}k}{s\log n}\right)\sum_{\ell=1}^{ \infty}\exp\left(-\dfrac{c^{\prime\prime}t_{\ell}^{1/2}}{s^{2}\log n}\right)$ $\displaystyle\leq$ $\displaystyle\;n\exp\left(-\dfrac{c^{\prime}k}{s\log n}\right)\left(O(1)+\int_ {0}^{\infty}\exp\left(-\dfrac{c^{\prime\prime}x}{s^{2}n^{3}\log n}\right){\rm d }x\right)$ $\displaystyle=$ $\displaystyle\;O\left(\exp\left(-\dfrac{c^{\prime}k}{s\log n}\right)\cdot s^{2 }n^{4}\log n\right)$ $\displaystyle=$ $\displaystyle\;o_{n}(1)\,,$ $$

where the last estimate holds if $k\geq C(\log n)^{5/2}$ for some sufficiently large $C>0$.
In conclusion, using the last display and Eq. ([49](https://arxiv.org/html/2503.08028v2#A4.E49)) we have shown that

$$ $\displaystyle\mathbb{P}\left(\exists\ell\geq 1:\left\|\eta_{s}\left(\dfrac{ \boldsymbol{w}_{t_{\ell}}}{\sqrt{t_{\ell}}}\right)\right\|_{\text{op}}\geq \dfrac{k}{2}+\dfrac{\sqrt{t_{\ell}}}{2s}\right)=o_{n}(1)\,.$ (53) $$

Now we control the in-between fluctuations. Noting that $\eta_{s}(\cdot)$ is a $1$-Lipschitz function, we have the following crude bound:

$$ $\displaystyle\max_{t_{\ell}\leq t\leq t_{\ell+1}}\left\|\eta_{s}\left(\dfrac{ \boldsymbol{w}_{t}}{\sqrt{t}}\right)-\eta_{s}\left(\dfrac{\boldsymbol{w}_{t_{ \ell}}}{\sqrt{t_{\ell}}}\right)\right\|_{\text{op}}\leq$ $\displaystyle\max_{t_{\ell}\leq t\leq t_{\ell+1}}\left\|\dfrac{\boldsymbol{w}_ {t}}{\sqrt{t}}-\dfrac{\boldsymbol{w}_{t_{\ell}}}{\sqrt{t_{\ell}}}\right\|_{F}$ $\displaystyle\leq$ $\displaystyle\dfrac{\max_{t_{\ell}\leq t\leq t_{\ell+1}}\|\boldsymbol{w}_{t}- \boldsymbol{w}_{t_{\ell}}\|_{F}}{\sqrt{t_{\ell}}}+\|\boldsymbol{w}_{t_{\ell}} \|_{F}\left(\dfrac{1}{\sqrt{t_{\ell}}}-\dfrac{1}{\sqrt{t_{\ell+1}}}\right)\,.$ $$

From Lemma [D.3](https://arxiv.org/html/2503.08028v2#A4.Thmclaim3), we obtain that

$$ $\displaystyle\mathbb{P}\left(\max_{t_{\ell}\leq t\leq t_{\ell+1}}\left\|\eta_{ s}(\boldsymbol{w}_{t}/\sqrt{t})-\eta_{s}(\boldsymbol{w}_{t_{\ell}}/\sqrt{t_{ \ell}})\right\|_{\text{op}}\geq 4n\cdot\sqrt{t_{\ell+1}-t_{\ell}}+4n\cdot\sqrt {t_{\ell}}\cdot\left(1-\sqrt{\dfrac{t_{\ell}}{t_{\ell+1}}}\right)\right)\leq 4 \,e^{-n^{2}t_{\ell}/4}\,.$ $$

By definition of $t_{\ell}$, simple algebra reveals that (we also use the fact that $n^{-1/2}\ll s^{-1}$):

$$ $4n\cdot\sqrt{t_{\ell+1}-t_{\ell}}+4n\cdot\sqrt{t_{\ell}}\cdot\left(1-\sqrt{ \dfrac{t_{\ell}}{t_{\ell+1}}}\right)\leq\\ \dfrac{\sqrt{t_{\ell}}}{2s}\,.$ $$

By union bound over $\ell\geq 1$,

$$ $\displaystyle\mathbb{P}\left(\exists\ell\geq 1:\max_{t_{\ell}\leq t\leq t_{ \ell+1}}\left\|\eta_{s}(\boldsymbol{w}_{t}/\sqrt{t})-\eta_{s}(\boldsymbol{w}_{ t_{\ell}}/\sqrt{t_{\ell}})\right\|_{\text{op}}\geq\dfrac{k}{2}+\dfrac{\sqrt{t_ {\ell}}}{2s}\right)$ $\displaystyle\leq$ $\displaystyle 4\sum_{\ell=1}^{\infty}\exp\left(-\dfrac{n^{2}}{8}-\dfrac{t_{ \ell}}{8}\right)$ $\displaystyle=$ $\displaystyle 4\exp(-n^{2}/8)\sum_{\ell=1}^{\infty}\exp\left(-\dfrac{1}{8} \left(1+\dfrac{\ell-1}{n^{3}}\right)^{2}\right)$ $\displaystyle\leq$ $\displaystyle 4\exp(-n^{2}/8)\left(O(1)+\int_{0}^{\infty}\exp\left(-\dfrac{x^{ 2}}{8n^{6}}\right)dx\right)=O(\exp(-n^{2}/8)n^{3})=o(1)\,.$ $$

Using this estimate together with Eq. ([53](https://arxiv.org/html/2503.08028v2#A4.E53)), we conclude that
with high probability the following holds simultaneously for all $t\geq 1$. Letting $\ell$ be largest such that $t_{\ell}\leq t$:

$$ $\left\|\eta_{s}(\boldsymbol{w}_{t}/\sqrt{t})\right\|_{\text{op}}\leq\left\| \eta_{s}(\boldsymbol{w}_{t_{\ell}}/\sqrt{t_{\ell}})\right\|_{\text{op}}+\left \|\eta_{s}(\boldsymbol{w}_{t}/\sqrt{t})-\eta_{s}(\boldsymbol{w}_{t_{\ell}}/ \sqrt{t_{\ell}})\right\|_{\text{op}}\leq k+\dfrac{\sqrt{t_{\ell}}}{s}\leq k+ \dfrac{\sqrt{t}}{s}\,,$ $$

and this finishes the proof.
∎

We remark that Assumption (C2) in the proof above is technically not needed, meaning that the additional noise stream $\boldsymbol{g}_{t}$ can in fact be discarded entirely: an appropriate thresholding of $\boldsymbol{v}_{t}$, the top eigenvector of $\eta_{s}(\boldsymbol{A}_{t,+})$, as in Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1), will also suffice to satisfy all conditions of Theorem [3](https://arxiv.org/html/2503.08028v2#Thmtheorem3), although $\boldsymbol{x}$ will not be recovered exactly; some $o(k)$ positions outside the support of $\boldsymbol{x}$ will also be chosen, at most. The reason for this is that the alignment $|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|=1-o_{n}(1)$ already, from a close inspection of our proof of Proposition [D.1](https://arxiv.org/html/2503.08028v2#A4.Thmclaim1). Regarding the proof of Proposition [D.4](https://arxiv.org/html/2503.08028v2#A4.Thmclaim4) above, one can easily realize that even if $\varepsilon=0$, it will go through without any modification. We choose to keep our formulation of Algorithm [2](https://arxiv.org/html/2503.08028v2#alg2) as faithful to the original work of Cai et al. [2017] as possible to discuss a variety of approaches, and leave this to the interested reader.

###### Proposition D.4 .

###### Proof of Proposition D.4 .

## Appendix E Proof of Proposition 2.1

We take the first row of $\hat{\boldsymbol{z}}_{1}$, and let $A=\{z_{11},\cdots,z_{1n}\}$. Let $r_{j}=\text{rank}(z_{1j})$ denote the rank of $z_{1j}$ with respect to the elements of $A$. Then since $z_{1j}\sim{\mathsf{N}}(0,1)$ across $j$, the collection of the first $k$ ranks $A_{k}=\{r_{1},\cdots,r_{k}\}$ constitutes a sample without replacement from $[n]$. Construct $\boldsymbol{v}$ a binary vector such that $v_{i}=1$ if and only if $i\in A_{k}$, and let $\boldsymbol{u}$
be a randomized-sign vector version of $(1/\sqrt{k})\boldsymbol{v}$. Let

$$ $\displaystyle\hat{\boldsymbol{m}}(\boldsymbol{y},t;\boldsymbol{g}_{1})= \boldsymbol{u}\boldsymbol{u}^{{\sf T}}=\boldsymbol{x}^{\prime}$ (54) $$

then it is clear that $\hat{\boldsymbol{m}}(\boldsymbol{y},t;\hat{\boldsymbol{z}}_{1})\sim\boldsymbol
{x}$ and is independent of $\boldsymbol{x}$ (as it is a function of only $\hat{\boldsymbol{z}}_{1}$). The identity from $(i)$ follows accordingly. To see that this error is clearly sub-optimal compared to polynomial time algorithms, observe that $\hat{\boldsymbol{m}}=\boldsymbol{0}$ is a polynomial time drift, which achieves error $1$ at every $t$.

Point $(ii)$ also follows immediately. Indeed, for every $\ell\geq 0$,

$$ $W_{1}(\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{\ell\Delta},\ell\Delta), \boldsymbol{x})=W_{1}(\boldsymbol{x}^{\prime},\boldsymbol{x})=0$ $$

Lastly, regarding point $(iii)$, note that since $\boldsymbol{x}^{\prime}$ is not dependent on $t$, we have, for every $\ell\geq 1$,

$$ $\hat{\boldsymbol{y}}_{\ell\Delta}=\hat{\boldsymbol{y}}_{(\ell-1)\Delta}+\Delta \boldsymbol{x}^{\prime}+\sqrt{\Delta}\hat{\boldsymbol{z}}_{\ell\Delta}$ $$

Simple induction gives

$$ $\hat{\boldsymbol{y}}_{\ell\Delta}=(\ell\Delta)\boldsymbol{x}^{\prime}+\sqrt{ \Delta}\sum_{j=1}^{\ell}\hat{\boldsymbol{z}}_{j\Delta}$ $$

We take the coupling of $(\hat{\boldsymbol{y}}_{\ell\Delta}/(\ell\Delta),\boldsymbol{x})$ such that $\boldsymbol{x}^{\prime}=\boldsymbol{x}$. Then by definition of the Wasserstein-$2$ metric,

$$ $W_{2}(\hat{\boldsymbol{y}}_{\ell\Delta}/(\ell\Delta),\boldsymbol{x})^{2}\leq \mathbb{E}\left[\left\|\dfrac{\sqrt{\Delta}\sum_{j=1}^{\ell}\hat{\boldsymbol{z }}_{j\Delta}}{\ell\Delta}\right\|^{2}\right]=\dfrac{n^{2}}{\ell\Delta}$ $$

It is clear that as $\ell\to\infty$, this quantity converges to 0 0. Hence we are done.

## Appendix F Proof of Proposition D.1

We conduct our analysis conditional on $\boldsymbol{u}$ (c.f. Equation ([1.2](https://arxiv.org/html/2503.08028v2#A17.EGx6))), and $S$ be the support of $\boldsymbol{u}$. Let $\boldsymbol{A}_{0}=\sqrt{t}\boldsymbol{u}\boldsymbol{u}^{{\sf T}}$
and notice that

$$ $\displaystyle\boldsymbol{A}_{t,+}=\boldsymbol{A}_{0}+\sigma_{+}\boldsymbol{Z} \,,\;\;\;\;\;\boldsymbol{A}_{t,-}=\boldsymbol{A}_{0}+\sigma_{-}\boldsymbol{W}\,,$ (55) $$

where $\sigma^{2}_{+}:=(1+\varepsilon)/2$, $\sigma^{2}_{-}:=(1+\varepsilon)/(2\varepsilon)$,
and $\boldsymbol{Z},\boldsymbol{W}\sim{\sf GOE}(n)$ are independent random matrices.

We have

$$ $\displaystyle\eta_{s}(\boldsymbol{A}_{t,+})=\boldsymbol{A}_{0}+\eta_{s}(\sigma _{+}\boldsymbol{Z})+\mathbb{E}[\boldsymbol{B}]+\big{(}\boldsymbol{B}-\mathbb{E }[\boldsymbol{B}]\big{)}\,,$ (56) $$

where

$$ $\displaystyle B_{ij}=\eta_{s}\Big{(}\sqrt{t}\cdot u_{i}u_{j}+\sigma_{+}Z_{ij} \Big{)}-\sqrt{t}\cdot u_{i}u_{j}-\eta_{s}(\sigma_{+}Z_{ij})\,.$ (57) $$

Our first order of business is to analyze $\mathbb{E}[\boldsymbol{B}]$. If $i\notin S$ or $j\notin S$, we have $\mathbb{E}[B_{ij}]=0$.
On the other hand, if $i,j\in S$, then

- (i)
Case 1: $u_{i}u_{j}=1/k$
In this case, we have $\mathbb{E}[B_{ij}]=-b_{0}+b_{1}\boldsymbol{1}_{i=j}$ where
(below $G\sim{\mathsf{N}}(0,1)$)
$\displaystyle b_{0}:=-\mathbb{E}\Big{\{}\eta_{s}\Big{(}\frac{\sqrt{t}}{k}+
\sigma_{+}G\Big{)}-\frac{\sqrt{t}}{k}\Big{\}}\,,\;\;\;b_{1}:=\mathbb{E}\Big{\{
}\eta_{s}\Big{(}\frac{\sqrt{t}}{k}+\sqrt{2}\sigma_{+}G\Big{)}-\eta_{s}\Big{(}
\frac{\sqrt{t}}{k}+\sigma_{+}G\Big{)}\Big{\}}\,.$
(58)
Recalling that $\sigma_{+}$ is bounded and bounded away from 0 0 (without loss of generality we can assume $\varepsilon<1/2$) and $s,\sqrt{t}/k-s$ grows with $n,k$, so that
$\eta_{s}(\sqrt{t}/k+Z_{ij})=\sqrt{t}/k+Z_{ij}-s$ with high probability; hence $|B_{ij}+s|=o_{P}(s)$ (as $Z_{ij}=o(s)$ with high probability). Noting that $|B_{ij}|\leq 2s$, we get $b_{0}=s(1+o(1))$ and
$b_{1}=o(s)$ (distribution on diagonal is different).
- (ii)
Case 2: $u_{i}u_{j}=-1/k$
By a similar reasoning, we have $\mathbb{E}[B_{ij}]=s(1+o(1))$.

We can thus rewrite

$$ $\displaystyle\eta_{s}(\boldsymbol{A}_{t,+})=(\sqrt{t}-kb_{0})\cdot\boldsymbol{ u}\boldsymbol{u}^{{\sf T}}+b_{1}\cdot{\boldsymbol{P}}_{S}+\eta_{s}(\sigma_{+} \boldsymbol{Z})+\big{(}\boldsymbol{B}-\mathbb{E}[\boldsymbol{B}]\big{)}\,,$ (59) $$

where $(\boldsymbol{P}_{S})_{ij}=1$ if $i=j\in S$ and $=0$ otherwise.

Next, we analyze the operator norm of $\eta_{s}(\sigma_{+}\boldsymbol{Z})$.
Let $\tilde{\boldsymbol{Z}}=(\tilde{Z}_{ij})_{i,j\leq n}$ be defined as

$$ $\displaystyle\tilde{Z}_{ij}=\eta_{s}(\sigma_{+}Z_{ij})\,\boldsymbol{1}(|\eta_{ s}(\sigma_{+}Z_{ij})|\leq\sqrt{2D\log n})\,.$ (60) $$

for some constant $C>0$ that we choose such that $\max_{ij\leq n}|Z_{ij}|\leq C\sqrt{\log n}$ with probability at least $1-n^{-D}$ ($D$ is defined in the statement of Proposition [D.1](https://arxiv.org/html/2503.08028v2#A4.Thmclaim1)). We have
$\tilde{\boldsymbol{Z}}=\eta_{s}(\boldsymbol{Z})$. By Bandeira and van Handel [2016], there exists an absolute constant $c>0$ such that for every $u>0$:

$$ $\displaystyle\mathbb{P}\left(\|\tilde{\boldsymbol{Z}}\|_{{\mbox{\rm\tiny op}}} \geq 4\sigma+u\right)\leq n\exp\left(-\frac{cu^{2}}{L^{2}}\right)\,,$ (61) $$

where

$$ $\displaystyle\sigma^{2}$ $\displaystyle=\max_{i\leq n}\sum_{j=1}^{n}\mathbb{E}[\eta_{s}(Z_{ij})^{2}]\,,$ (62) $\displaystyle L$ $\displaystyle=\max_{i,j\leq n}\|\tilde{Z}_{ij}\|_{\infty}\leq C\sqrt{\log n}\,.$ (63) $$

An immediate Gaussian calculation yields, for $i\neq j$:

$$ $\displaystyle\mathbb{E}[\eta_{s}(\sigma_{+}Z_{ij})^{2}]=$ $\displaystyle\int_{0}^{\infty}4z\cdot\mathbb{P}(\sigma_{+}Z_{ij}\geq z+s){\rm d }z\leq C_{1}e^{-s^{2}/(1+\varepsilon)}\,.$ (64) $$

for some constant $C_{1}>0$.

Proceeding analogously for $\eta_{s}(\sigma_{+}Z_{ii})$ and substituting in Eq. ([61](https://arxiv.org/html/2503.08028v2#A6.E61)),
we get $\sigma^{2}\leq 2C_{1}n\exp\{-s^{2}/(1+\varepsilon)\}$.
Applying Eq. ([61](https://arxiv.org/html/2503.08028v2#A6.E61)) there exists an absolute constant $C>0$ such that, with probability at least $1-2n^{-D}$,

$$ $\displaystyle\|\eta_{s}(\boldsymbol{Z})\|_{{\mbox{\rm\tiny op}}}$ $\displaystyle\leq C\Big{(}\sqrt{n}\exp\{-s^{2}/2(1+\varepsilon)\}\vee\sqrt{ \log n}\Big{)}\,$ (65) $\displaystyle\leq C\Big{(}k\vee\sqrt{\log n}\Big{)}\leq Ck\,.$ (66) $$

Lastly, consider $\boldsymbol{B}-\mathbb{E}[\boldsymbol{B}]$.
By Eq. ([57](https://arxiv.org/html/2503.08028v2#A6.E57)) we know that the entries of this matrix are independent
with mean 0 0 and bounded by $2s$, hence subgaussian. Further only a $k\times k$ submatrix
is nonzero, so that

$$ $\displaystyle\|\boldsymbol{B}-\mathbb{E}[\boldsymbol{B}]\|_{{\mbox{\rm\tiny op }}}\leq C_{1}\sqrt{k}s\,,$ (67) $$

with high probability (for instance, the operator norm tail bound above can be applied once more, which gives an error probability of at most $C\exp(-ck)$ for some absolute constant $C,c>0$).

Summarizing, we proved that

$$ $\displaystyle\eta_{s}(\boldsymbol{A}_{t,+})$ $\displaystyle=(\sqrt{t}-kb_{0})\cdot\boldsymbol{u}\boldsymbol{u}^{{\sf T}}+ \boldsymbol{\Delta}\,,$ (68) $\displaystyle\|\boldsymbol{\Delta}\|_{{\mbox{\rm\tiny op}}}$ $\displaystyle\leq C(k+\sqrt{k}s)\leq C^{\prime}k\,,$ (69) $$

where in the last step we used $k\gg(\log n)^{2}$

Recall that $\boldsymbol{v}_{t}$ denotes the top eigenvector of $\eta_{s}(\boldsymbol{A}_{t,+})$. By Davis-Kahan,

$$ $\displaystyle\min_{a\in\{+1,-1\}}\left\|\boldsymbol{v}_{t}-a\boldsymbol{u}\right\|$ $\displaystyle\leq\dfrac{Ck}{\sqrt{t}-kb_{0}}$ (70) $\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\dfrac{Ck}{\sqrt{t}-(1+ \varepsilon)ks}$ (71) $\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\varepsilon\,,$ (72) $$

where in $(a)$ we used the fact that $b_{0}=s+o(s)$ and in $(b)$ the fact that $t\geq(1+\delta)k^{2}\log(n/k^{2})$, whereby we can assume $\delta\geq C\varepsilon$ for $C$ a sufficiently large absolute constant.
Recalling the definition of the score $\hat{\boldsymbol{v}}_{t}$:

$$ $\hat{\boldsymbol{v}}_{t}=\boldsymbol{A}_{t,-}\boldsymbol{v}_{t}=\sqrt{t} \boldsymbol{u}\langle\boldsymbol{u},\boldsymbol{v}_{t}\rangle+\sigma_{-} \boldsymbol{W}\boldsymbol{v}_{t}$ $$

where we know that $\boldsymbol{G}=\boldsymbol{W}\boldsymbol{v}_{t}\sim{\mathsf{N}}(0,\boldsymbol{
I}_{n})$ by independence of $\boldsymbol{W}$ and $\boldsymbol{v}_{t}$.
Assuming to be definite that the sign of the eigenvector is chosen so that the last bound holds with $a=+1$, we get that $\langle\boldsymbol{u},\boldsymbol{v}_{t}\rangle\geq 1-\varepsilon^{2}$. We get that for every $i\in S$:

$$ $\displaystyle u_{i}>0$ $\displaystyle\Rightarrow\hat{v}_{t,j}\geq(1-\varepsilon^{2})\dfrac{\sqrt{t}}{k }-\sigma_{-}|G_{j}|\geq(1-\varepsilon^{2})\dfrac{\sqrt{t}}{k}-\dfrac{4}{ \varepsilon}\sqrt{\log n}$ (73) $\displaystyle u_{i}<0$ $\displaystyle\Rightarrow\hat{v}_{t,j}\leq-(1-\varepsilon^{2})\dfrac{\sqrt{t}}{ k}+\sigma_{-}|G_{j}|\leq-(1-\varepsilon^{2})\dfrac{\sqrt{t}}{k}+\dfrac{4}{ \varepsilon}\sqrt{\log n}$ (74) $$

where we use a union bound to get $|G_{j}|\leq 4\sqrt{\log n}$ for all $j\leq n$. Similarly, for all $i\notin S$,

$$ $|\hat{v}_{t,j}|\leq\sigma_{-}|G_{j}|\leq\dfrac{4}{\varepsilon}\sqrt{\log n}$ $$

These calculations reveal that: (i) the entries with the largest magnitudes are the elements of $S$, and (ii) if $u_{i}$ and $\hat{v}_{t,i}$ share the same sign for all $i\in S$. On this event, $\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|=0$.

Lastly, we claim that the top eigenvalue of $\eta_{s}(\boldsymbol{A}_{t,+})$ is larger than $k+\sqrt{t}/s$.
From triangle inequality applied to Eq. ([59](https://arxiv.org/html/2503.08028v2#A6.E59)),
we have

$$ $\displaystyle\lambda_{1}(\eta_{s}(\boldsymbol{A}_{t,+}))$ $\displaystyle\geq(\sqrt{t}-kb_{0})-b_{1}-\|\eta_{s}(\sigma_{+}\boldsymbol{Z}) \|_{{\mbox{\rm\tiny op}}}-\|\boldsymbol{B}-\mathbb{E}[\boldsymbol{B}]\|_{{ \mbox{\rm\tiny op}}}$ (75) $\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}(\sqrt{t}-(1+\varepsilon)ks)- Ck-C\sqrt{k}s$ (76) $\displaystyle\stackrel{{\scriptstyle(b)}}{{\geq}}\sqrt{t}-(1+C_{0}\varepsilon)ks$ (77) $\displaystyle\stackrel{{\scriptstyle(c)}}{{\geq}}\varepsilon ks\,.$ (78) $$

Here $(a)$ follows from Eqs. ([66](https://arxiv.org/html/2503.08028v2#A6.E66)) and ([67](https://arxiv.org/html/2503.08028v2#A6.E67)), $(b)$ because $k\gg 1$ and $s\gg 1$
and $(c)$ follows by taking $\delta\geq C\varepsilon$ for $C$ a sufficiently large absolute constant.
The claim follows because $ks\gg k$ and $ks\gg\sqrt{t}$.

## Appendix G Auxiliary lemmas for Section C

### G.1 Proof of Lemma C.1

We let $\boldsymbol{B}\sim{\mathsf{N}}(\boldsymbol{0},\boldsymbol{I}_{n^{2}})$ so that $\boldsymbol{W}=(\boldsymbol{B}+\boldsymbol{B}^{{\sf T}})/2$. For $\boldsymbol{v}\in\Omega_{n,k}$
we have $\langle\boldsymbol{v},\boldsymbol{W}\boldsymbol{v}\rangle=\langle\boldsymbol{v
},\boldsymbol{B}\boldsymbol{v}\rangle\sim{\mathsf{N}}(0,1)$. We thus have, by Gaussian tail bounds and a triangle inequality:

$$ $\mathbb{P}\left(|\langle\boldsymbol{v},\boldsymbol{W}\boldsymbol{v}\rangle| \geq C\sqrt{\log{\binom{n}{k}}}\right)\leq 2\exp\left(-\dfrac{C^{2}}{2}\log \binom{n}{k}\right)\,.$ $$

Taking the union bound over $\boldsymbol{v}\in\Omega_{n,k}$ gives the desired statement, since the cardinality of this set is $\binom{n}{k}2^{k}$.

### G.2 Proof of Lemma C.3

Using Lemma [C.2](https://arxiv.org/html/2503.08028v2#A3.Thmclaim2), we can take $x=n^{-1/4}$, say, and $\theta=\sqrt{1+\delta}$ for $\delta\geq n^{-c_{0}}$ for some small enough $c_{0}>0$, to get that

$$ $\mathbb{P}\left(\lambda_{1}(\boldsymbol{y})\leq\theta+1/\theta-n^{-1/4}-2/n \right)\leq C\exp(-cn^{1/3})\,,$ $$

for some absolute constants $C,c>0$.

We have the following identity, letting $\boldsymbol{W}\sim{\sf GOE}(n,1/n)$:

$$ $\displaystyle\langle\boldsymbol{u},\boldsymbol{v}_{1}\rangle^{2}=\dfrac{1}{ \theta^{2}\langle\boldsymbol{u},(\lambda_{1}(\boldsymbol{y})\boldsymbol{I}- \boldsymbol{W})^{-2}\boldsymbol{u}\rangle}\geq\dfrac{1}{\theta^{2}\cdot\|( \lambda_{1}(\boldsymbol{y})\boldsymbol{I}-\boldsymbol{W})^{-2}\|_{{\mbox{\rm \tiny op}}}}\,.$ $$

By standard Gaussian concentration, we know that, for any $\Delta>0$

$$ $\mathbb{P}\left(\|\boldsymbol{W}\|_{{\mbox{\rm\tiny op}}}\geq 2+\Delta\right) \leq C\exp(-cn\Delta^{2})\,.$ $$

In this inequality, we take

$$ $\Delta=\frac{1}{4}\big{(}\theta+1/\theta-n^{-1/4}-2/n-2\big{)}\,.$ $$

Note that with $\theta=\sqrt{1+\delta}$ and $\delta=o_{n}(1)$, we know that $\theta+1/\theta-2=\Theta(\delta^{2})$, so that $\Delta=\Theta(\delta^{2})$ if $\delta\geq n^{-c_{0}}$ with $c_{0}\leq 1/8$. Hence, by a union bound on the two concentration inequalities,

$$ $\mathbb{P}\left(\lambda_{\min}(\lambda_{1}(\boldsymbol{y})\boldsymbol{I}- \boldsymbol{W})\leq 2\Delta\right)\leq C\exp(-cn^{1/3})$ $$

and on the complement of this event, we know that

$$ $\langle\boldsymbol{v}_{1},\boldsymbol{u}\rangle^{2}\geq\dfrac{4\Delta^{2}}{ \theta^{2}}=\Theta(\delta^{4})$ $$

since $\theta=\Omega(1)$, and so $|\langle\boldsymbol{v}_{1},\boldsymbol{u}\rangle|=\Omega(\delta^{2})$.

## Appendix H Proof of Proposition B.1

In our proof, we will use the following elementary facts.

###### Fact H.1 .

For any deterministic unit vector $\boldsymbol{u}$, a unit vector $\boldsymbol{v}$ is uniformly random on the orthogonal subspace to $\boldsymbol{u}$ if and only if $\langle\boldsymbol{v},\boldsymbol{u}\rangle=0$ and $\boldsymbol{v}\stackrel{{\scriptstyle{\rm d}}}{{=}}\boldsymbol{Q}\boldsymbol{v}$ for every orthogonal matrix $\boldsymbol{Q}$ such that $\boldsymbol{Q}\boldsymbol{u}=\boldsymbol{u}$.

###### Fact H.2 .

Let $\boldsymbol{A}$ be a symmetric matrix, and $\boldsymbol{u}$ a unit vector. Denote $\boldsymbol{B}_{\alpha}=\alpha\boldsymbol{u}\boldsymbol{u}^{{\sf T}}+
\boldsymbol{A}$, and let $\boldsymbol{v}(\alpha)$ be a top eigenvector of $\boldsymbol{B}_{\alpha}$. Then $f(\alpha)=|\langle\boldsymbol{v}(\alpha),\boldsymbol{u}\rangle|$ is an increasing function of $\alpha>0$.

Let $\boldsymbol{u}\sim\operatorname{Unif}(\Omega_{n,k})$.
Recall that $\boldsymbol{A}_{t}=\sqrt{t}\boldsymbol{u}\boldsymbol{u}^{{\sf T}}+\sqrt{t}
\boldsymbol{W}$
where $\boldsymbol{W}\sim{\sf GOE}(n,1/2)$.

We conduct our analysis conditional on $\boldsymbol{u}$. Let $\boldsymbol{v}_{t}$ be a top eigenvector of $\boldsymbol{A}_{t}$.
For $t=(1+\delta)n/2$, $|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|\overset{a.s.}{\to}\sqrt{
\delta/(1+\delta)}$,
so that with high probability $|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|\geq\sqrt{\delta}/(2\sqrt{1+
\delta})$. If $t\geq(1+\delta)n/2$, we can use Fact [H.2](https://arxiv.org/html/2503.08028v2#A8.Thmclaim2) to obtain the same result. By choosing $\varepsilon$ such that $2\varepsilon<\sqrt{\delta/(1+\delta)}$, we know from standard concentration of the alignment (Lemma [B.6](https://arxiv.org/html/2503.08028v2#A2.Thmclaim6)) that $|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|\geq 2\varepsilon$ with probability at least $1-\exp(-cn)$ for some $c>0$ possibly dependent on $(\varepsilon,\delta)$.

By rotational invariance of $\boldsymbol{W}$,
$\boldsymbol{w}_{t}:=\dfrac{\boldsymbol{v}_{t}-\langle\boldsymbol{v}_{t},
\boldsymbol{u}\rangle\boldsymbol{u}}{\|\boldsymbol{v}_{t}-\langle\boldsymbol{v
}_{t},\boldsymbol{u}\rangle\boldsymbol{u}\|}$ is uniformly random on the orthogonal subspace to $\boldsymbol{u}$. hence, there exists
$\boldsymbol{g}\sim{\mathsf{N}}(\boldsymbol{0},\boldsymbol{I}_{n})$, such that

$$ $\boldsymbol{w}_{t}\sim\dfrac{(\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^ {{\sf T}})\boldsymbol{g}}{\|(\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{ {\sf T}})\boldsymbol{g}\|}\,.$ $$

Since $\|(\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{{\sf T}})\boldsymbol{g}\|
\sim\|\boldsymbol{g}^{\prime}\|$ for some $\boldsymbol{g}^{\prime}\sim{\mathsf{N}}(\boldsymbol{0},\boldsymbol{I}_{n-1})$,
we have, for some constant $c>0$,

$$ $\mathbb{P}\left(\|(\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{{\sf T}}) \boldsymbol{g}\|\leq\dfrac{\sqrt{n}}{2}\right)\leq\exp\left(-cn\right)\,.$ $$

Further, for every $1\leq i\leq n$, we know that

$$ $|((\boldsymbol{I}_{n}-\boldsymbol{u}\boldsymbol{u}^{{\sf T}})\boldsymbol{g})_{ i}|\leq|g_{i}|+\|\boldsymbol{u}\|_{\infty}\cdot|\langle\boldsymbol{u}, \boldsymbol{g}\rangle|\leq|g_{i}|+\dfrac{|\langle\boldsymbol{u},\boldsymbol{g} \rangle|}{\sqrt{k}}\,.$ $$

We next show that, with the claimed probability, only a few entries of $\boldsymbol{w}_{t}$ can have large magnitude. As a result,
less than $\ell$ entries of $\boldsymbol{u}$ can be estimated incorrectly (with $\ell=1$
if $k\ll n/\log n$).

Define $\ell=\lceil n\exp(-a_{n}\cdot n/k)\rceil\geq 1$ (with $a_{n}$
a sequence to be chosen later) and $g^{\text{abs}}_{(\ell)}$ as the $\ell$-th largest value among the $|g_{i}|$’s. We have

$$ $\mathbb{P}\left(|\langle\boldsymbol{u},\boldsymbol{g}\rangle|\geq\sqrt{na_{n}} \right)\leq\exp\left(-\dfrac{na_{n}}{2}\right)\,.$ $$

Furthermore, from a union bound, we get that

$$ $\displaystyle\mathbb{P}\left(g_{(\ell)}^{\text{abs}}\geq\dfrac{2\sqrt{na_{n}}} {\sqrt{k}}\right)$ $\displaystyle\leq{\binom{n}{\ell}}\cdot\exp\left(-\dfrac{2n\ell\cdot a_{n}}{k}\right)$ $\displaystyle\leq\left(\dfrac{en}{\ell}\right)^{\ell}\exp\left(-\dfrac{2n\ell \cdot a_{n}}{k}\right)$ $\displaystyle=\exp\left(-\dfrac{2n\ell\cdot a_{n}}{k}+\ell\log\dfrac{n}{\ell}+ \ell\right)\,.$ $$

By definition, we know that $\ell\geq\max\left\{n\exp(-a_{n}\cdot n/k),1\right\}$, so that

$$ $\dfrac{2na_{n}}{k}-\log\dfrac{n}{\ell}-1\geq\dfrac{2na_{n}}{k}-\min\left\{\log n ,a_{n}\cdot n/k\right\}-1\geq\dfrac{na_{n}}{2k}$ $$

as long as $na_{n}\gg k$. This means that

$$ $\mathbb{P}\left(g_{(\ell)}^{\text{abs}}\geq\dfrac{2\sqrt{na_{n}}}{\sqrt{k}} \right)\leq\exp\left(-\dfrac{na_{n}}{2k}\right)\,.$ $$

Define the set

$$ $\displaystyle{\mathcal{A}}_{n}(t):=\Big{\{}i\leq n:\,|w_{ti}|\geq 6\sqrt{\frac {a_{n}}{k}}\Big{\}}\,.$ $$

By the bounds above we have

$$ $\displaystyle\mathbb{P}\left(|{\mathcal{A}}_{n}(t)|\leq\ell-1\right)\geq 1-e^{ -cn}-e^{-na_{n}/2k}-e^{-na_{n}/2}\,.$ $$

Suppose that $|\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle|\geq 2\varepsilon$ also holds, and suppose without loss of generality that $\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle\geq 2\varepsilon$. Then, we have (as long as $6\sqrt{a_{n}}\leq(9/10)\varepsilon$)

$$ $\displaystyle i\in S,i\notin{\mathcal{A}}_{n}(t)\,u_{i}>0\Rightarrow v_{ti} \geq\dfrac{\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle}{\sqrt{k}}-\dfrac{6 \sqrt{a_{n}}}{\sqrt{k}}>\dfrac{\varepsilon}{\sqrt{k}}\Rightarrow i\in\hat{S}, \operatorname{sign}(v_{ti})>0,,$ $\displaystyle i\in S,i\notin{\mathcal{A}}_{n}(t),u_{i}<0\Rightarrow v_{ti}\leq -\dfrac{\langle\boldsymbol{v}_{t},\boldsymbol{u}\rangle}{\sqrt{k}}+\dfrac{6 \sqrt{a_{n}}}{\sqrt{k}}<-\dfrac{\varepsilon}{\sqrt{k}}\Rightarrow i\in\hat{S}, \operatorname{sign}(v_{ti})<0\,.$ $$

Analogously, for $i\notin S,i\notin{\mathcal{A}}_{n}(t)$, we have

$$ $|v_{ti}|\leq\dfrac{6\sqrt{a_{n}}}{\sqrt{k}}<\dfrac{\varepsilon}{\sqrt{k}}\,.$ $$

and we obtain that at most $\ell-1$ positions could be mis-identified.

Next, we show that the termination condition (Line 5, Algorithm [1](https://arxiv.org/html/2503.08028v2#alg1)) does not trigger for each $t\geq n^{2}$ (with high probability). We write

$$ $\boldsymbol{A}_{t}=\dfrac{\boldsymbol{y}_{t}+\boldsymbol{y}_{t}^{{\sf T}}}{2 \sqrt{t}}=\sqrt{t}\boldsymbol{u}\boldsymbol{u}^{{\sf T}}+\left(\dfrac{ \boldsymbol{B}_{t}+\boldsymbol{B}_{t}^{{\sf T}}}{2\sqrt{t}}\right)$ $$

From Weyl’s inequality:

$$ $\lambda_{1}(\boldsymbol{A}_{t})\geq\sqrt{t}-\left\|\dfrac{\boldsymbol{B}_{t}+ \boldsymbol{B}_{t}^{{\sf T}}}{2\sqrt{t}}\right\|_{\text{op}}$ $$

From standard operator norm results for GOE matrices (as $(\boldsymbol{B}_{t}+\boldsymbol{B}_{t}^{{\sf T}})/\sqrt{2t}\sim{\sf GOE}(n)$), we know that $\|(\boldsymbol{B}_{t}+\boldsymbol{B}_{t}^{{\sf T}})/(2\sqrt{t})\|_{\text{op}}
\leq 2\sqrt{n}$ with probability at least $1-\exp(-cn)$, for some $c>0$. Hence $\lambda_{1}(\boldsymbol{A}_{t})\geq\sqrt{t}-2\sqrt{n}>\sqrt{t}/2$ as $t\geq n^{2}\gg n$.

We obtain that

$$ $\mathbb{E}\left[\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)-\boldsymbol{x}\|^ {2}\right]=O\left(\dfrac{\ell-1}{k}+\exp\left(-\dfrac{na_{n}}{k}\right)\right) =O\left(\exp(-n\varepsilon^{2}/64k)\right)$ $$

where we picked $a_{n}>\varepsilon^{2}/64$ satisfying the bounds outlined above, namely $(i)$ $6\sqrt{a_{n}}\leq 0.9\varepsilon$, and $(ii)$ $na_{n}\gg k$. Notice that $na_{n}/k\gg\log(na_{n}/k)$ if $na_{n}\gg k$, and $\ell-1\leq n\cdot\exp(-a_{n}\cdot n/k)$.

###### Fact H.1 .

###### Fact H.2 .

## Appendix I Proof of Lemma D.3

###### Proof.

By the Markov Property, we know that $\max_{t_{\ell}\leq t\leq t_{\ell+1}}\|\boldsymbol{w}_{t}-\boldsymbol{w}_{t_{
\ell}}\|_{F}\overset{d}{=}\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{
w}_{t}\|_{F}$. By Gaussian concentration, we have

$$ $\mathbb{P}\left(\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{w}_{t}\|_{ F}-\mathbb{E}\left[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{w}_{t} \|_{F}\right]\geq x\right)\leq 2\exp\left(-\dfrac{x^{2}}{4(t_{\ell+1}-t_{\ell} )}\right)$ $$

This can be proven, e.g. by discretizing the interval $[0,t_{\ell+1}-t_{\ell}]$ into $r$ equal-length intervals and employing standard Gaussian concentration on vectors (then pushing $r\to\infty$). As the argument is standard, we omit the proof for brevity.

Now we bound $\mathbb{E}\left[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{w}_{t}\|_{
F}\right]$. We know that $\|\boldsymbol{w}_{t}\|_{F}$ is a non-negative submartingale, so that from Doob’s inequality:

$$ $\mathbb{E}\left[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{w}_{t}\|_{ F}^{2}\right]\leq 4\mathbb{E}[\|\boldsymbol{w}_{t_{\ell+1}-t_{\ell}}\|_{F}^{2} ]\leq 9(t_{\ell+1}-t_{\ell})n^{2}$ $$

so that from Cauchy-Schwarz, $\mathbb{E}\left[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{w}_{t}\|_{
F}\right]\leq 3\sqrt{t_{\ell+1}-t_{\ell}}n$. Hence as $t_{\ell}\geq 1$,

$$ $\mathbb{P}\left(\max_{t_{\ell}\leq t\leq t_{\ell+1}}\|\boldsymbol{w}_{t}- \boldsymbol{w}_{t_{\ell}}\|_{F}\geq 4\sqrt{(t_{\ell+1}-t_{\ell})t_{\ell}}\cdot n \right)\leq 2\exp\left(-\dfrac{n^{2}t_{\ell}}{4}\right)$ $$

The second tail bound follows immediately (at least, the proof would be analogous to the preceding display).
∎

###### Proof.

## Appendix J Proof of Lemma B.3

By Gaussian concentration, we have

$$ $\mathbb{P}\left(\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{W}_{t}\|_{ \text{op}}-\mathbb{E}\left[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol {W}_{t}\|_{\text{op}}\right]\geq x\right)\leq 2\exp\left(-\dfrac{x^{2}}{2(t_{ \ell+1}-t_{\ell})}\right)$ $$

This can be proven, e.g. by discretizing the interval $[0,t_{\ell+1}-t_{\ell}]$ into $r$ equal-length intervals and employing Gaussian concentration on vectors (then pushing $r\to\infty$). As the argument is standard, we omit the proof for brevity.

To evaluate $\mathbb{E}[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{W}_{t}\|_{\text
{op}}]$, we recognize that $\|\boldsymbol{W}_{t}\|_{\text{op}}$ is a submartingale, so that from Doob’s inequality:

$$ $\mathbb{E}\left[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{W}_{t}\|_{ \text{op}}^{2}\right]\leq 4\mathbb{E}\left[\|\boldsymbol{W}_{t_{\ell+1}-t_{ \ell}}\|_{\text{op}}^{2}\right]$ $$

Once again from Gaussian concentration,

$$ $\mathbb{P}\left(|\|\boldsymbol{W}_{1}\|_{\text{op}}-\mathbb{E}[\|\boldsymbol{W }_{1}\|_{\text{op}}]|\geq x\right)\leq 2\exp\left(\dfrac{-x^{2}}{2}\right)$ $$

so that $\mathbb{P}\left(|\|\boldsymbol{W}_{t_{\ell+1}-t_{\ell}}\|_{\text{op}}-\mathbb{
E}[\|\boldsymbol{W}_{t_{\ell+1}-t_{\ell}}\|_{\text{op}}]|\geq x\right)\leq 2
\exp\left(-x^{2}/(2(t_{\ell+1}-t_{\ell}))\right)$. Hence $\|\boldsymbol{W}_{t_{\ell+1}-t_{\ell}}\|_{\text{op}}$ is $(t_{\ell+1}-t_{\ell})$-subgaussian, implying that $\text{Var}(\|\boldsymbol{W}_{t_{\ell+1}-t_{\ell}}\|_{\text{op}})\leq 6(t_{\ell
+1}-t_{\ell})$. As $\mathbb{E}[\|\boldsymbol{W}_{t_{\ell+1}-t_{\ell}}\|_{\text{op}}]^{2}\sim 4(t_{
\ell+1}-t_{\ell})n$ (one can obtain this from the Bai-Yin Theorem along with sub-gaussianity, for instance), we get that $\mathbb{E}\left[\|\boldsymbol{W}_{t_{\ell+1}-t_{\ell}}\|^{2}\right]\leq 16(t_{
\ell+1}-t_{\ell})n$ eventually as $n$ gets large.

From Cauchy-Schwarz inequality, we get that

$$ $\mathbb{E}\left[\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{W}_{t}\|_{ \text{op}}\right]\leq 8\sqrt{t_{\ell+1}-t_{\ell}}\sqrt{n}$ $$

We conclude that

$$ $\mathbb{P}\left(\max_{0\leq t\leq t_{\ell+1}-t_{\ell}}\|\boldsymbol{W}_{t}\|_{ \text{op}}\geq 16\sqrt{(t_{\ell+1}-t_{\ell})n}\right)\leq 2\exp\left(-32n\right)$ $$

## Appendix K Proof of Lemma B.4

First, by orthogonal invariance of $\boldsymbol{W}_{t}$, we know that $\boldsymbol{v}_{t}$ is uniformly random over the unit sphere $\mathbb{S}^{n-1}$. We can write, using $\boldsymbol{g}\sim{\mathsf{N}}(0,\boldsymbol{I}_{n})$, the following representation

$$ $\boldsymbol{v}_{t}\sim\dfrac{\boldsymbol{g}}{\|\boldsymbol{g}\|}$ $$

As in the statement of the Lemma, we define the following set, for $\boldsymbol{v}\in\mathbb{R}^{n}$ and $C>0$:

$$ $A(\boldsymbol{v};C)=\left\{i:1\leq i\leq n,|v_{i}|\geq\dfrac{C\sqrt{\log(n/k)} }{\sqrt{n}}\right\}$ $$

As with the proof of Proposition [B.1](https://arxiv.org/html/2503.08028v2#A2.Thmclaim1), we first deal with the denominator $\|\boldsymbol{g}\|$: indeed, sub-exponential concentration gives us

$$ $\displaystyle\mathbb{P}\left(\sum_{j=1}^{n}g_{j}^{2}\leq\dfrac{n}{2}\right)$ $\displaystyle\leq 2\exp(-n/8)$ (79) $$

This leads us to define another set

$$ $B(\boldsymbol{g};C)=\left\{i:1\leq i\leq n,|g_{i}|\geq C\sqrt{\log(n/k)}\right\}$ $$

Let $p_{n}=\mathbb{P}(|g_{1}|\geq C\sqrt{\log(n/k)})$, then we have $|B(\boldsymbol{g};C)|\sim\text{Bin}(n,p_{n})$. From Gaussian tail bounds, we know that $p_{n}\leq(n/k)^{-C^{2}/2}$. We now use a Chernoff bound of the following form: for every $x\geq 4\mathbb{E}[X]$, where $X\sim\text{Bin}(n,p)$, then

$$ $\mathbb{P}\left(X\geq x\right)\leq\exp\left(-x/3\right)$ $$

It is clear that $np_{n}\ll k^{2}/n\leq\max\{k^{2}/n,\sqrt{k}\}$ when $C>2$, so that we have

$$ $\mathbb{P}\left(|B(\boldsymbol{g};C)|\geq\max\{\sqrt{k},k^{2}/n\}\right)\leq \exp\left(-\dfrac{1}{3}\max\{\sqrt{k},k^{2}/n\}\right)\leq\exp\left(-\dfrac{1} {3}n^{1/4}\right)$ $$

Therefore, with each fixed $t$, by union bound with probability at least $1-O(\exp(-\sqrt{n}))$, we have, for a possibly different $C>0$, $|A(\boldsymbol{v}_{t};C)|\leq\max\{\sqrt{k},k^{2}/n\}$. Our proof ends here, as $\max\{\sqrt{k},k^{2}/n\}\ll k/2$ for $\sqrt{n}\ll k\ll n$.

## Appendix L Proof of Lemma B.5

We know that

$$ $\displaystyle\boldsymbol{v}_{t}^{{\sf T}}\boldsymbol{W}_{t_{\ell}}\boldsymbol{ v}_{t}=$ $\displaystyle\boldsymbol{v}_{t}^{{\sf T}}\boldsymbol{W}_{t}\boldsymbol{v}_{t}- \boldsymbol{v}_{t}^{{\sf T}}(\boldsymbol{W}_{t}-\boldsymbol{W}_{t_{\ell}}) \boldsymbol{v}_{t}$ $\displaystyle=$ $\displaystyle\lambda_{1}(\boldsymbol{W}_{t})-\boldsymbol{v}_{t}^{{\sf T}}( \boldsymbol{W}_{t}-\boldsymbol{W}_{t_{\ell}})\boldsymbol{v}_{t}$ $\displaystyle=$ $\displaystyle\lambda_{1}(\boldsymbol{W}_{t_{\ell}})-\boldsymbol{v}_{t}^{{\sf T }}(\boldsymbol{W}_{t}-\boldsymbol{W}_{t_{\ell}})\boldsymbol{v}_{t}+(\lambda_{1 }(\boldsymbol{W}_{t})-\lambda_{1}(\boldsymbol{W}_{t_{\ell}}))$ $$

from which we obtain from Weyl’s inequality that

$$ $\sup_{t_{\ell}\leq t\leq t_{\ell+1}}\left|\boldsymbol{v}_{t}^{{\sf T}} \boldsymbol{W}_{t_{\ell}}\boldsymbol{v}_{t}-\lambda_{1}(\boldsymbol{W}_{t_{ \ell}})\right|\leq 2\sup_{t_{\ell}\leq t\leq t_{\ell+1}}\|\boldsymbol{W}_{t}- \boldsymbol{W}_{t_{\ell}}\|_{\text{op}}\leq 32\sqrt{(t_{\ell+1}-t_{\ell})n}$ $$

with probability at least $1-2\exp(-32n)$.

## Appendix M Proof of Lemma B.6

By Weyl’s inequality, $\boldsymbol{W}\mapsto\lambda_{1}(\boldsymbol{Y})$ (with $\boldsymbol{Y}=\theta\boldsymbol{v}\boldsymbol{v}^{{\sf T}}+\boldsymbol{W}$)
is a 1-Lipschitz function and therefore, by Borell inequality
(and Baik et al. [2005]), letting $\lambda_{*}(\theta):=\theta+1/\theta$,
for any $\varepsilon>0$,

$$ $\displaystyle\mathbb{P}\big{(}|\lambda_{1}(\boldsymbol{Y})-\lambda_{*}(\theta) |\geq\varepsilon\big{)}\leq 2e^{-n\varepsilon^{2}/4}\,.$ (80) $$

To prove concentration of $\langle\boldsymbol{v}_{1}(\boldsymbol{Y}),\boldsymbol{v}\rangle^{2}$, note that simple linear algebra yields

$$ $\displaystyle\frac{1}{\langle\boldsymbol{v}_{1}(\boldsymbol{Y}),\boldsymbol{v} \rangle^{2}}=\langle\boldsymbol{v},(\lambda_{1}(\boldsymbol{Y})\boldsymbol{I}- \boldsymbol{W})^{-2}\boldsymbol{v}\rangle=:F(\boldsymbol{W})\,.$ (81) $$

It is therefore sufficient to prove that $F(\boldsymbol{W})$ concentrates around a value that is bounded away from 0 0.
Fix $\varepsilon_{0}>0$ such that $2+3\varepsilon_{0}<\lambda_{*}(\theta)$ and define the event

$$ $\displaystyle\mathcal{E}:=\big{\{}\boldsymbol{W}:\,\|\boldsymbol{W}\|_{{\mbox{ \rm\tiny op}}}\leq 2+\varepsilon_{0}\,,\;|\lambda_{1}(\boldsymbol{Y})-\lambda_ {*}|\leq\varepsilon_{0}\big{\}}\,.$ (82) $$

By the Bai-Yin law and Gaussian concentration (plus the above concentration of $\lambda_{1}$), $\mathbb{P}(\mathcal{E})\geq 1-2e^{-c(\varepsilon_{0})n}$ for some $c(\varepsilon_{0})>0$.
Further, it is easy to check that $F(\boldsymbol{W})$ is Lipschitz on $\mathcal{E}$,
whence the concentration of $\langle\boldsymbol{u},\boldsymbol{v}_{1}(\boldsymbol{W})\rangle^{2}$ follows by another application of Borell inequality.

## Appendix N Proofs of reduction results

### N.1 Proof of Theorem 4

The algorithm consists in running the discretized diffusion ([4](https://arxiv.org/html/2503.08028v2#S1.E4))
with initialization $\hat{\boldsymbol{y}}_{t_{0}}=\boldsymbol{y}/\sigma^{2}$ at $t=t_{0}:=1/\sigma^{2}$. To avoid notational burden,
we will assume $(T-t_{0})/\Delta$ to be an integer. Let $\hat{\boldsymbol{y}}_{t_{0}}^{*}$ be generated by the
discretized diffusion with initialization at $\hat{\boldsymbol{y}}_{0}$ at $t=0$. Note that the
distribution of $\hat{\boldsymbol{y}}_{t_{0}}$ is the same as the one of $t_{0}\boldsymbol{x}+\sqrt{t}\boldsymbol{g}$ and hence by Assumption
$(b)$,
and Pinsker’s inequality

$$ $\displaystyle\|{\rm P}_{\hat{\boldsymbol{y}}_{t_{0}}}-{\rm P}_{\hat{ \boldsymbol{y}}^{*}_{t_{0}}}\|_{\mbox{\tiny\rm TV}}\leq\sqrt{\frac{1}{2}D_{ \mbox{\tiny\rm KL}}({\rm P}_{\hat{\boldsymbol{y}}^{*}_{t_{0}}}\|{\rm P}_{\hat{ \boldsymbol{y}}_{t_{0}}})}\leq\sqrt{\frac{1}{2}D_{\mbox{\tiny\rm KL}}( \overline{\rm P}_{\hat{\boldsymbol{y}}}^{T,\Delta}\|{\rm P}_{\boldsymbol{y}}^{ T})}\leq\sqrt{\frac{\varepsilon}{2}}\,.$ (83) $$

Hence $\hat{\boldsymbol{y}}_{t_{0}}$, $\hat{\boldsymbol{y}}^{*}_{t_{0}}$ can be coupled so that $\mathbb{P}(\hat{\boldsymbol{y}}_{t_{0}}\neq\hat{\boldsymbol{y}}^{*}_{t_{0}})
\leq\sqrt{\varepsilon/2}$.

We extend this to a coupling
of $(\hat{\boldsymbol{y}}^{*}_{t})_{t_{0}\leq t\leq T}$ and $(\hat{\boldsymbol{y}}_{t})_{t_{0}\leq t\leq T}$ in the obvious way:
we generate the two trajectories according to the discretized diffusion ([4](https://arxiv.org/html/2503.08028v2#S1.E4))
with the same randomness $\hat{\boldsymbol{z}}_{t}$.
Therefore $\mathbb{P}(\hat{\boldsymbol{y}}_{T}\neq\hat{\boldsymbol{y}}^{*}_{T})\leq\sqrt{
\varepsilon/2}$.
Another application of the assumption $D_{\mbox{\tiny\rm KL}}(\overline{\rm P}_{\hat{\boldsymbol{y}}}^{T,\Delta}\|{
\rm P}_{\boldsymbol{y}}^{T})\leq\varepsilon$
and Pinsker’s inequality yields $\mathbb{P}(\boldsymbol{y}_{T}\neq\hat{\boldsymbol{y}}^{*}_{T})\leq\sqrt{
\varepsilon/2}$, for $\boldsymbol{y}_{T}\stackrel{{\scriptstyle{\rm d}}}{{=}}T\boldsymbol{x}+\sqrt{T
}\boldsymbol{g}^{\prime}$ with $(\boldsymbol{x},\boldsymbol{g}^{\prime})\sim\mu\otimes{\mathsf{N}}(\boldsymbol
{0},\boldsymbol{I})$.
We conclude by triangle inequality $\mathbb{P}(\boldsymbol{y}_{T}\neq\hat{\boldsymbol{y}}_{T})\leq 2\sqrt{
\varepsilon/2}$, which coincides with the claim ([13](https://arxiv.org/html/2503.08028v2#S4.E13)).

Finally, Eq. ([14](https://arxiv.org/html/2503.08028v2#S4.E14)) follows by generating $N$ i.i.d. copies $\hat{\boldsymbol{x}}_{1},\dots,\hat{\boldsymbol{x}}_{N}$ using the above procedure, and letting $\hat{\boldsymbol{m}}(\boldsymbol{y})$ be their empirical average.

### N.2 Proof of Theorem 5

The algorithm consists in running the discretized diffusion ([4](https://arxiv.org/html/2503.08028v2#S1.E4))
with initialization $\hat{\boldsymbol{y}}_{t_{0}}=\boldsymbol{y}/\sigma^{2}$ at $t=t_{0}:=1/\sigma^{2}$. To avoid notational burden,
we will assume $(T-t_{0})/\Delta$ to be an integer. Let $\hat{\boldsymbol{y}}_{t_{0}}^{*}$ be generated by the
discretized diffusion with initialization at $\hat{\boldsymbol{y}}_{0}$ at $t=0$. Note that the
distribution of $\hat{\boldsymbol{y}}_{t_{0}}$ is the same as the one of $t_{0}\boldsymbol{x}+\sqrt{t}\boldsymbol{g}$ and hence by Assumption
$(b)$,

$$ $\displaystyle W_{1}({\rm P}_{\hat{\boldsymbol{y}}_{t_{0}}},{\rm P}_{\hat{ \boldsymbol{y}}^{*}_{t_{0}}})\leq W_{1}({\rm P}^{\hat{\boldsymbol{y}}}_{T, \Delta},{\rm P}^{\boldsymbol{y}}_{T})\leq\varepsilon\,.$ (84) $$

In other words there exists a coupling of $\hat{\boldsymbol{y}}^{*}_{t_{0}}$ and $\hat{\boldsymbol{y}}_{t_{0}}$ such that
$\mathbb{E}\|\hat{\boldsymbol{y}}^{*}_{t_{0}}-\hat{\boldsymbol{y}}_{t_{0}}\|_{2
}\leq\varepsilon$.

We extend this to a coupling
of $(\hat{\boldsymbol{y}}^{*}_{t})_{t_{0}\leq t\leq T}$ and $(\hat{\boldsymbol{y}}_{t})_{t_{0}\leq t\leq T}$ in the obvious way:
we generate the two trajectories according to the discretized diffusion ([4](https://arxiv.org/html/2503.08028v2#S1.E4))
with the same randomness $\hat{\boldsymbol{z}}_{t}$. A simple recursive argument (using the Lipschitz property of $\hat{\boldsymbol{m}}$,
in Assumption $(a)$)
then yields

$$ $\displaystyle\mathbb{E}\|\hat{\boldsymbol{y}}^{*}_{T}-\hat{\boldsymbol{y}}_{T} \|_{2}\leq\Big{(}1+L\Delta/d\Big{)}^{T/\Delta}\varepsilon\leq e^{LT/d} \varepsilon\,.$ (85) $$

(See for instance Montanari and Wu [2023] or Alaoui et al. [2023] for examples of this calculation.)
Let now $\boldsymbol{y}_{T}\stackrel{{\scriptstyle{\rm d}}}{{=}}T\boldsymbol{x}+\sqrt{T
}\boldsymbol{g}^{\prime}$ for $(\boldsymbol{x},\boldsymbol{g}^{\prime})\sim\mu\otimes{\mathsf{N}}(\boldsymbol
{0},\boldsymbol{I})$.
Another application of Assumption $(a)$ implies that this can be coupled to $\hat{\boldsymbol{y}}^{*}_{T}$ so that
$\mathbb{E}\|\boldsymbol{y}_{T}-\hat{\boldsymbol{y}}^{*}_{T}\|\leq\varepsilon$, and therefore

$$ $\displaystyle\mathbb{E}\|\hat{\boldsymbol{y}}_{T}-\boldsymbol{y}_{T}\|_{2}\leq 2 \,e^{LT/d}\varepsilon\,.$ (86) $$

As output, we return $\hat{\boldsymbol{x}}=\hat{\boldsymbol{y}}_{T}/T$. Using $\mathbb{E}\|\boldsymbol{y}_{T}-\boldsymbol{x}\|=\mathbb{E}\|\boldsymbol{g}\|/
\sqrt{T}$ and $T=\theta d$,

$$ $\displaystyle\mathbb{E}\|\boldsymbol{x}-\hat{\boldsymbol{x}}\|\leq 2e^{\theta L }\varepsilon+\frac{1}{\sqrt{\theta}}\,.$ (87) $$

Since the coupling has been constructed conditionally on $\boldsymbol{y}$, the claim
([15](https://arxiv.org/html/2503.08028v2#S4.E15)) follows.

Finally, Eq. ([14](https://arxiv.org/html/2503.08028v2#S4.E14)) follows by generating $N$ i.i.d. copies $\hat{\boldsymbol{x}}_{1},\dots,\hat{\boldsymbol{x}}_{N}$ using the above procedure, and letting $\hat{\boldsymbol{m}}(\boldsymbol{y})$ be their empirical average.

## Appendix O Proof of Lemma C.1

We let $\boldsymbol{B}\sim{\mathsf{N}}(\boldsymbol{0},\boldsymbol{I}_{n^{2}})$ so that $\boldsymbol{W}\sim(\boldsymbol{B}+\boldsymbol{B}^{\intercal})/2$. We consider a non-random vector $\boldsymbol{v}$ fitting the description, and note that

$$ $|\langle\boldsymbol{v},\boldsymbol{W}\boldsymbol{v}\rangle|\leq\dfrac{1}{2} \left\{|\langle\boldsymbol{v},\boldsymbol{B}\boldsymbol{v}\rangle|+|\langle \boldsymbol{v},\boldsymbol{B}^{\intercal}\boldsymbol{v}\rangle|\right\}$ $$

We know that $\langle\boldsymbol{v},\boldsymbol{B}\boldsymbol{v}\rangle,\langle\boldsymbol{v
},\boldsymbol{B}^{\intercal}\boldsymbol{v}\rangle\sim{\mathsf{N}}(0,1)$. We thus have, by Gaussian tail bounds and a triangle inequality:

$$ $\mathbb{P}\left(|\langle\boldsymbol{v},\boldsymbol{W}\boldsymbol{v}\rangle| \geq C\sqrt{\log{\binom{n}{k}}}\right)\leq 2\exp\left(-\dfrac{C^{2}}{2}\log \binom{n}{k}\right)$ $$

Union bounding over the set of all such vectors gives us the desired statement, as the cardinality of this set is $\binom{n}{k}2^{k}$.

## Appendix P Proof of Lemma C.3

Proof of Lemma [C.3](https://arxiv.org/html/2503.08028v2#A3.Thmclaim3): Using Lemma [C.2](https://arxiv.org/html/2503.08028v2#A3.Thmclaim2), we can take $x=n^{-1/4}$, say, and $\theta=\sqrt{1+\delta}$ for $\delta=o_{n}(1)$, to get that

$$ $\mathbb{P}\left(\lambda_{1}(\boldsymbol{y})\leq\theta+1/\theta-n^{-1/4}-2/n \right)\leq C\exp(-cn^{1/2})$ $$

for some absolute constants $C,c>0$.

We have the following identity, letting $\boldsymbol{W}\sim{\sf GOE}(n,1/n)$:

$$ $\displaystyle\langle\boldsymbol{u},\boldsymbol{v}_{1}\rangle^{2}=\dfrac{1}{ \theta^{2}\langle\boldsymbol{u},(\lambda_{1}(\boldsymbol{y})\boldsymbol{I}- \boldsymbol{W})^{-2}\boldsymbol{u}\rangle}\geq\dfrac{1}{\theta^{2}\cdot\|( \lambda_{1}(\boldsymbol{y})\boldsymbol{I}-\boldsymbol{W})^{-2}\|_{{\mbox{\rm \tiny op}}}}$ $$

By standard Gaussian concentration, we know that

$$ $\mathbb{P}\left(\|\boldsymbol{W}\|_{{\mbox{\rm\tiny op}}}\geq 2+x\right)\leq C \exp(-cnx^{2})$ $$

We take

$$ $x=\dfrac{\theta+1/\theta-n^{-1/4}-2/n-2}{4}$ $$

Note that with $\theta=\sqrt{1+\delta}$ and $\delta=o_{n}(1)$, we know that $\theta+1/\theta-2=\Theta(\delta^{2})$, so that $x=\Theta(\delta^{2})$ if $\delta\gg n^{-1/8}$. Furthermore we have, by Theorem 1 above,

$$ $\mathbb{P}\left(\lambda_{\min}(\lambda_{1}(\boldsymbol{y})\boldsymbol{I}- \boldsymbol{W})\leq 2x\right)\leq C\exp(-cn^{1/2})$ $$

and on the complement of this event, we know that

$$ $\langle\boldsymbol{v}_{1},\boldsymbol{u}\rangle^{2}\geq\dfrac{4x^{2}}{\theta^{ 2}}=\Theta(\delta^{4})$ $$

since $\theta=\Omega(1)$, and so $|\langle\boldsymbol{v}_{1},\boldsymbol{u}\rangle|=\Omega(\delta^{2})$. Hence we are done.

## Appendix Q Proof of Theorem 6

In order to simplify some of the formulas below we
center $\mu_{n,k}$. Namely, we redefine
$\mu_{n,k}$ to be the distribution of $\boldsymbol{x}=\boldsymbol{u}\boldsymbol{u}^{{\sf T}}-\mathbb{E}[\boldsymbol{u
}\boldsymbol{u}^{{\sf T}}]$
when $\boldsymbol{u}\sum\operatorname{Unif}(B_{n,k})$.

Throughout this proof, $C$ denotes a generic constant which depends on the constants in the assumptions, and is allowed to change from line to line.
We will write $\mathbb{E}_{n,k}$ for expectation under $\mu_{n,k}$
and $\mathbb{E}$ for expectation under $\overline{\mu}_{n,k}=\frac{1}{2}\mu_{n,k}+\frac{1}{2}\delta_{\boldsymbol{0}}$.
Further $\boldsymbol{y}_{t}=t\boldsymbol{x}+\boldsymbol{W}_{t}$, where the distribution of $\boldsymbol{x}$ is either $\mu_{n,k}$
or $\overline{\mu}_{n,k}$ as indicated.

The optimality of $\hat{\boldsymbol{m}}_{n}$ with respect to scalings $c\hat{\boldsymbol{m}}_{n}$ implies
by Pythagora’s theorem

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t )-\boldsymbol{x}\|^{2}\big{\}}=\mathbb{E}\{\|\boldsymbol{x}\|^{2}\}-\mathbb{E} \big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)\|^{2}\big{\}}\,,$ $$

whence, using $\mathbb{E}\{\|\boldsymbol{x}\|^{2}\}=1/2$ and the assumption ([17](https://arxiv.org/html/2503.08028v2#S5.E17)),
we obtain (for any polytime test $\phi$)

$$ $\displaystyle\varepsilon_{n}:=\sup_{t\leq(1-\delta)t_{\mbox{\tiny\rm alg}}} \Big{[}\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)\|^{2 }\big{\}}\vee\Big{|}\mathbb{E}\{\phi(\boldsymbol{y}_{t})\}-\mathbb{E}\{\phi( \boldsymbol{W}_{t})\}\Big{|}\Big{]}=o_{n}(n^{-1})\,.$ (88) $$

Since $\hat{\boldsymbol{m}}_{n}$ has no larger MSE than $\boldsymbol{P}_{{\sf B}}\hat{\boldsymbol{m}}_{n}$,
and $\|\boldsymbol{x}\|\leq 1$ almost surely, we can assume
without loss of generality $\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)\|\leq 1$ almost surely.
Hence Eq. ([88](https://arxiv.org/html/2503.08028v2#A17.E88)) implies

$$ $\displaystyle\varepsilon_{n}$ $\displaystyle\geq\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{ t},t)\|^{2}\big{\}}$ $\displaystyle=\frac{1}{2}\mathbb{E}_{n,k}\big{\{}\|\hat{\boldsymbol{m}}_{n}( \boldsymbol{y}_{t},t)\|^{2}\big{\}}+\frac{1}{2}\mathbb{E}\big{\{}\|\hat{ \boldsymbol{m}}_{n}(\boldsymbol{W}_{t},t)\|^{2}\big{\}}$ $\displaystyle\geq\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{W}_{ t},t)\|^{2}\big{\}}-\varepsilon_{n}\,,$ $$

and therefore

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{W}_{t},t )\|^{2}\big{\}}\leq 2\varepsilon_{n}\,.$ (89) $$

We write ${\sf KL}(\boldsymbol{X}\|\boldsymbol{Y})$ for the Kullback-Leiber distance between the laws of random vectors $\boldsymbol{X}$ and $\boldsymbol{Y}$, and ${\sf TV}(\boldsymbol{X},\boldsymbol{Y})$ for the corresponding total variation distance.
Recall that $\hat{\boldsymbol{y}}_{t}$ is defined by Eq. ([4](https://arxiv.org/html/2503.08028v2#S1.E4)), where will let $\sqrt{\Delta}\boldsymbol{z}_{t}=\boldsymbol{B}_{t+\Delta}-\boldsymbol{B}_{t}$ for a standard Brownian motion $(\boldsymbol{B}_{t}:t\geq 0)$.
By conditioning over $\sigma(\{\boldsymbol{B}_{s}:s\leq t\})$, we obtain

$$ $\displaystyle{\sf KL}(\boldsymbol{B}_{t+\Delta}\|\boldsymbol{y}_{t+\Delta}) \leq{\sf KL}(\boldsymbol{B}_{t}\|\boldsymbol{y}_{t})+C\mathbb{E}\{\|\hat{ \boldsymbol{m}}_{n}(\boldsymbol{B}_{t},t)\|^{2}\}\cdot\Delta\,.$ (90) $$

Therefore, by the above letting $t_{a}:=\max\{\ell\Delta:\;\ell\Delta\leq(1-\delta)t_{\mbox{\tiny\rm alg}}\}$

$$ $\displaystyle{\sf TV}(\boldsymbol{y}_{t_{a}},\boldsymbol{B}_{t_{a}})\leq\sqrt{ {\sf TV}(\boldsymbol{y}_{t_{a}},\boldsymbol{B}_{t_{a}})/2}\leq\sqrt{Cn \varepsilon_{n}}=o_{n}(1)\,.$ (91) $$

Next note that for $t_{b}:=\min\{\ell\Delta:\;\ell\Delta\geq(1+\delta)t_{\mbox{\tiny\rm alg}}\}$,

$$ $\displaystyle\hat{\boldsymbol{y}}_{t_{b}}=\hat{\boldsymbol{y}}_{t_{a}}+ \boldsymbol{B}_{t_{b}}-\boldsymbol{B}_{t_{a}}+\sum_{t\in{\mathbb{N}}\Delta\cap [t_{a},t_{b}]}\hat{\boldsymbol{m}}_{n}(\hat{\boldsymbol{y}}_{t},t)\Delta\,.$ (92) $$

Using Eq. ([91](https://arxiv.org/html/2503.08028v2#A17.E91)) we have that

$$ $\displaystyle\hat{\boldsymbol{y}}_{t_{b}}=\boldsymbol{m}_{0}+\boldsymbol{B}_{t _{b}}\,,\;\;\;{\mathbb{P}}(\|\boldsymbol{m}_{0}\|>2\delta_{n}n)=o_{n}(1)\,.$ (93) $$

Next we couple $(\hat{\boldsymbol{y}}_{t}:t\geq t_{b})$ to $(\hat{\boldsymbol{y}}_{t}^{0}:t\geq t_{b})$ defined by letting
$\hat{\boldsymbol{y}}^{0}_{t_{b}}=\boldsymbol{B}_{t_{b}}$ and,
for $t\in{\mathbb{N}}\Delta\cap[t_{b},\infty)$,

$$ $\displaystyle\hat{\boldsymbol{y}}^{0}_{t+\Delta}$ $\displaystyle=\hat{\boldsymbol{y}}^{0}_{t}+\hat{\boldsymbol{m}}_{n}(\hat{ \boldsymbol{y}}^{0}_{t},t)\Delta+\boldsymbol{B}_{t+\Delta}-\boldsymbol{B}_{t}\,.$ $$

By the assumed Lipschitz property of $\hat{\boldsymbol{m}}$ and Gromwall’s lemma:

$$ $\displaystyle\|\hat{\boldsymbol{y}}_{t}-\hat{\boldsymbol{y}}^{0}_{t}\|$ $\displaystyle\leq\prod_{t^{\prime}\in{\mathbb{N}}\Delta\cap[t_{b},t]}\Big{(}1+ \frac{C\Delta}{t^{\prime}}\Big{)}\cdot\|\boldsymbol{m}_{0}\|$ $\displaystyle\leq\Big{(}\frac{t}{t_{\mbox{\tiny\rm alg}}}\Big{)}^{C^{\prime}} \|\boldsymbol{m}_{0}\|\leq n^{1-c_{0}}\,,$ (94) $$

where the last inequality holds for some $c_{0}>0$ and all $t\leq n^{1+\delta_{0}}$, on the high probability
event $\|\boldsymbol{m}_{0}\|\leq 2\delta_{n}n$.

In order to finish the proof, we state and prove a useful lemma.

###### Lemma Q.1 .

Under the assumptions of Theorem [6](https://arxiv.org/html/2503.08028v2#Thmtheorem6), assume that $\delta_{n}$ vanishes
slowly enough. Then, for $t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}$,

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}_{n}(\boldsymbol{B}_{t},t )\|^{2}\big{\}}\leq Ce^{-(\sqrt{t}-\sqrt{t_{\mbox{\tiny\rm alg}}})^{4}/Cn}\,.$ (95) $$

###### Proof.

Let $\lambda_{1}(\boldsymbol{y}_{t})$ be the maximum eigenvalue of $(\boldsymbol{y}_{t}+\boldsymbol{y}_{T}^{{\sf T}})/\sqrt{2}$,
$\lambda_{*}(t):=\sqrt{2}(\sqrt{t}+\sqrt{t_{\mbox{\tiny\rm alg}}})^{2}$ and
$\phi(\boldsymbol{y}_{t}):=\boldsymbol{1}(\lambda_{1}(\boldsymbol{y}_{t})>
\lambda_{*}(t))$. Concentration
results about spiked GOE matrices imply, tor all $t\geq(1+\delta)t_{\mbox{\tiny\rm alg}}$,

$$ $\displaystyle{\mathbb{P}}_{n,k}\big{(}\phi(\boldsymbol{y}_{t})=0\big{)}\leq Ce ^{-n(\sqrt{t}-\sqrt{t_{\mbox{\tiny\rm alg}}})^{4}/C}\,,\;\;\;\;{\mathbb{P}} \big{(}\phi(\boldsymbol{B}_{t})=0\big{)}\leq C\,\exp\Big{\{}-\frac{1}{Cn}( \sqrt{t}-\sqrt{t_{\mbox{\tiny\rm alg}}})^{4}\Big{\}}\,.$ (96) $$

(To simplify notations, we omit the dependence of $\phi$ on $t$.)

By assumption, the MSE of $\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)$ is not larger than the one of $\hat{\boldsymbol{m}}_{n}(\boldsymbol{y}_{t},t)\phi(\boldsymbol{y}_{t})$.
Letting $\overline{\phi}(\boldsymbol{y}_{t}):=1-\phi(\boldsymbol{y}_{t})$:

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)- \boldsymbol{x}\|^{2}\big{\}}$ $\displaystyle\leq\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t )\phi(\boldsymbol{y}_{t})-\boldsymbol{x}\|^{2}\big{\}}$ $\displaystyle=\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)- \boldsymbol{x}\|^{2}\phi(\boldsymbol{y}_{t})\big{\}}+\mathbb{E}\big{\{}\| \boldsymbol{x}\|^{2}\overline{\phi}(\boldsymbol{y}_{t})\big{\}}$ $\displaystyle=\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)- \boldsymbol{x}\|^{2}\phi(\boldsymbol{y}_{t})\big{\}}+{\mathbb{P}}_{n,k}\big{(} \phi(\boldsymbol{y}_{t})=0\big{)}\,,$ (97) $$

whence

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)- \boldsymbol{x}\|^{2}\overline{\phi}(\boldsymbol{y}_{t})\big{\}}\leq{\mathbb{P} }_{n,k}\big{(}\phi(\boldsymbol{y}_{t})=0\big{)}\,.$ (98) $$

On the other hand

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t)- \boldsymbol{x}\|^{2}\overline{\phi}(\boldsymbol{y}_{t})\big{\}}$ $\displaystyle\geq\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{y}_{t},t )\|^{2}\boldsymbol{1}_{\boldsymbol{x}=0}\overline{\phi}(\boldsymbol{y}_{t}) \big{\}}$ $\displaystyle=\frac{1}{2}\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{ B}_{t},t)\|^{2}\overline{\phi}(\boldsymbol{W}_{t})\big{\}}$ $\displaystyle\geq\frac{1}{2}\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}( \boldsymbol{B}_{t},t)\|^{2}\big{\}}-\frac{1}{2}{\mathbb{P}}\big{(}\phi( \boldsymbol{W}_{t})=1\big{)}\,.$ (99) $$

Putting together Eqs. ([96](https://arxiv.org/html/2503.08028v2#A17.E96)), ([98](https://arxiv.org/html/2503.08028v2#A17.E98)), ([99](https://arxiv.org/html/2503.08028v2#A17.E99)),
we obtain (eventually adjusting the constant $C$)

$$ $\displaystyle\mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{B}_{t},t)\|^ {2}\big{\}}$ $\displaystyle\leq{\mathbb{P}}\big{(}\phi(\boldsymbol{B}_{t})=1\big{)}+2\,{ \mathbb{P}}_{n,k}\big{(}\phi(\boldsymbol{y}_{t})=0\big{)}$ $\displaystyle\leq C\exp\Big{\{}-\frac{1}{Cn}(\sqrt{t}-\sqrt{t_{\mbox{\tiny\rm alg }}})^{4}\Big{\}}\,.$ $$

∎

We are now in position to finish the proof of the theorem.
We couple the process $(\hat{\boldsymbol{y}}^{0}_{t}:\;t\geq t_{b})$ defined above with $(\boldsymbol{B}_{t}:t\geq t_{b})$
to get

$$ $\displaystyle{\sf KL}(\boldsymbol{B}_{t+\Delta}\|\hat{\boldsymbol{y}}^{0}_{t+ \Delta})$ $\displaystyle\leq{\sf KL}(\boldsymbol{B}_{t}\|\boldsymbol{y}^{0}_{t})+C\, \mathbb{E}\big{\{}\|\hat{\boldsymbol{m}}(\boldsymbol{B}_{t},t)\|^{2}\big{\}}\cdot\Delta$ $\displaystyle\leq{\sf KL}(\boldsymbol{B}_{t}\|\hat{\boldsymbol{y}}^{0}_{t})+C \Delta\exp\Big{\{}-\frac{n}{C}(\sqrt{t}-\sqrt{t_{\mbox{\tiny\rm alg}}})^{4} \Big{\}}\,.$ $$

Using ${\sf KL}(\boldsymbol{B}_{t_{b}}\|\hat{\boldsymbol{y}}^{0}_{t_{b}})=0$,
summing the last inequality over $t\geq t_{b}$, and applying Pinsker’s inequality we obtain
for $\delta_{n}\geq C^{\prime}(\log n/n)^{1/4}$ with $C^{\prime}$ a suitably large constant

$$ $\displaystyle\sup_{t\in{\mathbb{N}}\Delta\cap[t_{\mbox{\tiny\rm alg}}(1+\delta ),\infty)}{\sf TV}(\hat{\boldsymbol{y}}^{0}_{t},\boldsymbol{B}_{t})=o_{n}(1)\,.$ (100) $$

Putting together this bound and Eq. ([94](https://arxiv.org/html/2503.08028v2#A17.E94)) (which holds with high probability)
we obtain that

$$ $\displaystyle{\mathbb{P}}\big{(}\max_{t\in[t_{\mbox{\tiny\rm alg}}(1+\delta),n ^{1+\varepsilon_{0}}]}\|\hat{\boldsymbol{y}}_{t}-\boldsymbol{B}_{t}\|\geq n^{1 -c_{0}}\big{)}=o_{n}(1)\,.$ (101) $$

Using once more Lemma [Q.1](https://arxiv.org/html/2503.08028v2#A17.Thmclaim1) and the Lipschitz property of
$\hat{\boldsymbol{m}}_{n}$, we obtain that, for $\hat{\boldsymbol{x}}_{t}=\hat{\boldsymbol{m}}(\hat{\boldsymbol{y}}_{t},t)$,

$$ $\displaystyle{\mathbb{P}}\big{(}\max_{t\in[t_{\mbox{\tiny\rm alg}}(1+\delta),n ^{1+\varepsilon_{0}}]}\|\hat{\boldsymbol{m}}_{n}(\hat{\boldsymbol{y}}_{t},t)\| \geq n^{-c_{0}}\big{)}=o_{n}(1)\,,$ (102) $$

which implies the claim of the theorem.

###### Lemma Q.1 .

###### Proof.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. Would considering the 'average complexity' via diffusion over a class of easily-sampled distributions be more meaningful than focusing on 'bad' extreme counterexamples like the sparse submatrix prior?
2. Do your have conclusions about the training procedure in this paper, or does it imply that for the provided counterexample, the training of diffusion models will always fail due to its inherent computational hardness?
3. Is it possible this specific computational bottleneck is just an artifact of this particular formulation of diffusion sampling, or does it represent a fundamental hardness that alternative generative models like flow-matching models could not bypass?
soundness: 4
strengths: 1. This paper's core strength is its novel perspective, shifting the analysis of diffusion models from standard statistical limits to the overlooked statistical-computational gap.
2. The authors construct an excellent counterexample using exactly the sparse submatrix estimation to provide a distribution that is easy to sample from directly but computationally intractable for a diffusion model to learn.
3. Supportive simulations clearly validate the theory, showing that a computationally-bounded diffusion model fails to generate accurate samples from this "hard" distribution, just as the authors predict.
4. The paper's argument is mathematically rigorous, employing formal, well-established proofs.
summary: The paper's primary contribution is to rigorously prove that a computational bottleneck, not just a statistical one, limits the power of denoising diffusion models. The authors prove that not all easily-sampled probability distributions can be efficiently learned by these models. They construct a specific counterexample based on the sparse submatrix estimation problem, which is conjectured to be hard. This shows that even a nearly-perfect, efficiently-computed drift function will fail to generate accurate samples, thus identifying a fundamental gap between what is statistically possible and what is computationally feasible for this class of generative models.
weaknesses: 1. The paper's organization is not straightforward; the structure is relatively loose, which makes the argument hard to follow.
2. The planted clique counterexample feels specific, which makes it difficult to connect this theoretical computational bottleneck to the broad success of diffusion models on real-world data.

---

confidence: 2
contribution: 4
initial_rating: 8
presentation: 3
questions: * "l. However if such drifts exist, our results suggest that minimizing the score matching objective is not the right approach to find them (since the difference in value with bad drifts will be superpolynomially small)." (l.186) This is an interesting remark. There is some evidence that diffusion models fail at approximating the score. First, they are trained on the empirical measure hence the target distribution is the "optimal" target gives a fully memorizing denoiser.  Second, due to neural network approximations they often fail to correctly approximate the true denoiser. I would be keen to understand how the authors deal with that intrisic limitation of diffusion models and how this would affect their analysis. 

* In Proposition 2.1., I do not understand why we can have equality to 0 in (ii).

* Could you provide more details regarding the construction of the drift in Proposition 2.1? (the construction of a drift which is a bad approximation to the true denoiser but gives good diffusion sampling results). 

* Before Theorem 1 I think it would be worth providing just a few sentences regarding the validity of Assumption 2 and Assumption 1. At this stage it is very much unclear for the reader if those assumptions are easy or not to satisfy.
soundness: 4
strengths: * The paper is theoretically strong. I have to command the quality of the writing (apart from the introduction of diffusion models and some notation choice discussed in the "Weaknesses" section). Even though I am not a domain expert on stochastic localization and random matrix theory, the paper was pleasant to read and the ideas clearly exposed. 

* The results presented in the papers will be of great interest for the theoretical diffusion model community. Those results indeed fill a gap in the literature regarding what can be achieved with diffusion models.
summary: In this paper, the authors study the feasibility of sampling with diffusion models. More precisely, they are interested in answering the following question: "Are there distributions such that, there exist almost optimal denoisers that can be computed in polynomial time and for which diffusion sampling fail?" They answer this question by the affirmative. The paper proceeds as follows:

* In  Theorem 1, they show that for any arbitrary polynomial time denoiser which satisfies two generic assumptions, there exists a modified denoiser for which the diffusion sampling fails. In Corollary 3.1, they apply Theorem 1 to the case of sparse random matrix measures. Other sparse examples are considered 

*  While Theorem 1, shows that "good denoiser does not imply good sampling", in Theorem 2 they show that "good sampling implies good denoiser".

* In Theorem 3, they show that under an additional Lipschitz condition, the results can be strengthen to show that *any* polynomial time algorithm which is nearly optimal for the denoising task will lead to poor sampling quality.

The authors conclude the papers with toy numerical experiments.
weaknesses: * One of my main pain point with the paper is the choice of following the notation and convention of [1]. While I understand that it might be easier for the authors who seem to familiar with the techniques developed in [1] to build their theory on top of this framework it severely limits the adoption of the work, as it has to be translated back into "classical" diffusion model framework (see [2, 3, 4] for a (non-exhaustive) list of papers which adopt the classical diffusion model convention). The burden of the explanation and communication should be on the authors and not on the readers.

* In my opinion, the strongest result is Theorem 3 but there is little discussion on how realistic is the Lipschitz assumption the authors consider. I found Corollary 5.1. to be very hard to read. To come back to the Lipschitz justification I am not sure I understand the following sentence "We assume the Lipschitz constant to be C/t, because the input of the denoiser is yt = tx + Wt, and hence the two t-dependent factors cancel.)".

* Something that is unclear to me is the generality of the results presented in the paper. While Theorem 1 is very general it is not quite clear to me how to apply them beyond the sparse random matrix measure studied by the authors. In particular, I am wondering how specific those results are. To further specify my question: is (one of) the contribution(s) of the paper 1) the identification of a set of measures under which the assumptions stated in the general Theorems are valid thereby showing some limitations of diffusion models 2) the statement that most data distributions (even the ones considered in practice) fall under the category of the general Theorems thereby showing the much more stronger statement that diffusion models are flawed. It would be good to clarify  the scope of the paper in that respect. 

[1] Montanari (2023) -- Sampling, diffusions, and stochastic localization

[2] Conforti et al. (2024) -- KL Convergence Guarantees for Score diffusion models under minimal data assumptions

[3] Li et al. (2024) -- O(d/T) Convergence Theory for Diffusion Probabilistic Models under Minimal Assumptions

[4] Chen et al. (2022) -- Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: Alternative Objectives: Proposition 2.1 suggests that score matching is the root cause of the failure when ICGs are present. Does this motivate the development of alternative training objectives for diffusion models that do not strictly enforce denoising optimality, particularly at low SNR (small $t$)?

Optimization: Theorem 1 implies that the score matching loss landscape possessing near-optimal solutions (in the poly-time class) that are bad for sampling.  Do we know anything about whether standard optimization methods are likely to find these bad near-optimizers rather than potential good ones?
soundness: 4
strengths: Originality. The paper goes beyond the previous analysis of diffusion samplers devoted to KL/Wasserstein bounds or generalization arguments. It provides an interesting computational perspective. In particular, it identifies the information-computation gap as the key problem for correct diffusion sampling; it shows that even near-optimal score-matching drifts (among poly-time methods) can sample the wrong distribution and it establishes reductions from diffusion sampling to denoising.

Quality: Theorems 1--3 are stated with clear hypotheses, and Appendix proofs are thorough. The proofs are relatively easy to follow given the technicality of arguments.

Clarity: The paper is well-written. Assumptions specifying the ICG gap are stated clearly. The paper separates neatly the main conceptual arguments in the main text from the proofs. The figures nicely illustrate the theoretical claims.

Significance: This is a very interesting paper explaining when diffusion samplers fail for computational reasons—even when sampling from $\mu$ is itself easy. They also specify when positive results (which rely on accurate scores or favorable distribution classes) do not apply.
summary: This paper addresses a very interesting question: Can denoising diffusions efficiently sample from any distribution $\mu$ for which sampling is otherwise tractable?

The paper presents a rigorous negative answer by showing that diffusion samplers inherit the computational hardness of the underlying denoising task. The core mechanism of diffusion relies on learning the drift term $m(y,t)$, i.e. the Bayes-optimal denoiser. If this denoising problem exhibits an Information-Computation Gap (ICG), a regime where the statistically optimal error is low, but no polynomial-time algorithm can achieve it, then the authors show that diffusion sampling based on score-matching must fail.
The paper includes several key theoretical results.
Theorems 2 \& 5): Efficient diffusion sampling implies efficient near-Bayes-optimal estimation. The contrapositive confirms the central thesis: computational hardness in denoising implies diffusion sampling fails.
Theorem 1: The paper proves the existence of polynomial-time drifts that are super-polynomially close to the optimum (among poly-time algorithms) in the score-matching objective, yet produce samples maximally far from the target $\mu$ (in Wasserstein distance). So the standard training objective is insufficient and potentially misleading when an ICG exists.
Theorem 3: The paper shows that any poly-time, near-optimal drift satisfying a mild any \emph{Lipschitz} near-optimizer fails to sample correctly. Furthermore, \emph{Corollary~5.1} provides a concrete sufficient condition: a $C/t$-Lipschitz drift fails when $t \ge (1+\delta)\,t_{\mathrm{alg}}$.

These results are illustrated using the sparse PCA model where sampling is trivial but denoising is conjectured to be hard below a threshold $t_{\mathrm{alg}}$. Empirical results show that the diffusion process produces degenerate samples.

This work is very significant. It goes beyond analyzing statistical generalization or discretization errors, identifying fundamental limits of denoising diffusions.
weaknesses: Dependence on assumptions: The main results depend on Assumptions 1 and 2, which do not hold unless an ICG is present. This aligns with predictions for models like Sparse PCA, the underlying hardness remains conjectural in the general case.
It would help to map Assumptions 1 and 2 to specific hardness conjectures in the literature.

Robustness of Theorem 1: Th 1 shows that sampling can fail when the drift is $O(n^{-D})$ close to the optimum among poly-time algorithms. This is a fairly stringent definition of near-optimality. How robust are the results to this assumption? even an heuristic argument would be good. 

Construction of the Drift (Theorem 1) The drift $\hat{m}$ constructed to prove Theorem 1 is somewhat artificial. While this demonstrates the insufficiency of the score matching objective, it does not demonstrate that such deceptive drifts arise naturally during standard training procedures (e.g., gradient descent on neural networks).
        \item \emph{Actionable Insight:} Discuss whether standard optimization dynamics are likely to find these specific deceptive local minima, or if this is primarily an existential proof. Are there architectural choices or regularization techniques that might steer optimization away from these solutions?
    \end{itemize}

Lipschitz assumption. Th 3 establishes a negative result for Lipschitz drifts, specifically requiring a $C/t$-Lipschitz condition for $t \ge (1+\delta)t_{\mathrm{alg}}$. However, non-Lipschitz drifts might succeed.
 It would be good to clarify whether this Lipschitz condition is key to the failure of just a technical requirement. Are the typical architectures used in practice (e.g. your GNN in section 3) satisfy this specific $C/t$-Lipschitz property?
    
Applicability of the results. The analysis focuses on distributions where the ICG is well understood, e.g. Sparse PCA. How are these findings translate to standard distributions meet in physical systems; e.g. Lennard-Jones or \phi^4? 
It would be beneficial if the authors could discuss even informally the characteristics of  distributions that might lead to a ICG. 

Minor point: The authors follow a non-standard presentation of diffusion models (the one originating from stochastic localization). I know this is equivalent but i think it hurts the readability of the paper.


Decision: Accept (Poster)