Quantum enhanced rare event discovery and sampling Naixu Guo,1, ∗ Po-Wei Huang,2 Qisheng Wang,3, 4 Jayne Thompson,1, 5 Patrick Rebentrost,1, 6 Mile Gu,1, 7, † and Chengran Yang7, ‡ 1
arXiv:2606.06316v1 [quant-ph] 4 Jun 2026
Centre for Quantum Technologies, National University of Singapore, Singapore 117543, Singapore 2 Mathematical Institute, University of Oxford, Oxford OX2 6GG, United Kingdom 3 School of Computer Science, Shanghai Jiao Tong University, Shanghai 200240, China 4 School of Informatics, University of Edinburgh, Edinburgh EH8 9AB, United Kingdom 5 College of Computing and Data Science, Nanyang Technological University, Singapore 639798, Singapore 6 School of Computing, National University of Singapore, Singapore 117417, Singapore 7 Nanyang Quantum Hub, School of Physical and Mathematical Sciences, Nanyang Technological University, Singapore 637371, Singapore (Dated: June 5, 2026) Financial crashes, cascading failures in infrastructure, and critical errors in AI systems are frequently triggered by events that occur with extremely small probability. Efficiently discovering and sampling events with probability below a threshold is therefore of critical interest. Yet this task is highly non-trivial using existing classical or quantum methods. Being rare, such events require an immense sampling overhead to collect sufficient data samples. Moreover, because the rare events are not known in advance, they cannot be flagged for amplification using standard techniques. Here, we introduce a quantum algorithm for rare-event discovery and sampling without first learning which events are rare. The algorithm achieves the optimal quantum scaling with the rarity threshold. We further demonstrate that this can achieve a quadratic speedup for heavy-tailed systems whose tail has nonvanishing total mass, and translates into a robust polynomial speedup for stationary stochastic processes, with the exponent determined by its entropy-rate structure.
I.
INTRODUCTION
From edge-case scenarios that compromise AI safety [1, 2] to black-swan and seismic events [3, 4], rare events can have a disproportionate impact. The discovery of such events — together with the simulation of their potential consequences — is therefore critical for accurate risk assessment and for devising contingencies to avert catastrophe [5, 6]. Consequently, there is strong interest in efficient methods to sample only from the rare events of a distribution, namely those whose probabilities fall below some target threshold ∆. Yet the scarcity of such events renders them statistically difficult to model. Critically, we do not a priori possess a list of what outcomes are rare, nor knowledge of what effects they may have. This limits the applicability of existing quantum rare-event techniques, such as quantum Monte Carlo and amplitude amplification, which assume that the target events can be efficiently identified and flagged. Classical methods face similar challenges: importance sampling approaches [7, 8] rely on carefully designed biasing distributions whose effectiveness depends strongly on prior knowledge of the rare events themselves. Consequently, it is not immediately clear whether methods more efficient than classical Monte Carlo sampling exist. The standard approach is therefore first to estimate the probabilities of outcomes to determine which events are
∗ [email protected] † [email protected] ‡ [email protected]
rare, and only then to sample the process and post-select those events. Here, we introduce a quantum algorithm for sampling rare events without first knowing which events are rare. Given a probability distribution P through a quantum state preparation unitary UP , our algorithm produces samples from the distribution restricted to events whose probabilities lie below √ the threshold ∆. Its query complexity scales as O(1/ prare ∆), where prare is the probability that an ordinary sample from P is rare. We prove that the algorithm achieves optimal scaling with respect to ∆, offering a quadratic improvement over the optimal classical 1/∆ dependence. We further demonstrate speedups in heavy-tailed distributions and stochastic processes, where many individually unlikely events can together form a substantial tail. Beyond sampling, our algorithm can also synthesize a coherent quantum representation of the rare-event distribution that presents a crucial quantum resource for further downstream quantum advantage
II.
RESULTS
Rare event sampling. Consider a probability distribution P over N events accessible via sampling. For a given threshold ∆ > 0, we define the unknown but nonempty set of rare events as R = {xi | 0 <P P (xi ) ≤ ∆}, which has a total probability mass of prare = xi ∈R P (xi ). Our objective is to construct a procedure capable of generating samples from the conditional distribution ( P (xi ) if xi ∈ R PR (xi ) := prare . (1) 0 if xi ∈ /R
2 Quantum Sampler U
⟩
⟩
|0⟩ |0⟩ |0⟩ |0⟩
U
Spectral Embedding ⟩
Rare Event Amplifier Q-Rare-Event Sampler
0 0 0 0 [0 0 ]
0 0 0 0 [0 0 ]
}
Probability Thresholding
⟩
0 0 0 0 [0 0 ]
FIG. 1. Overview of quantum rare event sampling. Top: Starting from the initial sample state |P ⟩, the rare-event amplifier filters basis states with P (xi ) ≤ ∆ and amplifies the surviving subspace to prepare the target state |PR ⟩. Bottom: The amplifier is implemented by first √ converting amplitudes into a spectral representation (spectral embedding), and then applying an approximate threshold at ∆. This coherent pipeline enables sampling from the rare-event distribution without requiring prior identification of what events are rare.
A sharp cutoff at ∆ is unattainable with finite samples or bounded circuit depth. Because near-threshold events become indistinguishable from the true rare set, we relax the boundary by defining an ambiguous region of near-threshold events that may be mistakenly included: S∆ = {xi : ∆ < P (xi ) ≤ (1 + α)∆} for a small constant α > 0 (e.g., α = 0.001). This relaxation introduces an approximation error, measured P by the total ′varia′ tion distance Dtv (PR , PR ) = 12 xi |PR (xi ) − PR (xi )| between the ideal and achieved distributions. This error is controlled by the relative weight ζ = p∆ /prare , where P p∆ = xi ∈S∆ P (xi ). Our question is: how many samples from P (or applications of UP ) are required to construct a sampler O(ζ + ϵ)-close to PR for a tunable precision ϵ? To achieve this task quantum mechanically, we assume access to a quantum sampler (unitary) UP that prepares the quantum sample state |P ⟩ = UP |0⟩ =
N p X i=1
P (xi )|xi ⟩,
(2)
where |xi ⟩ denotes the computational basis state corresponding to the binary representation of event xi , analogous to other quantum sampling algorithms [9–12]. Our objective then translates to constructing a quantum circuit UPR that prepares the quantum rare event state |PR ⟩ = UPR |0⟩ = √
X p 1 P (xi )|xi ⟩. prare
(3)
xi ∈R
Measuring the state |P ⟩ or |PR ⟩ in the computational basis yields samples distributed according to P and PR respectively, matching the classical definitions. This correspondence allows us to compare classical and quantum algorithms based on the number of samples from P or applications of UP required to construct the rare event
Problem: Rare Event Sampling Input: Rare set: Goal: Cost:
Sampling access to P classically, or quantum access through UP . R = {x : 0 < P (x) ≤ ∆}. Construct a sampler for PR with total variation error O(ζ + ϵ). Number of samples from P , or calls to UP and UP† .
FIG. 2. Problem setting of rare event sampling. The task is to sample from the conditional distribution over events whose probabilities are at most ∆. The error term ζ captures the unavoidable ambiguity from near-threshold events.
sampler. Beyond sampling, the quantum circuit UPR can efficiently prepare many copies of |PR ⟩ for use in downstream quantum algorithms where such superpositions on critical. Our goal is to develop both classical and quantum algorithms for this task, characterize the optimal sample complexity, and identify the potential for quantum advantage. Further details are provided in Sec. D, and we summarize the problem setting as Fig. 2. Classical baseline. We first establish what can be achieved classically. If the set of rare events R were known a priori, one could directly apply rejection sampling: generate a sample xi from P and accept it if and only if xi ∈ R. Each accepted sample would then require O(1/prare ) samples from P on average. However, the rare set R is not known a priori. To decide whether to accept a sample xi , we must determine if P (xi ) ≤ ∆, which requires estimating outcome probabilities to precision O(α∆) ⊆ O(∆). By adapting
3 distribution-learning and estimation techniques from previous works [13, 14], we can construct a list of candidate rare events and use it for rejection sampling. Combining these steps yields the following result (see Sec. D for the detailed proof and construction). Result 1 (Classical rare event sampling). Given sampling access to a probability distribution P over N events, for any ϵ > 0, 1 1 N 1 1 O min log , 2 log + (4) ∆ ϵ ∆ ϵ prare samples from P suffice to construct a sampler for a distribution O(ζ + ϵ)-close to PR in total variation distance, where ζ = p∆ /prare . Here, the first term arises from the discovery of this list of rare events, while the second reflects the standard overhead of rejection sampling. As we will see, in many contexts, such as rare-event discovery in heavy-tailed distributions and stochastic processes, the former term dominates. This is because there can exist exponentially many rare events, ensuring that prare ≫ ∆. In such situations, our sampling cost minimally scales as Ω(1/∆). Our next result shows that there is way around this using any classical technique: Result 2 (Classical lower bound). For ∆ > 0, any classical algorithm requires Ω(1/∆) samples from P to construct a sampler for PR . The above lower bound is robust, and continues to hold when we construct a sampler for a distribution that is close to PR . The proof proceeds by reduction from the distribution distinguishing problem [15]. Details can be found in Sec. D. Quantum algorithm. We now turn to the quantum setting. A naive quantum approach would translate the classical strategy directly: use standard amplitude estimation [16] to determine P (xi ) and verify whether each outcome satisfies the rare-event criterion. This incurs a linear overhead in the sample space size N , negating the potential for a broad quantum advantage. Even multidimensional quantum amplitude estimation [14], while more efficient at identifying multiple outcomes, still requires O(1/∆) applications of UP to estimate probabilities with precision O(∆), which we prove to be suboptimal, and provides no scaling advantage. Our algorithm takes a different approach: it directly amplifies the rare components of the superposition state |P ⟩ without explicitly needing to first identify which events are rare. The conceptual strategy is illustrated in Fig. 1, and we formalize our contribution as follows (see Sec. D for details). Result 3 (Quantum rare event sampling). Given a quantum sampler UP for a probability distribution P over N events, for any ϵ > 0, 1 1 √ O log (5) ϵ prare ∆
applications of UP suffice to prepare and measure a quantum state, yielding an event from a distribution O(ζ + ϵ)-close to PR in total variation distance, where ζ = p∆ /prare . Our algorithm produces a quantum rare-event state, |PR ⟩, which is a proportional superposition of all possible rare events. Measuring |PR ⟩ yields a rare-event sample. The algorithm comprises three key components (detailed in Sec. D): 1. Rare-event amplification. We can prepare the target rare-event state |PR ⟩ as long as we can synthesis the unitary operator UΠR such that p UΠR |0⟩|P ⟩ = |0⟩ΠR |P ⟩ + 1 − prare |1⟩|⊥⟩, P where ΠR = xi ∈R |xi ⟩⟨xi | and |⊥⟩ is some arbitrary state. Measuring the ancilla and post-selecting √ on outcome 0 yields ΠR |P ⟩/ prare , which occurs with probability prare . To boost this success probability, we apply amplitude amplification [16], re√ quiring O(1/ prare ) applications of UP and UΠR . Thus, we reduce the problem to that of synthesizing UΠR . 2. Probability thresholding. To implement UΠR , we apply an approximate Heaviside step function to the of the Hermitian operator H = p P eigenvalues P (x )|x ⟩⟨x i i i |. The non-smooth Heaviside xi function is approximated by a polynomial of degree √ O(1/ ∆ log(1/ϵ)), where ϵ is the approximation error. This approximation has a finite√transition √ width of O(α ∆) around the threshold ∆, which corresponds directly to the ambiguous √ set S∆ defined earlier. This step requires O(1/ ∆ log(1/ϵ)) applications of a unitary UH that block encodes H. Thus the synthesis UH enables rare-event amplification. 3. Spectral embedding. We can synthesis UH using only access to UP . The key observation is that the amplitudes p of |P ⟩ equal the eigenvalues of H, both given by P (xi ). This correspondence allows us to extract the spectral information directly from the quantum sampler UP . By utilizing nonlinear amplitude transformation techniques [17–21], we convert the amplitude information in |P ⟩ into the eigenvalue structure of UH with only O(1) applications of UP . Combining these three advances yields a systematic way to prepare |PR ⟩ using Up that aligns with the state complexity. We complement this constructive upper bound with a matching lower bound, confirming that our algorithm achieves optimal scaling in ∆. The proof is established via reduction from the state and distribution discrimination problem [22, 23], showing the minimum resources required by any quantum algorithm for this task. Details can be found in Sec. D. Similar to the classical case, the lower bound also holds for sampling from an approximate distribution.
4
rare event original amplified
Quantum Classical
0.8
0.6
0.6 Error
Probability
0.8
0.4
0.4
0.2 0.2
0.066 0.044 0.022 0.0
0.0
11 11 10 11 01 11 00 11 11 10 10 10 01 10 00 10 11 01 10 01 01 01 00 01 11 00 10 00 01 00 00 00
0
500
1000
1500 2000 2500 3000 3500 Applications of (Quantum) Sampler
4000
4500
5000
Event
FIG. 3. Simulation results for a Dyson-Ising chain. Left: Original (red) versus amplified (blue) probabilities for a sequence of length L = 4. The dashed horizontal line represents the probability threshold ∆, and the shaded regions denote the target rare events. Right: Total variation distance (TVD) error versus the number of applications of the quantum sampler for a sequence of length L = 8 and Markov order χ = 3. The quantum algorithm converges more smoothly and rapidly than the classical baseline, which exhibits oscillatory behavior due to the intermittent sampling of rare events.
Result 4 (Optimality of quantum rare event sampling). √ For ∆ > 0, any quantum algorithm requires Ω(1/ ∆) applications of UP to construct a sampler for PR .
infrastructure failures [25], and seismic event catalogs [3]: each extreme realization may be exceptionally unlikely, yet the collection of such realizations can remain large enough to matter.
Importantly, the quantum algorithm produces more than a source of rare-event samples. Before the final measurement, it prepares the rare-event state |PR ⟩, a coherent representation of the full rare tail. This coherent form allows one to compare tails directly: for example, given rare-event states for two distributions on the same event space, a SWAP test estimates their rare-tail overlap and therefore whether the two models assign weight to similar extreme scenarios without needing to to compile a list of rare events. More broadly, rare-event state preparation upgrades rare-event sampling into a quantum primitive that can be used by downstream algorithms for analyzing distributional tails.
Application to Stochastic Processes. Consider a time-series governed by a stationary, ergodic stochastic processes described by a sequence of random variables X = . . . X0 , X1 , X2 . . .. Let it have some finite entropy rate H(X), defined as its average uncertainty per symbol 1 . The asymptotic equipartition property (AEP) tells us that almost all length-L trajectories lie in a typical set, each with probability approximately 2−LH(X) as we scale L [26]. Our goal is then to sample from only rare trajectories x0:L = x0 x1 · · · xL−1 - trajectories where
Application to Heavy Tail Distributions. Comparing Result 1 and Result 3, the scaling depends on the total probability mass is carried by the rare tail as a whole. When prare = Ω(1), √ the quantum algorithm saturates at the optimal O(1/ ∆) scaling, whereas the classical cost remains O(1/∆). The resulting quantum advantage is therefore quadratic. Heavy-tailed distributions provide a natural setting in which rare events can be individually unlikely but collectively important. In Methods C and Sec. F, we make this intuition concrete using rank-ordered power laws, where xi denotes the i-th most likely outcome and P (xi ) ∝ i−γ . The analysis gives a transition at γ = 1: for broad tails with 0 < γ ≤ 1, the rare tail can carry non-negligible total probability mass. In these regimes, the quantum algorithm reaches its ideal scaling, providing a quadratic improvement over the classical baseline. For steeper power laws with γ > 1, the rare tail loses total mass as ∆ decreases, so the additional amplification cost weakens the speedup below the ideal quadratic scaling. These models capture the central feature of fat-tailed risk in systems such as financial markets [24], cascading
P (x0:L ) ≤ ∆ = 2−αL
(6)
for some α ≥ H(X), such that their relatively likelihood compared to typical events dies off exponentially with L. In Sec. G 2, we show that this also setting in which we can determine exact - and often robust - quantum scaling advantage. Specifically, existing results in classical rare event generation tell us that such the number of trajectories such that P (x0:L = ∆) generally scales as 2s∆ L for some s∆ > 0 [27]. We can then prove that the total probability mass of the rare events prare , scales as ∼ ∆1−µ , where µ = s∆ /α. This algebraic relationship reduces the general quantum complexity of our protocol to Õ(∆(µ−2)/2 ). Near full quadratic quantum advantage can be achieved when µ is near to 1. This bound can be saturated. A particularly simple case is that of throwing L ≫ 1 coins with some bias p ≪ 1. In this setting, L length trajectories
1 The
entropy rate is 1 limL→∞ L H(X1 , X2 , . . . , XL ), entropy.
given where
H
by is
H(X) = the Shannon
5 in the typical set have probabilities scaling as 2−hp L where hp ≪ 1 is the binary entropy 2 . We then set rare-events threshold ∆ = 2−L . In this setting, s∆ = α = 1, and thus µ = 1. Numerical results. As an illustration, we apply our algorithm to sample rare events from the thermal distribution of a Dyson–Ising spin chain with finite interaction range χ. Its Hamiltonian is Hχ = −
χ XX i
Jk si si+k ,
(7)
k=1
where si ∈ {+1, −1} denotes the spin at site i, Jk is the coupling at separation k, and Jk = 0 for k > χ. When considering the spins sequentially along the chain, the thermal distribution at each temperature T induces a stochastic process over spin sequences. For a finite interaction range χ, the conditional distribution of the next spin depends only on the previous χ spins, so the resulting process has Markov order χ. We set α = 2H(X), χ = 3, Jk = 2−k , T = 0.8. We then apply our algorithm for a target ∆ ≈ 4.32 × 10−3 . In this setting, we find µ ≈ 0.866, signaling almost full quadratic quantum advantage. The results are shown in Fig. 3. The left panel shows that after 600 applications of the quantum sampler, the algorithm selectively amplifies the rare events. The right panel plots the total variation distance (TVD) as a function of the application count of the quantum sampler. The TVD of the quantum algorithm decreases smoothly and more rapidly than the classical baseline. The classical curve also exhibits oscillatory behavior, as rare events are observed only intermittently. Small changes in the total sample count could distort the probability outcome if the number of rare event samples does not change.
III.
DISCUSSION
We have developed a quantum algorithm for discovering and sampling rare events of a given distribution P , defined as those occurring with probability below a threshold ∆. The algorithm prepares a quantum state superposing all such rare events weighted by their √ likelihood, requiring a general query complexity of Õ(1/ prare ∆). When the aggregate rare-event mass prare is non-vanishing, this com√ plexity saturates the optimal Θ(1/ ∆) scaling limit. For heavy-tailed distributions, this yields a quadratic speedup over all classical methods. Meanwhile, for stochastic processes, we can derive a general closed-form solution for prare and thus can analytically determine the anticipated quantum speed-up - which should be near quadratic in physically motivated settings
2 That is h
p = −p log p − (1 − p) log(1 − p)
A key feature of our algorithm is that it produces a quantum rare-event state |PR ⟩, opening pathways beyond simply sampling rare events. Combined with the quantum SWAP test, for instance, this enables efficient measurement of whether the rare-event profiles of two distributions closely match, a quantity of significant interest in extreme value theory [28], financial tail risk analysis [24], and AI safety [2]. More broadly, the ability to coherently manipulate rare-event states may serve as a primitive for quantum machine-learning methods that focus on distributional tails. Realizing these theoretical benefits in practice requires coherent access to a quantum sampler and fault-tolerant quantum primitives. Despite these demanding hardware requirements, our algorithm interfaces naturally with existing quantum simulation subroutines. For stochastic processes in particular, recurrent quantum circuits provide systematic methods to construct the required quantum samplers [9, 12, 29–31], and these circuits can require drastically less memory than their classical counterparts [32]. Indeed, such memory advantages have already been demonstrated in the context of rare-event sampling for high-temperature spin chains [27]. By leveraging these constructions as subroutines, our framework enables simultaneous memory and computational advantages for rare-event analysis of complex systems. Ultimately, our work demonstrates that a rigorous, optimal quantum advantage in rare-event analysis persists even when the target events cannot be identified in advance.
IV. A.
METHODS
Classical rare event sampling
We briefly describe the steps in the classical algorithm. We need to estimate and construct a set of rare events R with the given sampling access. There are two ways to achieve this: one is to utilize the Monte Carlo estimation with precision O(∆), and the other is to learn the distribution with ℓ∞ -distance O(∆). The complexity of the Monte Carlo method is O(∆−1 log N ) with a constant success probability. As shown in Ref. [13], the complexity of ℓ∞ -tomography is Θ(∆−2 ). Therefore, this step requires O(min{∆−1 log N, ∆−2 }) samples. Based on this constructed set R, we further implement rejection sampling to achieve the sampler. For each sample xi obtained from the distribution, we check whether xi ∈ R or not. If yes, it is our output of the rare event. Otherwise, we keep resampling. This step takes O(1/prare ) samples to give a single rare event.
B.
Quantum rare event sampling
Here, we provide more details about the probability thresholding and spectral embedding components of our quantum algorithm.
6 Probability thresholding. The construction requires classifying each basis state |xi ⟩ as rare or non-rare based on its probability P (xi ). This is equivalent to applying a filtering function: X p ΠR = f P (xi ) |xi ⟩⟨xi |, xi
where f√ is a Heaviside step function satisfying f (p) = 1 if p ≤ ∆ and f (p) = 0 otherwise. To implement this p filtering, we observe that the amplitudes P (xi ) of |P ⟩ are precisely the eigenvalues of the Hermitian operator Xp H= P (xi )|xi ⟩⟨xi |. xi
The problem thus reduces to applying f to the eigenvalues of H, a task known as eigenvalue transformation. To implement this eigenvalue transformation on a quantum computer, we use quantum signal processing [33, 34], which applies polynomial functions to eigenvalues. We therefore approximate the Heaviside function with a polynomial. This √ introduces a finite transition width around the threshold ∆, which is a fundamental limitation since polynomials cannot implement perfectly sharp transitions. The transition width corresponds to√the ambiguous set S∆ . A polynomial of degree d = O(1/ ∆ log(1/ϵ)) achieves approximation error ϵ. Using quantum signal processing techniques, UΠR can then be constructed with O(d) applications of a unitary UH that embeds H. Spectral embedding. It remains to construct UH from the quantum sampler UP . Noticing that the spectral information of H is already encodedpin UP : the first column of UP contains the amplitudes P (xi ), which are exactly the eigenvalues of H: p ⟨xi |UP |0⟩ = P (xi ) = λi (H). The challenge is to convert this amplitude information into an eigenvalue encoding coherently. Using nonlinear amplitude transformation techniques [17, 18], which leverage quantum singular value transformation and coherent amplitude manipulation, we construct UH with only O(1) applications of UP and UP† .
total probability mass is prare =
1 ZN,γ
N X
k −γ ,
k=k∆
where k∆ := min{k : P (xk ) ≤ ∆}. To expose the dependence on the support size, we set N (∆) = Θ(∆−q ), and restrict to the nontrivial regime P (xN ) ≤ ∆ < P (x1 ) so that rare events exist, but not all events are rare. As shown in Sec. F, the rare-tail mass changes behavior at γ = 1. For broad tails with 0 < γ < 1, the nontrivial support-growth regime is 1 ≤ q < 1/(1 − γ). In this regime, prare = Ω(1). The quantum algorithm reaches SQ (∆) = O(∆−1/2 ). For polynomial-size support, the classical identification cost scales as SC (∆) = O ∆−1 log(1/∆) . Thus broad power-law tails yield a quadratic leadingpower improvement. At the critical exponent γ = 1, the same ideal quantum scaling holds when the support grows faster than 1/∆, namely when q > 1. In that case, 1 prare → 1 − , q and again SQ (∆) = O(∆−1/2 ). At the boundary q = 1, the rare mass decays only logarithmically, so the scaling is only logarithmically worse than the ideal case. For steeper power laws with γ > 1, rare events exist when q ≥ 1/γ. In this regime, prare = Θ ∆(γ−1)/γ , and hence SQ (∆) = O ∆−1+1/(2γ) . The algorithm still gives a polynomial improvement over the classical baseline, but the speedup is weaker than the ideal quadratic regime because the rare tail itself has vanishing total probability mass. ACKNOWLEDGMENTS
C.
Power-law tails
We illustrate the role of the aggregate rare-event mass using rank-frequency power-law distributions. Consider P (xk ) =
k −γ , ZN,γ
where xk is the k-th most likely event, γ > 0, and ZN,γ = PN −γ . Since P (xk ) decreases with k, the rare-event j=1 j set R∆ = {xk : P (xk ) ≤ ∆} is a tail in rank space. Its
NG, JT, PR, and MG are supported by the National Research Foundation, Singapore, through the National Quantum Office, hosted in A*STAR, under its Centre for Quantum Technologies Funding Initiative (S24Q2d0009). NG also acknowledges support through the Research Excellence Scholarship from SandboxAQ. PWH acknowledges support from the Engineering and Physical Sciences Research Council (EPSRC) Doctoral Training Partnership (EP/W524311/1), with a CASE Conversion Studentship in collaboration with Quantum Motion. PWH
7 further acknowledges support from the Ministry of Education, Taiwan, for a Government Scholarship to Study Abroad (GSSA) and St. Catherine’s College, University of Oxford, for an Alan Tayler Scholarship. QW acknowledges the support of startup funding from Shanghai Jiao Tong University and the Engineering and Physical Sciences Research Council under grant EP/X026167/1. MG acknowledges the support of the National Research Foundation of Singapore through the NRF Investigatorship Program (Award No. NRFNRFI09-0010). CY is funded by Schmidt Sciences, LLC. CQT acknowledges funding from OCBC via a joint NUS-OCBC research project.
CODE AND DATA AVAILABILITY
The source code and data for numerics can be found at https://github.com/georgepwhuang/rare_event. AUTHOR CONTRIBUTIONS
NG and CY conceived the original idea. PWH conducted the numerical experiments. QW provided the proof idea for the lower bounds. JT and MG provided valuable feedback on the manuscript. PR provided input on the theory. All authors contributed to the theoretical analysis and the writing of the paper. COMPETING INTERESTS
All authors are inventors on a patent application related to this work.
[1] D. Hendrycks and K. Gimpel, A baseline for detecting misclassified and out-of-distribution examples in neural networks, in International Conference on Learning Representations (2017). [2] D. Amodei, C. Olah, J. Steinhardt, P. Christiano, J. Schulman, and D. Mané, Concrete problems in AI safety (2016), arXiv:1606.06565 [cs.AI]. [3] T. Utsu, Representation and analysis of the earthquake size distribution: A historical review and some new approaches, Pure Appl. Geophys. 155, 509 (1999). [4] C. Shyalika, R. Wickramarachchi, and A. P. Sheth, A comprehensive survey on rare event prediction, ACM Comput. Surv. 57 (2024). [5] P. Embrechts, C. Klüppelberg, and T. Mikosch, Modelling Extremal Events (Springer Berlin Heidelberg, 1997). [6] D. Sornette and G. Ouillon, Dragon-kings: Mechanisms, statistical methods and empirical evidence, Eur. Phys. J. Spec. Top. 205, 1–26 (2012). [7] S. T. Tokdar and R. E. Kass, Importance sampling: a review, WIREs Comput. Stat. 2, 54–60 (2009). [8] M. Balesdent, J. Morio, and J. Marzat, Kriging-based adaptive importance sampling algorithms for rare event estimation, Struct. Saf. 44, 1–10 (2013). [9] D. Aharonov and A. Ta-Shma, Adiabatic quantum state generation and statistical zero knowledge, in Proceedings of the Thirty-Fifth Annual ACM Symposium on Theory of Computing, STOC ’03 (Association for Computing Machinery, New York, NY, USA, 2003) p. 20–29. [10] A. Gilyén and T. Li, Distributional property testing in a quantum world, in 11th Innovations in Theoretical Computer Science Conference (ITCS 2020), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 151, edited by T. Vidick (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2020) pp. 25:1–25:19. [11] C. Zoufal, A. Lucchi, and S. Woerner, Quantum generative adversarial networks for learning and loading random distributions, npj Quantum Inf. 5 (2019).
[12] K. Temme and P. Wocjan, Quantized Markov chain couplings that prepare qsamples, Quantum 9, 1951 (2025). [13] B. Waggoner, ℓp testing and learning of discrete distributions, in Proceedings of the 2015 Conference on Innovations in Theoretical Computer Science, ITCS ’15 (Association for Computing Machinery, New York, NY, USA, 2015) p. 347–356. [14] J. van Apeldoorn, Quantum probability oracles & multidimensional amplitude estimation, in 16th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2021), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 197, edited by M.-H. Hsieh (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2021) pp. 9:1–9:11. [15] Z. Bar-Yossef, The Complexity of Massive Data Set Computations, Ph.D. thesis, University of California, Berkeley (2002). [16] G. Brassard, P. Høyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, in Quantum computation and information, Contemporary Mathematics, Vol. 305 (American Mathematical Society, Providence, RI, USA, 2002) pp. 53–74. [17] N. Guo, K. Mitarai, and K. Fujii, Nonlinear transformation of complex amplitudes via quantum singular value transformation, Phys. Rev. Res. 6, 043227 (2024). [18] A. G. Rattew and P. Rebentrost, Non-linear transformations of quantum amplitudes: Exponential improvement, generalization, and applications (2023), arXiv:2309.09839 [quant-ph]. [19] N. Guo, Z. Yu, M. Choi, Y. Han, A. Agrawal, K. Nakaji, A. Aspuru-Guzik, and P. Rebentrost, Quantum transformer: Accelerating model inference via quantum linear algebra (2024), arXiv:2402.16714 [quant-ph]. [20] A. G. Rattew, P.-W. Huang, N. Guo, L. Pira, and P. Rebentrost, Accelerating inference for multilayer neural networks with quantum computers, in The Fourteenth International Conference on Learning Representations
8 (2026). [21] Y. Du, X. Wang, N. Guo, Z. Yu, Y. Qian, K. Zhang, M. Hsieh, P. Rebentrost, and D. Tao, A Gentle Introduction to Quantum Machine Learning, Artificial Intelligence (R0) (Springer Nature Singapore, 2025). [22] C. H. Bennett, E. Bernstein, G. Brassard, and U. Vazirani, Strengths and weaknesses of quantum computing, SIAM J. Comput. 26, 1510 (1997). [23] A. Belovs, Quantum algorithms for classical probability distributions, in 27th Annual European Symposium on Algorithms (ESA 2019), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 144, edited by M. A. Bender, O. Svensson, and G. Herman (Schloss Dagstuhl – LeibnizZentrum für Informatik, Dagstuhl, Germany, 2019) pp. 16:1–16:11. [24] T. Bollerslev, V. Todorov, and S. Z. Li, Jump tails, extreme dependencies, and the distribution of stock returns, J. Econom. 172, 307–324 (2013). [25] I. Dobson, B. A. Carreras, V. E. Lynch, and D. E. Newman, Complex systems analysis of series of blackouts: Cascading failure, critical points, and self-organization, Chaos 17 (2007). [26] C. E. Shannon, A mathematical theory of communication, Bell Syst. Tech. J. 27, 379 (1948). [27] C. Aghamohammadi, S. P. Loomis, J. R. Mahoney, and J. P. Crutchfield, Extreme quantum memory advantage for rare-event sampling, Phys. Rev. X 8, 011025 (2018). [28] J. Blanchet, F. He, and K. Murthy, On distributionally robust extreme value analysis, Extremes 23, 317–347 (2020). [29] F. C. Binder, J. Thompson, and M. Gu, Practical unitary simulator for non-Markovian complex processes, Phys. Rev. Lett. 120, 240502 (2018). [30] C. Yang, F. C. Binder, V. Narasimhachar, and M. Gu, Matrix product states for quantum stochastic modeling, Phys. Rev. Lett. 121, 260602 (2018). [31] D. Layden, R. Sweke, V. Havlı́ček, A. Chowdhury, and K. Neklyudov, Wavefunction flows: Efficient quantum simulation of continuous flow models (2025), arXiv:2510.08462 [quant-ph]. [32] C. Yang, M. Florido-Llinàs, M. Gu, and T. J. Elliott, Dimension reduction in quantum sampling of stochastic processes, npj Quantum Inf. 11 (2025). [33] G. H. Low and I. L. Chuang, Optimal Hamiltonian simulation by quantum signal processing, Phys. Rev. Lett. 118, 010501 (2017). [34] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics, in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) pp. 193–204. [35] S. Chakraborty, A. Gilyén, and S. Jeffery, The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation, in 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 132, edited by C. Baier, I. Chatzigiannakis, P. Flocchini, and S. Leonardi (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2019) pp. 33:1–33:14. [36] C. W. Helstrom, Detection theory and quantum mechanics, Inf. Control 10, 254–291 (1967).
[37] A. Holevo, Statistical decision theory for quantum systems, J. Multivar. Anal. 3, 337–394 (1973). [38] J. Watrous, The Theory of Quantum Information (Cambridge University Press, 2018). [39] L. Lin and Y. Tong, Near-optimal ground state preparation, Quantum 4, 372 (2020). [40] Y. Dong, X. Meng, K. B. Whaley, and L. Lin, Efficient phase-factor evaluation in quantum signal processing, Phys. Rev. A 103, 042419 (2021). [41] C. Durr and P. Høyer, A quantum algorithm for finding the minimum (1999), arXiv:quant-ph/9607014 [quant-ph]. [42] J. van Apeldoorn and A. Gilyén, Improvements in quantum SDP-solving with applications, in 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 132, edited by C. Baier, I. Chatzigiannakis, P. Flocchini, and S. Leonardi (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2019) pp. 99:1–99:15. [43] A. Ahuja and S. Kapoor, A quantum algorithm for finding the maximum (1999), arXiv:quant-ph/9911082 [quantph]. [44] Y. Dong, L. Lin, and Y. Tong, Ground-state preparation and energy estimation on early fault-tolerant quantum computers via quantum eigenvalue transformation of unitary matrices, PRX Quantum 3, 040305 (2022). [45] H. Buhrman, S. Gharibian, Z. Landau, F. Le Gall, N. Schuch, and S. Tamaki, Beating the natural Grover bound for low-energy estimation and state preparation, Phys. Rev. Lett. 135, 030601 (2025). [46] L. Rabiner and B. Juang, An introduction to hidden Markov models, IEEE ASSP Mag. 3, 4 (1986). [47] P. H. Algoet and T. M. Cover, A sandwich proof of the Shannon–McMillan–Breiman theorem, Ann. Probab. 16, 899 (1988). [48] C. Aghamohammadi and J. P. Crutchfield, Minimum memory for generating rare events, Phys. Rev. E 95, 032101 (2017). [49] R. Chao, D. Ding, A. Gilyen, C. Huang, and M. Szegedy, Finding angles for quantum signal processing with machine precision (2020), arXiv:2003.02831 [quant-ph]. [50] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, Grand unification of quantum algorithms, PRX Quantum 2, 040203 (2021). [51] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019). [52] J. Wang, Y. Dong, and L. Lin, On the energy landscape of symmetric quantum signal processing, Quantum 6, 850 (2022).
9
Appendices for “Quantum enhanced rare event discovery and sampling” CONTENTS
A. Notations
9
B. Problem definition of rare event sampling
10
C. Preliminary to quantum linear algebra
11
D. Quantum and classical algorithm for rare event sampling 1. Classical algorithm 2. Quantum algorithm
12 12 16
E. Discussion 1. Polynomial approximation 2. Related works
21 21 22
F. Phase transition of quantum advantages in power-law tails 1. Rank-frequency power-law model 2. Evaluation of the rare-event mass 3. Summary of quantum advantages and interpretations
23 23 24 29
G. Application to stochastic process 1. Stochastic processes and the Asymptotic Equipartition Property 2. Thermodynamic mapping for rare events 3. Quantum rare event sampling for stochastic processes 4. Quantum stochastic modeling
29 29 31 32 33
H. Numerical results 1. Setup 2. Perturbed coin model
34 34 35 Appendix A: Notations
In this appendix, we collect the notation used throughout the paper. When we need to distinguish them, we reserve Ω for a generic finite sample space and X for the alphabet of a stochastic process. A generic finite distribution is written as P, Q : Ω → [0, 1], with Ω = {x1 , . . . , xN } and X X P (x) = 1, Q(x) = 1. (A.1) x∈Ω
x∈Ω
We write E[·] for expectation. The total variation distance (TVD) between two probability distributions P and Q on the same sample space is Dtv (P, Q) :=
1X |P (x) − Q(x)|. 2
(A.2)
x∈Ω
For quantum states we use Dirac notation: |ψ⟩ denotes a state vector, ⟨ψ| its conjugate transpose, and ⟨ψ|ϕ⟩ the inner product between |ψ⟩ and |ϕ⟩. A unitary U satisfies U † U = U U † = I. For vectors, ∥ · ∥2 denotes the Euclidean norm; for matrices, ∥ · ∥ denotes the operator norm. For any subset A ⊆ Ω, we write X ΠA := |x⟩⟨x| (A.3) x∈A
for the projector onto the computational-basis subspace indexed by A.
10 Given a probability distribution P over Ω = {x1 , . . . , xN }, a quantum sampler for P is a unitary UP satisfying UP |0⟩ =
N p X P (xi )|xi ⟩.
(A.4)
i=1
e to suppress polylogarithmic factors. We use standard asymptotic notation O(·), Ω(·), and Θ(·). We write O(·) We denote a discrete-time stochastic process by X = {Xt }t∈Z ,
(A.5)
where each Xt takes values in a finite alphabet X of size |X | < ∞. Appendix B: Problem definition of rare event sampling
In this section, we formally define the rare event sampling problem. We consider a general setting where two probability distributions, P and Q, are defined over the same sample space. One may understand the problem as “Give me samples from Q, but only if they are unlikely under P .” We start by formally defining the access model. Definition S1 (Classical and quantum sampler). Let P be a probability distribution over the sample space of N elements {x1 , x2 . . . , xN }, where P (xi ) denotes the probability of xi . We say we have access to a (classical) sampler for P if we can sample xi according to the distribution P . We define a unitary operator UP as a quantum sampler for P if UP |0⟩ =
N p X P (xi )|xi ⟩.
(B.1)
i=1
Rare events are defined with respect to the distribution P as elements with probabilities smaller than the threshold ∆ > 0. Our goal is to generate samples of these events according to the distribution Q. This distinction is important for applications involving a change of measure. For instance, if Q is the uniform distribution, we sample rare events uniformly; if Q = P , we sample them proportional to their original probabilities. To facilitate the lower bound proof, we distinguish between “yes” and “no” instances, where the “no” instance corresponds to the case where the total probability of rare events is negligible. We introduce parameters p1 and p2 to separate these cases. For simplicity in later sections, we will implicitly assume p2 = 0 and fix p1 to a suitable value without loss of generality. Furthermore, we define the ambiguous region as the set of events with probabilities in the interval [∆, 3∆/2]. We note that the upper bound 3∆/2 is chosen for concreteness; this can be generalized to (1 + α)∆ for any constant α > 0 without affecting the validity of our results. Problem S1 (General rare event sampling). Assume we are given access to samplers for probability distributions P and Q over the same sample space of N events {xi }, respectively. Let ∆ > 0 and ϵ > 0 be given parameters. We P define the set of rare events as R = {xi : P (xi ) ≤ ∆} and their total probability under Q as prare = xi ∈R Q(xi ). P Additionally, we define the boundary region S∆ = {xi : ∆ < P (xi ) < 3∆/2} and its mass q∆ = xi ∈S∆ Q(xi ). Given p1 , p2 ≥ 0 such that p1 − p2 > ϵ + q∆ , the task is to perform the following with success probability at least 2/3: 1. If prare ≥ p1 , then construct a sampler for a distribution that is Θ(ϵ + ζ)-close to QR in total variation distance, where ζ = q∆ /prare and ( Q(xi )/prare , xi ∈ R, QR = (B.2) 0, xi ∈ / R. 2. If prare ≤ p2 , then output “Impossible”. The rare-event sampling problem considered in the main text is recovered as a special case of Problem S1. In the generalized formulation, the distribution P defines which events are rare, while the distribution Q defines how these rare events should be sampled. The main text focuses on the natural self-sampling case, where the same distribution plays both roles: Q = P . Under this specialization, the rare-event set becomes R∆ := {x : 0 < P (x) ≤ ∆},
(B.3)
11 P and its total probability mass is prare := P (R∆ ) = x∈R∆ P (x). The target distribution QR in Problem S1 then reduces to the conditional rare-event distribution ( P (x)/prare , x ∈ R∆ , PR (x) = (B.4) 0, x∈ / R∆ . Similarly, the ambiguous region becomes S∆ := {x : ∆ < P (x) ≤ 3∆/2}, with mass p∆ := P (S∆ ), so that the unavoidable near-threshold error parameter is ζ := p∆ /prare . The parameters p1 and p2 are promise parameters used only in the generalized decision version of the problem. They allow the algorithm to distinguish between the case where rare events have sufficient total mass and the case where rare-event sampling is impossible. The main-text problem is the sampling-only, yes-instance regime, where the rare set is assumed to have nonzero mass. Equivalently, if one wants to embed the main-text problem into Problem S1, one may take p2 = 0, p1 = pmin ≤ prare , where pmin > 0 is a promised lower bound. The additional promise condition in Problem S1, p1 − p2 > ϵ + q∆ , is needed only if the algorithm is required to certify the “Impossible” case. When the “Impossible” branch is omitted, as in the main text, the algorithm and its complexity are naturally stated directly in terms of the actual rare-event mass prare . Problem S2 (Rare event sampling). Let P be a probability distribution over a finite sample space Ω = {x1 , . . . , xN }, and let ∆ > 0 be a probability threshold. Define the rare-event set as Eq. (B.3), and we assume that the rare-event mass is nonzero prare > 0. The target rare-event distribution is the conditional distribution ( P (x)/prare , x ∈ R∆ , PR (x) := (B.5) 0, x∈ / R∆ . Define the ambiguous region S∆ := {x ∈ Ω : ∆ < P (x) ≤ (1+α)∆}, with α > 0 a constant and total mass p∆ := P (S∆ ). We write ζ := p∆ /prare for the relative near-threshold ambiguity. Given a precision parameter ϵ > 0, the task is to e construct a sampler whose output distribution PR satisfies Dtv PeR , PR = O(ϵ + ζ). Appendix C: Preliminary to quantum linear algebra
In this section, we introduce quantum linear algebra, including block encoding, amplitude encoding, and quantum singular value transformation. Block encoding is an encoding scheme that embeds a general matrix A into a block of a unitary matrix. Without loss of generality, we focus on the top left block of the unitary matrix U , that is, (⟨0a | ⊗ In )U (|0a ⟩ ⊗ In ). Here, a is the number of required ancillary qubits. n
n
Definition S2 (Block encoding [34, 35]). We say a unitary UA is an (αBE , a, ϵ)-encoding of matrix A ∈ C2 ×2 if ∥A − αBE (⟨0a | ⊗ In )UA (|0a ⟩ ⊗ In )∥ ≤ ϵ.
(C.1)
In our work, we focus on a specific type of block encoding, so-called the amplitude block encoding, the unitary matrix containing a diagonal matrix A filled with the amplitude of a given quantum state, namely, A = diag(ψ1 , . . . , ψN ). PN Given access to a quantum circuit Uψ that prepares the quantum state |ψ⟩ = j=1 ψj |j⟩, the amplitude block encoding can be efficiently constructed. Since we focus on the stochastic process, without loss of generality, we consider amplitudes as real values. Theorem S1 (Amplitude block encoding [17, 18]). Given access to a n-qubit quantum circuit Uψ that prepares an PN n-qubit state Uψ : |0⟩ → |ψ⟩ = j=1 ψj |j⟩, where {ψj } are real and ∥ψ∥2 = 1, one can construct an (1, n + 2, 0)encoding of the diagonal matrix A = diag(ψ1 , . . . , ψN ) with O(n) circuit depth and O(1) queries to controlled-Uψ and controlled-Uψ† . Once given a block encoding UA of a Hermitian matrix A, we can construct a unitary operator that is a block encoding of applying certain polynomials g to the matrix A. QSVT(Φ, UA ) =
g(A) · · ·
(C.2)
where Φ is a set of phase angles required to implement the polynomial g. In general, that is, when A is not a Hermitian matrix, g(A) is a matrix that is obtained by applying the polynomial g to the singular values of A. The circuit depth scales with the degree of the polynomial. More formally, we have the following result.
12 Lemma S1 (Polynomial eigenvalue transformation [34]). Given U that is an (α, a, ϵ)-encoding of a Hermitian matrix A, and a real d-degree polynomial g(x) with |g(x)| ≤ 21 for x ∈ [−1, 1]. Let δ > 0, one can prepare a p (1, a + n + 4, 4d ϵ/α + δ)-encoding of g(A/α) by using O(d) queries to UA and O(d(a + 1)) one- and two-qubit quantum gates. The description of the quantum circuit can be computed classically in time O(poly(d, log(1/δ))). Therefore, we can implement many useful functions on matrices by finding their good polynomial approximations. Here, “good” means that the degree of the polynomial scales logarithmically with the error in the interval [−1, 1]. In general, especially for functions that are not smooth in this interval, such polynomial approximations do not exist. However, a good approximation may still be found on the subset of [−1, 1]. In the following, we list some known results for the sign and rectangle functions. The rectangle function f (x) is defined as ( 1, x ∈ [−t, t] ft (x) := (C.3) 0, x ∈ [−1, −t) ∪ (t, 1]. Lemma S2 (Polynomial approximation of the rectangle function [34]). For any error ϵ ∈ (0, 1/2), there exist an even d-degree polynomial gd ∈ R[x] that approximates the rectangular function ft (x), such that 0 ≤ gd (x) ≤ 1, ∀x ∈ [−1, 1], and |ft (x) − gd (x)| ≤ ϵ,
∀x ∈ / [−t − Γ, −t + Γ] ∪ [t − Γ, t + Γ].
(C.4)
where Γ ∈ (0, 1/2) denotes the size of ambiguity range, and the degree d scales as O Γ1 log (1/ϵ) . Appendix D: Quantum and classical algorithm for rare event sampling
Here, we provide the classical and quantum algorithm for Problem S1. We also establish query complexity lower bounds for both cases, demonstrating that our quantum algorithm is optimal with respect to the dominant factor defining the rare event. Furthermore, we show a quadratic separation between the classical and quantum complexities. 1.
Classical algorithm
Theorem S2 (Classical algorithm for general rare event sampling). There exists a classical algorithm that solves Problem S1 using N 1 1 1 log O min , 2 log (D.1) ∆ ϵ ∆ ϵ samples from P and O(max{1/(p1 − p2 − q∆ )2 , 1/prare }) samples from Q. Proof. We first construct a set of candidate rare events, which will be used to distinguish between the two cases and subsequently for rejection sampling. We now show that the candidate rare set can be constructed using O(∆−1 log(N/δ)) samples from P . The point is that we do not need additive ℓ∞ -estimation of the whole distribution. It suffices to perform threshold classification, allowing arbitrary behavior in the ambiguous region. Draw T independent samples from P , and let T 1X PbT (x) = 1{Xt = x} T t=1
be the empirical frequency of event x. Define the candidate set 4∆ b b R := x : PT (x) ≤ . 3
(D.2)
(D.3)
First, fix x with P (x) ≤ ∆. Let Nx = T PbT (x). Then Nx is binomial with mean T P (x). Since the event PbT (x) > 4∆/3 is increasing in P (x), its probability is maximized over P (x) ≤ ∆ at P (x) = ∆. Therefore, by the multiplicative Chernoff bound, 4∆ 4T ∆ T∆ b Pr PT (x) > ≤ Pr Bin(T, ∆) > ≤ exp − . (D.4) 3 3 27
13 Thus a truly rare event is falsely excluded with probability at most exp(−T ∆/27). Second, fix x with P (x) ≥ 3∆/2. The event PbT (x) ≤ 4∆/3 is a lower-tail event, and its probability is maximized when P (x) = 3∆/2. By the Chernoff bound, we have 4∆ 3∆ 4T ∆ T∆ b Pr PT (x) ≤ ≤ Pr Bin T, ≤ ≤ exp − . (D.5) 3 2 3 108 Therefore, a non-rare event with P (x) ≥ 3∆/2 is falsely included with probability at most exp(−T ∆/108). Taking a union bound over all N events, the probability that any decisive event is misclassified is at most T∆ T∆ T∆ N exp − + N exp − ≤ 2N exp − . 27 108 108
(D.6)
b ⊆ R∆ ∪ S∆ with probability at least 1 − δ, where Thus, choosing T = O (1/∆ log(N/δ)) ensures that R∆ ⊆ R R∆ = {x : P (x) ≤ ∆} and S∆ = {x : ∆ < P (x) ≤ 3∆/2}. An alternate approach to this is to first learn the distribution directly without the two-level sampling to precision ∆/6 in the ℓ∞ distance with probability 1 − δ with sample complexity Θ( ∆12 log( 1δ )) by the algorithm in Ref. [13], b such that R∆ ⊆ R b ⊆ R∆ ∪ S∆ with probability 1 − δ by taking which is also able to give us the set R 1 1 log T =O . (D.7) ∆2 δ Combining these two results, it suffices to take N 1 1 1 log , 2 log . T = O min ∆ δ ∆ δ
(D.8)
Now we consider distinguishing between the two cases prare ≥ p1 and prare ≤ p2 . On the good event of the b = R∆ ∪ Se∆ , where Se∆ ⊆ S∆ . threshold-classification step, the constructed candidate set has the form R e Let qe∆ := Q(S∆ ), we have 0 ≤ qe∆ ≤ q∆ . If prare ≥ p1 , then b ≥ prare ≥ p1 . Q(R)
(D.9)
b ≤ p2 + qe∆ ≤ p2 + q∆ . Q(R)
(D.10)
If prare ≤ p2 , then Thus the gap between the two cases is at least g := p1 − p2 − q∆ .
(D.11)
b to constant accuracy relative to this gap requires By the promise condition, g > 0. Estimating Q(R) 1 1 =O O g2 (p1 − p2 − q∆ )2
(D.12)
samples from Q. Equivalently, since the promise gives p1 − p2 − q∆ > ϵ, one may state this step as requiring O(1/ϵ2 ) samples from Q. b The accepted distribution is If the first case holds, we then use rejection sampling based on R. Q(x) , x ∈ R ∪ Se , ∆ ∆ e Q(x) = prare + qe∆ (D.13) 0, otherwise. The acceptance probability is prare + qe∆ , so one accepted sample requires O(1/(prare + qe∆ )) ∈ O(1/prare ) samples from Q. e and the ideal rare-event distribution QR satisfies The total variation distance between Q X X 1 Q(x) Q(x) Q(x) e QR ) = Dtv (Q, − + 2 prare + qe∆ prare prare + qe∆ x∈R∆
x∈Se∆
14 qe∆ prare + qe∆ q∆ ≤ = ζ. prare
=
(D.14)
Including the failure probability δ of the threshold-classification step gives total error at most ζ + δ. Taking δ = O(ϵ), the final total variation error is O(ϵ + ζ). In total, the number of samples from P is 1 1 N 1 O min log , 2 log , (D.15) ∆ ϵ ∆ ϵ and the number of samples from Q is O max
1
1 , (p1 − p2 − q∆ )2 prare
.
(D.16)
Following this, one can immediately provide the classical algorithm for rare event sampling. Corollary S1 (Classical algorithm for rare event sampling). There exists a classical algorithm that solves Problem S2 using 1 N 1 1 1 O min log , 2 log + ∆ ϵ ∆ ϵ prare samples from P . Proof. The result follows from the construction in Theorem S2 specialized to the self-sampling case Q = P . In Problem S2, there is no “Impossible” branch, so the decision step distinguishing prare ≥ p1 from prare ≤ p2 is unnecessary. Using the threshold-classification step from the proof of Theorem S2, with probability at least 1 − δ one constructs a b satisfying R∆ ⊆ R b ⊆ R∆ ∪ S∆ using candidate set R 1 N 1 1 O min log , 2 log (D.17) ∆ δ ∆ δ samples from P . b = R∆ ∪ Se∆ for some Se∆ ⊆ S∆ . Rejection sampling from P using R b therefore accepts with On this good event, R b e probability P (R) = prare + pe∆ , where pe∆ := P (S∆ ) ≤ p∆ . Hence one accepted sample requires O(1/prare ) samples from P . The same TVD calculation as in the proof of Theorem S2 gives Dtv (PeR , PR ) ≤
p∆ pe∆ ≤ = ζ. prare + pe∆ prare
(D.18)
Including the failure probability of the candidate-set construction gives error at most ζ + δ. Taking δ = O(ϵ) yields total variation error O(ϵ + ζ) and gives the stated sample complexity. Now we prove the classical lower bound, where the proof is achieved by reduction from a standard distribution distinguishing problem. Lemma S3 (Lower bound for distribution distinguishing [15]). For any two probability distributions P1 and P2 over the same sample space, any classical algorithm that distinguishes P1 from P2 requires Ω(1/d2H (P1 , P2 )) samples, where s 2 p 1 X p dH (P1 , P2 ) := P1 (i) − P2 (i) (D.19) 2 i is the Hellinger distance between P1 and P2 .
15 Theorem S3 (Classical lower bound in ∆). Any classical algorithm that solves Problem S1 or Problem S2 with 0 < ∆ < 14 requires Ω(1/∆) samples drawn from P . Proof. We prove the lower bound using a single pair of yes instances. This avoids relying on the “Impossible” branch and therefore applies both to the general rare-event sampling problem and to the self-sampling version. Fix a constant output accuracy parameter η < 1/10. Consider the following two probability distributions on the sample space {1, 2, 3}: ∆ 5∆ P0 = , 2∆, 1 − , (D.20) 2 2 ∆ 5∆ P1 = 2∆, , 1 − . (D.21) 2 2 For 0 < ∆ < 1/4, both are valid probability distributions, and 1−
3∆ 5∆ > . 2 2
(D.22)
Hence, for both distributions, the boundary region S∆ = {x : ∆ ≤ P (x) ≤ 3∆/2} is empty. We first compute the distance between the two input distributions. Their squared Hellinger distance is r !2 r !2 √ √ ∆ ∆ 1 2∆ − − 2∆ + d2H (P0 , P1 ) = 2 2 2 =
√
r 2∆ −
∆ 2
!2 =
∆ = Θ(∆). 2
(D.23)
By Lemma S3, any classical algorithm that distinguishes P0 from P1 with constant success probability requires Ω(1/∆) samples. We now show that any rare-event sampling algorithm would distinguish these two distributions. In the general problem Problem S1, choose the fixed distribution 1 1 1 Q= , ,0 , p1 = , p2 = 0. (D.24) 2 2 3 Since the boundary region is empty, we have q∆ = 0. The promise condition is satisfied for any sufficiently small constant precision parameter, because p1 − p2 = 1/3 > 0 = q∆ . If the unknown distribution is P0 , then the rare set is R∆ (P0 ) = {1}, and hence the target rare-event distribution is the delta function δ1 . Moreover, prare = Q(R∆ (P0 )) = Q({1}) =
1 ≥ p1 . 2
(D.25)
If the unknown distribution is P1 , then the rare set is R∆ (P1 ) = {2}, and hence the target rare-event distribution is the delta function δ2 . Similarly, prare = Q(R∆ (P1 )) = Q({2}) =
1 ≥ p1 . 2
(D.26)
Thus both inputs are yes instances of the general problem, but their correct output samplers are different. The same hard pair also applies to the self-sampling problem, where Q = P . Indeed, if P = P0 , then the rare set is again {1} and the conditional rare-event distribution is P0,R = δ1 . If P = P1 , then the rare set is again {2} and the conditional rare-event distribution is P1,R = δ2 . In both cases the rare-event mass is nonzero, i.e., prare = ∆/2 > 0. Therefore the pair P0 , P1 consists of valid yes instances also for the sampling-only self-sampling problem. Now suppose that there exists a classical algorithm A that solves either the general problem above or the self-sampling problem using S samples from the unknown distribution P . We use A to distinguish whether the unknown input is P0 or P1 . Run A on the corresponding rare-event sampling instance, then draw one sample from the sampler produced by A. If the sample is 1, output P0 . If the sample is 2, output P1 . In all other cases, output either distribution arbitrarily. Conditioned on A succeeding, its output distribution is within total variation distance at most η from the correct rare-event distribution. Therefore, if P = P0 , the produced sampler outputs 1 with probability at least 1 − η. If P = P1 , the produced sampler outputs 2 with probability at least 1 − η. Since η < 1/10 and A succeeds with constant
16 probability, this gives a distinguisher for P0 and P1 with constant success probability. By the Hellinger-distance lower bound above, such a distinguisher requires 1 S=Ω . (D.27) ∆ Therefore any classical algorithm solving the rare-event sampling problem requires Ω(1/∆) samples from P . The conclusion holds both for the general problem and for the sampling-only self-sampling problem. Note that the above proof naturally generalizes to the case where we can only sample from the rare distribution approximately. Using the hard instance in the proof of Theorem S3, PR = (1, 0, 0). Performing the rare event sampling algorithm on distributions P1 and P2 , we output failure and distribution (1, 0, 0), respectively. Even if we can only obtain the final output to precision 1/3 in total variation distance, we can still distinguish P1 and P2 by taking a constant number of samples from the output distribution.
2.
Quantum algorithm
Theorem S4 (Quantum algorithm for rare event sampling). There exists a quantum algorithm that solves Problem S1 using 1 1 1 1 √ log O max √ (D.28) √ ,√ p1 − p2 prare ϵ ∆ n o 1√ √ 1 queries to UP and O max √p1 − queries to UQ with O(n) ancilla qubits. Here, UP and UQ denote the p2 , prare quantum samplers for P and Q, respectively. † Proof. By Theorem S1, we can use O(1) p times of controlled-UP and controlled-UP to construct an (1, n + 2, 0)-encoding UA of a diagonal matrix A = diag({ p(xi )}). The next step is to apply the rectangular function to the diagonal matrix A of the encoding UA , which approximately constructs a projection ΠR into the space filled with rare events. Note that in the quantum setting, we operate on the square roots of the probabilities rather than the probabilities √ themselves. To ensure the ambiguity range is contained in S∆ , for quantum we choose the ambiguity range as 0.2 ∆ such that √ √ ( ∆ + 0.2 ∆)2 = 1.44∆ < 1.5∆. (D.29)
Lemma S2 guarantees that the rectangular function can error ϵ by an even polynomial P with √ √ be approximated to √ degree O( √1∆ log(1/ϵ)), with t = ( ∆ + ξ)/2 = 1.1 ∆ and Γ = ξ/2 = 0.1 ∆. Lemma S1 claims that such an even polynomial can be applied to the diagonal matrix A by using O( √1∆ log(1/ϵ)) queries to controlled-UA and controlled-UA† . With these, we can approximately construct a block encoding of the operator ΠR =
X xi ∈R
|xi ⟩⟨xi |.
(D.30)
Then, we apply the block encoding UΠR to the initial state prepared by UQ , that is, UΠR (I ⊗ UQ )|0⟩|0⟩ = UΠR |0⟩
Xp xi
q(xi )|xi ⟩ ≈ |0⟩
X p xi ∈R
e q(xi )|xi ⟩ + |1⟩|⊥⟩,
(D.31)
e is an arbitrary unnormalized state. For simplicity, write |ϕ⟩ = UQ |0⟩. The probability of measuring the where |⊥⟩ ancilla qubits as |0⟩ is ∥Π̃R |ϕ⟩∥2 ≥ ∥ΠR |ϕ⟩∥2 (1 − ϵ)2 = prare (1 − ϵ)2 ,
(D.32)
where Π̃R is the actual constructed approximated projector. We also have ∥Π̃R |ϕ⟩∥2 ≤ ∥ΠR |ϕ⟩∥2 + q∆ = prare + q∆ + ϵ2 .
(D.33)
17 To decide whether we are in the first or the second case as required in Problem S1, it suffices to achieve this by measuring ∥Π˜R |ϕ⟩∥2 with amplitude estimation [16]. To separate these two cases, we need to make sure p1 (1 − ϵ)2 ≥ p2 + q∆ + ϵ2 .
(D.34)
p1 (1 − ϵ)2 = p1 (1 − 2ϵ + ϵ2 ) ≥ p1 (1 − 2ϵ),
(D.35)
Note that
it suffices to take p1 − p2 > 2ϵ + q∆ ≥ 2ϵp1 + q∆ . One can further rescale 2ϵ to ϵ. Note that Eq. (D.31) is also a state preparation unitary. One can ignore the second register and only consider the first register, which can be written as q Uprep |0⟩ = ∥P̃rare |ϕ⟩∥|0⟩ + 1 − ∥P̃rare |ϕ⟩∥2 |1⟩. (D.36) √ √ √ √ p + p By using Theorem S1 and implementing the approximate rectangular function f (x− 1 2 2 ), with Γ′ = ( p1 − p2 )/2, we can boost the gap between p1 and p2 to Ω(1). Then we can use amplitude amplification [16] with precision O(1) to distinguish the first case and the second case with probability at least 2/3. To distinguish whether we are in the first 1 √ √1 1 √ 1√ case or in the second case uses O( √p1 − p2 ∆ log ϵ ) queries to UP and O( p1 − p2 ) queries to UQ . Now we consider if we are in the first case. By further using the amplitude amplification O( √p1rare ) times [16], one can prepare the target state
√
X p 1 Q(xi )|xi ⟩. prare
(D.37)
xi ∈R
Combining the QSVT procedure, the overall queries to UP is O( √∆p1 log 1ϵ ) and queries to UQ is O( √p1rare ). rare 2 2 e R |ϕ⟩ = P ⟨xi |Π e R |ϕ⟩ be the actual post-selection probability. The resulting post-selected distribuLet r̃ := Π i tion is e R (xi ) = Q
e R |ϕ⟩ ⟨xi |Π
2
.
r̃
(D.38)
Let Sun := Ω \ (R ∪ S∆ ) denote the non-rare region. By the polynomial approximation guarantee, for xi ∈ R, e R |ϕ⟩ (1 − ϵ)2 Q(xi ) ≤ ⟨xi |Π
2
≤ Q(xi ).
(D.39)
For xi ∈ S∆ , we use the trivial bound e R |ϕ⟩ ⟨xi |Π
2
≤ Q(xi ),
(D.40)
≤ ϵ2 Q(xi ).
(D.41)
and for xi ∈ Sun , e R |ϕ⟩ ⟨xi |Π
2
Therefore, (1 − ϵ)2 prare ≤ r̃ ≤ prare + q∆ + ϵ2 .
(D.42)
|r̃ − prare | ≤ 1 − (1 − ϵ)2 prare + q∆ + ϵ2 .
(D.43)
Moreover,
For the rare region, we first note that e R |ϕ⟩ ⟨xi |Π X Q(xi ) − prare r̃
xi ∈R
2
18 2
≤
X xi ∈R
e R |ϕ⟩ X Q(xi ) − ⟨xi |Π 1 + Q(xi ) − . prare r̃ r̃ 1
(D.44)
xi ∈R
e R |ϕ⟩ The second term is nonnegative because ⟨xi |Π
2
≤ Q(xi ) for xi ∈ R.
eR : We evaluate the total variation distance between the target distribution QR and the actual distribution Q X eR ) = 1 e R (xi ) Dtv (QR , Q QR (xi ) − Q 2 x i 2 e ⟨x | Π |ϕ⟩ 2 X X X i R 1 Q(xi ) Q(xi ) ϵ Q(xi ) ≤ + − + 2 prare r̃ r̃ r̃ xi ∈S∆
xi ∈R
≤
1X 2
xi ∈R
Q(xi )
1 1 + − prare r̃
xi ∈Sun
e R |ϕ⟩ X Q(xi ) − ⟨xi |Π xi ∈R
r̃
2
+
q∆ ϵ + r̃ r̃
) 1 − (1 − ϵ)2 prare |r̃ − prare | q∆ ϵ2 + + + r̃ r̃ r̃ r̃ ( ) 1 2 1 − (1 − ϵ)2 prare + 2q∆ + 2ϵ2 ≤ 2 r̃ 2 1 − (1 − ϵ) prare + q∆ + ϵ2 ≤ (1 − ϵ)2 prare 2 2ϵ − ϵ ζ ϵ2 = + + , (1 − ϵ)2 (1 − ϵ)2 (1 − ϵ)2 prare 1 ≤ 2
2
(
(D.45)
where ζ = q∆ /prare . We choose the polynomial approximation error such that ϵ ≤ p1 . In the yes case of Problem S1, prare ≥ p1 . Hence ϵ2 prare
≤
ϵ2 ≤ ϵ. p1
(D.46)
Therefore, assuming ϵ ≤ 1/2 without loss of generality, e R ) = O(ϵ + ζ). Dtv (QR , Q
(D.47)
If we are in the second case, we can directly output “impossible”. Combining these two cases, one needs 1 1 1 1 √ , log (D.48) O max √ √ √ p1 − p2 prare ϵ ∆ n o 1√ √ 1 queries to UP , and O max √p1 − queries to UQ . p2 , prare Following this, we can immediately write the case for rare event sampling. Corollary S2 (Quantum algorithm for rare event sampling). There exists a quantum algorithm that solves Problem S2 using 1 1 O √ log (D.49) ϵ prare ∆ queries to UP and UP† with O(n) ancilla qubits. Here, UP denote the quantum samplers for P . Proof. This is the self-sampling case Q = P of Theorem S4, with the “Impossible” decision branch omitted. In Problem S2, the rare-event mass prare is assumed to be nonzero, and the goal is only to prepare a sampler close to the conditional rare-event distribution.
19
Quantum algorithm for rare event sampling c-UP /UP† ∆ ϵ
Inputs: – UP generates the quantum sample state |P ⟩. – Threshold for rare events. – Error parameter.
Outputs: |PR ⟩ – A quantum rare-event sample state that superposes all rare events. Algorithm: 1. By Theorem S1, use c-UP /UP† a constant number of times to construct the diagonal block encoding of amplitudes UA . √ 2. Implement a polynomial approximation of the rectangle function f (x) via Lemma S1, using O(1/ ∆ log(1/ϵ)) applications of UA . 3. Feed the quantum sample state |P ⟩ into the constructed quantum circuit, and use amplitude amplification to boost the √ success probability and prepare the target state |PR ⟩. This further uses O(1/ prare ) applications of the constructed √ quantum circuit. In total, the algorithm uses O(1/ ∆prare log(1/ϵ)) applications of c-UP /UP† .
FIG. S1. Description of the quantum rare-event sampling algorithm, with the special case P = Q.
The thresholding construction from Theorem S4 uses 1 1 O √ log ϵ ∆
(D.50)
queries to UP and UP† to implement the approximate rare-event projector. Applying this projector to the quantum √ sample state UP |0⟩ succeeds with probability O(prare ). Amplitude amplification therefore adds a factor O(1/ prare ). Hence the total query complexity is 1 1 O √ log . (D.51) ϵ prare ∆ The total variation error bound is the same as in the proof of Theorem S4, with Q = P , giving O(ϵ + ζ). In the following, we provide the lower bound proof for the quantum case. Theorem S5 (Quantum lower bound in ∆). For any 0 < ∆ < 1/4, any quantum algorithm that solves Problem S1 with quantum samplers UP , UQ , parameters p1 = 2/3, p2 = 1/3, and 0 < ϵ < 1/3, requires 1 Ω √ ∆ queries to UP and UP† , even when Q is fixed and known. Proof. We prove the lower bound by reducing a simple quantum-sampler distinguishing task to rare-event sampling. Consider two probability distributions P0 and P1 over the sample space Ω = {0, 1}: P0 (1) = ∆/2 and P0 (0) = 1 − ∆/2; P1 (1) = 2∆ and P1 (0) = 1 − 2∆. Since 0 < ∆ < 1/4, these are valid probability distributions. Also we have P0 (1) < ∆, P1 (1) > 3∆ 2 , and both P0 (0) and P1 (0) are larger than 3∆/2. Therefore the boundary region S∆ is empty for both instances. For b ∈ {0, 1}, let the corresponding quantum samplers be the one-qubit rotations p p Ub |0⟩ = 1 − Pb (1)|0⟩ + Pb (1)|1⟩. (D.52) Equivalently, write Ub |0⟩ = cos θb |0⟩ + sin θb |1⟩,
(D.53)
20 √ 2∆. For 0 < ∆ < 1/4, this gives |θ1 − θ0 | = Θ( ∆). We first show that √ distinguishing these two quantum samplers requires Ω(1/ ∆) queries. Choose the canonical rotation implementation cos θb − sin θb Ub = . (D.54) sin θb cos θb
where θ0 = arcsin
q
∆ 2
and θ1 = arcsin
√
Then √ θ1 − θ 0 = O(|θ1 − θ0 |) = O( ∆). (D.55) 2 √ The same bound holds for the inverse unitaries ∥U0† − U1† ∥ = O( ∆). Consider an arbitrary quantum algorithm making T queries to the unknown sampler and its inverse, interleaved (T ) (T ) with arbitrary known unitaries. Let |Ψ0 ⟩ and |Ψ1 ⟩ denote the final pure states of the algorithm when the oracle is U0 and U1 , respectively. By the standard hybrid argument for quantum query algorithms [22], replacing the oracle calls one by one gives n o √ (T ) (T ) (D.56) |Ψ0 ⟩ − |Ψ1 ⟩ ≤ T · max ∥U0 − U1 ∥, ∥U0† − U1† ∥ = O(T ∆). ∥U0 − U1 ∥ = 2 sin
On the other hand, if the algorithm distinguishes the two cases with success probability at least 2/3, then the final states must have constant trace distance. Indeed, by the Helstrom–Holevo theorem [36–38], the optimal success probability for distinguishing two equally likely states ρ0 , ρ1 is psucc =
1 1 + ∥ρ0 − ρ1 ∥1 . 2 4 (T )
(D.57) (T )
Thus psucc ≥ 2/3 implies ∥ρ0 − ρ1 ∥1 ≥ 2/3. For pure states ρb = |Ψb ⟩⟨Ψb |, we have r 2 1 (T ) (T ) (T ) (T ) ∥ρ0 − ρ1 ∥1 = 1 − ⟨Ψ0 |Ψ1 ⟩ ≤ |Ψ0 ⟩ − |Ψ1 ⟩ . 2
(D.58)
Therefore, ∥ρ0 − ρ1 ∥1 ≥ 2/3 implies (T )
(T )
|Ψ0 ⟩ − |Ψ1 ⟩ ≥
1 = Ω(1). 3
√ Combining this with the hybrid bound yields O(T ∆) = Ω(1), and therefore 1 T =Ω √ . ∆
(D.59)
(D.60)
We now show that any rare-event sampling algorithm would distinguish U0 from U1 . Set Q = (0, 1), p1 = 2/3, p2 = 1/3, and 0 < ϵ < 1/3. Since the boundary region is empty for both instances, we have q∆ = 0. Hence the promise condition p1 − p2 > ϵ + q∆ is satisfied. Suppose there exists a quantum algorithm A that solves Problem S1 using S queries to UP and UP† . We use A to distinguish whether the unknown sampler is U0 or U1 . Run A with quantum sampling access to the unknown P , and with the fixed known sampler UQ for Q = (0, 1). If A outputs “Impossible”, we output P1 . Otherwise, we output P0 . If P = P0 , then outcome 1 is rare because P0 (1) = ∆ 2 < ∆. Since Q = (0, 1), the Q-mass of the rare set is prare = Q({1}) = 1 ≥ p1 . Therefore A must construct a sampler for the rare-event distribution with probability at least 2/3. If P = P1 , then no outcome is rare: P1 (1) = 2∆ > 3∆/2, while P1 (0) > 3∆/2. Hence prare = Q(∅) = 0 ≤ p2 . By the definition of Problem S1, A must output “Impossible” with probability at least 2/3. Thus A distinguishes the two quantum samplers U0√and U1 with constant success probability using S queries to UP and UP† . Since this distinguishing task requires Ω(1/ ∆) queries, we conclude that 1 S=Ω √ . ∆
21 Corollary S3 (Quantum lower bound for self-sampling rare-event sampling). Fix a constant α > 0 in the definition S∆ = {x : ∆ < P (x) ≤ (1 + α)∆}. For all sufficiently small ∆ > 0, any quantum algorithm that solves Problem S2 with constant success probability and constant accuracy requires 1 Ω √ ∆ queries to UP and UP† . This holds even for distributions with empty ambiguous region. Proof. Choose any constant c > 1 + α, and assume ∆ < 1/(c + 3/2 + α). Consider two distributions on Ω = {1, 2, 3}:
∆ 1 , c∆, 1 − c + ∆ , 2 2 ∆ 1 P1 = c∆, , 1 − c + ∆ . 2 2
P0 =
(D.61) (D.62)
For both distributions, the third probability is larger than (1 + α)∆, and c∆ > (1 + α)∆. Hence the ambiguous region is empty. Under P0 , the rare set is {1} and the conditional rare-event distribution is δ1 . Under P1 , the rare set is {2} and the conditional rare-event distribution is δ2 . The squared Hellinger distance between the two input distributions is r !2 r !2 √ √ 1 ∆ ∆ − c∆ + (D.63) c∆ − d2H (P0 , P1 ) = 2 2 2 =
√
1 c− √ 2
2 ∆ = Θ(∆).
(D.64)
By the quantum distribution-distinguishing lower bound of Ref. [23], distinguishing quantum samplers for P0 and P1 requires 1 1 =Ω √ Ω (D.65) dH (P0 , P1 ) ∆ queries. It remains to reduce √ this distinguishing task to rare-event sampling. Suppose there were an algorithm A solving Problem S2 with o(1/ ∆) queries. Given an unknown sampler for either P0 or P1 , run A and draw one sample from the sampler it outputs. If the sample is 1, output P0 ; if the sample is 2, output P1 ; otherwise output arbitrarily. Conditioned on the success of A, the output distribution is within constant total variation distance of either δ1 or δ2 . Thus this procedure distinguishes P0 from P1 with constant success probability. This contradicts the Hellinger-distance √ lower bound above. Therefore any quantum algorithm solving Problem S2 requires Ω(1/ ∆) queries.
Appendix E: Discussion 1.
Polynomial approximation
Here, we provide some discussions of the polynomial approximation used in the quantum algorithm. In this paper, we focus on the stochastic process, which implies that all quantum amplitudes (corresponding to the square root of the classical distribution) are real and positive. To achieve the rare event-related tasks, we need a good polynomial approximation for the function ( √ 1 for all x ∈ [0, ∆], √ f (x) = (E.1) 0 for all x ∈ ( ∆, 1], and |f (x)| ≤ 1 for x ∈ [−1, 0). The shifted Heaviside function is a good candidate. There are two ways for approximating the Heaviside function: one is achieved by the error function [33], which is an approximation of the sign function, and another is achieved by taking the integral of the filtering function [39], which is an approximation of the delta function.
22 Response function and error for symmetric QSP protocol 1.0
0.5
Matrix element
target polynomial target function Im[Poly(a)]
0.0 −0.5
0.0 −0.5 −1.0
10−2
10−2
10
−6
QSP/Poly error
10−10
QSP/target error
10−14 10−18 −1.00
−0.75
−0.50
−0.25
0.00 0.25 Input signal
0.50
0.75
1.00
target polynomial target function Im[Poly(a)]
0.5
−1.0
Absolute error
Absolute error
Matrix element
1.0
Response function and error for symmetric QSP protocol
10−6
QSP/Poly error
10−10
QSP/target error
10−14 10−18 −1.00
−0.75
−0.50
−0.25
0.00 0.25 Input signal
0.50
0.75
1.00
FIG. S2. Polynomial approximation for the sign function (left) and rectangle function (right).
Note that for the sign function case, one needs to add a constant 1. Asymptotically, the sign function approximation is better than the second method. However, from practical performance, we notice that directly approximating the rectangle function ( √ √ 1 for all x ∈ [− ∆, ∆], √ √ (E.2) g(x) = 0 for all x ∈ [−1, − ∆) ∪ ( ∆, 1], via optimization-based angle finding methods works better, which is somewhat surprising, as one usually constructs the rectangle function by taking linear combinations of the Heaviside function in theoretical proofs. The reason is that in practice, we use an optimization-based method to find the polynomial approximation. Especially for the quasi-Newton method [40], which is currently the state-of-the-art method to find many polynomial approximations, one can see from Fig. S2 that the overall approximation for the rectangle function is much better than the sign function. This difference in approximation error leads to the result that when we use the sign function in the numerical experiment, it requires a much higher polynomial degree to compress the largest amplitudes than the rectangle function case. Further, the rectangle function needs fewer ancilla qubits than the sign function when we implement it via QSVT in the block-encoding setting. The reason is that to implement the shifted sign function, one needs to make a shift and take a linear combination with the identity, as the shifted sign function is neither odd nor even and cannot be directly implementable via QSVT phase angles. This requires two ancilla qubits, and the coefficient from LCU will make amplitudes smaller. For the rectangle function, as we can directly find the polynomial approximation in practice, it does not require additional ancilla qubits. Therefore, in this paper, we choose to use the rectangle function. The above discussions do not affect much for the asymptotic results, which are described in the theorems. However, for practical considerations like the numerical experiment, these do make a difference. 2.
Related works
The rare event sampling problem can be decomposed into two distinct tasks: first, identifying the rare events, and second, sampling from them. Our quantum algorithm achieves efficiency by performing these steps coherently, thereby avoiding the need to explicitly read out or measure the rare events. In contrast, relying on the standard amplitude estimation algorithm [16] would incur a linear overhead in the size of the sample space, as it necessitates estimating the probability of each event individually. Even utilizing multidimensional quantum amplitude estimation [14] proves suboptimal, as it requires O(1/∆) queries to estimate all probabilities to the necessary precision of O(∆). Further, to achieve a quadratic speedup with respect to prare during the sampling phase, one must effectively implement a projector onto the rare event subspace. Without our coherent approach, constructing such a projector would require complex quantum arithmetic circuits. However, the trade-off for this coherent implementation is that our complexity √ √ scales with the product of O(1/ ∆) and O(1/ prare ). In the worst-case scenario where there exists only a single rare event (that is, prare = O(∆)), our complexity converges to that of multidimensional quantum amplitude estimation.
23 Related problems include quantum minimum and maximum finding. The original quantum minimum finding algorithm, described in Ref. [41], is based on Grover’s search algorithm. This approach has been generalized to the two-phase quantum search in Ref. [42]. These methods assume that one can prepare the state X 1 p |x⟩|p(x)⟩, |X L | x∈X L
(E.3)
and one has access to an oracle O such that ( −|p(x)⟩ O|p(x)⟩ = |p(x)⟩
if p(x) ≤ ∆, otherwise,
(E.4)
where |p(x)⟩ corresponds to the binary representation of p(x). To prepare the state, a standard way is to assume another oracle Ox : |x⟩|0⟩ → |x⟩|p(x)⟩,
(E.5)
and implement this oracle onto the uniform superposition state. Constructing these oracles may pose significant implementation challenges. Similar limitations apply to the quantum maximum finding algorithm [43]. While recent work [18] addresses maximum finding using a state preparation input model, their approach relies on restrictive assumptions regarding the largest amplitude and spectral gap (the difference between the largest and second-largest amplitudes). Our approach draws conceptual inspiration from quantum ground state preparation (or low-energy state preparation) [39, 44, 45]. In this analogy, our target output corresponds to a superposition of low-energy eigenstates, defined as those events with probability below the threshold ∆. To utilize this insight, however, one must construct an effective Hamiltonian derived from the given state preparation unitary, which is highly nontrivial. Furthermore, we explicitly prove the optimality of our quantum algorithm and demonstrate a separation between classical and quantum complexities. Our proof technique for these bounds is distinct from the method in Ref. [39]. Appendix F: Phase transition of quantum advantages in power-law tails
In this section, we provide a rigorous discussion to determine in a concrete model family, when the rare-event mass prare stays macroscopic and when it vanishes. This is the key quantity controlling whether the quantum algorithm attains its ideal Θ(∆−1/2 ) scaling, or instead suffers an additional penalty from amplitude amplification. Throughout this section, we fix P = Q, ϵ = Θ(1) and focus only on the asymptotic dependence on ∆. Under this convention, the result as stated in Theorem S2, Theorem S4 simplifies to 1 1 1 SC (∆) = O min log N, 2 + , (F.1) ∆ ∆ prare 1 SQ (∆) = O √ . (F.2) ∆prare When N = poly(1/∆), one has log N = Θ(log(1/∆)), and therefore 1 1 1 SC (∆) = O log + . ∆ ∆ prare
(F.3)
Thus, once N is polynomial in 1/∆, the asymptotic competition between classical and quantum algorithms is determined entirely by the scaling of prare . 1.
Rank-frequency power-law model
We consider a rank-ordered family of distributions over N events, P (xk ) =
k −γ , ZN,γ
(F.4)
24 PN where γ ≥ 0 is the tail exponent and ZN,γ := j=1 j −γ is the generalized harmonic number (also the partition function). The events are ordered so that P (x1 ) ≥ P (x2 ) ≥ · · · ≥ P (xN ). Since the probabilities are monotone decreasing, the rare-event set is always a tail in rank space. For γ > 0, define k∆ := min{k ∈ {1, . . . , N } : P (xk ) ≤ ∆}.
(F.5)
Following this, the rare-event set is R∆ = {xk : k ≥ k∆ }, and the threshold rank satisfies P (xk∆ ) ≤ ∆ < P (xk∆ −1 ).
(F.6)
The total rare-event mass is prare =
N X
P (xk ) =
k=k∆
N X
1 ZN,γ
k −γ .
(F.7)
k=k∆
Because k 7→ k −γ is monotone, integral comparison implies that the sum and integral have the same asymptotic scaling (a standard consequence of the integral test, also known as the Maclaurin–Cauchy test), especially when N is large: ! Z N 1 prare = Θ x−γ dx . (F.8) ZN,γ k∆ To make the role of the state-space size explicit, we parameterize N (∆) = Θ(∆−q ),
q > 0.
(F.9)
We also require that for sufficiently small ∆, P (xN ) ≤ ∆ < P (x1 ).
(F.10)
The left inequality ensures that the rare set is nonempty, while the right inequality ensures that not every state is rare. 2.
Evaluation of the rare-event mass
We now evaluate prare in the three tail regimes γ > 1, γ = 1, and 0 < γ < 1. The case γ = 0 is the uniform case, and one cannot reasonably distinguish which event is rare. Regime 1: γ > 1. In this regime, the generalized harmonic number converges: ZN,γ = ζ(γ) + O(N 1−γ ) = Θ(1),
(F.11)
where ζ(·) is the zeta function. Therefore, the threshold rank obeys k∆ = Θ(∆−1/γ ). The nonemptiness condition P (xN ) ≤ ∆ is equivalent to N −γ ≲ ∆, ζ(γ) + o(1)
(F.12)
which, under N = Θ(∆−q ), requires q ≥ 1/γ (up to constant factors at the boundary). In this regime, Eq. (F.8) yields ! ! Z N 1−γ 1−γ k − N ∆ . (F.13) prare = Θ x−γ dx = Θ γ−1 k∆ Whenever N ≳ k∆ , the lower limit dominates, so 1−γ prare = Θ k∆ = Θ ∆(γ−1)/γ ,
(F.14)
implying that the rare-event mass vanishes polynomially as ∆ → 0. Although the distribution has a long tail in cardinality, that tail carries too little total probability. The rare subspace therefore becomes increasingly difficult to amplify quantumly.
25 Substituting Eq. (F.14) into Eq. (F.3), we obtain SC (∆) = O
1 1 log ∆ ∆
,
(F.15)
because 1
= Θ ∆−(γ−1)/γ = o(∆−1 ).
prare
(F.16)
The quantum complexity is
1 SQ (∆) = Θ √ = Θ ∆−1/2 ∆−(γ−1)/(2γ) ∆prare 1 = Θ ∆−1+ 2γ .
(F.17)
1 < − 12 , the exponent in Eq. (F.17) is more negative than −1/2. In other words, Since γ > 1, we have −1 + 2γ the quantum algorithm is strictly worse than the ideal Θ(∆−1/2 ) scaling because it pays an additional amplification penalty coming from the shrinking mass of the rare subspace. Therefore, the speedup is strictly sub-quadratic. Regime 2: γ = 1. At the critical exponent γ = 1, the normalization diverges logarithmically:
ZN,1 =
N X 1 j=1
= log N + γE + O(1/N ) = Θ(log N ),
j
where γE is the Euler’s constant. The threshold rank is therefore 1 . k∆ = Θ ∆ log N
(F.18)
(F.19)
Using Eq. (F.8), prare = Θ
1 log N
Z N
dx k∆ x
!
=Θ
log(N/k∆ ) log N
.
(F.20)
Now substitute N = Θ(∆−q ). Then 1 + O(1), ∆
(F.21)
1 1 − log q log + O(1). ∆ ∆
(F.22)
log N = q log while Eq. (F.19) implies log k∆ = log Hence log
N = log N − log k∆ k∆ 1 1 = (q − 1) log + log log + O(1). ∆ ∆
(F.23)
Substituting into Eq. (F.20) gives prare = Θ
(q − 1) log(1/∆) + log log(1/∆) q log(1/∆)
This formula reveals three distinct subregimes.
.
(F.24)
26 (a) q > 1. In this case, prare → 1 −
1 ∈ (0, 1). q
(F.25)
Thus, the rare-event mass remains macroscopic. The classical and quantum complexities become 1 1 SC (∆) = O log , ∆ ∆ 1 . SQ (∆) = Θ √ ∆
(F.26) (F.27)
Therefore S (∆) pQ =O SC (∆)
!
1 p
→ 0.
log(1/∆)
(b) q = 1. This is the boundary case N = Θ(1/∆). Eq. (F.24) simplifies to log log(1/∆) prare = Θ . log(1/∆)
(F.28)
(F.29)
Thus, the rare-event mass is not constant, but it vanishes only logarithmically slowly. The quantum complexity becomes s ! log(1/∆) 1 . (F.30) SQ (∆) = Θ √ ∆ log log(1/∆) Meanwhile, the classical complexity remains SC (∆) = O
1 1 log ∆ ∆
.
(F.31)
This is a weaker boundary version of the phase transition. (c) q < 1. In this case, P (xN ) ≫ ∆, so the rare set is empty. This regime is outside the problem setting.
From the calculation, we see that in the regime of γ = 1, the logarithmic divergence of the normalization exactly compensates the decay of the probabilities. This is the critical point at which the tail stops losing polynomial mass. If N grows slightly faster than 1/∆, the rare-event set already carries a constant fraction of the total probability, and the quantum algorithm reaches its ideal Θ(∆−1/2 ) scaling.
Regime 3: 0 < γ < 1. We now turn to the most favorable heavy-tail regime for the quantum algorithm. For 0 < γ < 1, the normalization diverges polynomially with N : ZN,γ =
N X j=1
j −γ =
N 1−γ 1 + o(1) . 1−γ
(F.32)
Therefore, the rank-frequency law takes the asymptotic form P (xk ) =
k −γ = (1 − γ)N −(1−γ) k −γ 1 + o(1) . ZN,γ
(F.33)
Since P (xk ) is monotone decreasing in k, define the threshold rank by k∆ := min{k ∈ {1, . . . , N } : P (xk ) ≤ ∆}.
(F.34)
Whenever k∆ ≤ N , the rare-event set is precisely R∆ = {xk : k ≥ k∆ }.
(F.35)
27 Using Eq. (F.33), the condition P (xk ) ≤ ∆ is equivalent, to leading order, to (1 − γ)N −(1−γ) k −γ ≲ ∆,
(F.36)
which gives k∆ =
1/γ
1−γ ∆N 1−γ
1 + o(1) .
(F.37)
1 + o(1) .
(F.38)
Dividing by N , we obtain the especially useful ratio k∆ = N
1−γ ∆N
1/γ
This expression makes the nondegeneracy conditions transparent. To have a genuine rare-event problem, we must require both: 1. the rare set is non-empty, and 2. the rare set does not coincide with the entire support. The first requirement is that at least one state lies below the threshold. Because P (xk ) decreases with k, this is equivalent to demanding that the least probable state be rare: P (xN ) ≤ ∆.
(F.39)
N −γ 1−γ 1 + o(1) . = ZN,γ N
(F.40)
Using Eq. (F.33) at k = N , we find P (xN ) =
Hence, the rare set is non-empty only if (1 − γ)/N ≲ ∆, equivalently, N ≳ (1 − γ)/∆. If we parameterize the support size as N (∆) = Θ(∆−q ),
(F.41)
then it implies, at the exponent level, q ≥ 1. This lower bound is therefore not an extra assumption; it is exactly the condition ensuring that the rare set is not empty. More precisely: • if q < 1, then P (xN ) ≫ ∆, so no state is rare asymptotically; • if q = 1, writing N ∼ c/∆, non-emptiness requires c ≥ 1 − γ; • if q > 1, then P (xN ) ≪ ∆, so the rare set is automatically non-empty. The second requirement is that not every state be rare. Since P (x1 ) is the largest probability, this means P (x1 ) > ∆. Using Eq. (F.33) at k = 1, we obtain P (x1 ) =
1 ZN,γ
=
1−γ 1 + o(1) . 1−γ N
(F.42)
Thus, to keep the problem nontrivial, we need 1−γ ≳ ∆. N 1−γ
(F.43)
Under the scaling ansatz N (∆) = Θ(∆−q ), this becomes q<
1 1−γ
(F.44)
at the exponent level. Indeed, if q > 1/(1 − γ), then even the most probable state satisfies P (x1 ) ≪ ∆, so all states are rare asymptotically. The boundary case q = 1/(1 − γ) depends on prefactors and is excluded here for simplicity.
28 Combining these two requirements, the genuinely nondegenerate regime is 1≤q<
1 , 1−γ
(F.45)
with the lower boundary q = 1 requiring the prefactor condition N ∼ c/∆ and c > 1 − γ if one wants a strictly positive limiting rare mass. We now evaluate the aggregate rare-event mass. By integral comparison, prare = =
N X
1 ZN,γ
k −γ
k=k∆ 1−γ N 1−γ − k∆ 1−γ
1 ZN,γ
! + o(1).
(F.46)
Using Eq. (F.32), this simplifies to prare = 1 −
k∆ N
1−γ + o(1).
(F.47)
Finally, substituting Eq. (F.38) gives prare = 1 −
1−γ ∆N
1−γ γ
1−γ c
1−γ γ
+ o(1).
(F.48)
∈ (0, 1).
(F.49)
There are two relevant nondegenerate subcases: (a) q = 1. Write N ∼ c/∆ with c > 1 − γ. Then prare → 1 −
Thus, even at the minimal support growth required for non-emptiness, the rare tail already carries a strictly positive constant fraction of the total probability mass, provided c > 1 − γ. (b) 1 < q < 1/(1 − γ). Here ∆N → ∞, so Eq. (F.48) gives prare → 1.
(F.50)
In other words, as soon as the support grows faster than 1/∆ but still remains below the full-degeneracy threshold, almost all of the probability mass lies in the rare tail. Therefore, throughout the entire nondegenerate regime in Eq. (F.45), prare = Ω(1).
(F.51)
This is the decisive structural property of the 0 < γ < 1 regime. Although each individual rare event has a probability at most ∆, the set of all such events retains a macroscopic fraction of the total mass. The tail is therefore heavy not only in cardinality, but also in aggregate probability. As a result, the amplitude-amplification stage of the quantum algorithm does not suffer any asymptotic penalty from a vanishing target subspace. Substituting prare = Ω(1) into the general complexity bounds yields 1 1 SC = O log , (F.52) ∆ ∆ 1 SQ = O √ . (F.53) ∆ Hence the quantum algorithm retains its ideal ∆−1/2 scaling, while the classical method still pays the logarithmic identification cost.
29 TABLE SI. Asymptotic regimes for rank-frequency power laws P (xk ) ∝ k−γ with N (∆) = Θ(∆−q ). Tail regime
Condition on q
γ>1
q ≥ 1/γ
γ=1
q>1
γ=1
q=1
0<γ<1
1 ≤ q < 1/(1 − γ)
Ω(1)
Quantum cost SQ (∆) 1 −1+ 2γ Θ ∆ Θ ∆−1/2 q Θ ∆−1/2 loglog(1/∆) log(1/∆) Θ ∆−1/2
γ=0
–
Degenerate
Degenerate
Rare mass prare Θ ∆(γ−1)/γ 1 − 1q + o(1) log(1/∆) Θ loglog(1/∆)
Advantage Subquadratic Ideal quadratic Nearly quadratic Ideal quadratic Degenerate
Here SQ (∆) suppresses logarithmic factors in log(1/ϵ). In the nondegenerate regimes, the classical cost satisfies SC (∆) = O ∆−1 log(1/∆) and SC (∆) = Ω(∆−1 ). The transition occurs at γ = 1.
3.
Summary of quantum advantages and interpretations
The preceding case analysis shows that the critical quantity is not merely the number of rare states, but the total probability mass carried by those states. The quantum algorithm always pays the thresholding cost Θ(∆−1/2 ), but it −1/2 pays an additional amplification factor prare . Thus: • if prare vanishes polynomially, as it does for γ > 1, the speedup is degraded to sub-quadratic; • if prare = Ω(1), as it does generically for 0 < γ < 1 and for γ = 1 with q > 1, the quantum complexity saturates at Θ(∆−1/2 ); • at the boundary (γ, q) = (1, 1), the rare mass shrinks only logarithmically, quantum algorithm does not achieve the ideal Θ(∆−1/2 ) form. This identifies a sharp phase transition at γ = 1 within the power-law family of Eq. (F.4), and we list all the results in Table SI. Above the critical exponent, the rare tail is too light and prare collapses. At and below the critical exponent, the tail remains macroscopically populated, and the amplitude-amplification penalty disappears. The argument proves a phase transition within the rank-frequency power-law model. We do not claim that γ ≤ 1 is the only possible way for an arbitrary distribution family to satisfy prare = Ω(1) together with N = poly(1/∆). For the Zipfian power-law families analyzed here, the threshold γ = 1 is where the asymptotic behavior changes. Appendix G: Application to stochastic process 1.
Stochastic processes and the Asymptotic Equipartition Property
A discrete stochastic process generates a random output Xt at each time step t, which takes value xt from an alphabet X of finite alphabet size. A consecutive L output sequence xt:t+L := xt , xt+1 , · · · xt+L−1 is governed by the joint probability distribution Pr(xt:t+L ). Stochastic processes can be effectively simulated using recurrent models, which comprise a memory system and an output mechanism [46]. At time step t, given a memory state Si , the model generates an output xt and updates its memory to state Sj according to the joint probability P (Sj , xt |Si ). Following this update, the model recursively applies these transition rules to generate the subsequent output sequence. Definition S3 (Joint entropy). Given a collection of random variables X1 , . . . , XL , with support on Ω1 , . . . , ΩL respectively, we define the joint entropy of the collection of discrete random variables as X X H(X1 , . . . , XL ) = − ··· P (x1 , . . . , xL ) log P (x1 , . . . , xL ). (G.1) x1 ∈X
xL ∈X
For discrete-time stochastic processes, the entropy rate is defined asymptotically. However, note that the entropy rate may not exist for an arbitrary stochastic process. Definition S4 (Entropy rate of stochastic process). For a discrete-valued, discrete-time stochastic process X, the entropy rate is defined as 1 H(X0:L ). L→∞ L
H(X) = lim
(G.2)
30 Here, we focus on the so-called stationary distribution, where any consecutive L outputs are statistically invariant with respect to time, P (x0:L ) = P (xt:t+L )
∀t, L.
(G.3)
For a stationary process, it can be shown that the entropy rate exists and can be quantified via the conditional entropy. Lemma S4. For a stationary stochastic process, the entropy rate exists and is equal to H(X) = lim H(XL−1 |XL−2 , . . . , X0 ). L→∞
(G.4)
The following theorem states that for a well-conditioned stochastic process, almost all sequences we expect to see a probability p(x0:L ) ≈ 2−LH(X) . A stationary stochastic process is ergodic if the time average of the output sequence converges to the ensemble average (or space average), namely, L
1X Xt = E[X]. L→∞ L t=1 lim
(G.5)
Lemma S5 (Asymptotic Equipartition Property (AEP) [47]). Assume X is a stationary and ergodic process, the probability of a sequence x0:L converges to the entropy rate of the stochastic process, that is, lim −
L→∞
1 log P (x0:L ) = H(X). L
(G.6)
This can be understood as an information-theoretic analog of the law of large numbers. We define the typical sequence to have a tighter bound for these sequences. Definition S5 (ϵt -typical sequence). For ϵt > 0, we say a L-length sequence x0:L is ϵt -typical if −
1 log P (x0:L ) − H(X) ≤ ϵt . L
We define typical set T (ϵt , L) as the set of all ϵt -typical L-length sequences, that is, 1 T (ϵt , L) := x0:L : − log P (x0:L ) − H(X) ≤ ϵt . L
(G.7)
(G.8)
The AEP theorem implies the following theorem about the typical set: Lemma S6 (Properties of typical sequence). Let ϵt > 0 be fixed. For any δt > 0, for sufficiently large L, the probability that a sequence is ϵt -typical is at least 1 − δ, namely, P (x0:L ∈ T (ϵt , L)) ≥ 1 − δt ,
(G.9)
(1 − δt )2L(H(X)−ϵt ) ≤ |T (ϵt , L)| ≤ 2L(H(X)+ϵt ) .
(G.10)
and the size of the typical set satisfies
An important observation is that the size of the typical set is exponentially smaller than the whole possible set, unless the stochastic process is nearly uniformly distributed. This can be seen from evaluating the ratio between the size of the typical set and the size of all events. |T (ϵt , L)| 2L(H(X)+ϵt ) 2L(H(X)+ϵt ) ≤ = = 2L(H(X)+ϵt −log |X |) . |X |L |X |L 2L log |X |
(G.11)
A ε-machine is a special type of Hidden Markov Model that consists of several hidden states, which serve as the − := x memory states Si that carry information from the past. The past is defined as a semi-infinite sequence ← x −∞:0 . ← − An epsilon machine maps the past x to one of the hidden states Si . And at each time step, the machine generates an output x and updates its memory state from Si to Sj with probability Pr(x, Sj |Si ). Such a ε-machine is also called a unifilar hidden Markov model whose output xt and the hidden state St determine the next hidden state St+1 uniquely.
31 2.
Thermodynamic mapping for rare events
By definition, rare events are events that are less likely to happen than typical events. The AEP theorem, per Lemma S6, claims that all of the typical sequences happen with probability roughly 2−LH(X) . Definition S6 (Rare event for stochastic process). For a stochastic process X with entropy rate H(X), we define the set of rare events as −KLH(X) RL } := {x0:L |P (x0:L ) ≤ ∆}, K = {x0:L |P (x0:L ) ≤ 2
(G.12)
where K > 1 is some large constant, and ∆ = 2−LKH(X) . To verify whether the event is rare or not, there is a small range of uncertainty. Definition S7 (Unsure event for stochastic process). For a stochastic process X with entropy rate H(X), we define the set of unsure events as Sξ = {x0:L |∆ ≤ Pr(x0:L ) ≤ ∆ + 2ξ},
(G.13)
where ∆ = 2−LKH(X) , and 0 < ξ < ∆. To bound the ratio of rare events is to notice that the entropy rate is similar to the “energy density”, inspired by the concept of Boltzmann weight in statistical physics. By changing the “temperature”, it is possible to construct a mapping, such that the atypical set of the original process is then transformed into the typical set of the new process [48]. Analogously to the energy, we write u = KH(X). We define a set log2 Pr(x0:L ) X =u . (G.14) Λu,L = x0:L − L Lemma S7 (β-mapping [48]). For the stationary stochastic process X, there exists a new stochastic process Xβ such that X
β lim ΛX u,L = lim Λuβ ,L ,
L→∞
L→∞
(G.15)
where uβ = βu + log2 λ, and β, λ > 0. For any u, there exists a β such that uβ = H(Xβ ). With this lemma in hand, we can bound the size of the rare event set Sξ . Theorem S6 (Bounds for the rare event). For sufficiently large L, the probability prare rare event set can be lower bounded by prare = Ω(∆1−µ )
(G.16)
H(X )
β where µ = KH(X) < 1.
Proof. According to Lemma S7, given a u, we can find a β such that X
β lim ΛX u,L = lim Λuβ ,L , uβ = H(Xβ ).
L→∞
L→∞
(G.17)
We can define an alternative set ΛX,ϵ u,L =
log2 Pr(x0:L ) x0:L | − − u| ≤ ϵ . L
(G.18)
For any tolerant error ϵ2 , we can find a tolerant error ϵ1 and sequence length L, such that the typical set Tβ (ϵ1 , L) of process X satisfies 2 Tβ (ϵ1 , L) ⊆ ΛX,ϵ u,L
(G.19)
Let u = KH(X) + η, where η is a very small shift. We choose ϵ2 < η. Then we can obtain L 2 ΛX,ϵ u,L ⊆ RK
(G.20)
32 Then Tβ (ϵ1 , L) ⊆ RL K. By Lemma S6, we can have the following lower bound −L(KH(X)+η+ϵ2 ) 2 P (ΛX,ϵ |Tβ (ϵ1 , L)| u,L ) ≥ 2 −L(KH(X)+η+ϵ2 )
≥2
1−µ
= Ω(∆
)
(1 − δ)2
L(H(Xβ )−ϵ1 )
(G.21) (G.22) (G.23)
H(X )
β < 1. where µ = KH(X) 2 L Since ΛX,ϵ u,L ⊆ RK , we can have the lower bound of the probability of rare event set,
1−µ 2 prare ≥ P (ΛX,ϵ ) u,L ) = Ω(∆
3.
(G.24)
Quantum rare event sampling for stochastic processes
In the following, we detail the applications of our algorithm to stochastic processes, which have additional properties provided by the AEP theorem and β-mapping. This corresponds to a special case of Theorem S4 when Q = P and p1 − p2 = O(1). Given sampling access to a stochastic process, our goal is to sample rare events whose probability is significantly smaller than that of the typical events. This goal is similar to the low-energy state preparation in quantum physics. One may regard a given distribution p P (xi ) as the energy of a physical system’s state xi . Consequently, the physical P system has the Hamiltonian H = P (xi )|xi ⟩⟨xi |. Then, finding the low-energy states of the Hamiltonian, whose i √ energy is less than ∆, is equivalent to finding the rare events of the stochastic process. Corollary S4 (Quantum rare event sampling for stochastic process). Consider a stochastic process X satisfying the asymptotic equipartition property and the β-mapping. Let the rare event set RL K be defined as Definition S6 and prare be the sum of probabilities of all rare events. Let the ambiguity be ξ, and the ratio of the probability of unsure events to rare events be ζ. Given access to a sampler UL for this stochastic process, one can construct a quantum sampler for a distribution that is O(ϵ + ζ)-close to the distribution ( P (xi )/prare , xi ∈ RL K, PR = (G.25) 0, xi ∈ / RL K, by using O(
q
H(Xβ ) ∆µ 1 ξ 2 log ϵ ) times of UL , where µ = KH(X) < 1.
Proof. The algorithm follows mostly the same as Theorem S4. To make sure the ambiguity range is ξ, for quantum we need to set ξ ′ such that √ ( ∆ + 2ξ ′ )2 = ∆ + 2ξ. (G.26) √ √ √ It suffices to take ξ ′ = − ∆+2 ∆+2ξ . Since ξ is small compared to ∆, we further have ξ ′ = Θ(ξ/ ∆). Lemma S2 guarantees that the rectangular function can be approximated to error ϵ by an even polynomial P with degree O( ξ1′ log(1/ϵ)). Lemma S1 claims that such an even polynomial can be applied to the diagonal matrix A by using
O( ξ1′ log(1/ϵ)) queries to controlled-UL and controlled-UL† . √ We can set t = ∆ + ξ ′ , Γ = ξ ′ in Lemma S2, then construct a construct a (1, n + 3, ϵ)-encoding UΠR of a reflection operator by using O( ξ1′ log( 1ϵ )) times of controlled-U and controlled-U † . The unsure set in Lemma S2 is Sξ = {x0:L |∆ ≤ Pr(x0:L ) ≤ ∆ + 2ξ}. Let ΠR denote the projector into all rare events and Πun denote the projector into the unknown area, i.e., X X ΠR = |x0:L ⟩⟨x0:L |, Πun = |x0:L ⟩⟨x0:L |. x0:L ∈RL K
x0:L ∈Sξ
(G.27)
(G.28)
33 Then, we apply the block encoding UΠR to the initial state Xp Xp e P (x0:L )|x0:L ⟩ ≈ |0⟩ΠR P (x0:L ) + |1⟩|⊥⟩ UΠR |0⟩ x0:L
x0:L
= |0⟩
X x0:L ∈RL K
p e P (x0:L )|x0:L ⟩ + |1⟩|⊥⟩,
(G.29)
e is an unnormalized state. The probability of measuring the ancilla qubits as |0⟩ is where |⊥⟩ ∥Π̃R |ϕ⟩∥2 ≥ ∥ΠR |ϕ⟩∥2 (1 − ϵ)2 = prare (1 − ϵ)2 .
(G.30)
By further using the amplitude amplification O( √p1rare ) times, one can prepare the target state √
1 prare
X x0:L ∈RL K
p P (x0:L )|x0:L ⟩.
(G.31)
Combining the QSVT procedure requiring O(1/ξ ′ log 1ϵ ), the overall queries to the controlled-U and controlled-U † is √
log 1ϵ ). Since prare = Ω(∆1−µ ) where µ < 1, the final complexity is O( ξ√p∆ rare O
1 ∆µ/2 log ξ ϵ
Following exactly the same error analysis as in Theorem S4, one can show the total variation distance is O(ϵ + ζ). √ In most cases ξ = O(∆), and we obtain the complexity as O( ∆µ−2 log 1ϵ ). Under the same conditions, following 1 Theorem S2, the classical complexity would be O( ∆ log 1ϵ ), which yields at least a polynomial speedup. 4.
Quantum stochastic modeling
Given a stochastic process governed by distribution Pr(x0:L ), we can always find a unifilar hidden Markov model, called ε-machine. The ε-machine consists of a set of hidden memory states Si with transition probability Pr(Sj , x|Si ). At each time step, the memory state updates to Sj from Si and outputs x with probability Pr(Sj , x|Si ). The ε-machine is unifilar, meaning that the current memory state Si and output x determine the next memory state Sj . From the ε-machine, we can construct a quantum model. The first step is to construct the MPS state q Ax,ij = Pr(Sj , x|Si ). (G.32) Then, we can evaluate its left canonical form Ãx such that X Æx Ãx = I.
(G.33)
x
Ãx forms a set of Kraus operators that can be used to construct the unitary operator U such that X U |i⟩|0⟩ = Ãx |i⟩|x⟩.
(G.34)
x
The quantum model will also map the classical state Si to a quantum state |σi ⟩, which are usually non-orthogonal to each other. The memory states |σi ⟩ can also be computed from MPS states Ax,ij . It is proven that for certain stochastic processes, the memory dimension dq of the quantum model, which is defined to be the dimension of the Hilbert space that the quantum states lie in, can be significantly less than the logarithm of the number of classical states Si [30]. Given an initial state |σi ⟩, we can couple the memory state with L output systems initialized in state |0⟩⊗L , which results in Xp UL |σi ⟩|0⟩⊗L = Pr(x0:L |Si )|σj ⟩|x0:L ⟩, (G.35) x0:L
34 where |σj ⟩ denotes the subsequent memory state determined by the initial state |σi ⟩ and the output sequence x0:L . Meanwhile, the decoupling procedure decouples the memory state from the output systems by applying the reverse unitary Ur Assume the stochastic process has Markov order χ, which means |σi ⟩ is determined by the last χi bits. Then we can construct the reverse unitary Ur which is a χ-qubit control unitary such that X Ur = |x0:χ ⟩⟨x0:χ | ⊗ Uσ†i (G.36) x0:χ
where Uσi |0⟩ = |σi ⟩. Applying Ur will map the memory state back to |0⟩. Therefore, combining UL and Ur will give the desired quantum state ! Xp ⊗L Ur UL |σi ⟩|0⟩ = |0⟩ Pr(x0:L |Si )|x0:L ⟩ . (G.37) x0:L
Appendix H: Numerical results 1.
Setup
We simulate the results of our algorithm using numpy by directly obtaining the matrix representation of the block-encoding unitaries and operating on them. For simulations on small-scale systems, we first construct a (1, (L + χ) log |X | + 2, 0)-block encoding where diagonal encodes the amplitudes of the output quantum state of the recurrent quantum circuit via the protocol as listed in the main paper. Note that the additional qubits stem in the block encoding when compared to Theorem S1 from the fact that the memory qubit is reverted to the zero state, and serves as an ancilla in the block encoding. The amplitudes are then accordingly transformed via our algorithm to produce the required block encoding given an extra ancilla qubit. For larger systems, however, constructing the full (1, (L + χ) log |X | + 2, 0)-block encoding requires an additional polynomial cost for classical statevector simulators due to the much larger statevector from the excess ancilla qubits. We instead classically construct a (1, χ log |X |, 0)-block encoding unitary as follows: √ D I − D2 U (D) = √ (H.1) I − D2 −D where D is a diagonal square matrix that encodes the amplitudes of the output state of the recurrent quantum circuit on its diagonal. While an efficient implementation of such a block encoding on a quantum circuit is an open question, for simulation purposes, it serves the same effect as the diagonal block encoding we construct in our paper, given that the PREP† · SEL · PREP methodology used to obtain the diagonal block encoding in [18] produces a Grover-like reflection operator. This smaller-scale construction can be used to simulate much larger systems with much less time and memory consumption than a full implementation of the block encoding. To obtain the phase angles in the QSVT algorithm that approximates the thresholding function, we use the optimization method introduced in [40] and implemented in pyqsp [49, 50] package. Instead of the constructive methods used to construct angles typically used in proofs for QSVT [34, 51], the phase angles that correspond to the approximation of a polynomial of fixed degree d are found by optimization methods that minimize the difference between the found polynomial and the target function. In particular, the optimizer we utilize finds symmetrical phase angles via an iterative quasi-Newton method [40]. Similar research using optimization methods to find phase angles for QSP and QSVT includes [49, 52]. To apply the angles from QSP to those used in QSVT, one has to first convert the angles from the Wx convention used in QSP to the reflection convention of block encoding construction used in QSVT. We modify the phase angles as follows, where ϕ′k are the angles obtained for the Wx convention, and ϕk are the angles for the reflection convention. π dπ ′ ϕ0 + 4 − 2 k = 0 π ′ ϕk = ϕk + 2 (H.2) 1≤k ≤d−1. ϕ′ + π k=d d 4 Given that the obtained angles from QSP/QSVT implements g ∈ C[x] such that f (x) = Re(g(x)), to implement the block encoding of f (D), we need to implement f (x) = 12 (g(x) + g ∗ (x)) where the phase angles of g ∗ (x) are obtained by negating the phase angles obtained for g(x). While this can be achieved by simply taking the sum of two QSVT circuits by LCU, one can alternatively modify the QSVT circuit such that a Hadamard gate is added to the ancilla
35 ⇥d |0i
H
RZ
RZ
H
n+2
|0i
0
DiagBE
n
0
FIG. S3. Circuit for implementation of a matrix polynomial. 1 − p, 0
1 − p, 1 p, 1
s0
s1 p, 0
FIG. S4. Perturbed coin dynamics. State s0 denotes heads and state s1 denotes tails. At each step, the coin flips with probability p and otherwise remains in the current state.
qubit for the phase-controlled rotation gates before and after the QSVT circuit, as shown in Fig. S3, as the additional input of the |1⟩ state would negate the phase angles, as shown in [50]. For classical simulations, we first generate the original probability distribution, then sample from the distribution. From the samples, we compute the sample distribution and truncate samples whose occurrences in the sample distribution fall above the threshold – thus finding the set of rare events that can be later used for importance/rejection sampling. Note that when one increases classical queries, most queries would produce non-rare events, which would increase the TVD until the increased queries would allow for an extra sample of rare events, and the TVD drops again. This results in a near-periodic fluctuation of the TVD in the classical sampling method that is absent in the quantum method. Our simulation results show that our algorithm can be used to effectively identify rare events, suppress non-rare events, and linearly boost the probabilities of the rare events, given that we can approximate the threshold function to a sufficient closeness.
2.
Perturbed coin model
We apply our example to a simple Markovian process that can be generated by a perturbed coin. At each time step, the coin flips with probability p, otherwise remains at its original state, as shown in Fig. S4. As time evolves, one may
rare event original amplified
Quantum Classical
0.8
0.6
0.6 Error
Probability
0.8
0.4
0.4
0.2 0.2
0.057 0.038 0.019 0.0
0.0
11
10
11
01
11
00
11
11
11
10
10
01
10
00
10
11
10
10
01
01
01
00
01
11
01
10
00
01
00
00
00
00
0
500
1000
1500 2000 2500 3000 3500 Applications of (Quantum) Sampler
4000
4500
5000
Event
FIG. S5. Simulation results for a perturbed coin. The left panel shows the output of the quantum algorithm on a smaller-sized version of the p-coin simulation with L = 4. The new distribution of rare events (blue) amplifies rare events while non-rare events are suppressed compared to the original (red). The right panel shows the convergence of the TVD error for L = 8 as the number of queries increases. The quantum algorithm converges faster and more smoothly compared to its classical counterpart.
36 observe a sequence of heads (denoted as 0) and tails (denoted as 1). The classical perturbed coin contains only two classical memory states: heads (denoted as s0 ) and tails (denoted as s1 ), which carry information from the past. Quantum models use the quantum memory states {|σ0 ⟩, |σ1 ⟩} to carry the information from past rather than the classical states [29, 30]. These quantum states can be expressed as p √ |σ0 ⟩ = 1 − p|0⟩ + p|1⟩, (H.3) p √ |σ1 ⟩ = p|0⟩ + 1 − p|1⟩, |σ0 ⟩ and |σ1 ⟩ are non-orthogonal to each other, which leads to the memory advantage over the classical model. The unitary operator U implements the transition between the quantum memory states. At each time step, the quantum model implements the coupling unitary operator U to couple the memory state |σi ⟩ and the output register, initialized in the vacuum state |0⟩. Since the perturbed coin is a Markov process, the output 0/1 determines the next memory state to be |σ0/1 ⟩. The unitary operator U satisfies A0 = ⟨0|U |0⟩ = |σ0 ⟩⟨0|, A1 = ⟨1|U |0⟩ = |σ1 ⟩⟨1|.
(H.4)
The unitary operator U can be constructed from the Kraus operators Ax . Specifically, we choose the first qubit to be the memory qubit and the second qubit to be the output qubit. The matrix form of U is √ 1−p ∗ 0 ∗ √ ∗ p ∗ 0 U = √ (H.5) p ∗ √ 0 ∗ 0 ∗ 1−p ∗ where ∗ denotes the missing entries. The unitary operator can then be obtained by using the Gram-Schmidt procedure to obtain the missing entries. We encode the probability distribution in the amplitudes of a quantum state and block-encode the distribution as a diagonal block-encoding. We set sequence length L = 8, rareness parameter K = 2.2, and the probability of the coin to be p = 0.1. The history is set to be S0 . We evaluate the TVD distance of the sampling distribution. To compare with classical Monte Carlo methods of simulation, we fix the number of degrees d for polynomial approximation for the √ √ quantum simulations, which would amount to d⌈1/ prare ⌉ queries after ⌈1/ prare ⌉ rounds of amplitude amplification. We disregard the final post-selection of the quantum case and directly take the probability distribution that would be generated from the quantum state if post-selection were 0 after amplification. In the classical case, we match the √ number of queries d⌈1/ prare ⌉ to identify rare events, and further, to compete on fair grounds, we do not perform rejection sampling in the importance sampling phase, but only account for the runtime to identify the rare events. Further implementation details can be found in Sec. H. We plot the TVD distance of our quantum algorithm together with that of the Monte Carlo method in Fig. S5.