Learning the structure of any Hamiltonian
from minimal assumptions

The key technique enabling our learning algorithms is the pseudo-Choi state, a concept introduced by Castaneda and Wiebe [CW23]. They define a version with a reference control qubit, but it is straightforward to consider a referenceless version as well. We will use both variants in our algorithms.

This state draws analogies to the usual Choi state $(\mathcal{E}\otimes\mathbb{I})(|\Omega\rangle\!\langle\Omega|)$ of a quantum channel $\mathcal{E}$. However, since $H$ does not necessarily induce a quantum channel, pseudo-Choi states do not exhibit channel–state duality in the sense of the Choi–Jamiołkowski isomorphism. Nonetheless, it serves as a powerful learning resource: [CW23] shows that Clifford classical shadows [HKP20] can be used to estimate all $m$ parameters of $H$ with $\mathcal{O}(\Delta^{2}\log(m)/\epsilon^{2})$ copies of $|\Psi_{\mathrm{ref}}(H/\Delta)\rangle$. Here, $\Delta\geq 2\|H\|$ is a normalization factor necessary for the preparation of this state. This measurement scheme is therefore able to simultaneously learn a large number of observables with no restriction on their locality. We describe this protocol in more detail in Section 3.

To prepare such states, [CW23] calls upon a subroutine to block encode the Hamiltonian from its time evolution [LC17, GSLW19]. This involves block encoding the operator $\sin(H/\Delta)=\frac{1}{2}(e^{\mathrm{i}H/\Delta}+e^{-\mathrm{i}H/\Delta})$ as a linear combination of unitaries (LCU) [CW12], followed by invoking the QSVT [GSLW19] to apply a polynomial approximation of the $\arcsin$ function onto that block. The normalization $\|H/\Delta\|\leq 1/2$ ensures that the spectrum being transformed lies in the appropriate domain of the $\arcsin$ approximation. This technique is very accurate: the polynomial only needs degree $d=\mathcal{O}(\log(1/\varepsilon_{\mathrm{arcsin}}))$ to get spectral error $\varepsilon_{\mathrm{arcsin}}$, and choosing $\varepsilon_{\mathrm{arcsin}}=\Theta(\epsilon/\Delta)$ suffices for the learning problem. Then from $\mathcal{O}(d)$ queries to the controlled time evolution with resolution $|t|=\Theta(1/\Delta)$, one can construct the block encoding of $H/\Delta$. Finally, the pseudo-Choi state $|\Psi_{\mathrm{ref}}(H/\Delta)\rangle$ is prepared probabilistically from the block encoding by measuring some ancilla flag qubits; this succeeds with probability at least $1/2$, so a constant number of repetitions per copy suffices.

The limitations of this approach are two-fold. First, it does not supply an efficient way to learn the terms if they are a priori unknown. Second, the state-preparation procedure requires controlled forms of the time evolution, $\texttt{ctrl}(e^{-\mathrm{i}Ht})$ and $\texttt{ctrl}(e^{+\mathrm{i}Ht})$. Access to these resources is unnatural in most learning frameworks, and there are no-go conditions for controlling black-box unitaries [AFCB14, GST24]. Furthermore, in the time-forward access model one cannot directly query the (controlled) inverses either. In order to make the pseudo-Choi state method applicable to more conventional learning scenarios, we seek a way to lift these limitations.

To get our main results Theorems 1.5 and 1.6, we show how to prepare pseudo-Choi states without either controlled time evolutions or (in the second access model) time reversal. The guarantees on our state-preparation algorithms are as follows.

Observe that the success probability for $|\Psi_{\mathrm{ref}}(H/\Delta)\rangle$ is always at least $1/2$. Meanwhile for $|\Psi(H)\rangle$, the dependence of $\|\lambda\|^{2}$ will cancel nicely with its use case (which we describe later). The demand for $\mathcal{O}(\epsilon/\Delta)$ accuracy is a consequence of the fact that $|\Psi_{\mathrm{ref}}(H/\Delta)\rangle$ encodes the normalized parameters $\lambda_{a}/\Delta$. In order to get the desired in accuracy in learning $\lambda_{a}$, we need the state-preparation error to be $\leq c\epsilon/\Delta$ for a small constant $c\ll 1$. Below, we describe the techniques that make these constructions possible.

