Robust and Efficient Hamiltonian Learning

Wenjun Yu
QICI, Department of Computer Science, The University of Hong Kong, Pokfulam Road, Hong Kong, China
Center on Frontiers of Computing Studies, Peking University, Beijing 100871, China

Jinzhao Sun
Clarendon Laboratory, University of Oxford, Oxford OX1 3PU, UK
Quantum Advantage Research, Beijing 100080, China

Zeyao Han
School of Physics, Peking University, Beijing 100871, China

Xiao Yuan
Center on Frontiers of Computing Studies, Peking University, Beijing 100871, China

With the fast development of quantum technology, the sizes of both digital and analog quantum systems increase drastically.
In order to have better control and understanding of the quantum hardware, an important task is to characterize the interaction, i.e., to learn the Hamiltonian, which determines both static and dynamic properties of the system.
Conventional Hamiltonian learning methods either require costly process tomography or adopt impractical assumptions, such as prior information on the Hamiltonian structure and the ground or thermal states of the system.
In this work, we present a robust and efficient Hamiltonian learning method that circumvents these limitations based only on mild assumptions.
The proposed method can efficiently learn any Hamiltonian that is sparse on the Pauli basis using only short-time dynamics and local operations without any information on the Hamiltonian or preparing any eigenstates or thermal states.
The method has a scalable complexity and a vanishing failure probability regarding the qubit number. Meanwhile, it performs robustly given the presence of state preparation and measurement errors and resiliently against a certain amount of circuit and shot noise.

We numerically test the scaling and the estimation accuracy of the method for transverse
field Ising Hamiltonian with random interaction strengths and molecular Hamiltonians, both with varying sizes and manually added noise.
All these results verify the robustness and efficacy of the method, paving the way for a systematic understanding of the dynamics of large quantum systems.

§ INTRODUCTION

Recently, we have witnessed a rapid development of quantum technologies with drastically increasing sizes of fully or partially controllable quantum systems [Arute et al., 2019, Neill et al., 2021, Ebadi et al., 2021, Wu et al., 2021, Gong et al., 2021, Zhong et al., 2020, Mi et al., 2021, Kandala et al., 2017, Kjaergaard et al., 2020, Zhang et al., 2017].

Following the law of quantum mechanics, the quantum devices allow for controllable quantum operations and offer opportunities to naturally probe quantum many-body behaviors in a programmable or analog way.

To further upgrade the quantum technology or even to realize a universal quantum computer, an essential requirement is to better understand the underlying quantum mechanisms of the systems [Georgescu et al., 2014, Altman et al., 2021, Hauke et al., 2012].