We remove the need for controlled time evolutions through a technique called controlization. While various approaches for controlization have been proposed [Jan02, Zho+11, FDDB14, NSM15, DNSM19, CBW24], we employ the one described by Odake, Kristjánsson, Taranto, and Murao [OKTM23] due to its simplicity and universal applicability. Let $H_{0}=H-\operatorname{tr}(H)\mathbb{I}_{n}/2^{n}$ be the traceless part of $H$. The subroutine generates a mixed-unitary channel $\mathcal{E}(t)$ which approximates $\texttt{ctrl}(e^{-\mathrm{i}H_{0}t})$ to any target diamond-norm distance $\gamma$. In terms of evolution time, controlization has the nice feature that it exhibits zero overhead: to approximate $\texttt{ctrl}(e^{-\mathrm{i}H_{0}t})$, we query $e^{-\mathrm{i}H\delta t}$ for a total of $N=t/\delta t$ times. The downside is that each query requires a small time resolution of $t_{\mathrm{min}}=t/N$, so our main task is bounding what $N$ suffices to achieve the total target error.

The controlization technique of [OKTM23] is based on Pauli twirling and randomized product formulas. It is well know that averaging over the Pauli group gives the identity

This can be extended by appending a control qubit and replacing each Pauli gate by $\texttt{ctrl}_{0}(P)$. In this case the identity becomes

Exponentiating this twirled Hamiltonian gives, up to some global phase, the unitary $\texttt{ctrl}_{1}(e^{-\mathrm{i}H_{0}t})$. We then show that this idea can be easily extended to synthesize multi-controlled unitaries $\texttt{ctrl}_{b}(e^{-\mathrm{i}H_{0}t})$, where $b\in\{0,1\}^{k}$ is any $k$-qubit control string.

In order to evolve by the twirled Hamiltonian, Trotter’s product formula tells us that $e^{-\mathrm{i}\sum_{j}h_{j}}\approx\prod_{j}e^{-\mathrm{i}h_{j}}$. It is enticing to simply concatenate the unitaries $e^{-\mathrm{i}h_{j}}=\texttt{ctrl}_{0}(P^{\dagger})e^{-\mathrm{i}Ht/4^{n}}\texttt{ctrl}_{0}(P)$, but clearly this is inefficient; even ignoring the Trotter error, there are exponentially many terms. Fortunately, Campbell’s qDRIFT protocol [Cam19] was designed with precisely this concern in mind. Rather than include each term in the simulation, qDRIFT constructs random product formulas wherein terms are selected with probability proportional to their norm. In our context, all terms have the same weight, so we just sample from $\mathcal{P}_{n}$ uniformly. By treating this ensemble of random product formulas as a mixed-unitary channel (wherein we draw a new formula for each instance), the analysis of [Cam19] (and tightened in [CHKT21]) shows that qDRIFT sequences of length $N=\mathcal{O}(t^{2}\|H\|^{2}/\gamma)$ suffice to get diamond-norm distance at most $\gamma$. Since we only need applications running for time $t=\mathcal{O}(1/\Delta)=\mathcal{O}(1/\|H\|)$, this simplifies to $N=\mathcal{O}(1/\gamma)$.

To determine what $\gamma$ we need to maintain our learning guarantees, note that our approximation prepares a mixed state from which we probabilistically prepare the pseudo-Choi states. A careful error analysis shows that $\mathcal{O}(\epsilon/\Delta)$ error in trace distance of the mixed states (hence in diamond-norm distance of the channels) suffices. Each block-encoding circuit makes $Q_{1}=\mathcal{O}(\log(\Delta/\epsilon))$ queries to $\texttt{ctrl}(e^{\pm\mathrm{i}H/\Delta})$, so choosing $\gamma=\mathcal{O}(\epsilon/(Q_{1}\Delta))=\widetilde{\mathcal{O}}(\epsilon/\Delta)$ nets the claim of Theorem 1.9. (For the success probabilities, observe that $\Pr(\text{Preparing }|\Psi_{\mathrm{ref}}(H/\Delta)\rangle)\geq 1/2$ always holds because of the reference term, while $\Pr(\text{Preparing }|\Psi(H)\rangle)$ is proportional to the norm of the block-encoded matrix.)

The construction above removes the need for $\texttt{ctrl}(e^{\pm\mathrm{i}Ht})$, but still requires negative-time queries in order to implement the QSVT (and even the LCU for $\sin(H/\Delta)$). To develop a completely time-reversal-free approach, we construct an entirely different quantum circuit for block encoding $H/\Delta$. The basic idea is simple: we appeal to the well-known power series for the (matrix) logarithm [Hal15].

For any square complex matrix $X$ such that $\|X\|<\log 2$, the following identity holds:

Letting $X=-\mathrm{i}H/\Delta$, this formula gives us a relationship between the Hamiltonian and a linear combination of its time evolution operators which only run forward in time. Block encoding the LCU still requires controlled unitaries, but by our prior discussion we now have a technique to approximate them. Note that our normalization requirement is slightly relaxed compared to the QSVT approach, but for sake of clarity let us impose $\Delta\geq 2\|H\|$ as before.³³3We can actually relax the QSVT normalization further: one can allow $\|H/\Delta\|\leq 1-\delta$ for any $\delta\in(0,1)$ by increasing the polynomial degree to $d=\mathcal{O}(\log(\Delta/\epsilon)/\delta)$ [GSLW19, Lemma 70].

Technically, Eq. 7 is not written in the form of an LCU. It is also a series with an infinite number of terms. Our analysis is then to determine where we should truncate the series to maintain sufficient accuracy, and how to express the truncated sum as a proper LCU. The former is straightforward: since $\|X\|\leq\frac{1}{2}$, we have $\|e^{X}-\mathbb{I}\|\leq\frac{1}{2}$. Let $L_{K}$ be the partial sum of Eq. 7 up to $k\leq K$ terms. An application of the triangle inequality gets $\|X-L_{K}\|\leq 2^{-(K+1)}$, which means that $K=\lceil\log_{2}(\frac{1}{2\gamma})\rceil$ ensures a truncation error of at most $\gamma$.

To express $L_{K}$ as an LCU, we expand the terms $(e^{-\mathrm{i}H/\Delta}-\mathbb{I})^{k}$ using the binomial theorem. With the help of combinatorial identities, we can derive analytical expressions for coefficients $c_{j}$ such that $L_{K}=\sum_{j=0}^{K}c_{j}e^{-\mathrm{i}(H/\Delta)j}$. Techniques for implementing LCUs are well established [CW12, GSLW19], requiring controlled instances of the unitaries on only $\lceil\log_{2}(K+1)\rceil=\mathcal{O}(\log\log(1/\gamma))$ control qubits. The bottleneck of any LCU implementation is the size of its coefficient $\ell_{1}$-norm, $\Lambda=\sum_{j}|c_{j}|$. Using the analytical form of the $c_{j}$’s, we can compute an essentially tight bound of $\Lambda=\mathcal{O}(2^{K})$. Recalling that we need to take error below $\gamma=\mathcal{O}(\epsilon/\Delta)$, we immediately get $\Lambda=\mathcal{O}(\Delta^{2}/\epsilon^{2})$ and Theorem 1.10 follows. For learning the Hamiltonian parameters, our LCU now prepares the state $|\Psi_{\mathrm{ref}}(H/(\Lambda\Delta))\rangle$ (still with at least $1/2$ probability), which requires us to refine the learning precision by $1/\Lambda$. In turn, this raises the copy complexity to $\Lambda^{2}\cdot\mathcal{O}(\Delta^{2}\log(m)/\epsilon^{2})=\widetilde{\mathcal{O}}(\Delta^{4}/\epsilon^{4})$. Each copy costs $\widetilde{\mathcal{O}}(1/\Delta)$ evolution time to produce, leading to the claimed result in Theorem 1.6.