According to the Schrödinger equation, the static and dynamical behaviors of quantum systems are described by their Hamiltonian.
Successful detection of the Hamiltonian, i.e., the so-called Hamiltonian learning [Granade et al., 2012, Wiebe et al., 2014, Wiebe et al., 2014, Krastanov et al., 2019, Wang et al., 2017, O'Brien et al., 2021, Haah et al., 2021] task, could guide us to have a better understanding and hence control of analog quantum simulators [Lloyd, 1996, Das and Chakrabarti, 2008, Friedenauer et al., 2008, Aspuru-Guzik and Walther, 2012, Georgescu et al., 2014, Hauke et al., 2012], digital quantum processors [Bernstein and Vazirani, 1997, Shor, 1995, Grover, 1996, Harrow et al., 2009, Krantz et al., 2019, Kjaergaard et al., 2020, Trabesinger, 2012], quantum sensors [Pirandola et al., 2018, Degen et al., 2017, Boss et al., 2017], etc.

A straightforward choice of learning Hamiltonian is via quantum process tomography [Chuang and Nielsen, 1997, Altepeter et al., 2003, Leung, 2003, Merkel et al., 2013, Rahimi-Keshari et al., 2011].

However, since the dimension of a many-body Hamiltonian grows exponentially with the system size, complete tomography of the Hamiltonian is costly and formidable in real experiments [Mohseni et al., 2008, Baldwin et al., 2014].

Here we briefly review some other efficient learning methods and describe their limitations.
Wiebe et al. [Wiebe et al., 2014, Wiebe et al., 2014] firstly use a trustworthy quantum simulator to inverse the unknown evolution and iteratively approximate the exact inverse.
However, since the scaling of complexity is case-specific, the order of magnitude thereof could be prohibitively large for a certain Hamiltonian family.
Besides, the accurate simulation of some arbitrary inverse Hamiltonian evolution remains hard for experimental realizations in NISQ era [Preskill, 2018].
To avoid the usage of expensive and inaccessible quantum resources, some protocols focus on using states that contain information on the Hamiltonia, including steady states and thermal states [Bairey et al., 2019, Evans et al., 2019, Anshu et al., 2021, Qi and Ranard, 2019].
These methods are scalable with the system size by measuring states with desired properties to infer the Hamiltonian.
However, they rely on specific states containing the information of the target Hamiltonian, which are generally hard to prepare and are vulnerable to measurement errors.

There are also learning proposals that get rid of the requirements of preparing some special initial states.
For example, Li et al. [Li et al., 2020] construct equations based on the energy conservation of arbitrary states during evolution.
Given the structure of the Hamiltonian, this work considers the effects of measurement errors and gives a bound on the measurement errors based on the gap between the smallest two singular values of the equation matrix.
Notably, this protocol is steady under both state preparation and measurement (SPAM) errors.
This steadiness, nevertheless, relies heavily on the model assumptions of the measurement errors and Hamiltonian structure (ansatz) errors as well as the eigenstate thermalization hypothesis, and it is hard to give a firm claim for the real-world cases based on these constraints.
Zubida et al. [Zubida et al., 2021] estimate the Hamiltonian parameters by constructing equations from Ehrenfest theorem about the derivatives of observables involved in the structure of the Hamiltonian.
The number of measurements scales with precision quartically $\Omega(1/\epsilon^{4})$ due to the trade-off between statistical and systematic errors, but its accuracy would be affected by state preparation and measurement errors.
Hangleiter et al. [Hangleiter et al., 2021] measure state properties by evolving it to identify the nearest-neighbor coupled Hamiltonian in superconducting systems.
Since the structure and system nature are known, they show that measurements contain the eigenvalues of the Hamiltonian.
By fixing a time evolution as a base, this method can be robust to state preparation errors.
Nevertheless, all of the above three methods need a precise prior about the structure of the Hamiltonian, rendering them unfavorable for tasks of learning an unknown system from scratch or reliably verifying Hamiltonian simulation with extra noise terms.

In this work, we report an efficient Hamiltonian learning algorithm that circumvents the above requirements and assumptions. Our method only utilizes short-time Hamiltonian evolution — getting rid of the requirement of preparing ground or thermal states, assumes sparse Hamiltonians — avoiding the second assumption without specifying the Hamiltonian structure [The sparsity is defined with respect to the decomposition of the Hamiltonian in the Pauli matrix basis. Note that realistic Hamiltonian, including molecules, lattice models, quantum field theories, etc., all satisfy the sparsity requirement.], and exploits ideas from randomized benchmarking (RB) [Flammia and Wallman, 2020] — being robust against effects from SPAM errors.
The method requires only local input states, local Pauli measurements, and local Pauli operations inserted in the Hamiltonian dynamics, which is accessible for various experimental platforms, such as superconducting qubits [Neill et al., 2021, Kjaergaard et al., 2020, Krantz et al., 2019, Mezzacapo et al., 2014, García-Álvarez et al., 2015, Asaad et al., 2016, Weber et al., 2017], trapped ions [Zhang et al., 2017, Häffner et al., 2005, Blatt and Roos, 2012], Rydberg atoms [Saffman, 2016, Ebadi et al., 2021], nuclear magnetic resonance [Cory et al., 1998, O'Brien et al., 2021], etc.

After our paper was posted to the arXiv, there come out various interesting proposals toward efficiently learning unknown Hamiltonians [Franca et al., 2022, Gu et al., 2022, Wilde et al., 2022, Huang et al., 2023].
Wilde et al. [Wilde et al., 2022] recruit the first-order Trotter to simulate Hamiltonian with prior parameters and use the effective fidelity between the simulated and evolved states to generate the loss function.
By imposing constraints on the loss function, they also show the accuracy guarantee of the optimization solution.
Huang et al. [Huang et al., 2023] randomly implement Pauli gates to project the target Hamiltonian onto some stabilizer spaces in the Pauli basis, so the eigenspace of the new Hamiltonian is easily determined.

They achieve the Heisenberg-limited scaling complexity by introducing the phase estimation method as a subroutine and assuming the bounded-degree structure.
We also expect to find more and more interesting proposals for learning an unknown Hamiltonian with fewer assumptions about the Hamiltonian itself.

For the first time, we introduce a provably efficient and robust method to detect an arbitrary Hamiltonian operator containing sparse decomposition parameters over Pauli basis.
In order to be SPAM error robust, our protocol employs the exponential fitting idea from randomized benchmarking to estimate the absolute values of parameters and uses simplified process tomography equations to fulfill the need of sign estimation.
The sign estimation remains robust as it only returns discrete values.
Based on this protocol, we give a provable guarantee in Section <ref> to the estimation results sketched from theoretical bounds of all subroutines.
Given the sparsity $s$ and qubit number $n$, we achieve a measurement complexity as $\order{{sn}/{\epsilon^4}}$ where $\epsilon$ depicts the desired precision.
We simulate the protocol on multiple systems with various Hamiltonian structures in Section <ref> to verify and monitor the scaling properties, including an 8-qubit molecular Hamiltonian with over 100 nonzero Pauli terms.
All these results conclude that our Hamiltonian learning protocol on $n$-qubit systems is scalable and robust. It is also comparably practical since it does not require the correct Hamiltonian structure or the preparation of eigenstates or thermal states and only needs single-qubit operations.

During the estimation of the absolute values (which we denote by stage 1), the protocol estimates the Pauli fidelity of the Hamiltonian evolution channel and transforms these values to be Pauli error rates, which are equal to squares of the decomposition parameters.
We construct the fidelity estimator by implementing the random circuits [Flammia and Wallman, 2020] in an interleaved manner to estimate the extraneous Hamiltonian channel.
We need to assume that the circuit noise performs benignly as in Assumption <ref> and <ref>, and the exponential decay and fittings achieve SPAM-robust estimation.
Suppose Hamiltonian operators to be sparse as in A1 and A2, we adopt the sparse transforming idea [Harper et al., 2021] to reduce the exponential calculation to an efficient sampling.
We also relieve the original assumption about the unbiased Gaussian noise to biased noise.
In the second stage, we introduce an optimization problem inspired by the process tomography protocol [Shabani et al., 2011] to estimate the sign information.
With the assumption about the SPAM errors in Assumption <ref>, the sign estimator can be SPAM robust as it tolerates small disturbances.

The remainder of this paper is organized as follows.
We provide the basic problem setting and intuition in <Ref>.
We present an overview and the main ideas of our protocol in <Ref> with a rigorous analysis of the validity of the protocol in <Ref>.
The numerical simulations for the Hamiltonian learning protocol are exhibited in <Ref>.
We conclude the results in <Ref>.

§ BACKGROUND

In this section, we introduce the main setting of the Hamiltonian learning problem and some basics thereof.

We then summarize the assumptions about the Hamiltonian in our protocol.
We also introduce the relationship between the time evolution channel and the Hamiltonian operator as a tool to extract the estimation.
In order to better control the SPAM errors and circuit noise, we propose some restrictions to keep them performing benignly.
The assumptions and constraints help to learn the information of the Hamiltonian operator robustly and efficiently.

§.§ Problem Setting

The Hamiltonian learning problem aims to estimate the unknown underlying Hamiltonian operator of a quantum system.

Consider an $n$-qubit Hermitian Hamiltonian $H$ with dimension $d\times d$ and $d=2^n$,

which can be decomposed on the Pauli operator basis as

\begin{gather}\label{eq:decomp}
H=\sum_{P_\alpha \in{\sf P}^n}s_{\alpha}P_{\alpha},
\end{gather}

with real coefficients $s_\alpha$.
Here the ${\sf P}^n=\{I,X,Y,Z\}^{\otimes n}$ is the refined Pauli group consisting of $4^n$ Pauli matrices.
We label each element therein by a $2n$-bit string as $\alpha$ for $P_\alpha$, and Appendix <ref> includes the detailed definitions of this labeling and other Pauli-related concepts.
We denote the vector ${\bf s}\in\mathbb{R}^{d^2}$ by the decomposition parameter, and each element can be obtained by $s_\alpha=\frac{1}{d}\Tr(HP_\alpha)$ since ${\sf P}^n$ forms a basis under the Hilbert Schmidt norm.
Consequently, the Hamiltonian learning problem is treated as the reconstruction of the decomposition parameters ${\bf s}$.

In the literature, several Hamiltonian learning methods have been proposed to efficiently learn the decomposition parameters ${\bf s}$ [Bairey et al., 2019, Hou et al., 2020, Garrison and Grover, 2018, Anshu et al., 2021, Qi and Ranard, 2019, Hangleiter et al., 2021] by imposing

* a certain Hamiltonian structure, such as nearest-neighbor or local interactions

* the ability to prepare an eigenstate or thermal state of the Hamiltonian.

Yet, these constraints may fail in practical scenarios.

For the first requirement, a realistic system may have different types of interactions, and it would be impractical to assume a specific Hamiltonian structure. For example, we may not have full prior information on the Hamiltonian structure for some partially controllable systems [Altman et al., 2021]. Meanwhile, even if we know the Hamiltonian of an ideal quantum system, the imperfection of the system may lead to unknown and complex interactions [Hauke et al., 2012, Altman et al., 2021, Zhang et al., 2017], such as the multiqubit crosstalk of superconducting qubits [Krantz et al., 2019, Kjaergaard et al., 2020, Takita et al., 2017, Kandala et al., 2019, Sun et al., 2021].
Indeed, the dimension of the vector $\bf s$ increases exponentially with the number of qubits, and a specific structure of the Hamiltonian will restrict the problem size to an exponentially smaller subspace. In this work, we release the assumption to avoid prior information on the Hamiltonian. We notice that while the Hamiltonian structure may be unknown, most real-world physical systems, such as superconducting qubits, cold atoms, etc., generally have sparse decomposition parameters ${\bf s}$.
For convenience, we denote the set of nonzero parameters by ${\bf s}^\star$.
In this work, our protocol and analysis mainly focus on the Hamiltonian with sparse nonzero parameters, and we make the following assumptions to determine the decomposition parameters ${\bf s}$ of the target Hamiltonian.

Given Hamiltonian $H$, we have

$\bf A1$ (Random sparse support) The support of the nonzero decomposition parameters ${\bf s}^\star$ is chosen uniformly at random from all subsets with size $s$ of the group ${\sf P}^n$. Here $s$ is polynomial with the qubit number, $n$.

$\bf A2$ (Spectral gap) The smallest nonzero parameter of ${\bf s}^\star$ is larger than a constant $\sqrt{\epsilon_0}$.

$\bf A1$ ensures the sparsity of the Pauli decomposition of the target Hamiltonian, which is much more general compared to the structure assumptions, i.e., indicating which Pauli term has a nonzero coefficient.
Although this assumption supposes uniformly random support, this is due to the concern of theoretical analysis and proofs.
In the numerical results, we will show the feasibility of our protocol when learning sparse Hamiltonians with different structures, including a comparably complicated molecular Hamiltonian.

The assumed lower bound of the spectral gap in $\bf A2$ helps to distinguish the nonzero parameters from overwhelming noise effects.
These two assumptions help to narrow the estimation to the sparse parameters ${\bf s}^\star$, which makes the execution efficient regarding the number of qubits.

For the second requirement, preparing the eigenstate or thermal state of an arbitrary Hamiltonian is proven to be QMA-hard [Kempe et al., 2006], and this eigen- or thermal state preparation also seems impractical. Instead of solving those computationally hard problems, a more natural way is to consider the time evolution of the Hamiltonian, which is accessible for both analog and digital quantum computers [Georgescu et al., 2014, Trabesinger, 2012, Childs et al., 2018, Low and Chuang, 2019, Sun et al., 2021]. Using the Schrödinger equation,

the time-evolved state for an initial state $\rho$ is represented as a unitary

\begin{gather}
\mathcal{H}_t(\rho)=U_H(t)\rho U_H(t)^\dagger,
\end{gather}

where $U_H(t)=e^{-iHt}$ corresponds to Hamiltonian $H$ and time $t$.
Here and henceforth, we use $\mathcal{H}_t$ to denote the Hamiltonian evolution channel with time $t$.
For a short enough time, we can approximate the Hamiltonian evolution channel by expanding and truncating the infinite series as

\begin{align}\label{eq:Hamiltonian}
=\rho+it\sum_{\alpha\in{\sf P}^n}s_\alpha(\rho P_\alpha-P_\alpha\rho)\notag\\
+&t^2\sum_{\alpha,\beta\in{\sf P}^n}s_\alpha s_\beta\left[ P_\alpha\rho P_\beta-\frac{1}{2}(P_\alpha P_\beta\rho+\rho P_\alpha P_\beta)\right]+o(t^3).
\end{align}

For simplicity, here and after in this paper we will use the index $\alpha\in{\sf P}^n$ to indicate the corresponding operator $P_\alpha\in{\sf P}^n$ when it would not incur ambiguity.
We note that the first- and second-order terms of the evolution channel faithfully store the decomposition parameters of the Hamiltonian operator, which offers the possibility for extracting the Hamiltonian information from the dynamics.
Together, we will utilize the short-time evolution under the target unknown Hamiltonian and present a Hamiltonian learning protocol for any unknown Hamiltonian with sparse decomposition parameters.

§.§ Noise Models

In way to construct a robust Hamiltonian learning protocol, a crucial obstacle remaining is the unknown SPAM errors and circuit noise.
Due to the limitation of the quality and controlling quantum systems, these errors and noise are ubiquitous in NISQ devices.
The imperfect state preparation and readout noise will undermine the availability of measurement outcomes, which renders the estimation untrustworthy.
To make the protocol SPAM robust, we first need  to restrict the strength of these errors; otherwise, the unbounded errors will always be able to destroy all the effective information from measurements.
For theoretical convenience, we will hold the following assumption to quantify the device's SPAM errors.

We denote the noisy state by $\tilde{\rho}$ and the noisy measurement operator by $\tilde{M}$.
The noisy device satisfies the following constraints for every pair of positive Pauli eigenstates and Pauli measurements,

\begin{gather}
\abs{\Tr(\tilde{M}\mathcal{H}(\tilde{\rho}))-\Tr(M\mathcal{H}(\rho))}\leq\tau,
\end{gather}

where $\mathcal{H}$ here represents the Hamiltonian evolution with all possible short evolving $t$ used in the implementation.
The positive bound $\tau$ must be in the region of $\left[0,O\left(t_1\sqrt{\frac{\epsilon_0}{s}}\right)\right]$ where $t_1$ is a typical unit time used in the regression during the sign estimation.

In order to eliminate these depicted SPAM effects, we will show shortly that the protocol would choose special circuit components, especially random Pauli gates, and implement them in a interleaved way to build a SPAM-robust estimation.
In this sense, we need to clarify the form and amplitude of the introduced gate noise in a rigorous way.
Following the idea and convention in randomized benchmarking protocols, [Flammia and Wallman, 2020, Helsen et al., 2022, Helsen et al., 2019] we impose time-independent and Markovian quantum noise.
Furthermore, for every noisy implementation of the Pauli gate, we assume that there exists some certain noise channel $\Lambda$ such that

\begin{gather}\label{eq:independent}
\tilde{\mathcal{P}}=\Lambda\circ\mathcal{P},\ \forall\,P\in {\sf P}^n,
\end{gather}

where $\mathcal{P}$ is the implementation of the unitary $P$.
We denote the noise effect satisfies Eq. (<ref>) by gate-independent noise.
Therefore, we have the following model for circuit noise.

The noise from the implementation of quantum circuits is time-independent, Markovian, and gate-independent.

Even though only gate-independent noise allows the RB circuit to construct the standard twirling structure, Wallman [Wallman, 2018] proved that a slight gate-dependent disturbance of the gate-independent noise still keeps the twirling results and causes limited harm to the final analysis of RB.
This might give us a clue about the possibility to relax this gate-independent requirement.

An additional assumption is ad-hoc for the execution of the interleaved protocol.
Since the implementation of random Pauli gates carries some noise inevitably, our protocol needs to calibrate and compensate for these effects to extract the net information of the Hamiltonian.
However, in order to eliminate them, the noise is supposed to be diagonal in the superoperator representation, i.e., a stochastic channel.

The noise $\Lambda$ of the circuit implementation for Pauli gates is a stochastic channel (Pauli channel), and it is close to the identity according to the operator norm as

\begin{gather}
\|\Lambda-\mathcal{I}\|\leq\frac{1}{3}.
\end{gather}

Note that Pauli gates are actually implemented from single-qubit pulses in realization.
From some current experimental verification and detection results [Harper et al., 2020, Barends et al., 2014, Somoroff et al., 2021], the implementation of local operations can be close to their ideal version up to high precision.
Therefore, Assumption <ref> is reasonable.
In fact, the constant $\frac{1}{3}$ can be further enlarged to an arbitrary constant less than 1, which still makes sure that all the Pauli fidelity terms of $\Lambda$ are larger than 0.

An overview of the Hamiltonian learning execution. The first stage can be viewed as a Pauli error rate estimator along with an independent Pauli fidelity estimator. In circuit (a), we introduce the kernel of the random circuits to get the desired fidelity terms, where $\tilde{P}$ denotes the noisy Pauli gate. By choosing the random matrix $M$, the index $\ell$, and the offset $d$, the protocol accordingly queries the fidelity terms $f_{M'\ell+d}$, where $M'$ is the shuffled $M$ as discussed in Section <ref>.
According to the query, the procedure then prepares input states followed by the cascading structures including the noisy Pauli gates and Hamiltonian evolution.
The sub-subscripts represent the corresponding single-qubit parts of Pauli gates.
By measuring circuits with different depths, it acquires the fidelity terms with bounded noise and combines them to be the bin sets as in Eq. (<ref>).
The formal illustration of the circuits can be viewed in Section <ref>.
In subfigure (b), the procedure uses multiple subsampling groups and bin sets aligned therein to detect the type of each bin set and peel to generate bins with single nonzero Pauli error rates. For example, in this two-qubit case, the procedure finds the first bin set in the first group to be single-ton and extracts the nonzero term and index, $p_{XX}$, with the aid of different signs. For the second bin set, the procedure gets stuck due to the summation of multiple terms. It then switches to the second group to detect the first bin set with result $p_{IZ}$. By learning $p_{IZ}$ is nonzero, the procedure uses it to peel that stuck case and gets $p_{YX}$.
Then the procedure uses estimated error rates to peel the second bin set of group two and certifies it as a zero-ton bin.
This peeling can be viewed more detailedly in Section <ref>.
The second stage aims to complete the Hamiltonian learning by estimating signs of decomposition parameters as illustrated in Section <ref>. In circuit (c), the procedure inputs random local Pauli eigenstates and measures them with Pauli operators after evolving under Hamiltonian, where $\tilde{C}$ and $\tilde{C'}$ denote the noisy local Clifford gates for state preparation and measurements which are not conjugate. By constructing the relationship between the measurement outcomes and the unknown parameters, we can solve them and recover the sign information.

§ PROTOCOL OVERVIEW

In this section, we illustrate the main ingredients of our protocol to estimate the unknown Hamiltonian operator.
The basic idea is to extract the information of the sparse parameter vector $\bf s^\star$ sequentially from the second- and first-order terms of the series expansion of Hamiltonian evolution in Eq. (<ref>).
In this sense, our protocol is naturally analog to a quantum channel detection scheme.

We briefly summarize our protocol as follows. We also refer to Figure <ref> for a schematic summary.

* At the first stage, the procedure first  extracts the Pauli fidelity with Pauli inputs and measurements in a similar vein to RB as in Figure <ref> (a), which is robust against SPAM errors.
By linear regression with $t^2$, the results correspond to the second-order terms in the evolution channel.
Then we show a procedure to efficiently transform the extracted fidelity to the squares of $\bf s^\star$, and an example is illustrated in Figure <ref> (b).

* At the second stage, we construct and solve the process equations to estimate the remaining sign information of $\bf s^\star$ from the first-order terms as in Figure <ref> (c).
With knowledge of the support position about $\bf s^\star$ from the first stage and the nature that small noise on results will not cause any sign flips, the sign estimator can be executed robustly and efficiently for both quantum and classical complexity.

One may argue that a straightforward strategy is to directly consider the first-order expansion of Eq. (<ref>) and extract $\bf s^\star$ with certain input states and measurements.
However, the first-order term corresponds to coherent (off-diagonal) information of the channel, which is generally sensitive to SPAM errors. Indeed, most quantum channel detection schemes that are robust to SPAM errors, such as randomized benchmarking [Magesan, 2012, Dankert et al., 2009, Erhard et al., 2019, Helsen et al., 2019], can only extract incoherent (diagonal) information of channels.
Moreover, as in Eq. (<ref>), the overhead  for  classical computing to construct the process equations from state preparations and measurements would be exponentially large without the information of the sparse support of the Hamiltonian effect since there are exponential  items on the Pauli basis.
Therefore, we choose this composite algorithm to circumvent these obstacles.

We discuss in detail the two stages in the remainder of this section and refer to the next section and appendix for an informal summary and detailed proofs of the protocol, respectively.

§.§ Stage 1. Absolute value estimation

The first stage is to estimate the squares (absolute values) of the decomposition parameters, which come from the second-order effects of the Hamiltonian evolution.

Considering the $\chi$-representation using the Pauli basis, the Hamiltonian evolution channel can be decomposed as

\begin{align}\label{eq:processmatrix}
\mathcal{H}_t(\rho)=\sum_{\alpha,\beta\in{\sf P}^n}\chi_{\alpha\beta}P_\alpha\rho P_\beta,
\end{align}

where $p_\alpha\coloneqq\chi_{\alpha\alpha}$ is denoted by the Pauli error rate.
Comparing Eq. (<ref>) with Eq. (<ref>), the relationship between the parameters in $\bf s$ and Pauli error rates $\{p_\alpha\}_{\alpha\in{\sf P}^n}$ is

\begin{gather}\label{eq:chi}
\begin{cases}
s_\alpha^2 t^2+o(t^3)& \text{$\alpha\neq 0$}\\
1-t^2\sum_{\beta\in{\sf P}^n,\beta\neq0}s_\beta^2+o(t^3)& \text{$\alpha=0$}
\end{cases}.
\end{gather}

Since the square of decomposition parameter $s_\alpha^2$ is proportional to the Pauli error rate $p_\alpha$ with a small remainder from small $t$, we can approximately obtain every $|s_\alpha|$ given the corresponding Pauli error rate $p_\alpha$.
Note that the sparsity of $|s_\alpha|$ also implies the sparsity of $p_\alpha$.
Therefore, the basic idea is to exploit the recently proposed Pauli estimation protocol to learn all the sparse Pauli error rates [Harper et al., 2021] from Pauli fidelity (which is defined later).
Nevertheless, there are two challenges if we directly adopt existing Pauli estimation protocols.

* First, the Hamiltonian channel also contains coherent parts, while we only need the incoherent information here.

* Second, Pauli error rates of a short-time Hamiltonian channel are only approximately sparse to the second-order expansion. It is actually dense if we consider infinite expansions.

Fortunately, we can circumvent these in our estimation.
For the first one, our circuit would apply the Pauli twirling to the Hamiltonian evolution channel to convert it into a Pauli (diagonal) channel.
In brief, Pauli twirling annihilates the off-diagonal elements of a channel and keeps all the Pauli (diagonal) information the same.
Denote a Pauli gate by $\mathcal{P}_\alpha$, and the evolution channel $\mathcal{H}_t$ after Pauli twirling must be

\begin{gather}\label{eq:PauliTwirling}
\mathcal{H}_t^{{\sf P}^n}(\rho)\coloneqq\frac{1}{\abs{{\sf P}^n}}\sum_{\alpha\in{\sf P}^n}\mathcal{P}_\alpha^\dagger\circ\mathcal{H}_t\circ\mathcal{P}_\alpha(\rho)=\sum_{\alpha\in{\sf P}^n}p_\alpha P_\alpha \rho P_\alpha.
\end{gather}

For the second one, we apply regression to the results with different small $t$ to suppress the effects induced by higher-order terms in the evolution.

We further prove that the protocol is noise-resilient and efficient for large-scale Hamiltonian evolution channels.

An illustrative diagram of the gate-cascading estimator of the Pauli fidelity of the Hamiltonian evolution channel. This estimator performs robustly against both SPAM errors and the implementation noise from Pauli gates. (a) In order to rule out the effects of the implementation noise $\Lambda$, we recruit the circuit in [Flammia and Wallman, 2020] as a reference circuit. The procedure inputs the stabilizer states of stabilizer group $G$ in the Pauli group and measures the Pauli syndromes using local Clifford gates $C$ and $C^\dagger$.
Random Pauli gates $P$ twirl noise channels into the diagonal version with only fidelity left.
The sub-subscripts of the locally tensored gates represent the corresponding single-qubit parts.
Therefore, the procedure varies the gate length and gets a decaying factor by tracking the overall syndromes from both Pauli gates and the syndrome measurements, which helps to estimate the fidelity information of $\Lambda$. (b) Since we have access to the short-time evolution of the target Hamiltonian, we implement it in the estimation circuit. Similarly to the previous circuit, we estimate the fidelity of the composite channel $\mathcal{H}_t\circ\Lambda$. Under Assumption <ref>, a ratio estimator contributes to the extraction of the net fidelity information of the Hamiltonian channel from the composite results from both circuits.

§.§.§ Pauli fidelity estimation

In this subsection, we focus on getting data from measuring quantum circuits consisting of the Hamiltonian evolution.
Particularly, the protocol  extracts Pauli fidelity as a type of Pauli information for further analysis and extraction of the target parameters.

Inputting and measuring the same operator $P_x$, the outcome corresponds to a Pauli fidelity of an arbitrary channel $\mathcal{E}$ as defined in the following

\begin{equation}\label{eq:Paulifidelity}
f_x \coloneqq \frac{1}{2^n}\tr[P_x \mathcal{E}(P_x)] = \sum_{\alpha\in{\sf P}^n}(-1)^{\langle x,\alpha\rangle_p}p_\alpha,\ \forall\, x\in{\sf P}^n,
\end{equation}

where $\langle x,\alpha\rangle_p\in\{0,1\}$ is a useful operation named Pauli inner product.
Depending on whether $P_\alpha$ and $P_x$ commute or not, $\langle x,\alpha\rangle_p$ returns 0 or 1.
We give the mathematical definition of this inner product using bit strings of indices in Appendix <ref>.
Eq. (<ref>) and Eq. (<ref>) depict the major tool and the main target of the estimation in this step, respectively.

The fidelity estimator contains two circuits to extract the Pauli fidelity of the target Hamiltonian channel, as shown in Figure <ref>. The idea of the circuit (a), as well as the Pauli estimation, is proposed in [Flammia and Wallman, 2020].
We introduce an additional circuit (b) to rule out the noise effects.
Therefore, each round of execution consists of implementing these two circuits.

As shown in circuit (a) of Figure <ref>, we prepare a stabilizer state using a local Clifford gate $C$ and $\ket{0}^{\otimes n}$, where the state is the simultaneous positive eigenstate of all Pauli operators in a stabilizer group $G$ that includes $P_x$ given the target fidelity is $f_x$.
Here the group can be chosen arbitrarily given it contains $P_x$.
The concepts related to the stabilizer group in the Pauli group are reviewed in Appendix <ref>.
We then implement random Pauli gates, and
these gates twirl the noise channel $\Lambda$ into a Pauli channel $\Lambda^{\sf{P}^n}$ which is a diagonal matrix in the superoperator representation with elements to be all Pauli fidelity.
In this case, the circuit is equal to the combination of diagonal channel from the twirled noise channel.
We finally measure the syndrome $\beta\in A_G=\faktor{{\sf P}^n}{G}$ of the state, i.e., do $C^\dagger$ to the state and measure it on $Z^{\otimes n}$ basis.
A brief introduction of the syndrome measurement is also in Appendix <ref>.

The estimator sums all the indices of the sampled Pauli gates as $a$ and decides the corresponding syndrome or the decomposition of $a$ in the quotient group: $\alpha\in A_G$.
It is found that the distribution probability of the overall error syndrome $\alpha+\beta$ is the linear combination of exponentially decaying fidelity terms with signs depending on the inner products of fidelity indices and the overall syndrome.
We can further isolate every single decaying by adding an additional sign for each overall syndrome and taking an expectation over these syndromes as shown in Lemma 6 of [Flammia and Wallman, 2020].
By varying $m$, we can extract the desired fidelity term $f_x$ of $\Lambda$ robustly as SPAM errors do not affect the decay rates.
An advantage of considering the stabilizer formalism is that we can simultaneously estimate all the fidelity in one stabilizer group within one round of experiments by different post-processing.
While it would be hidden in the asymptotic complexity, this advantage can expedite real-world processing considerably.

In circuit (b) of Figure <ref>, we interleave Pauli gates and the additional time evolution and keep the state preparation and measurements the same. The target channel now becomes $\mathcal{H}_t\circ\Lambda$ with the same analysis as the previous circuit. We can, therefore, get the Pauli fidelity $f_{\mathcal{H}_t\circ\Lambda}$.
In the remainder of this paper, we will omit the subscript when it is clear of which channel the fidelity is.

With the executions of both circuits, we can eliminate the Pauli gate noise $\Lambda$ and extract the net information of the Pauli fidelity of $\mathcal{H}_t$ solely.
As in assumption <ref>, we assume that the implementation noise $\Lambda$ is itself a Pauli channel, and a simple ratio estimator helps to detect the corresponding Pauli fidelity of the Hamiltonian evolution $\mathcal{H}_t$.
Therefore, we get a noise-robust estimation of the Pauli fidelities of the Hamiltonian evolution channel $\mathcal{H}_t$.

Observing Eq. (<ref>) and Eq. (<ref>), we realize that fidelity terms of the evolution channel also contain higher-order effects of the evolving time.
In order to extract the net second-order coefficients that directly contain the information of decomposition parameters, our procedure applies an ordinary linear regression for multiple pairs of Pauli fidelity of evolution channels and squares of evolving times.
This helps to estimate the second-order Pauli fidelity, $\{f_x^{(2)}\}$, of the Hamiltonian evolution channel, which is the second-order coefficient of the ordinary fidelity of the evolution,

\begin{align}\label{eq:fidelityexpan}
f_{\mathcal{H}_t,x}=\frac{1}{2^n}\tr[P_x \mathcal{H}_t^{{\sf P}^n}(P_x)]=f_x^{(0)}+f_x^{(2)}t^2+o(t^2).
\end{align}

This fidelity denotes the second-order coefficient of the ordinary Pauli fidelity, and we have the following definition from Eq. (<ref>) and Eq. (<ref>)

\begin{align}
f_x^{(2)}\coloneqq\sum_{\alpha\in{\sf P}^n}s^2_\alpha[(-1)^{\langle\alpha,x\rangle}-1],
\end{align}

which is time-independent.

We summarize the whole fidelity estimator in Algorithm <ref>.

Fidelity Estimator: FE($\sf G,X$, $l,t_0$)
[1]
Input: a target set $\sf X$ belonging to a stabilizer group $\sf G$, a positive integer $l$, and the unit time length $t_0$
Initialize the identity array $I$ with size $|\sf{X}|$ to be all 1
$\hat{f}_\Lambda\leftarrow I-\textbf{FEstimator}({\sf G,X},l)$
$t\leftarrow i\cdot t_0$
$\hat{f}_{\mathcal{H}_t\circ\Lambda,i}\leftarrow I-\textbf{HFEstimator}({\sf G,X},l,t)$
$\hat{f}_{\mathcal{H}_t,i}\leftarrow \hat{f}_{\mathcal{H}_t\circ\Lambda,i}/\hat{f}_\Lambda$ratio estimator
Do ordinary linear regression for every fidelity in $\sf X$ from $(\hat{f}_{\mathcal{H}_{it_0}},(i\cdot t_0)^2)$ pairs and get $\hat{f}^{(2)}$
the second-order fidelity array $\hat{f}^{(2)}$

The subroutines, FEstimator and HFEstimator, denote the implementation of circuit $(a)$ and $(b)$ in Figure <ref> with post-processing, respectively.
The algorithm uses the former estimation as a reference and implements the latter for multiple evolving times to both eliminate the Pauli noise channel and extract the second-order terms.
For instance, here, we choose five evenly-increasing evolving times for Hamiltonian evolutions.
Dividing the composite fidelity by the reference fidelity of noise channel $\Lambda$, we can estimate the fidelity of the Hamiltonian channel with the corresponding evolving time, as shown in Lemma <ref>.
We denote this as the ratio estimator.
We will elaborate on these subroutines in detail in Appendix <ref>.

For this second-order fidelity estimation, we hold a guarantee about the accuracy and complexity of the results.
In this theoretical statement, we assign the ordinary least square as our regression method for convenience.

Let $\sf G$ be a stabilizer group in ${\sf P}^n$ and ${\sf X}\subset{\sf G}$,
and suppose Assumption <ref> and <ref> hold.
For any small $\epsilon,\delta>0$, the following happens with probability $1-\delta$:
Run Algorithm <ref> $\textbf{FE}({\sf X,G},l,t_0)$ with $t_0$ satisfies $\|\mathcal{H}_{5t_0}-\mathcal{I}\|<\frac{1}{4}$ and  $l = \frac{2}{\epsilon^2} \log\bigl(\frac{2|\kappa|}{\delta}\bigr)$ where $\kappa$ is the set of variant sequence lengths. regression, and the resulting regression array $\hat{f}^{(2)}$ satisfies

\begin{gather*}
\abs{\hat{f}^{(2)}_x-f^{(2)}_x}\leq\frac{\sigma O(\epsilon)(f_{\mathcal{H}_t,x}r_{\Lambda,x}+r_{\mathcal{H}_t\circ\Lambda,x})}{\left(f_{\Lambda,x}-O(\epsilon)r_{\Lambda,x}\right)t_0^2}+o(t_0^2),
\end{gather*}

where the upper bound term adopts the largest bound term referred from Lemma <ref> among multiple evolution with different $t$, and $\sigma$ is a constant regrading the regression method.

The formal proof can be viewed in Appendix <ref>.
The main idea is to combine the accuracy bound from the subroutines of this algorithm.
By tracking the processes of the ratio estimator and the regression, the accuracy of the second-order fidelity can be determined from those inherited bounds.

Note that in this analysis, we choose the ordinary least square regression with an instance of five evenly increasing evolving times, and the parameter to represent the noise propagation during fitting is denoted by $\sigma$.
In order to minimize the overall noise in this phase, we can adjust the unit time $t_0$ to reconcile these two noise sources.
Specifically, under this fitting model, we need to choose $t_0\sim\sqrt[4]{ O(\sigma\epsilon r_\Lambda)}$.

Algorithm <ref> serves as a supplier of second-order Pauli fidelity for the later bin construction as stated in Algorithm <ref>, which is a key structure in the transform from fidelity to error rates in Sec. <ref>.
The bin construction specifies the fidelity terms needed with indices included in the subset ${\sf X}\subset {\sf G}$ since we do not need the whole estimation as explained in the next subsection.
Algorithm <ref> will return the set of desired second-order fidelity $\sf X$ in stabilizer group $\sf G$.
Even though this estimator can estimate all the fidelity in a stabilizer group using only one round of measurements, we still count each query separately regarding the worst-case complexity.

§.§.§ From Pauli fidelity to Pauli error rates

In the last subsection, we have constructed a second-order fidelity estimator.
According to the previous discussion, the rest of this stage is to use this subroutine to get the desired values, namely, to transform the second-order fidelity to squares of decomposition parameters.
We first define the second-order Pauli error rate $p_{\alpha}^{(2)}$ which stands for the second-order coefficient of $p_\alpha$.
According to Eq. (<ref>), we know the specific representation of $p^{(2)}$:

\begin{gather}\label{eq:chi_wot}
\begin{cases}
s_\alpha^2& \text{$\alpha\neq 0$}\\
-\sum_{\beta\in{\sf P}^n,\beta\neq0}s_\beta^2& \text{$\alpha=0$}
\end{cases}.
\end{gather}

Therefore, the target is reduced to implement the transform from second-order fidelity to second-order Pauli error rates.

We know from the definition that the following Walsh-Hadamard transform with Pauli inner product depicts the relation between Pauli fidelity and Pauli error rates.
In the remainder of this paper, we use Walsh-Hadamard transform to denote the Pauli variant when it causes no ambiguity.

\begin{align}\label{eq:PauliWHT}
\begin{split}
p_\alpha=&\frac{1}{4^n}\sum_{x\in{\sf P}^n}(-1)^{\langle x,\alpha\rangle_p}f_x,\ \forall\, \alpha\in{\sf P}^n,\\
f_x=&\sum_{\alpha\in{\sf P}^n}(-1)^{\langle x,\alpha\rangle_p}p_\alpha,\ \forall\, x\in{\sf P}^n.
\end{split}
\end{align}

Since second-order fidelity and Pauli error rates are linear components of the ordinary fidelity and Pauli error rates, they can also be linked by the Walsh-Hadamard transform.
Even though the second-order error rates are sparse according to Assumption $\bf A1$, the second-order fidelity is generally dense.
A faithful implementation of Walsh-Hadamard transform thus requires the summation of an exponential number of terms.
Instead of directly following Eq. (<ref>) with an exponential number of summations, our protocol recruits the idea for Pauli properties estimation [Harper et al., 2021] to efficiently obtain the sparse Pauli error information.

Along with the idea firstly introduced in Ref. [Harper et al., 2021], our protocol utilizes a special structure named bin.
A bin $U$ is a linear combination of a subset of noisy second-order fidelity terms, as defined in Eq. (<ref>), which also indicates the value of a set of second-order Pauli error rates.

The procedure first determines a size index $b< 2n$.
Then each bin is constructed from $2^b$ randomly sampled second-order fidelity terms.

In order to implement the random sampling, the procedure employs several $2n\times b$ random binary matrices $M$ and transforms them to corresponding $M'$, where $M'\coloneqq J_nM$ and $J_n\coloneqq I_n\otimes X$.
By enumerating all $b$-bit string $\ell$, the random matrices $M'$ help to disperse them to $2n$-bit which are valid Pauli indices.
This conversion from $M$ to $M'$ aims to keep the dispersion in the form of Eq. (<ref>).
Furthermore, the procedure recruits additional $2n$-bit vectors $d$ as offsets.
The detailed construction can be viewed in Algorithm <ref>.
After this construction, we observe that each bin consists of $2^{2n-b}$ second-order Pauli error rates and that different bins have no overlap for a given $M$ (or $M'$) as follows,

\begin{align}\label{eq:hashfunction}
U[j]&\coloneqq\frac{1}{2^b}\sum_{\ell\in\mathbb{F}^b_2}(-1)^{\langle\ell,j\rangle}\tilde{f}^{(2)}_{M'\ell+d}=\sum_{\alpha:M^T\alpha=j}(-1)^{\langle \alpha,d\rangle_p}\tilde{p}^{(2)}_\alpha.
\end{align}

This equation can be checked in Lemma <ref>.
Here $\tilde{f}^{(2)}_{M'\ell+d}$ denotes the noisy Pauli fidelity after the linear regression, and $\tilde{p}^{(2)}_\alpha$ is a noisy second-order Pauli error rate with noise comes from the noisy fidelity.
Note that each bin $U[j]$ selects all the error rates of which indices satisfy $M^T\alpha=j$.
The bin construction actually disperses the whole set of error rates into bins with different indices.
This is similar to a hash function that uses $\alpha$ to uniquely decide the bin index $j$ that contains $\tilde{p}^{(2)}_\alpha$.
Therefore, Given a matrix $M$ (or $M'$), all bins form a partition over the second-order Pauli error rates.
We denote every sampling matrix $M$ (or $M'$) a subsampling group as indices it disperses form a subgroup in $\mathbb{F}_2^{2n}$.

Based on the sparse support assumption and a large enough subsampling size $b\sim O(\log s)$, the procedure constructs bins that are likely to contain only one or few nonzero Pauli error rates.
Therefore, we only need polynomial number of fidelity terms.
This detection of the number of nonzero terms in every bin can be achieved with the aid of offset strings $d$.
Intuitively, we construct a bin $U[j]$ with multiple offsets $d$, and we denote these bins with the same bin index $j$ and different offsets to belong to a bin set.
Notably, all bins in a single bin set contain the same set of second-order error rates with different signs according to Eq. (<ref>).
Suppose a bin set only contains one nonzero term $p^{(2)}_\alpha$ and some small noise, we would find absolute values of bins in this bin set remain close.
Similarly, absolute values that are all close to zero indicate that bins contain all zero Pauli error rates, while significant but diverse absolute values imply that there are multiple nonzero rates in a bin set.
Therefore, we can classify bin sets by zero-ton, single-ton, and multi-ton bins.
This bin detector is formalized in Algorithm <ref> of Appendix <ref>.

When we decide the bin is a single-ton, the averaged absolute value would be regarded as the estimation of that second-order Pauli error rate.
Furthermore, the Pauli index can be learned from signs of bins with different offsets.
According to Eq. (<ref>), when the noise is not large enough to flip the sign, there comes an equation to describe the linkage between a bin's sign and the Pauli index of the nonzero term,

\begin{gather}
\sgn{U}=\langle \alpha,d\rangle_p,
\end{gather}

where the sign function is defined as follows,

\begin{equation}
\sgn{x} = \begin{cases}
0 & \text{if } x \ge 0,\\
1 & \text{if } x < 0.
\end{cases}
\end{equation}

By choosing multiple offsets $d$, the procedure can solve the index $\alpha$.

However, it is hard to deal with a multi-ton bin since it is a summation of multiple error rates, and we can hardly extract the information of any single rates directly.
In order to avoid this stuck case, the procedure uses multiple sampling matrices $M$ and constructs different groups of bins.
In each group, the partition of Pauli error rates varies from others.
In this case, if the procedure finds a multi-ton bin, it is supposed to find some single-ton bins in another group with overlapped Pauli error rates.

By peeling out the overlapping error rates detected from single-ton bins, the multi-ton bin contains fewer nonzero error rates and is degraded to single-ton eventually.
We denote this  by the peeling process, and we show an example of the peeling part of a 2-qubit system in Figure <ref>.
Here in Algorithm <ref>, we give a concise overview of how the peeling process works as mentioned above in this subsection and temporarily remove all the technical ingredients.
Shortly, this algorithm invokes Algorithm <ref> to initialize several groups of random bins as above-explained.
With these bins, the peeling process enumerates all of them and uses Algorithm <ref> to detect their nature, including their types, and possible values and indices.
We can then record the desired second-order Pauli rates when finding single-ton bins.

Peeling Decoder (Informal)
[1]
${\tt Input}:$ Number of sampling matrices $C$, number of bins $2^b$
${\tt Initialize}:$ $\mathcal{P}\gets$ empty list to store the sparse indices and second-order Pauli error rates ($\alpha,{p}^{(2)}_{\alpha}$)
Run Algorithm <ref> to construct sampling matrices $\{\mathbf{M}_c\}$ and bins $\{\mathbf{U}_c[j]\}$
we can find new single-ton bins

$c\in C, j\in \mathbb{F}_2^b$
Detect bins $\mathbf{U}_c[j]$ by Algorithm <ref>
the bin is single-ton
record ($\alpha,{p}^{(2)}_{\alpha}$) in $\mathcal{P}$
other groups $c'\neq c$
find bins $\mathbf{U}_{c'}[j']$ contains this ${p}^{(2)}_{\alpha}$
subtract ${p}^{(2)}_{\alpha}$ from $\mathbf{U}_{c'}[j']$
the bin is not single-ton
continue to next bin
Return: the stored indices and second-order Pauli error rates $\mathcal{P}$

We will further give an exhaustive illustration of this peeling process in Algorithm <ref> in Appendix <ref>.

It is worth noting that since the procedure needs to peel out nonzero rates and generate new single-ton bins, it triggers noise propagation among bins during peeling.
In Appendix <ref>, we analyze the behavior of the whole transformation with noisy fidelity terms gained from Algorithm <ref>.
We show that an index array $\mathbf{T}$ indicates the size of variances for noise in each bin.
Moreover, since the noise of fidelity terms from Algorithm <ref> is not necessarily unbiased, we also consider the bias effects in the noise propagation during the peeling process, which is more complicated than the ideal version of transform introduced in [Harper et al., 2021].
Combining the noise analysis and the mechanisms of the sampling and peeling, as well as Proposition <ref>, we can claim the following.

Suppose the Assumption A1 & A2, <ref>, and <ref> hold.
Execute Algorithm <ref> with $t_0$ satisfies $\|\mathcal{H}_{5t_0}-\mathcal{I}\|<\frac{1}{4}$ and $l=\frac{2}{\epsilon^4}\log\left(\frac{4ns|\kappa|}{\delta}\right)$ sequences with a set $\kappa$ of variant sequence lengths, $B=2^b=\max\left\{O(s),O\left(\frac{\epsilon^4}{t_0^4\epsilon_0^2}\right)\right\}$, $C=O(1)$, the unit time length $t_0$, and offsets $\mathbf{D}$ with $P=O(n)$ for each subsampling group.
The transformer will estimate all Pauli error rates $\hat{\mathbf{p}}$
with accurate support information and error bounds to be $\|\hat{\mathbf{p}}-\mathbf{p}\|_\infty\leq O\left(\frac{\epsilon^2}{t_0^2\sqrt{s}}\right)$.
Therefore, the absolute values of these nonzero decomposition parameters can be estimated by $\|\abs{\hat{s}^\star}-\abs{s^\star}\|_\infty\leq O\left(\frac{\epsilon}{t_0\sqrt[4]{s}}\right)$
The estimation works successfully with probability at least $1-\delta-O\left(\frac{1}{s}\right)$.

The proof relies on several lemmas stated in Appendix <ref>.
The proof starts by fetching the error bounds from Proposition <ref> and analyzing the noise in every bin.
By proving Lemma <ref> and Lemma <ref>, it shows that sampling bins serve as hash functions and that the bin detector works accurately with a given level of noise.
Based on these, the index $T$ records the noise size and helps the procedure to adjust the noise expectation.
Therefore, we can prove the success of this peeling process.
The formal proof can be viewed in Appendix <ref>

In summary, the first stage can detect the second-order Pauli error rates of the Hamiltonian evolution channel or, equivalently, the squared decomposition parameters of the Hamiltonian.
The protocol first uses random gate twirling to detect the Pauli fidelity of the Hamiltonian channel, followed by regression to eliminate evolving times.
It then constructs bins from several sets of second-order fidelity as the estimator of all nonzero second-orderPauli error rates.
The detailed protocol with formalized subroutines is illustrated in Appendix <ref>.
The validity and efficiency of estimation for the whole set of squared decomposition parameters ${\bf s}^\star$ is proven rigorously in Appendix <ref>.
We also discuss the possible benefits the estimation gets if some prior knowledge about the underlying structure of the Hamiltonian is available.
The knowledge might be collected from the understanding of the experimental devices or the properties of the tasks running on the devices.
We can expect a more accurate and reliable result with this kind of aid, which is illustrated in Appendix <ref>.

§.§ Stage 2. Sign Estimation

After the first stage, we can detect all nonzero second-order Pauli error rates of the Hamiltonian channel.
Nevertheless, these terms only carry the square values of parameters as stated in Eq. (<ref>).
As a complement, the protocol for this stage is to estimate all the sign information.
The idea of our sign estimator is based on constructing equations as in quantum process tomography [Altepeter et al., 2003, Mohseni et al., 2008, Poyatos et al., 1997] and utilizing the information gained from the first stage.

Suppose the procedure chooses $m$ states and measurement settings $\{(\rho_k,M_k)\}_{k=1}^m$.
As shown in circuit (c) of Figure <ref>, the procedure measures these states on the Pauli operators after evolving under the target Hamiltonian.
For each pair, the procedure collects multiple measurement results of the evolved state with different evolving times $\{t\}$.
By the regression over $\{t\}$, the procedure focuses on the first-order coefficients of the measurements.
According to the expansion in Eq. (<ref>), the first-order measurement for a pair of SPAM settings is depicted by the process equation as
\begin{gather}\label{eq:linearordereq}
\Tr(\mathcal{H}^{(1)}_t(\rho)M)=i\sum_{\alpha\in{\sf P}^n}s_\alpha \Tr(\rho [P_\alpha, M]),
\end{gather}

where $\mathcal{H}^{(1)}_t$ is the first-order expansion of the Hamiltonian evolution channel, and $[P_\alpha,M]$ is the commutator.
After the procedure computes all these first-order coefficients from measurements, we collect a series of observations for different SPAM settings
where each corresponds to a linear combination of trace results for commutators as in Eq. (<ref>).
This reduces the problem of estimating Hamiltonian parameters ${\bf s}$ to solving the following process equations,

\begin{gather}\label{eq:LEQ}
\mathcal{M}=\Phi \cdot {\bf s},
\end{gather}

where $\Phi$ is the matrix consisting of $\{i \Tr(\rho_k [P_\alpha, M_k])\}_{k,\alpha}$.

Generally, the size of unknown parameters $\bf s$ scales exponentially with $n$, rendering the solution intractable.
In our execution, the sparsity assumption guarantees that there are only sparse nonzero variables to be solved.
Besides, from the first stage, our procedure has extracted the square values of decomposition parameters, which also includes the precise sparse support information.
The size of linear equations is reduced to $s$ by shrinking the target parameter from $\bf s$ to ${\bf s}^\star$ and reducing the coefficient matrix.
Therefore, we have the following equations to illustrate the polynomial-size problem,

\begin{gather}
\mathcal{M}=\Phi \cdot {\bf s}^\star.
\end{gather}

The problem will be more knotty when SPAM errors are taken into account.
Suppose measurements of the evolved states for the chosen observable are corrupted by additive noise.
The noise will be transmitted in the linear regression stage, and the effects of noise will be enlarged due to the short evolving time $t$.
In this case, the first-order part of measurements $\mathcal{M}$ will also be corrupted by the noise as $\hat{\mathcal{M}}=\mathcal{M}+\omega$.

To suppress these effects, by choosing the coefficient matrix $\Phi$ with random local Pauli eigenstates $\rho$ and Pauli measurements $M$, we can construct a coefficient matrix which limits the noise propagation during solving the unknown parameters.
Regarding the noise and biased observation vector $\hat{\mathcal{M}}$, the original equations will be derived to a new optimization problem,

\begin{align}\label{eq:optimize}
\begin{split}
\min_x \|x\|_1&\\
||\Phi x- \hat{\mathcal{M}}||_2&\leq \epsilon,
\end{split}
\end{align}

where $\epsilon$ is set to avoid the stuck case when the equation has no solution, and $x$ is an unknown vector with dimension $s$.
It is easy to see that when $\epsilon\geq\|\omega\|_2$, the true parameter ${\bf s}^\star$ lies in the feasible set.
The algorithm for sign estimation is shown in the following with a proper choice of $\epsilon\geq\|\omega\|_2$ according to Lemma <ref>.

Sign Estimation: SE($\{\abs{p_\alpha},\alpha\},m,t_1$)
[1]
Input: The nonzero error rates support of $\mathcal{P}$ from Algorithm <ref>, a positive integer $m$, and the unit time length $t_1$
Prepare $m$ random local Pauli eigenstates and $m$ random Pauli operators for measurements $(\rho_{\gamma_i},P_{\beta_i})_{i=1}^m$
Choose evolving times $\{i\cdot t_1\}$ which are multiples of $t_1$ and construct the Hamiltonian evolution channel $\{\mathcal{H}_t\}$ for evolving time
Measure $\{\Tr(\mathcal{H}_t(\rho_\gamma)P_\beta)\}$ and calculate the coefficient matrix $\Phi$ for Pauli indices belong to $\mathcal{P}$
Process ordinary linear regression on the measurement outcomes along the time $\{i\cdot t_1\}$ and get observations $\hat{\mathcal{M}}$
Calculate $\sigma\coloneqq\frac{\sum\abs{i\cdot t_1-\bar{t}}}{\sum(i\cdot t_1-\bar{t})^2}\cdot t_1$ as the regression coefficient
Choose $\epsilon=\sqrt{m}(\sigma\cdot\tau/t_1+o(t_1))$ for the optimization problem in Eq. (<ref>)
Solve the organized optimization problem with solution $x^\star$
solution $x^\star$

By sampling the input states and measurement operators, the $\Phi$ matrix possesses approximate restricted isometry property (RIP) as defined in Definition <ref>.
This property, by definition, makes sure that the norm of an arbitrary vector after the matrix product is close to the norm of the original vector.
We prove this RIP of the randomly sampled process matrix for all possible decomposition parameters with a high probability according to concentration inequalities in Lemma <ref>, which we refer to as the approximate RIP.
Based on the definition of RIP, the noise of resulting nonzero parameters can be bounded by the chosen relaxation $\epsilon$.
Since the sign estimator aims at extracting the discrete $\pm1$ information, it tolerates a reasonable amount of noise that would not cause a sign flip in every solved parameter.
Therefore, the noise is not destructive for the estimator, given that we use multiple rounds of equations to minimize the statistical effects.

Suppose the absolute value estimator perfectly returns the support information of decompose parameters of the Hamiltonian, and Assumption A1 & A2 and <ref> hold.
Run Algorithm <ref> with the support information of decomposition parameters and the unit time length $t_0$.
By setting $m$ to be $O(s)\leq m\leq O\left(\frac{t_1^2\epsilon_0}{\tau^2}\right) $ The solution of Algorithm <ref> contains perfect sign information of all nonzero decomposition parameters with probability at least $1-\mathrm{e}^{-O(m-s)}$.

The full proof is shown in Appendix <ref>.

As there exists an upper bound of the $m$ due to the accumulation of the SPAM errors, we cannot directly use larger $m$ to significantly suppress the sign errors. On the other hand, we can choose multiple blocks of equations to estimate signs for multiple times. Based on this, we can further reduce the failure probability by an additional majority vote mechanism when we hold multiple sets of signs.

In this stage, we propose a protocol to extract the sign information of the Hamiltonian channel by constructing equations from random state preparation and measurements.
With the aid of absolute values from the former stage, the size of this problem is reduced to linear with the sparsity.
Based on the randomness, the equations will be robust against the measurement noise.
Therefore, this efficient and robust estimation fulfills the need for sign estimation.

§ MAIN THEOREM

Executing the two-stage protocol illustrated above, we can estimate the whole information of the nonzero decomposition parameters ${\bf s}^\star$, namely, we can recover the Hamiltonian operator faithfully from this estimation as in Algorithm <ref>.

Hamiltonian Estimator
[1]
Input: An oracle to control the unknown Hamiltonian $H$, positive integer $l$, $b$, $C$, $P$ and $m$, and the unit time lengths $t_0$ and $t_1$
Input: Offsets $\mathbf{D}$ for subsampling
Call Algorithm <ref> with $l$, $b$, $C$, $P$, $t_0$ and offsets $\mathbf{D}$ to construct bin sets with random subsampling groups and estimate all the second-order error rates $\mathcal{P}$
$\text{Supp}({\bf s}^\star)\leftarrow\{\alpha\,|\,p_\alpha\neq0\}$
Call Algorithm <ref> with $\text{Supp}({\bf s}^\star)$, $m$ and $t_1$ to estimate sign information $\sgn{{\bf s}^\star}$ nonzero parameters
${\bf s}^\star\leftarrow\{\sqrt{p_\alpha}\cdot\sgn{s_\alpha}\,|\,p_\alpha\in\mathcal{P}\}$
${\bf s}^\star$

We combine the guarantee of our protocol from the guarantees of the two stages since our main protocol is a composite algorithm.
Our guarantee, therefore, does rely on some mild assumptions, which are either inherited from foregoing subroutines or from requirements to bond these stages in a self-consistent way.

We summarize the results in Theorem <ref>, which states that our algorithm can estimate the decomposition parameters of an unknown sparse Hamiltonian with high accuracy and vanishing failure probability.
It shows that our protocol is provably robust against circuit noise and SPAM errors under mild assumptions, which can be viewed from the guarantees of the fidelity estimator and the sign estimator separately.
This algorithm is also provably scalable for both quantum measurement complexity and classical computational complexity.
Moreover, the quantum complexity claimed above is for the worst case, and this complexity can be further reduced by merging all queries of qubitwise commutative Pauli fidelity, which relies on the choices of stabilizer groups elaborated in Appendix <ref>.
This merging can usually reduce the needed rounds of fidelity detectors by one or two orders of magnitude.
Therefore, our protocol provides an implementable method for Hamiltonian learning.

Suppose Assumption A1 & A2, <ref>, <ref>, and <ref> hold. Run Algorithm <ref> with $l=\frac{2}{\epsilon^4}\log\left(\frac{4ns|\kappa|}{\delta}\right)$ sequences for each length where $\kappa$ is the set of variant sequence lengths, $B=2^b=\max\left\{O(s),O\left(\frac{\epsilon^4}{t_0^4\epsilon_0^2}\right)\right\}$, $C=O(1)$, $m=O(s)$, unit times $t_0$ satisfies $\|\mathcal{H}_{5t_0}-\mathcal{I}\|<\frac{1}{4}$ and  and $t_1$, and the offsets $D$ with size $P=O(n)$ for each subsampling group.
The Hamiltonian estimator will return all nonzero decomposition parameters with the perfect support estimation and $\|\hat{\bf s}^\star-{\bf s}^\star\|_\infty\leq O\left(\frac{\epsilon}{t_0\sqrt[4]{s}}\right)$, which succeeds with probability at least $1-\delta-O\left(\frac{1}{s}\right)$.
The circuit measurement complexity of this execution is $\tilde{O}\left(\frac{sn}{\epsilon^4}\right)$ where $\tilde{O}$ notation ignores the logarithmic terms.
The post-processing complexity is $\text{poly}(n,s)$, where poly denotes that the scaling is polynomial with the elements.

This statement can be regarded as a  combination of Proposition <ref>, <ref>, and <ref>.
During the combination, we need first to consider the overall precision resulting from subsequently invoking subroutines.
According to the target overall precision, we shall calculate the precision for each subroutine and then count and sum up the required complexities for all these subroutines.
The rigorous counting and derivation are given in Appendix <ref>.

It is true that we can choose other fitting models and strategies for extracting the corresponding information in both stages, and we leave a detailed study of the strategy for future studies.
In the execution of the algorithm, we need two inputs, $t_0$ and $t_1$, to serve as the fitting times of stages 1 and 2, respectively.
Nevertheless, there is a natural trade-off toward the size of those evolving times $t_0$ and $t_1$.
If the evolving times are too short, this could suppress systematic errors from high-order terms in the expansion, while it would enlarge the errors from imperfect observations and measurements by fitting.
If the times are too long, this will lead to large systematic errors.
Therefore, the balanced choice of $t_0$ and $t_1$ can be viewed in Remark <ref> and <ref>, respectively.
To further improve the performance, a possible modification is to consider those higher-order terms in the fitting model to mitigate the systematic errors.
In appendix <ref>, we exhibit different trials of executions by fitting higher-order terms and taking the desired coefficients.
This can vastly reduce the amplitude of the systematic errors and allow us to implement a longer time evolution.
Moreover, the idea from [Yuan et al., 2016] can help to construct an exponential suppression of the high-order terms with the number of data points.
In general, we expect the choices of evolving times to satisfy that the effects of fluctuation errors from observations dominate.

§ NUMERICAL RESULTS

In this section, we exhibit several numerical results to show the validity and accuracy of the protocol for learning an unknown Hamiltonian with sparse decomposition parameters on the Pauli basis.
These numerical results work as strong evidence that verifies the foregoing intuitively illustrated protocol.
They also reveal some critical properties of the Hamiltonian learning method in the realistic implementation as we discuss shortly.

In order to simplify the numerical simulation without modifying the core components of the protocol, we make some mild assumptions in the simulation, which would not weaken the results arguably.
Among all of these simulations, instead of running the cascading circuits to extract the fidelity information as in <Ref>, we assume there exists an oracle to return measurement outcomes of input states under the evolution with the input time of the unknown Hamiltonian, which helps to directly get the Pauli fidelity terms of evolution.
In this sense, we count each desired fidelity term as one query, and the number of calls for the oracle represents the total measurements.
This is based on the understanding that the complexity of the first stage is always dominant in the total number of circuit measurements, which makes the tracking of fidelity query a faithful indicator of the total quantum complexity.

As for the simulation of noise effects, we add zero-mean Gaussian noise $\mathcal{N}(0,10^{-3})$ to the queried Hamiltonian fidelity terms to imitate the fluctuation noise from the fidelity estimator.
Along the fitting process and the sparse transformation, the noise effects would be exhibited in the absolute value estimation.
Moreover, we add the Gaussian noise $\mathcal{N}(0,10^{-3})$ to the measurement outcomes in the sign estimation as a representation of the SPAM effects when constructing the process equations.
These errors will be propagated by the fitting and the optimization stated in stage 2, and the resulting signs reflect overall noise effects.
This makes it possible to observe the noise-resilience of our learning protocol from the accuracy of each part in the simulation.

Estimation of the 6-qubit random transverse field Ising model with different choices of the bin numbers, $b$. In the box plot, we exhibit the statistical properties of errors from the outcomes using different $b$.
In each case, the data is gained from the estimation of 50 random TFIM Hamiltonians with 10 rounds each.
For every box, the orange line denotes the median of the error distribution, while the upper and lower bounds represent the $75\%$ and $25\%$, respectively.
By increasing $b$, the procedure estimates all the decomposition parameters with smaller errors. A sharp vanishing over the distribution of the estimation error is witnessed when $b$ grows up to $4$.
In the inset, we use the violin plot to show a complete view of the distribution of the query numbers and recovery errors.
Every violin along its axis represents a distribution of the corresponding value, and the widths everywhere illustrate the probability density.
There is also a sharp shrink of the error bar while the querying number is growing along with the increasing $b$.

Numerical results for the transverse field Ising model and molecular Hamiltonians. In (a), we show the distribution of the results for random TFIM Hamiltonian with different system sizes varying from 1 to 7.
For each $n$, the procedure runs for 50 random TFIM Hamiltonian operators with 10 rounds for each to detect the average reconstruction results.
The parameter $b$ is chosen by witnessing the threshold behaviors by varying $b$ as a trial process. For example, we choose $b=5$ for 6-qubit random TFIM Hamiltonian learning, and all others can be viewed in the appendix.
In (b), we show the ability to recover different molecular Hamiltonians, H$_2$ (4 qubits), H$_3$ (6 qubits), and LiH (6 qubits).
For each molecular system, we execute our method for 500 independent rounds to gather these statistical properties.
The bin parameters are also chosen according to the threshold behaviors.
In (c), we show an average reconstruction among 20 rounds of the 7-qubit TFIM Hamiltonian, and we use the error bars to display the standard deviations of parameters. The terms are sorted by the absolute values of decomposition parameters.
The labels of the $x$-axis indicate both the Pauli indices as well as the signs of the estimation (the signs are all corrected estimated, as one can see from (a)).
In (d), the procedure estimates the Hamiltonian for the H$_4$ (8 qubits) molecule, and we exhibit the 50 dominant terms of both ideal and reconstructed parameters among 184 nonzero terms.
This result is an average of 8 rounds, and the error bars are constructed according to standard deviations.
In each round, the simulation uses 5 blocks to vote for the final estimation of signs to suppress the failure probability of sign flips.
Similar to (c), this plot also uses $x$-axis coordinates to denote the Pauli indices and signs, and the estimation shows great agreement with the ideal parameters.

First, we study how different choices of the number of bins $b$ affect the efficacy of the protocol.

We consider the 6-qubit transverse field Ising model (TFIM) Hamiltonian with random interaction parameters,

\begin{gather}\label{eq:randtfim}
\end{gather}

For different choices of $b$, we run the protocol for 50 randomly chosen Hamiltonian operators, where each of them is estimated for 10 independent rounds.
As defined in <Ref>, the parameter $b$ reflects the number of bin sets in every subsampling group.
With a larger $b$, the procedure hashes Pauli error rates more dispersively, meaning that every bin contains fewer underlying Pauli error rates as well as nonzero error rates.
Since the peeling process relies on the sparsity of nonzero terms in every bin, procedures with larger $b$ would generate more accurate estimations, where the reconstruction error is depicted by the relative 1-norm distance,

\begin{gather}\label{Eq:relativeerror}
\end{gather}

On the other hand, a larger $b$ is equivalent to querying more fidelity terms in the subsampling step, which leads to a trade-off of choosing a proper $b$.
In Figure <ref>, we can view in the box-plot that the recovery errors decay vastly when the procedure increases $b$ from $3$ to 4.
The error is even smaller when it continues to increase $b$.
However, as shown in the inset of Figure <ref>, the execution overhead also increases with a $b$.
Therefore, a wise choice of $b$ can be explored as the threshold under which the error is significantly large, where we choose it as $b=5$ in our following simulation.
This remains consistent with the previous analysis about the peeling step, where the procedure expects the bin number to be larger than the sparsity.

However, we may not have access to the error in practice for an unknown quantum system.
The problem then rises as to how to decide a proper $b$ in practice.
That relies on adaptively increasing $b$ from $b_0$. to choose the proper bin size $b$.
For every $b$, we run the protocol multiple times to first observe the peeling process to find out whether some multi-ton bins get stuck after the peeling process.
If this happens with significant frequencies, it indicates that the current $b$ is not enough.
We can further observe the convergence of results from different rounds.
As shown in Figures <ref> and <ref>, the variance between different rounds of estimation will decrease sharply when $b$ is large enough i.e., $b\geq\lceil\log s\rceil$, and we denote the first $b$ to observe this decaying as the threshold.
Since the execution complexity scales as $O(n2^b/\epsilon^4)$ when we choose $b$ arbitrarily, the overall adaptive trial of $b$ to $b=\lceil\log s\rceil$ only increases the complexity $O(ns/\epsilon^4)$ by a constant factor.

Next, we study the scaling of the protocol for different Hamiltonians with varying system sizes.

We first consider the TFIM with different numbers of qubits.
As in Figure <ref>(a), we show a complete view of the distribution of outcomes in different $n$.
We choose $b$ for each $n$ according to the abovementioned threshold idea.
The errors are also measured by the relative 1-norm distance, and the major parts concentrate on small errors in each case.
While an interesting observation is that the lower bound of these errors grows with the system size, this is mainly due to the regression error as we execute the ordinary least-square linear regression in this simulation.
Even though we only consider the second-order effects here, We further analyze the effects of different fitting settings in the following and in Appendix <ref>, which suggests that considering the higher even-order terms in the fitting can significantly improve the estimation.

The TFIM model still has a simple structure with nearest-neighbour 1D interactions.
We then consider more complicated examples for benchmark, i.e., molecular Hamiltonians, which have been demonstrated to be potential candidates for quantum computing and quantum simulation [Troiani et al., 2005, Santini et al., 2011].
We encode the fermionic Hamiltonian with qubits, which generally has much more complicated Pauli decomposition parameters that are not local but have sophisticated multiqubit interaction.

We note that though the nonzero terms are sparse compared to the exponential size of Pauli decomposition, the nonzero support is still quite large, which grows quartically with the system size. If we could measure each Pauli decomposition parameter to accuracy $\varepsilon$, the reconstruction error defined in Eq. (<ref>) also grows quartically. In order to show a system-size independent error rate,  we recruit the average 1-norm distance,

\begin{gather}
\end{gather}

to depict the reconstruct errors in the molecular cases.

In Figure <ref>(d), we display the distribution of the estimation of several molecular Hamiltonians only considering the zero- and second-order terms during fitting.
Generally, both the quantum complexity increase when system sizes grow, while the average reconstruction errors remain comparably constant.

In the remaining two subplots, we use the bar plot to exhibit the estimation of every parameter.
In Figure <ref>(c), we exhibit the average estimation of 50 independent rounds of executions for a 7-qubit TFIM Hamiltonian with random parameters.
Furthermore, we use the higher-order fitting considering the fourth-order terms of the Hamiltonian expansion to extract a much more precise estimation.
The higher-order fitting will require an additional one or two times query complexity, while it brings a greatly precise estimation.
Therefore, as can be viewed in (c), the distances between pairs of original values and reconstructed values are generally close to zero.
In Figure <ref>(d), we estimate the Hamiltonian of the hydrogen chain $H_4$.
We show the 50 dominant terms with both the ideal and reconstructed decomposition parameters [There are $184$ nonzero terms in the $H_4$ Hamiltonian in total.].
Since the number of nonzero terms is great and they vary from large terms to very small terms, we also consider the fourth-order terms' effects during fitting to improve the estimation of absolute values.
In order to provide the perfect sign estimation, we add a further majority voting from 5 blocks of resulting signs when deciding the final signs.
In this average estimation of 8 independent rounds, random noise causes little harm to the reconstruction, and our procedure returns precise estimations of the most significant terms among the nonzero parameters.
On the other hand, the error may grow when the underlying value turns small since the procedure cannot divide the noise effects from the small parameters accurately.
This result sheds light on how the algorithm would execute on Hamiltonian with much more nonzero decomposition parameters.

We refer to Appendix <ref> for more discussions on the implementation and numerical results. The codes to generate these simulations are shared in [Han and Yu, 2021].

§ CONCLUSION

In this work, we have proposed a universal, robust, and efficient Hamiltonian learning protocol. The protocol is universal in the sense that it estimates any $n$-qubit Hamiltonians with sparse Pauli decomposition parameters without prior knowledge of the Hamiltonian structure. For example, the protocol can estimate Hamiltonians consisting of global interactions, or namely those high-weight Pauli terms.
Next, our protocol is robust since it estimates the Hamiltonian noise-resiliently against SPAM errors and certain amounts of the circuit and shot noise.
The protocol is also efficient with polynomially scaling classical and quantum complexities.
Besides, the protocol directly exploits the time evolution of the Hamiltonian instead of requiring eigenstates or thermal states of the system, which make it reliable in the real-world implementation.
All these results have been rigorously proven in our theoretical analysis in Appendix as well as numerically verified for different sizes of transverse field Ising models and different electronic Hamiltonians in the main text.

Our work provides a scalable and more practically-minded Hamiltonian learning method compared to the previous protocols based either on explicit structures or informative states.
This protocol could be applied for calibrating quantum computing hardware, studying the interaction structures of intricate many-body quantum systems, detecting interesting many-body physics phenomena, etc.
As our method exhibits a direct linkage between the fidelity estimation and Hamiltonian learning, we can expect that development in the fidelity estimation also improve the learning of Hamiltonian.
For example, the proposal in [Chen et al., 2022] implies that entanglement states can help to reduce the overheads for fidelity estimation, which can also be employed in the Hamiltonian learning.
Our method only exploits the incoherent part of the Hamiltonian time evolution process to estimate the absolute value of the decomposition parameters, and one may expect to make use of the coherent part (the first-order expansion) to more efficiently extract the Hamiltonian information. Nevertheless, as we have discussed previously, this would require a new SPAM-robust process detection mechanism, which also deserves its own investigation.
Another point full of practical interest is to combine the fitting method that holds a clear trade-off between the higher-order effects and the random noise with our protocol.
As inspired by Figure <ref> and <ref>, a wise choice of fitting method can significantly improve the estimation, while we still need to balance the effects of the overall errors.
Moreover, a rigorous analysis of more general non-sparse Hamiltonian learning remains open.
we are interested in the case when the parameters of the unknown Hamiltonian distribute unevenly rather than sparsely.
As a rough idea, compressed sensing might help in this circumstance.

§ ACKNOWLEDGEMENT

This work is supported by the National Natural Science Foundation of China Grant No. 12175003.
The numerics is supported by High-performance Computing Platform of Peking University.

[Arute et al., 2019]
Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C. Bardin, Rami
Barends, Rupak Biswas, Sergio Boixo, Fernando G. S. L. Brandao, David A.
Buell, Brian Burkett, Yu Chen, Zijun Chen, Ben Chiaro, Roberto Collins,
William Courtney, Andrew Dunsworth, Edward Farhi, Brooks Foxen, Austin
Fowler, Craig Gidney, Marissa Giustina, Rob Graff, Keith Guerin, Steve
Habegger, Matthew P. Harrigan, Michael J. Hartmann, Alan Ho, Markus Hoffmann,
Trent Huang, Travis S. Humble, Sergei V. Isakov, Evan Jeffrey, Zhang Jiang,
Dvir Kafri, Kostyantyn Kechedzhi, Julian Kelly, Paul V. Klimov, Sergey Knysh,
Alexander Korotkov, Fedor Kostritsa, David Landhuis, Mike Lindmark, Erik
Lucero, Dmitry Lyakh, Salvatore Mandrà, Jarrod R. McClean, Matthew
McEwen, Anthony Megrant, Xiao Mi, Kristel Michielsen, Masoud Mohseni, Josh
Mutus, Ofer Naaman, Matthew Neeley, Charles Neill, Murphy Yuezhen Niu, Eric
Ostby, Andre Petukhov, John C. Platt, Chris Quintana, Eleanor G. Rieffel,
Pedram Roushan, Nicholas C. Rubin, Daniel Sank, Kevin J. Satzinger, Vadim
Smelyanskiy, Kevin J. Sung, Matthew D. Trevithick, Amit Vainsencher, Benjamin
Villalonga, Theodore White, Z. Jamie Yao, Ping Yeh, Adam Zalcman, Hartmut
Neven, and John M. Martinis.
Quantum supremacy using a programmable superconducting processor.
Nature, 5740 (7779):0 505–510, October 2019.

[Neill et al., 2021]
C Neill, T McCourt, X Mi, Z Jiang, MY Niu, W Mruczkiewicz, I Aleiner, F Arute,
K Arya, J Atalaya, et al.
Accurately computing the electronic properties of a quantum ring.
Nature, 5940 (7864):0 508–512, 2021.

[Ebadi et al., 2021]
Sepehr Ebadi, Tout T Wang, Harry Levine, Alexander Keesling, Giulia Semeghini,
Ahmed Omran, Dolev Bluvstein, Rhine Samajdar, Hannes Pichler, Wen Wei Ho,
et al.
Quantum phases of matter on a 256-atom programmable quantum
Nature, 5950 (7866):0 227–232, 2021.

[Wu et al., 2021]
Yulin Wu, Wan-Su Bao, Sirui Cao, Fusheng Chen, Ming-Cheng Chen, Xiawei Chen,
Tung-Hsun Chung, Hui Deng, Yajie Du, Daojin Fan, et al.
Strong quantum computational advantage using a superconducting
quantum processor.
Physical review letters, 1270 (18):0 180501,

[Gong et al., 2021]
Ming Gong, Shiyu Wang, Chen Zha, Ming-Cheng Chen, He-Liang Huang, Yulin Wu,
Qingling Zhu, Youwei Zhao, Shaowei Li, Shaojun Guo, et al.
Quantum walks on a programmable two-dimensional 62-qubit
superconducting processor.
Science, 3720 (6545):0 948–952, 2021.

[Zhong et al., 2020]
Han-Sen Zhong, Hui Wang, Yu-Hao Deng, Ming-Cheng Chen, Li-Chao Peng, Yi-Han
Luo, Jian Qin, Dian Wu, Xing Ding, Yi Hu, et al.
Quantum computational advantage using photons.
Science, 3700 (6523):0 1460–1463, 2020.

[Mi et al., 2021]
Xiao Mi, Matteo Ippoliti, Chris Quintana, Ami Greene, Zijun Chen, Jonathan
Gross, Frank Arute, Kunal Arya, Juan Atalaya, Ryan Babbush, et al.
Time-crystalline eigenstate order on a quantum processor.
Nature, pages 1–1, 2021.

[Kandala et al., 2017]
Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink,
Jerry M Chow, and Jay M Gambetta.
Hardware-efficient variational quantum eigensolver for small
molecules and quantum magnets.
Nature, 5490 (7671):0 242–246, 2017.

[Kjaergaard et al., 2020]
Morten Kjaergaard, Mollie E. Schwartz, Jochen Braumüller, Philip Krantz, Joel
I.-J. Wang, Simon Gustavsson, and William D. Oliver.
Superconducting qubits: Current state of play.
Annual Review of Condensed Matter Physics, 110
(1):0 369–395, 2020.

[Zhang et al., 2017]
Jiehang Zhang, Guido Pagano, Paul W Hess, Antonis Kyprianidis, Patrick Becker,
Harvey Kaplan, Alexey V Gorshkov, Z-X Gong, and Christopher Monroe.
Observation of a many-body dynamical phase transition with a 53-qubit
quantum simulator.
Nature, 5510 (7682):0 601, 2017.

[Georgescu et al., 2014]
Iulia M Georgescu, Sahel Ashhab, and Franco Nori.
Quantum simulation.
Reviews of Modern Physics, 860 (1):0 153,

[Altman et al., 2021]
Ehud Altman, Kenneth R Brown, Giuseppe Carleo, Lincoln D Carr, Eugene Demler,
Cheng Chin, Brian DeMarco, Sophia E Economou, Mark A Eriksson, Kai-Mei C Fu,
et al.
Quantum simulators: Architectures and opportunities.
PRX Quantum, 20 (1):0 017003, Feb 2021.

[Hauke et al., 2012]
Philipp Hauke, Fernando M Cucchietti, Luca Tagliacozzo, Ivan Deutsch, and
Maciej Lewenstein.
Can one trust quantum simulators?
Reports on Progress in Physics, 750 (8):0
082401, 2012.

[Granade et al., 2012]
Christopher E Granade, Christopher Ferrie, Nathan Wiebe, and David G Cory.
Robust online hamiltonian learning.
New Journal of Physics, 140 (10):0 103013,

[Wiebe et al., 2014]
Nathan Wiebe, Christopher Granade, Christopher Ferrie, and David Cory.
Quantum hamiltonian learning using imperfect quantum resources.
Physical Review A, 890 (4):0 042314,

[Wiebe et al., 2014]
Nathan Wiebe, Christopher Granade, Christopher Ferrie, and David G Cory.
Hamiltonian learning and certification using quantum resources.
Physical review letters, 1120 (19):0 190501,

[Krastanov et al., 2019]
Stefan Krastanov, Sisi Zhou, Steven T Flammia, and Liang Jiang.
Stochastic estimation of dynamical variables.
Quantum Science and Technology, 40 (3):0
035003, 2019.

[Wang et al., 2017]
Jianwei Wang, Stefano Paesani, Raffaele Santagati, Sebastian Knauer, Antonio A
Gentile, Nathan Wiebe, Maurangelo Petruzzella, Jeremy L O’Brien, John G
Rarity, Anthony Laing, et al.
Experimental quantum hamiltonian learning.
Nature Physics, 130 (6):0 551–555, 2017.

[O'Brien et al., 2021]
Thomas E O'Brien, Lev B Ioffe, Yuan Su, David Fushman, Hartmut Neven, Ryan
Babbush, and Vadim Smelyanskiy.
Quantum computation of molecular structure using data from
challenging-to-classically-simulate nuclear magnetic resonance experiments.
arXiv preprint arXiv:2109.02163, 2021.

[Haah et al., 2021]
Jeongwan Haah, Robin Kothari, and Ewin Tang.
Optimal learning of quantum hamiltonians from high-temperature gibbs
arXiv preprint arXiv:2108.04842, 2021.

[Lloyd, 1996]
Seth Lloyd.
Universal quantum simulators.
Science, pages 1073–1078, 1996.

[Das and Chakrabarti, 2008]
Arnab Das and Bikas K Chakrabarti.
Colloquium: Quantum annealing and analog quantum computation.
Reviews of Modern Physics, 800 (3):0 1061,

[Friedenauer et al., 2008]
Axel Friedenauer, Hector Schmitz, Jan Tibor Glueckert, Diego Porras, and Tobias
Simulating a quantum magnet with trapped ions.
Nature Physics, 40 (10):0 757–761, 2008.

[Aspuru-Guzik and Walther, 2012]
Alán Aspuru-Guzik and Philip Walther.
Photonic quantum simulators.
Nature physics, 80 (4):0 285–291, 2012.

[Bernstein and Vazirani, 1997]
Ethan Bernstein and Umesh Vazirani.
Quantum complexity theory.
SIAM Journal on computing, 260 (5):0
1411–1473, 1997.

[Shor, 1995]
Peter W. Shor.
Scheme for reducing decoherence in quantum computer memory.
Phys. Rev. A, 520 (4):0 R2493–R2496, Oct

[Grover, 1996]
Lov K Grover.
A fast quantum mechanical algorithm for database search.
In Proceedings of the twenty-eighth annual ACM symposium on
Theory of computing, pages 212–219, 1996.

[Harrow et al., 2009]
Aram W Harrow, Avinatan Hassidim, and Seth Lloyd.
Quantum algorithm for linear systems of equations.
Physical review letters, 1030 (15):0 150502,

[Krantz et al., 2019]
Philip Krantz, Morten Kjaergaard, Fei Yan, Terry P Orlando, Simon Gustavsson,
and William D Oliver.
A quantum engineer's guide to superconducting qubits.
Applied Physics Reviews, 60 (2):0 021318,

[Trabesinger, 2012]
Andreas Trabesinger.
Quantum simulation.
Nature Physics, 80 (4):0 263–263, 2012.

[Pirandola et al., 2018]
Stefano Pirandola, B Roy Bardhan, Tobias Gehring, Christian Weedbrook, and Seth
Advances in photonic quantum sensing.
Nature Photonics, 120 (12):0 724–733, 2018.

[Degen et al., 2017]
Christian L Degen, F Reinhard, and Paola Cappellaro.
Quantum sensing.
Reviews of modern physics, 890 (3):0 035002,

[Boss et al., 2017]
Jens M Boss, KS Cujia, Jonathan Zopes, and Christian L Degen.
Quantum sensing with arbitrary frequency resolution.
Science, 3560 (6340):0 837–840, 2017.

[Chuang and Nielsen, 1997]
Isaac L. Chuang and M. A. Nielsen.
Prescription for experimental determination of the dynamics of a
quantum black box.
J. Mod. Opt., 440 (11-12):0 2455–2467, 1997.

[Altepeter et al., 2003]
J. B. Altepeter, D. Branning, E. Jeffrey, T. C. Wei, P. G. Kwiat, R. T. Thew,
J. L. O'Brien, M. A. Nielsen, and A. G. White.
Ancilla-assisted quantum process tomography.
Phys. Rev. Lett., 90:0 193601, May 2003.

[Leung, 2003]
Debbie W Leung.
Choi’s proof as a recipe for quantum process tomography.
Journal of Mathematical Physics, 440 (2):0
528–533, 2003.

[Merkel et al., 2013]
Seth T. Merkel, Jay M. Gambetta, John A. Smolin, S. Poletto, A. D.
Córcoles, B. R. Johnson, Colm A. Ryan, and M. Steffen.
Self-consistent quantum process tomography.
Phys. Rev. A, 87:0 062119, 2013.

[Rahimi-Keshari et al., 2011]
Saleh Rahimi-Keshari, Artur Scherer, Ady Mann, Ali T Rezakhani, AI Lvovsky, and
Barry C Sanders.
Quantum process tomography with coherent states.
New Journal of Physics, 130 (1):0 013006,

[Mohseni et al., 2008]
Masoud Mohseni, Ali T Rezakhani, and Daniel A Lidar.
Quantum-process tomography: Resource analysis of different
Physical Review A, 770 (3):0 032322,

[Baldwin et al., 2014]
Charles H Baldwin, Amir Kalev, and Ivan H Deutsch.
Quantum process tomography of unitary and near-unitary maps.
Physical Review A, 900 (1):0 012110, 2014.

[Preskill, 2018]
John Preskill.
Quantum Computing in the NISQ era and beyond.
Quantum, 2:0 79, August 2018.
ISSN 2521-327X.
URL <https://doi.org/10.22331/q-2018-08-06-79>.

[Bairey et al., 2019]
Eyal Bairey, Itai Arad, and Netanel H. Lindner.
Learning a local hamiltonian from local measurements.
Phys. Rev. Lett., 1220 (2):0 020504, Jan 2019.
ISSN 1079-7114.
URL <http://dx.doi.org/10.1103/PhysRevLett.122.020504>.

[Evans et al., 2019]
Tim J Evans, Robin Harper, and Steven T Flammia.
Scalable bayesian hamiltonian learning.
arXiv preprint arXiv:1912.07636, 2019.

[Anshu et al., 2021]
Anurag Anshu, Srinivasan Arunachalam, Tomotaka Kuwahara, and Mehdi
Sample-efficient learning of interacting quantum systems.
Nature Physics, 170 (8):0 931–935, 2021.

[Qi and Ranard, 2019]
Xiao-Liang Qi and Daniel Ranard.
Determining a local hamiltonian from a single eigenstate.
Quantum, 3:0 159, 2019.

[Li et al., 2020]
Zhi Li, Liujun Zou, and Timothy H. Hsieh.
Hamiltonian tomography via quantum quench.
Phys. Rev. Lett., 124:0 160502, Apr 2020.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.124.160502>.

[Zubida et al., 2021]
Assaf Zubida, Elad Yitzhaki, Netanel H Lindner, and Eyal Bairey.
Optimal short-time measurements for hamiltonian learning.
arXiv preprint arXiv:2108.08824, 2021.

[Hangleiter et al., 2021]
Dominik Hangleiter, Ingo Roth, Jens Eisert, and Pedram Roushan.
Precise hamiltonian identification of a superconducting quantum
arXiv preprint arXiv:2108.08319, 2021.

[Flammia and Wallman, 2020]
Steven T. Flammia and Joel J. Wallman.
Efficient estimation of pauli channels.
ACM Transactions on Quantum Computing, 10 (1), dec
ISSN 2643-6809.
URL <https://doi.org/10.1145/3408039>.

[Mezzacapo et al., 2014]
Antonio Mezzacapo, U Las Heras, JS Pedernales, L DiCarlo, E Solano, and
L Lamata.
Digital quantum rabi and dicke models in superconducting circuits.
Scientific reports, 40 (1):0 1–4, 2014.

[García-Álvarez et al., 2015]
L García-Álvarez, J Casanova, A Mezzacapo, IL Egusquiza, L Lamata,
G Romero, and E Solano.
Fermion-fermion scattering in quantum field theory with
superconducting circuits.
Physical review letters, 1140 (7):0 070502,

[Asaad et al., 2016]
Serwan Asaad, Christian Dickel, Nathan K Langford, Stefano Poletto, Alessandro
Bruno, Michiel Adriaan Rol, Duije Deurloo, and Leonardo DiCarlo.
Independent, extensible control of same-frequency superconducting
qubits by selective broadcasting.
npj Quantum Information, 20 (1):0 1–7, 2016.

[Weber et al., 2017]
Steven J Weber, Gabriel O Samach, David Hover, Simon Gustavsson, David K Kim,
Alexander Melville, Danna Rosenberg, Adam P Sears, Fei Yan, Jonilyn L Yoder,
et al.
Coherent coupled qubits for quantum annealing.
Physical Review Applied, 80 (1):0 014004,

[Häffner et al., 2005]
Hartmut Häffner, Wolfgang Hänsel, CF Roos, Jan Benhelm, Michael
Chwalla, Timo Körber, UD Rapol, Mark Riebe, PO Schmidt, Christoph Becher,
et al.
Scalable multiparticle entanglement of trapped ions.
Nature, 4380 (7068):0 643–646, 2005.

[Blatt and Roos, 2012]
Rainer Blatt and Christian F Roos.
Quantum simulations with trapped ions.
Nature Physics, 80 (4):0 277, 2012.

[Saffman, 2016]
Mark Saffman.
Quantum computing with atomic qubits and rydberg interactions:
progress and challenges.
Journal of Physics B: Atomic, Molecular and Optical Physics,
490 (20):0 202001, 2016.

[Cory et al., 1998]
David G Cory, Mark D Price, and Timothy F Havel.
Nuclear magnetic resonance spectroscopy: An experimentally accessible
paradigm for quantum computing.
Physica D: Nonlinear Phenomena, 1200 (1-2):0
82–101, 1998.

[Franca et al., 2022]
Daniel Stilck Franca, Liubov A Markovich, VV Dobrovitski, Albert H Werner, and
Johannes Borregaard.
Efficient and robust estimation of many-qubit hamiltonians.
arXiv preprint arXiv:2205.09567, 2022.

[Gu et al., 2022]
Andi Gu, Lukasz Cincio, and Patrick J Coles.
Practical black box hamiltonian learning.
arXiv preprint arXiv:2206.15464, 2022.

[Wilde et al., 2022]
Frederik Wilde, Augustine Kshetrimayum, Ingo Roth, Dominik Hangleiter, Ryan
Sweke, and Jens Eisert.
Scalably learning quantum many-body hamiltonians from dynamical data.
arXiv preprint arXiv:2209.14328, 2022.

[Huang et al., 2023]
Hsin-Yuan Huang, Yu Tong, Di Fang, and Yuan Su.
Learning many-body hamiltonians with heisenberg-limited scaling.
Phys. Rev. Lett., 130:0 200403, May 2023.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.130.200403>.

[Harper et al., 2021]
Robin Harper, Wenjun Yu, and Steven T. Flammia.
Fast estimation of sparse quantum noise.
PRX Quantum, 2:0 010322, Feb 2021.
URL <https://link.aps.org/doi/10.1103/PRXQuantum.2.010322>.

[Shabani et al., 2011]
A. Shabani, R. L. Kosut, M. Mohseni, H. Rabitz, M. A. Broome, M. P. Almeida,
A. Fedrizzi, and A. G. White.
Efficient measurement of quantum dynamics via compressive sensing.
Phys. Rev. Lett., 1060 (10):0 100401, Mar

[Hou et al., 2020]
Shi-Yao Hou, Ningping Cao, Sirui Lu, Yi Shen, Yiu-Tung Poon, and Bei Zeng.
Determining system hamiltonian from eigenstate measurements without
correlation functions.
New Journal of Physics, 220 (8):0 083088,

[Garrison and Grover, 2018]
James R. Garrison and Tarun Grover.
Does a single eigenstate encode the full hamiltonian?
Phys. Rev. X, 8:0 021026, Apr 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevX.8.021026>.

[Takita et al., 2017]
Maika Takita, Andrew W Cross, AD Córcoles, Jerry M Chow, and Jay M
Experimental demonstration of fault-tolerant state preparation with
superconducting qubits.
Physical review letters, 1190 (18):0 180501,

[Kandala et al., 2019]
Abhinav Kandala, Kristan Temme, Antonio D Córcoles, Antonio Mezzacapo,
Jerry M Chow, and Jay M Gambetta.
Error mitigation extends the computational reach of a noisy quantum
Nature, 5670 (7749):0 491, 2019.

[Sun et al., 2021]
Jinzhao Sun, Xiao Yuan, Takahiro Tsunoda, Vlatko Vedral, Simon C. Benjamin, and
Suguru Endo.
Mitigating realistic noise in practical noisy intermediate-scale
quantum devices.
Phys. Rev. Applied, 15:0 034026, Mar
URL <https://link.aps.org/doi/10.1103/PhysRevApplied.15.034026>.

[Kempe et al., 2006]
Julia Kempe, Alexei Kitaev, and Oded Regev.
The complexity of the local hamiltonian problem.
Siam journal on computing, 350 (5):0
1070–1097, 2006.

[Childs et al., 2018]
Andrew M Childs, Dmitri Maslov, Yunseong Nam, Neil J Ross, and Yuan Su.
Toward the first quantum simulation with quantum speedup.
Proceedings of the National Academy of Sciences, 1150
(38):0 9456–9461, 2018.

[Low and Chuang, 2019]
Guang Hao Low and Isaac L Chuang.
Hamiltonian simulation by qubitization.
Quantum, 3:0 163, 2019.

[Sun et al., 2021]
Jinzhao Sun, Suguru Endo, Huiping Lin, Patrick Hayden, Vlatko Vedral, and Xiao
Perturbative quantum simulation.
arXiv preprint arXiv:2106.05938, 2021b.

[Helsen et al., 2022]
J. Helsen, I. Roth, E. Onorati, A.H. Werner, and J. Eisert.
General framework for randomized benchmarking.
PRX Quantum, 3:0 020357, Jun 2022.

[Helsen et al., 2019]
Jonas Helsen, Joel J Wallman, Steven T Flammia, and Stephanie Wehner.
Multiqubit randomized benchmarking using few samples.
Physical Review A, 1000 (3):0 032304,

[Wallman, 2018]
Joel Wallman.
Randomized benchmarking with gate-dependent noise.
Quantum, 2:0 47, 2018.

[Harper et al., 2020]
Robin Harper, Steven T Flammia, and Joel J Wallman.
Efficient learning of quantum noise.
Nature Physics, 160 (12):0 1184–1188, 2020.

[Barends et al., 2014]
R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White,
J. Mutus, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro,
A. Dunsworth, C. Neill, P. O'Malley, P. Roushan, A. Vainsencher, J. Wenner,
A. N. Korotkov, A. N. Cleland, and John M. Martinis.
Superconducting quantum circuits at the surface code threshold for
fault tolerance.
Nature, 5080 (7497):0 500–503, apr 2014.

[Somoroff et al., 2021]
Aaron Somoroff, Quentin Ficheux, Raymond A Mencia, Haonan Xiong, Roman V
Kuzmin, and Vladimir E Manucharyan.
Millisecond coherence in a superconducting qubit.
arXiv preprint arXiv:2103.08578, 2021.

[Magesan, 2012]
Easwar Magesan.
Characterizing Noise in Quantum Systems.
PhD thesis, University of Waterloo, Waterloo, Ontario, Canada, 2012.

[Dankert et al., 2009]
Christoph Dankert, Richard Cleve, Joseph Emerson, and Etera Livine.
Exact and approximate unitary 2-designs and their application to
fidelity estimation.
Phys. Rev. A, 80:0 012304, Jul 2009.

[Erhard et al., 2019]
Alexander Erhard, Joel J. Wallman, Lukas Postler, Michael Meth, Roman Stricker,
Esteban A. Martinez, Philipp Schindler, Thomas Monz, Joseph Emerson, and
Rainer Blatt.
Characterizing large-scale quantum computers via cycle benchmarking.
Nat. Commun., 100 (1):0 5347, November 2019.
URL <https://doi.org/10.1038/s41467-019-13068-7>.

[Helsen et al., 2019]
Jonas Helsen, Xiao Xue, Lieven MK Vandersypen, and Stephanie Wehner.
A new class of efficient randomized benchmarking protocols.
npj Quantum Inf., 5:0 71, 2019b.

[Mohseni et al., 2008]
M. Mohseni, A. T. Rezakhani, and D. A. Lidar.
Quantum-process tomography: Resource analysis of different
Phys. Rev. A, 77:0 032322, Mar 2008b.
URL <https://link.aps.org/doi/10.1103/PhysRevA.77.032322>.

[Poyatos et al., 1997]
J. F. Poyatos, J. I. Cirac, and P. Zoller.
Complete characterization of a quantum process: The two-bit quantum
Phys. Rev. Lett., 78:0 390–393, Jan 1997.
URL <https://link.aps.org/doi/10.1103/PhysRevLett.78.390>.

[Yuan et al., 2016]
Xiao Yuan, Zhen Zhang, Norbert Lütkenhaus, and Xiongfeng Ma.
Simulating single photons with realistic photon sources.
Physical Review A, 940 (6):0 062305, 2016.

[Troiani et al., 2005]
Filippo Troiani, Alberto Ghirri, Marco Affronte, S Carretta, P Santini,
G Amoretti, S Piligkos, G Timco, and REP Winpenny.
Molecular engineering of antiferromagnetic rings for quantum
Physical review letters, 940 (20):0 207208,

[Santini et al., 2011]
P Santini, S Carretta, F Troiani, and G Amoretti.
Molecular nanomagnets as quantum simulators.
Physical review letters, 1070 (23):0 230502,

[Han and Yu, 2021]
Zeyao Han and Wenjun Yu.
Github: HamiltonianLearning.
<https://github.com/zyHan2077/HamiltonianLearning>, December

[Chen et al., 2022]
Senrui Chen, Sisi Zhou, Alireza Seif, and Liang Jiang.
Quantum advantages for pauli channel estimation.
Phys. Rev. A, 105:0 032435, Mar 2022.
URL <https://link.aps.org/doi/10.1103/PhysRevA.105.032435>.

[Li et al., 2014]
Xiao Li, Joseph Kurata Bradley, Sameer Pawar, and Kannan Ramchandran.
The SPRIGHT algorithm for robust sparse Hadamard transforms.
In 2014 IEEE International Symposium on Information Theory,
Honolulu, Hawaii, USA, June 2014. IEEE.

[Candes and Tao, 2005]
E. J. Candes and T. Tao.
Decoding by linear programming.
IEEE Transactions on Information Theory, 510
(12):0 4203–4215, 2005.

[Candes et al., 2006]
Emmanuel J Candes, Justin K Romberg, and Terence Tao.
Stable signal recovery from incomplete and inaccurate measurements.
Communications on Pure and Applied Mathematics: A Journal
Issued by the Courant Institute of Mathematical Sciences, 590
(8):0 1207–1223, 2006.

[Baraniuk et al., 2008]
Richard Baraniuk, Mark Davenport, Ronald DeVore, and Michael Wakin.
A simple proof of the restricted isometry property for random
Constructive Approximation, 280 (3):0
253–263, 2008.

[Candès, 2008]
Emmanuel J. Candès.
The restricted isometry property and its implications for compressed
Comptes Rendus Mathematique, 3460 (9):0
589–592, 2008.
ISSN 1631-073X.

§ PRELIMINARIES

§.§ Pauli Concept

In this section, we will introduce the commonly-used Pauli group and the label to represent every element.
We will start from the single-qubit case.

For the single-qubit Pauli group $\mathbb{P}$, it is generated by three Pauli operators,

\begin{gather}
\mathbb{P}=\langle X,Y,Z\rangle.
\end{gather}

The multiplication among these generators will cause some extra phase terms for some operators, so this group contains several elements that are not Hermitian.
On the other hand, all the Pauli elements serve either as observables or gates.
Therefore, we only take the refined Pauli group $\sf P$ with all Hermitian elements into account.
In group $\sf{P}$, we treat elements with only differences on phase the same, which is equivalent to making a quotient over the phase element

\begin{gather}
\sf{P}=\faktor{\mathbb{P}}{\langle iI\rangle}.
\end{gather}

Moreover, for every integer $n$, the refined $n$-qubit Pauli group $\sf{P}^n$ is constructed from tensor products of $\sf{P}$, and we have

\begin{gather}
\sf{P}^n=\bigotimes_{i=1}^n\sf{P}.
\end{gather}

In order to specify every element among $4^n$ in the refined $n$-qubit Pauli group, a canonical way will be using a $2n$-bit string to label every Pauli operator in the refined group.
We can start to illustrate this convention from the case of the single-qubit Pauli group $\sf P$.
This refined group contains four elements which are the three Pauli matrices and an identity operator.
As for the label, we refer to a bit string $\alpha$ with a length of two,

$I$      $X$      $Y$      $Z$

$\alpha$     00     01     10     11

Therefore, if we want to label an $n$-qubit Pauli element, we recruit $2n$-bit strings to record the Pauli element in each qubit.
For a $2n$-bit label $\alpha$, the $2i-1^{\text{th}}$ and $2i^{\text{th}}$ bits, $\alpha_{2i-1}\alpha_{2i}$, characterize the corresponding single Pauli operator of $P_\alpha$ on the $i$th qubit.
For example, if we want to represent the operator $XZIX$, then its label must be $\alpha=01110001$.

Since there exists an isomorphism between the $2n$-bit labels and the Pauli elements, we need to define the corresponding operations that use the label to express some properties of Pauli operators.
The basic operation is the addition.
We define the addition of Pauli labels following the multiplication of Pauli operators.
Regardless of phases, the summation is equal to the multiplication of every pair of Pauli operators so that

\begin{gather}
P_\alpha\cdot P_\beta=P_{\alpha+\beta},\ \ \forall\,P_\alpha,P_\beta\in \sf{P}^n.
\end{gather}

Note that we will use $\alpha\in\sf{P}^n$ for short of denoting $P_\alpha\in\sf{P}^n$ when it would not incur any confusion.
We next introduce a new inner product between labels of two Pauli operators.
For two $n$-qubit Pauli labels $\alpha$ and $\beta$, we define the Pauli inner product as follows

\begin{gather}
\langle\alpha,\beta\rangle_p\coloneqq\sum_{i=1}^n(\alpha_{2i-1}\cdot\beta_{2i}+\beta_{2i-1}\cdot\alpha_{2i})\mod{2},
\end{gather}

where we use the subscript to distinguish the Pauli inner product from the ordinary inner product.
This operation actually describes the commutation of every pair of Pauli operators.
By calculating the qubitwise commutation and accumulating them, the resulting Boolean value indicates the commutation property as

\begin{align}
[P_\alpha,P_\beta]=\left(1-(-1)^{\langle\alpha,\beta\rangle_p}\right)P_\alpha P_\beta.
\end{align}

§.§ Pauli Stabilizer Group

There exists a kind of special subgroup in the refined Pauli group, namely the stabilizer group.

A stabilizer group $S$ is defined to be a subgroup of the refined Pauli group ${\sf P}^n$ such that $\langle a,b\rangle_p=0$ for all $P_a,P_b\in S$.

It can be directly inferred from the definition that each stabilizer group is a normal subgroup of the refined Pauli group.
Therefore, we can find further the quotient group $A_S\coloneqq \faktor{{\sf P}^n}{S}$.
By group theory, we know that for an arbitrary $P_c\in{\sf P}^n$, there exists a unique pair of $e_c\in A_S$ and $s_c\in S$ so that $P_c=e_c\cdot s_c$.
Since only $e_c$ contributes to the syndrome measurement after implementing $P_c$ to a stabilizer state, we denote $e_c$ by the error syndrome of $P_c$.

It is easy to prove that the maximum size of a stabilizer group is $2^n$ for ${\sf P}^n$.
Therefore, the size of the quotient group of a maximum stabilizer group is $2^n$.
Consider a maximum stabilizer group generated by $n$ generator $G\coloneqq\langle g_1,\cdots,g_n\rangle$, and for any Pauli operator $P_c=e_c\cdot s_c$, we can get the full commuting relation between $P_c$ and $G$ by the commuting relation between $P_c$ and set $\{g_1,\cdots,g_n\}$, which a $n$-bit string can represent.
Therefore, each element in $A_S$ maps onto a certain $n$-bit string, and this serves the basics of the stabilizer algebra.
Consequently, if we implement a simultaneous measurement of all the generators referred to as the syndrome measurement, the results turn out to be the bit string representation of the error syndrome of the equivalent Pauli gate to the whole circuit.

During the analysis of the fidelity estimator in Sec. <ref>, the protocol tracks and adds the syndromes from sampled Pauli gates and the final measurement.
The summation is processed by the corresponding bit strings by binary summation.
And the result is the overall syndrome of the circuit.

§.§ Pauli Channel Characterizations

The commonly used Pauli channel can be regarded as a stochastic quantum channel.
The is due to the definition of Pauli channels,

\begin{gather}\label{eq:Paulichannel}
\mathcal{E}^{\sf P}(\cdot)\coloneqq\sum_{\alpha\in{\sf P}^n}p_\alpha P_\alpha(\cdot)P_\alpha,
\end{gather}

where coefficients $\{p_\alpha\}$ are Pauli error rates.
By the CPTP conditions for every physical channel, we have the following constraints over these Pauli error rates,

\begin{gather}
\sum_{\alpha\in{\sf P}^n}p_\alpha=1,\ \ p_\alpha\geq0\ \forall\, \alpha\in{\sf P}^n.
\end{gather}

Therefore, the set $\{p_\alpha\}$ constitutes a probability distribution over $4^n$ possible Pauli operations.
Consequently, a Pauli channel with the form as Eq. (<ref>) can be interpreted as applying a Pauli gate $P_\alpha(\cdot)P_\alpha$ with probability $p_\alpha$.

Moreover, there is a second set of $4^n$ parameters that can fully characterize a Pauli channel $\mathcal{E}^{\sf P}$, which is defined in <Ref> as Pauli fidelity,

\begin{gather}\label{eq:Paulifidelity2}
f_\alpha\coloneqq\frac{1}{2^n}\Tr(P_\alpha\mathcal{E}^{\sf P}(P_\alpha)),\ \ \forall\,\alpha\in{\sf P}^n.
\end{gather}

Note we can also deduce this value from the definition of Pauli channels.
From Eq. (<ref>), every Pauli operator is an eigenoperator of an arbitrary Pauli channel, while the eigenvalue of $P_\alpha$ is just $f_\alpha$.
Thereby, some also refer values defined in Eq. (<ref>) as Pauli eigenvalues.

§ PAULI PROPERTIES ESTIMATION

In this work, we estimate an arbitrary Hamiltonian operator with a sparse decomposition vector efficiently.
As stated in <Ref>, the decomposition vector $\bf s$ can faithfully describe a Hamiltonian operator.
Therefore, the problem of Hamiltonian estimation is, in principle, a scalable task.
According to Eq. (<ref>), the Hamiltonian channel contains the full information of the vector.
We choose the Hamiltonian evolution as the target, so the question is deduced to estimate a Hamiltonian channel.
In this section, we introduce the fidelity estimator in Appendix <ref>. The estimator applies to the reference circuits and the circuit with the Hamiltonian evolution sequentially to extract the net Pauli fidelity of the Hamiltonian channel.
Moreover, we recruit the sparse Pauli error rates estimator in Appendix <ref> to detect the error rates of the Hamiltonian channel, which are directly linked with the decomposition parameters.

§.§ Pauli Fidelity Estimation

Here we first introduce the procedure and theoretical guarantees presented in [Flammia and Wallman, 2020], and it will be recruited as a subroutine for the estimation of the Hamiltonian channel's Pauli fidelity.
This algorithm was first proposed to estimate the implementation quality of all Pauli gates.
To quantify this quality, we regard the practical implementation of a Pauli gate $P$ as $\tilde{P}=\Lambda\circ P$.
Therefore, Algorithm <ref> estimates the Pauli fidelity of this noise channel $\Lambda$ over a set of Pauli indices, $\sf X$, which belongs to a stabilizer subgroup $\sf G$ in the whole Pauli group $\sf{P}^n$.

For the implementation, we choose the input as a stabilizer state of stabilizer subgroup $\sf{G}$, and the final measurement will be a syndrome measurement of the group $\sf{G}$.
Note both the input states and the measurements contain some noise.
Hence the protocol must be noise-resilient of SPAM errors.
The execution starts with a series of random Pauli gates.
Intuitively, it first sets the length index $m$ as 0 as a reference, and $m$ increases exponentially from 1.
We denote this sequence of $m$ by $\kappa$.
The procedure then applies $m+1$ layers of uniformly random Pauli gates to the input state for each $m\in\kappa$.
Since the randomly twirled channel will be a diagonal channel in the superoperator representation, the expected channel of this circuit equals the composition of $m$ identically diagonal channels with some SPAM error effects, where the diagonal channel is equal to the diagonal part of the noise channel, $\Lambda$, namely, the Pauli fidelity terms of $\Lambda$.
By employing the stabilizer states and syndrome measurements, the results bring some diagonal information about the noise channel tarnished by SPAM errors.
More specifically, the results will be an exponential decay with the length index $m$ with factors to be the diagonal terms.
When a certain measurement result decays down to one-third of the reference, the algorithm then collects all results along $m\in\kappa$ to fit and get the diagonal information.
The Algorithm <ref> is shown in the following, where $l$ denotes the number of random sequences for every  length needed to construct the random estimator.

FEstimator(${\sf G,X},l$)
[1]
Input: a target set $\sf X$ belonging to a stabilizer group $\sf{G}$, and a positive integer $l$
Initialize an array $\hat{r}_\Lambda$ with size $|\sf X|$ to be all -1; Initialize an array $\bm{V}$ with size $|\sf X|$ to be all 0
Input states will be positive stabilizer states $\tilde{\rho}_{\sf G}$; The measurement will be syndrome measurement $M_{\sf G}$ of $\sf G$
Apply a random Pauli gate $P_{\alpha}\in\sf{P}^n$ to state $\tilde{\rho}_{\sf{G}}$, and perform measurement $M_{\sf{G}}$ with result $b$
$x\in \sf{X}$
$\bm{V}[x]+=(-1)^{\langle x,\alpha+b\rangle_p}$
$\hat{V} \leftarrow \frac{1}{l}\bm{V}$; $\bm{V}\leftarrow \bm{0}$; $m \leftarrow 1$
$\exists x \in \sf{X}$ such that $\hat{r}_\Lambda[x]= -1$
Apply $m+1$ random Pauli gates $\{P_{\alpha_0},\cdots,P_{\alpha_m}\}$ to state $\tilde{\rho}_{\sf G}$, and perform measurement $M_{\sf G}$ with result $b$
$x\in \sf{X}$
$\bm{V}[x]+=(-1)^{\langle x,\sum_{i=0}^m \alpha_i+b\rangle_p}$
$\hat{W} \leftarrow \frac{1}{l}\bm{V}$; $\bm{V}\leftarrow \bm{0}$
$x\in \sf{X}$ such that $\hat{r}_\Lambda[x]=-1$
$v \leftarrow \hat{V}[x]$ and $w \leftarrow \hat{W}[x]$
If $w \leq v / 3$ and $w,v>0$, $\hat{r}_\Lambda[x] \leftarrow 1 - \bigl(w/v\bigr)^{1/m}$
Else If $w\le0$ or $v\le0$, $\hat{r}_\Lambda[x] \leftarrow 1$
the array $\hat{r}_\Lambda $

Even though the above algorithm works for a subset of a stabilizer group, we can simply generalize this method to detect a general subset $\sf{X}$ in the Pauli group.
To achieve that, we employ the concept of stabilizer covering.
For a set $\sf{X}$ of Pauli operators, its stabilizer covering is a collection of stabilizer groups, of which the union set covers $\sf{X}$.
Therefore, we can divide the target set exclusively into multiple stabilizer groups and execute the algorithm for each piece to estimate the whole set.

Let $\sf G$ be a stabilizer group in ${\sf P}^n$ and ${\sf X}\subset{\sf G}$.
Suppose the noise $\Lambda$ satisfies the assumption <ref> and $\|\Lambda-\mathcal{I}\|<\frac{1}{2}$, and the following holds with probability $1-\delta$ by choosing small $\epsilon,\delta>0$:
Running ${\bf FEstimator}({\sf G,X},l)$ with $l = \frac{2}{\epsilon^2} \log\bigl(\frac{2|\sf{X}||\kappa|}{\delta}\bigr)$ where $\kappa$ is the set of exponentially increasing sequence lengths.
The procedure uses at most $O\left(\log{\frac{1}{\Delta_{\sf{X}}}}\right)$ sequence lengths with largest length $m$ to be $O\bigl( \tfrac{1}{\Delta_{\sf{X}}}\bigr)$, where $\Delta_{\sf{X}}$ denotes the smallest residual $r_\Lambda\coloneqq1-f_\Lambda$ among the set $\sf{X}$ regarding Pauli fidelity $f_\Lambda$ of noise channel $\Lambda$.
Moreover, the output $\hat{r}_\Lambda$ satisfies

\begin{align*}
|\hat{r}_{\Lambda,x} - r_{\Lambda,x}| \leq O(\epsilon) r_{\Lambda,x}\ \ \forall x\in\sf{X}.
\end{align*}

Beyond this standard Pauli fidelity estimation for the implementation noise of Pauli gates, we want further to detect the fidelity information of the Hamiltonian evolution.
The idea of this detection is similar to that of the standard procedure for the noise channel.
We will use random Pauli gates to twirl the desired Hamiltonian channel.

For convenience, we denote the extended algorithm by $\textbf{HFEstimator}({\sf G,X},l,t)$ where the additional parameter $t$ represents the evolving time of the Hamiltonian channel.
The major difference between this protocol and the original one is that the new procedure will estimate the composite channel of the implementation noise and the Hamiltonian evolution, which is $\mathcal{H}_t\circ \Lambda$.
After implementing a random Pauli gate, the new procedure appends the Hamiltonian channel $\mathcal{H}_t$ to the circuit, as shown in Algorithm <ref>.
In summary, the modified procedure will implement $m+1$ random gates along with the channel $\mathcal{H}_t$, alternatively.
The circuit is illustrated in Figure <ref>(b).

HFEstimator(${\sf G,X},l,t$)
[1]
Input: a target set $\sf X$ belonging to a stabilizer group $\sf{G}$, a positive integer $l$, and the evolving time $t$
Initialize an array $\hat{r}_{\mathcal{H}\circ\Lambda}$ with size $|\sf X|$ to be all -1; Initialize an array $\bm{V}$ with size $|\sf X|$ to be all 0
Input states will be positive stabilizer states $\tilde{\rho}_{\sf G}$; The measurement will be syndrome measurement $M_{\sf G}$ of $\sf G$
Apply a random Pauli gate $P_{\alpha}\in\sf{P}^n$ and Hamiltonian evolution $\mathcal{H}_t$ to state $\tilde{\rho}_{\sf{G}}$; Perform measurement $M_{\sf{G}}$ with result $b$
$x\in \sf{X}$
$\bm{V}[x]+=(-1)^{\langle x,\alpha+b\rangle_p}$
$\hat{V} \leftarrow \frac{1}{l}\bm{V}$; $\bm{V}\leftarrow \bm{0}$; $m \leftarrow 1$
$\exists x \in \sf{X}$ such that $\hat{r}_{\mathcal{H}\circ\Lambda}[x]= -1$
Apply $m+1$ random Pauli gates $\{P_{\alpha_0},\cdots,P_{\alpha_m}\}$ interleaved with $\mathcal{H}_t$ to $\tilde{\rho}_{\sf G}$; Perform measurement $M_{\sf G}$ with result $b$
$x\in \sf{X}$
$\bm{V}[x]+=(-1)^{\langle x,\sum_{i=0}^m \alpha_i+b\rangle_p}$
$\hat{W} \leftarrow \frac{1}{l}\bm{V}$; $\bm{V}\leftarrow \bm{0}$
$x\in \sf{X}$ such that $\hat{r}_{\mathcal{H}\circ\Lambda}[x]=-1$
$v \leftarrow \hat{V}[x]$ and $w \leftarrow \hat{W}[x]$
If $w \leq v / 3$ and $w,v>0$, $\hat{r}_{\mathcal{H}\circ\Lambda}[x] \leftarrow 1 - \bigl(w/v\bigr)^{1/m}$
Else If $w\le0$ or $v\le0$, $\hat{r}_{\mathcal{H}\circ\Lambda}[x] \leftarrow 1$
the array $\hat{r}_{\mathcal{H}\circ\Lambda} $

Similar to Algorithm <ref>, this algorithm actually replaces the noise channel $\Lambda$ by the composite channel $\mathcal{H}_t\circ\Lambda$.
Therefore, the resulting estimation from Algorithm <ref> keeps very similar guarantees to Proposition <ref>.

Let $\sf{X}\subseteq\sf{P}^n$ and $\sf G$ be a stabilizer group containing $\sf X$.
Suppose the noise $\Lambda$ satisfies assumption <ref> and assumption <ref>, and the following holds with probability $1-\delta$ by choosing small $\epsilon,\delta>0$:
Running ${\bf HFEstimator}({\sf G,X},l,t)$ with $t$ satisfies $\|\mathcal{H}_t-\mathcal{I}\|<\frac{1}{4}$ and $l = \frac{2}{\epsilon^2} \log\bigl(\frac{2|\sf{X}||\kappa|}{\delta}\bigr)$ where $\kappa$ is the set of exponentially increasing sequence lengths. The procedure uses at most $O\left(\log{\frac{1}{\Delta_{\sf{X}}}}\right)$ sequence lengths with largest length $m$ to be $O\bigl( \tfrac{1}{\Delta_{\sf{X}}}\bigr)$, where $\Delta_{\sf{X}}$ denotes the smallest residual $r_{\mathcal{H}\circ\Lambda}\coloneqq1-f_{\mathcal{H}\circ\Lambda}$ among the set $\sf{X}$ regarding Pauli fidelity $f_{\mathcal{H}\circ\Lambda}$ of the composite channel $\mathcal{H}_t\circ\Lambda$.
Moreover, the output $\hat{r}_{\mathcal{H}\circ\Lambda}$ satisfies

\begin{align*}
|\hat{r}_{\mathcal{H}\circ\Lambda,x} - r_{\mathcal{H}\circ\Lambda,x}| \leq O(\epsilon) r_{\mathcal{H}\circ\Lambda,x}\ \ \forall x\in\sf{X}.
\end{align*}

Since the Hamiltonian channel is defined to be $\mathcal{H}_t(\rho)=e^{-iHt}\rho e^{iHt}$, the channel is close to the identity channel when the evolving time $t$ is small.
With the assumption stated above, the composite channel $\mathcal{H}_t\circ\Lambda$ is close to the identity under the Pauli twirling.
Specifically, we have

\begin{align}\label{eq:closetoi}
\|(\mathcal{H}_t\circ\Lambda)^{{\sf P}^n}-\mathcal{I}\|=&\max_{\alpha\in{\sf P}^n}\abs{1-f_{\mathcal{H}_t,\alpha}f_{\Lambda,\alpha}}\notag\\
\leq&1-\min_{\alpha\in{\sf P}^n}f_{\mathcal{H}_t,\alpha}\min_{\beta\in{\sf P}^n}f_{\Lambda,\beta}\notag\\
\end{align}

where the first inequality comes from the fact that all Pauli fidelity is upper bounded by 1, and the last inequality comes from the requirement of $t$ and the assumption <ref>.
Therefore, according to the Proposition <ref>, HFEstimator can estimate the composite channel precisely as stated.

With these two estimators, the full detection of the Hamiltonian's Pauli fidelity can be executed as in Algorithm <ref>.
The procedure first employs FEstimator to estimate information on the implementation noise of Pauli gates.
We call this detection the reference circuit and denote the detected fidelity by $f_\Lambda$.
Since we will also implement Pauli gates in the following circuits, this reference knowledge helps to rule out the effect of this noise.
In the second step, the procedure invokes HFEstimator to estimate circuits with the unknown Hamiltonian channel.
By definition, this subroutine returns the fidelity $f_{\mathcal{H}\circ\Lambda}$ about the composite of implementation noise and Hamiltonian channels.
By Assumption <ref>, the implementation noise is the Pauli channel,  or namely, it is a diagonal matrix.
Thus, with a direct ratio estimator, the procedure eliminates the effects of noise and keeps the net estimation of Pauli fidelity of the Hamiltonian channel.
We denote this fidelity by $f_{\mathcal{H}_t}$ with residual $r_{\mathcal{H}_t}$.

Let $\sf{X}\subseteq\sf{P}^n$ and $\sf G$ be a stabilizer group containing $\sf X$.
For any small $\epsilon,\delta>0$ and noise under Assumption <ref> and <ref>, the following holds with probability $1-\delta$:
Run Algorithm <ref> with $t_0$ satisfies $\|\mathcal{H}_{5t_0}-\mathcal{I}\|<\frac{1}{4}$ and $l = \frac{2}{\epsilon^2} \log\bigl(\frac{2|\kappa|}{\delta}\bigr)$ where $\kappa$ is the set of variant sequence lengths.
The ratio estimator can estimate the Pauli fidelity information of the desired $\mathcal{H}_t$ satisfying that

\begin{gather*}
\abs{\hat{f}_{\mathcal{H}_t,x}-f_{\mathcal{H}_t,x}}\leq\frac{O(\epsilon)(f_{\mathcal{H}_t,x}r_{\Lambda,x}+r_{\mathcal{H}_t\circ\Lambda,x})}{f_{\Lambda,x}-O(\epsilon)r_{\Lambda,x}},\ \forall\,x\in\sf{X}.
\end{gather*}

The guarantees from Proposition <ref> and Corollary <ref> can help to provide this accuracy bound.
In the first round, the procedure estimates the fidelity information of the native noise $\Lambda$ by the stated complexity up to an accuracy

\begin{gather*}
\abs{\hat{f}_{\Lambda,x}-f_{\Lambda,x}}\leq O(\epsilon) r_{\Lambda,x}.
\end{gather*}

Similarly, in the second round, the procedure executes the circuit including the composite $\mathcal{H}_t\circ\Lambda$ and returns another accuracy with a very similar complexity

\begin{gather*}
\abs{\hat{f}_{\mathcal{H}_t\circ\Lambda,x}-f_{\mathcal{H}_t\circ\Lambda,x}}\leq O(\epsilon) r_{\mathcal{H}_t\circ\Lambda,x}.
\end{gather*}

According to Assumption <ref>, we can show the fidelity is $f_{\mathcal{H}_t\circ\Lambda,x}=f_{\Lambda,x}f_{\mathcal{H}_t,x}$.
As $\Lambda$ is itself a Pauli channel, we have

\begin{align}
f_{\mathcal{H}_t\circ\Lambda,x}=&\frac{1}{2^n}\Tr(P_i(\mathcal{H}_t\circ\Lambda)^{\sf P}(P_i))=\frac{1}{2^n}\Tr(P_i\mathcal{H}_t\circ\Lambda(P_i))\notag\\
=&\frac{f_{\Lambda,x}}{2^n}\Tr(P_i\mathcal{H}_t(P_i))=\frac{f_{\Lambda,x}}{2^n}\Tr(P_i\mathcal{H}_t^{\sf P}(P_i))\notag\\
\end{align}

Thus we can calculate the estimation by the ratio estimator.
The upper and lower bounds of this estimation are therefore calculated as follows,

\begin{align}
\frac{\hat{f}_{\mathcal{H}_t\circ\Lambda,x}}{\hat{f}_{\Lambda,x}}-f_{\mathcal{H}_t,x}&\leq\frac{f_{\Lambda,x}f_{\mathcal{H}_t,x}+O(\epsilon) r_{\mathcal{H}_t\circ\Lambda,x}}{f_{\Lambda,x}-O(\epsilon) r_{\Lambda,x}}-f_{\mathcal{H}_t,x}\notag\\
\end{align}

For the lower bound, we have a similar way of quantifying that

\begin{align}
f_{\mathcal{H}_t,x}-\frac{\hat{f}_{\mathcal{H}_t\circ\Lambda,x}}{\hat{f}_{\Lambda,x}}&\leq f_{\mathcal{H}_t,x}-\frac{f_{\Lambda,x}f_{\mathcal{H}_t,x}-O(\epsilon) r_{\mathcal{H}_t\circ\Lambda,x}}{f_{\Lambda,x}+O(\epsilon) r_{\Lambda,x}}\notag\\
\end{align}

In the statement, we claim the errors of estimations will be bounded by a multiplicative accuracy.
Although there are different components involved in the relative part, this is still precise since the residual is small and all the fidelity terms are less than 1.
As we have introduced in the main text, the local gate noise is usually very close to the identity channel [Harper et al., 2020, Barends et al., 2014, Somoroff et al., 2021], which means $f_{\Lambda}$ is close to 1, and $r_{\Lambda}$ is negligible.
Likewise, the Hamiltonian channels are assigned with short evolving times in our learning method.
Thus, these channels are also well-behaved with large $f_{\mathcal{H}_t}$ and vanishing $r_{\mathcal{H}_t\circ\Lambda}$.
Therefore, our estimation will be of a multiplicative precision, and this error is much less than the predetermined $\epsilon$.

In this analysis, we do not rule out the effects of evolving times, so this bound of estimation may vary from different Hamiltonian and evolving times.
To resolve this dependence, the procedure FE increases the evolving times evenly among multiple calls of HFEstimator.
With a regression over these multiple $t$, the procedure can estimate the second-order fidelity $f^{(2)}$ as stated in Eq. (<ref>).
With the bounded noise proved in Lemma <ref>, the regression parameter can be estimated with high accuracy.
Let $\sf{X}\subseteq\sf{P}^n$ and $\sf G$ be a stabilizer group containing $\sf X$, and suppose Assumption <ref> and <ref> hold.
For any small $\epsilon,\delta>0$, the following happens with probability $1-\delta$:
Run Algorithm <ref> $\textbf{FE}({\sf X,G},l,t_0)$ with $t_0$ satisfies $\|\mathcal{H}_{5t_0}-\mathcal{I}\|<\frac{1}{4}$ and  $l = \frac{2}{\epsilon^2} \log\bigl(\frac{2|\kappa|}{\delta}\bigr)$ where $\kappa$ is the set of variant sequence lengths. regression, and the resulting regression array $\hat{f}^{(2)}$ satisfies

\begin{gather*}
\abs{\hat{f}^{(2)}_x-f^{(2)}_x}\leq\frac{\sigma O(\epsilon)(f_{\mathcal{H}_t,x}r_{\Lambda,x}+r_{\mathcal{H}_t\circ\Lambda,x})}{\left(f_{\Lambda,x}-O(\epsilon)r_{\Lambda,x}\right)t_0^2}+o(t_0^2),
\end{gather*}

where the upper bound term adopts the largest bound term referred from Lemma <ref> among multiple evolution with different $t$, and $\sigma$ is a constant regrading the regression method.

According to the ordinary least squares, the parameters $\hat{f^{(2)}}$ can be calculated from the noisy observations and ideal predictors,

\begin{gather}\label{eq:OLS}
\hat{f}^{(2)}_x\coloneqq\frac{\sum_{j=1}^k(t_j^2-\bar{t^2})\hat{f}_{\mathcal{H}_{t_j},x}}{\sum_{j=1}^k(t_j^2-\bar{t^2})^2},
\end{gather}

where $k$ denotes the number of observations in the regression.
According to Lemma <ref>, the errors of estimation of $f_{\mathcal{H}_t}$ are also bounded by the predetermined precision $\epsilon$ multiplying the size of residual.
In Eq. (<ref>), the estimator is a linear summation.
Thus we choose the largest bias among these fidelity terms with different evolving times $t$ and use it to replace all other bounds,

\begin{gather}
\abs{\hat{f}^{(2)}_x-f^{(2)}_x}\leq\frac{\sum_{j=1}^k\abs{t_j^2-\bar{t^2}}}{\sum_{j=1}^k(t_j^2-\bar{t^2})^2}\cdot\max_j\{\Delta f_{\mathcal{H}_{t_j},x}\}+o(t^2).
\end{gather}