So far we have discussed the utility of the pseudo-Choi state with reference, as introduced by [CW23]. Let us now turn to the referenceless version as in Eq. 4. We invoke it in order to accomplish the structure identification task. The reference was essential to learn values of the Hamiltonian coefficients $\lambda_{a}$, but it serves no purpose when we only care about detecting which coefficients are large relative to the others. By performing Bell measurements [Mon17] on copies of $|\Psi(H)\rangle$, we can learn which terms in $H$ have weight greater than $\epsilon$.

Observe that while we need $N=\widetilde{\mathcal{O}}(\|\lambda\|^{2}/\epsilon^{2})$ copies, the probability of producing a single copy is $p=\Omega(\|\lambda\|^{2}/\Delta^{2})$ (time-reversal) or $p=\Omega(\epsilon^{2}\|\lambda\|^{2}/\Delta^{4})$ (time-forward) from $\widetilde{\mathcal{O}}(1/\Delta)$ evolution time. It is a standard fact that $\mathcal{O}(N/p)$ Bernoulli trials are sufficient to get $N$ successes. Hence, $\widetilde{\mathcal{O}}(\Delta/\epsilon^{2})$ (time-reversal) or $\widetilde{\mathcal{O}}(\Delta^{3}/\epsilon^{4})$ (time-forward) total time evolution is used to learn the structure of $H$. Unlike prior approaches to estimating Pauli spectra of unitary operators [MO10, Ang23, Gut24, BLMT24a], our algorithm can in principle be applied to any $n$-qubit operator $A$, as long as one can consistently generate copies of its pseudo-Choi state $|\Psi(A)\rangle$. We are also not restricted by any locality considerations, and the total number of queries to $e^{-\mathrm{i}Ht}$ essentially only depends on $\epsilon/\|H\|$ (i.e., the target error in “natural” units of the Hamiltonian’s norm).

Algorithm 1 describes the core idea of our structure-identifying algorithm. For each copy of $|\Psi(H)\rangle$, we measure in the $2n$-qubit Bell basis, which can be described as $\{(P\otimes\mathbb{I}_{n})|\Omega\rangle:P\in\mathcal{P}_{n}\}$. In turn, this induces a distribution governed by $\Pr(E_{a})=|\lambda_{a}|^{2}/\|\lambda\|^{2}$. Sampling from this distribution allows us to learn the Pauli terms in $H$ that have relatively high weight.

To analyze the cost of this algorithm, we need to bound the number of measurements such that we observe at least all the terms with $|\lambda_{a}|>\epsilon$. This is a variant of the coupon collector problem, which is well studied [Shi07]. Our setup has two modifications from the classic version of the problem. First, we have a nonuniform distribution, which bottlenecks the sample complexity based on the least likely coupon. Second, we only care to collect low-rarity coupons (i.e., those such that $|\lambda_{a}|>\epsilon$). Hence the rarest coupon we need to collect has at least $\Pr(E_{a})>\epsilon^{2}/\|\lambda\|^{2}$. This condition allows us to bound the tail of the cumulative distribution where, after $N$ i.i.d. draws of the coupons (Pauli terms), we still have not observed at least one of each type. We show that this probability is at most $\delta$ whenever $N=\mathcal{O}(\|\lambda\|^{2}\log(m/\delta)/\epsilon^{2})$.

Besides this main argument, there are some smaller technical details to take care of. First is the fact that we do not prepare $|\Psi(H)\rangle$ exactly, but rather an approximation $|\Psi(\widetilde{H})\rangle$. Suppose that $\|H-\widetilde{H}\|\leq\varepsilon$ for some small $\varepsilon$, which we can tune by constructing a more precise block encoding. Then, we show that $S_{\epsilon}(H)\subseteq S_{\epsilon+\varepsilon}(\widetilde{H})$. While we may pick up terms with $|\lambda_{a}|\leq\epsilon$ (this is possible even if sampling from the exact pseudo-Choi state), it is fine to learn a superset. This is because the next stage predicts all $\lambda_{a}$ using classical shadows, which only has an explicit logarithmic cost in the number of parameters. Any parameters deemed too small can simply be discarded. Clearly since we only take a polynomial number of samples, this superset has bounded size. But we can prove something even stronger: $S_{\epsilon+\varepsilon}(\widetilde{H})\subseteq\{E_{1},\ldots,E_{m}\}$. That is, even if the approximation $\widetilde{H}$ has support that leaks onto other Pauli terms, its $(\epsilon+\varepsilon)$-support never contains those erroneous terms.

Finally, we remark that the controlization approximation affects this algorithm more severely than the parameter estimation stage. We find that the Hamiltonian coefficients become perturbed on the order of $\Delta\sqrt{Q_{1}\gamma}$, where $Q_{1}$ is the number of queries to the time evolution for one block encoding and $\gamma$ is the controlization error per query. In order to ensure the correctness of the Bell sampling procedure, we need to adjust the $\epsilon$ threshold of Theorem 1.11 and choose a sufficiently small $\gamma$. The upshot of this analysis is a time resolution of $t_{\mathrm{min}}=\widetilde{\Omega}(\epsilon/\Delta^{3})$ under the time-reversal access model, and $t_{\mathrm{min}}=\widetilde{\Omega}(\epsilon^{4}/\Delta^{5})$ under the conventional access model, for this stage of the algorithm. As is inherent to controlization, $t_{\mathrm{total}}$ and $N_{\mathrm{exp}}$ are unaffected.

Having described an algorithm that learns at the standard quantum limit, we consider how to the enhance the procedure to reach the Heisenberg limit. Using the pseudo-Choi state approach, this is ultimately achieved by a technique called uniform spectral amplification [LC17], a block-encoding generalization of amplitude amplification. As such, we will need access to time evolution inverses, so the following discussion is only viable under the time-reversal model.

To achieve Heisenberg-limited scaling, we use the precision bootstrapping framework, an idea that has enjoyed many fruitful successes in the dynamical learning paradigm [KLY15, HKOT23, DOS23, Zha+23, BLMT24a]. The first use of the strategy for many-body Hamiltonian learning is attributed to Dutkiewicz, O’Brien, and Schuster [DOS23], inspired by the feedback control scheme of Kura and Ueda [KU18] for qudit learning. Its broader applicability was subsequently appreciated by [BLMT24a]. We briefly review the idea here.

Let $\operatorname{\mathcal{A}}:(H,\epsilon,\delta)\mapsto\widehat{\lambda}\in\mathbb{R}^{m}$ be some base learning algorithm with the guarantee that $\Pr(\|\lambda-\widehat{\lambda}\|_{\infty}\leq\epsilon)\geq 1-\delta$. The idea is to bootstrap applications of this base algorithm, with constant precision $\epsilon_{0}=1/2$, to the Heisenberg limit $t_{\mathrm{total}}\sim 1/\epsilon$ by adaptively adjusting the Hamiltonian being learned via quantum control. Consider a previously learned estimate $\widehat{H}$ and form a “residual Hamiltonian” $R=H-\widehat{H}$. Assuming that the current estimate is only $\eta$-good, we can amplify $R$ by a factor of $1/\eta$ to boost the residual errors up to a constant. Then if we learn $R/\eta$ using $\operatorname{\mathcal{A}}$ with merely constant precision $\epsilon_{0}=1/2$, we can improve our estimate to $\eta/2$ error. Recursing on this process for multiple rounds $j=0,1,\ldots,T$ for some $T$, it is clear that each round produces an estimate with error $\eta_{j}=2^{-j}$. Hence $T=\lfloor\log_{2}(1/\epsilon)\rfloor$ suffices to get the target accuracy $\epsilon$. Taking some care with the acceptable amount of failure probability establishes the correctness of the algorithm. This overall strategy is outlined in Algorithm 2.

To achieve Heisenberg-limited scaling, the process of boosting the residual by $1/\eta$ needs to take $\mathcal{O}(\eta)$ time, otherwise the advantage is lost. In the conventional time-evolution setting, this holds because applying $e^{-\mathrm{i}R/\eta}$ takes time $1/\eta$ by definition. In [DOS23], they consider continuous quantum control, wherein the experimenter can implement $e^{-\mathrm{i}(H-\widehat{H})t}$ directly, while [BLMT24a] assumes only discrete control, requiring them to approximate the residual evolution via Trotterization. In our algorithm, we prepare pseudo-Choi states from a block encoding of the Hamiltonian. Thus the bootstrapping strategy requires us to block encode the residual Hamiltonian, which comes with unique advantages and challenges. The main advantage is that we can block encode $-\widehat{H}$ exactly as an LCU [CW12], thereby avoiding the complications of Trotter error in the discrete control model.

On the other hand, boosting the block-encoded residual by $1/\eta$ is more challenging. To do this, we invoke Low and Chuang’s uniform spectral amplification technique [LC17],⁴⁴4A generalization for amplifying the singular values of a not-necessarily-Hermitian matrix [GSLW19, Theorem 30] also accomplishes this task. which essentially allows us to map $R\mapsto R/\eta$, up to some controllable error $\varepsilon_{\mathrm{amp}}$, as long as $\|R/\eta\|\leq 1$. This costs $\mathcal{O}(\log(1/\varepsilon_{\mathrm{amp}})/\eta)$ queries to the block encoding of $R=H-\widehat{H}$. Meanwhile, recall that block encoding $H$ itself costs $\mathcal{O}(\log(1/\varepsilon_{\mathrm{arcsin}}))$ queries to $\texttt{ctrl}(e^{\pm\mathrm{i}H/\Delta})$. A careful error analysis shows that we need to take $\varepsilon_{\mathrm{arcsin}}=\Theta(\eta/\Delta)$ (on the other hand, $\varepsilon_{\mathrm{amp}}=\Theta(\epsilon_{0}/\Delta)$ suffices). Thus in terms of $\eta$, the cost of each round of the bootstrap scales as $\log(1/\eta)/\eta$, giving us nearly Heisenberg-limited scaling up to the factor of $\log(1/\eta)\leq\log(1/\epsilon)$.

Note that this argument applies to the structure identification stage of our algorithm as well. For intuition here, consider the first round $j=0$, wherein we identify all the terms which have $|\lambda_{a}|>1/2$ and estimate them, giving us $\widehat{H}_{1}$. In the next round $j=1$, we synthesize $H-\widehat{H}$ and amplify by $2$. This splits the Hamiltonian terms into two sets: those in $S_{1/2}(H)$, and those in the complement $S_{1/2}(H)^{\mathsf{C}}$. We only care to identify the elements of $S_{1/2}(H)^{\mathsf{C}}$; after amplification, all the such terms larger than $1/4$ become larger than $1/2$, which again can be identified with just a constant $1/2$ precision. Iterating in this way, we can deduce the identity of all terms larger than $\epsilon$ with only $\lfloor\log_{2}(1/\epsilon)\rfloor$ rounds. Combining this with the parameter estimation procedure gives us the claim in Theorem 1.5.

One final detail to point out is the normalization of the residual Hamiltonian. To guarantee that $\|R/\eta\|\leq 1$, we make the particular choice $\Delta=2m$. This is because our error metric is with respect to the vector-norm of the coefficients, which typically has a loose conversion from the operator norm:

That said, there exist Hamiltonians which saturate this inequality, so in the worst case we cannot expect to tighten it. While we use the $\ell_{\infty}$-norm as our metric, converting to different norms simply incurs factors of $m$ in the sample complexity.