Note that the higher-order errors come from the systematic errors in the estimation of $f^{(2)}$ as we  only consider the second-order terms of Pauli fidelity.
Furthermore, $o(t^2)$ is due to the observation that only even number order terms contribute to the fidelity.
In the execution of FE, we implement multiple evolutions with $t$ varies from $t_0$ up to $5t_0$ evenly.
Therefore, we have the following constant

\begin{gather}\label{eq:linearreg}
\frac{\sigma}{t_0^2}=\frac{\sum_{i=1}^5\abs{t_i^2-\bar{t^2}}}{\sum_{i=1}^5(t_i^2-\bar{t^2})^2}\approx\frac{1}{10t_0^2}.
\end{gather}

Combining the above two equations, we can get the claimed bound.

§.§ Sparse Pauli Error Rates

In this section, we focus on the post-processing of the fidelity information detected in the last part.
Even though the whole set of fidelity can fully characterize an unknown Pauli channel, it costs an exponential-size calculation to transform the fidelity to Pauli error rates, as stated in Eq. (<ref>).
These error rates imply the absolute values of our desired parameters according to Eq. (<ref>).
Fortunately, the proposal in [Harper et al., 2021] offers an efficient method to complete this transformation given the target channel has sparse nonzero Pauli error rates.
We will illustrate this method with some reorganized algorithms and guarantees.

The protocol first constructs bins to collect the linear combinations of certain fidelity terms.

Subsampling and Bin Construction: $\textbf{BC}(\{\mathbf{d}_{c;t}\},l,b,C,P,t_0)$
[1]
Number of sampling groups $C$, number of observations P, offsets $\{\mathbf{d}_{c;t}\}_{CP}$, positive integers $l$ and $b$, and the unit time length $t_0$  a $4^n$ array $\hat{f}^{(2)}$ with all -1
Random subsampling matrices $\mathbf{M}_c\in\mathbb{F}_2^{2n\times b}$ for $c \in [C]$.
$\mathbf{M}_c'\leftarrow J_n \mathbf{M}_c$ $\forall\,c\in[C]$.
$c\in[C]$, $t\in[P]$, and $\ell\in\mathbb{F}_2^{b}$
$k \leftarrow \mathbf{M}_c' \ell+\mathbf{d}_{c;t}$
Query Algorithm <ref> with $(G_k,\{P_k\},l,t_0)$
$B \leftarrow 2^{b}$
$c\in[C]$ and $t\in[P]$

$U_{c;t}[j] \leftarrow\frac{1}{B} \sum_{\ell\in\mathbb{F}_2^{b}} (-1)^{\langle j,\ell\rangle}\hat{f}^{(2)}_{\mathbf{M}_c' \ell+\mathbf{d}_{c;t}}$
Return $U_{c;t}[j]$

As shown in Algorithm <ref>, the procedure chooses fidelity according to the random indices and queries Algorithm <ref> for the fidelity value.
According to the Pauli index defined in Appendix <ref>, the stabilizer group $G_k$ is generated by $\langle P_{\phi(\overline{k_1k_2},1)},P_{\phi(\overline{k_1k_2},1)},\cdots,P_{\phi(\overline{k_1k_2},n)}\rangle$, where result of $\phi: \mathbb{F}_2^2\times \{1,\cdots,n\}\rightarrow\mathbb{F}_2^{2n}$ is 0-extended to $2n$-bit string:

\begin{gather}
\phi(\alpha,i)=\begin{cases}
\overline{\alpha}\times 2^{2(n-i)} & \alpha\neq00\\
\overline{01}\times2^{2(n-i)} & \alpha=00.
\end{cases}
\end{gather}

Intuitively, that means we extend the original Pauli to 1-weight Pauli on every qubit as a generator and use $X$ when encountering $I$ on some qubits.

$\mathbf{M}$ here is some random Boolean matrix, and the procedure reshuffles the fidelity by this matrix.
In order to make bins represented in the same form of a hash function, the procedure further introduces the additional $J_n\coloneqq  I_n\otimes X$ to modify random matrices.
Vector $\mathbf{d}$ works as offsets to create variant phases.
Therefore, every bin is labeled by the matrix index $c$, the offset index $t$, and the intrinsic bit string $j$ with length $b$.
We denote the bins with the same $j$ and $c$ to be in a single bin set.
Therefore, there are $B=2^b$ bin sets regarding a random matrix.
With the summation and Walsh-Hadamard Transform, the bin set $\{U_{c;t}[j]\}_t$ actually constructs a partition over the complete Pauli error rates as stated in Lemma <ref>.

The $B$-point WHT subsampled bin coefficients with index $j\in\mathbb{F}_2^{b}$ can be written as:

\begin{align}\label{eq:U_cp}
{U}_{c;t}[j] &= \sum_{\alpha:\,\mathbf{M}_c^T\alpha = j} p_\alpha(-1)^{\langle\mathbf{d}_{c;t},\alpha\rangle_p}+{W}_{c;t}[j], \,\forall t\in[P].
\end{align}

Moreover, the sampling error is as follows

\begin{align*}
{W}_{c;t}[j] &= \sum_{\alpha:\,\mathbf{M}_c^T\alpha = j} {W}_\alpha(-1)^{\langle\mathbf{d}_{c;t},\alpha\rangle_p},
\end{align*}

where ${W}_\alpha$ is the noise of the Pauli error rate $p_\alpha$.

From Lemma <ref>, matrix $\mathbf{M}$ determines the hash function to partition the whole set of Pauli error rates into these bin sets.
We thus denote these matrices by subsampling matrices.
Since each bin set is a partitioned piece and there are only sparse nonzero error rates, it will be likely that every bin set contains only a single nonzero error rate, which allows us to learn the value directly.
To achieve this, the procedure must recruit bin sets number $B=O(s)$.

The procedure also constructs methods to determine the precise status of every bin set to verify we can acquire the desired error rate directly.
The following Algorithm <ref> gives a detailed illustration of the construction of this detector.
We can shed the most light on the case of noiseless processing.
According to Lemma <ref>, an offset vector decides the sign for every term in a bin.
Because every bin set employs $P$ different offsets, if all the terms are zero, the absolute values of these bins would be all zero.
In a similar manner, when there exists one nonzero term, all the absolute values would be the same.
In contrast, if these absolute values vary from different bins, the procedure claims there are several nonzero terms in a bin set.
For convenience, we denote the above three types of bins by zero-ton, single-ton, and multi-ton bins, respectively.
Obviously, bins in a bin set remain in the same type, so we may also use these types to label a whole bin set without leading to ambiguity.

Besides, this detector aims to extract both the value and the index of the nonzero error rate when it is the only nonzero term in a bin set.
The idea of index estimation relies on the choice of offsets $\mathbf{d}$.
Ideally, if there is one nonzero term without noise, the procedure would choose the set of offsets $\{d_{c,1},\cdots d_{c,P}\}$ to contain all the linear independent unit vectors in the binary space of length $2n$.
Therefore, the sign of each of the bins suggests the binary value of $\alpha$ on the corresponding bit according to Lemma <ref>.
As for the value of the corresponding error rate, the procedure arbitrarily fetches a bin in the bin set and eliminates the effect of the random sign to get the rate since we have known the index of this rate.

Suppose there exists some noise during the construction of bins.
The idea of distinguishing a bin set among the aforementioned three types keeps the same.
The procedure checks bins with the first $P_1$ random offsets and determines the type of this bin set.
Given a single-ton bin set with the error rate $p_\alpha$, we have the following

\begin{gather}
\sgn{\mathbf{U}_c[j]}=\mathbf{D}_c^T\alpha+W,
\end{gather}

where $\mathbf{D}_c$ is offset matrix from $P=P_1+P_2$ $\mathbf{d}$s, and $W$ represents the effects from noise.
In order to protect the information about the index $\alpha$ of the only nonzero term against noise, the procedure employs some linear error-correcting codes and implements each offset $\bf d$ among the last $P_2$ as a column of the code generating matrix.
Therefore, the signs of the last $P_2$ bins serve as an element of codewords that remains resilient against the noise from corrupting the signs.
In the detection, the procedure thus decodes signs back to index $\alpha$ if the code distance is large enough.
Given the index of the nonzero term, the procedure will eliminate signs in the first $P_1$ bins and average over these bins.
This helps to suppress the statistical fluctuations in these bins and get a much more precise estimation of the nonzero error rate.

Bin Detector: $\mathrm{BD}(\mathbf{U}_c,\mathbf{D}_c,T)$
[1]
Input: bin $\mathbf{U}_c$, offsets $\mathbf{D}_c$ and the number $T$ to indicate error size;
Parameter: real numbers $\gamma_1,\gamma_2\in(0,1)$;
$\frac{1}{P_1}\sum_{t=0}^{P_1-1}U_{c;t}^2\leq T(1+\gamma_1)\nu^2$
$\widehat{\mathfrak{B}} \leftarrow $
Return ($\widehat{\mathfrak{B}}$, nil, nil)    zero-ton verification
$\widehat{\alpha}\leftarrow \texttt{Decode}([\sgn {U_{P_1}},\cdots,\sgn{U_{P-1}}]^T)$
$\widehat{p}_{\widehat{\alpha}}\leftarrow\frac{1}{P_1}\sum_{t=0}^{P_1-1}(-1)^{\langle\mathbf{d}_{c;t},\widehat{\alpha}\rangle_p}U_{c;t}$ single-ton search
$\frac{1}{P_1}\sum_{t=0}^{P_1-1}(U_{c;t}-(-1)^{\langle\mathbf{d}_{c;t},\widehat{\alpha}\rangle_p}\widehat{p}_{\widehat{\alpha}})^2\leq T(1+\gamma_2)\nu^2$
$\widehat{\mathfrak{B}} \leftarrow$
Return ($\widehat{\mathfrak{B}}$, $\widehat{\alpha}$, $\widehat{p}_{\widehat{\alpha}}$) single-ton verification
$\widehat{\mathfrak{B}} \leftarrow$
Return ($\widehat{\mathfrak{B}}$,nil, nil)

However, the procedure cannot extract all the nonzero error rates in this way since it is highly possible that there are some bins including several nonzero terms.
The protocol implements the peeling among different bin sets using the single-ton bin set with the support of the bin detector.
As shown in Algorithm <ref>, the procedure will fetch the single-ton bin set to peel out the error rate of all other bin sets in different sampling groups which contain the same error rate.
By constantly peeling, the procedure generates more artificial single-ton bins, and the peeling will keep until there are no remaining nonzero error rates.
In the following, we give the full version of the peeling process with detailed parameters and calculations.
The procedure also maintains an array $\mathbf{T}$ to keep track of noise propagated among bins, which helps the bin detector return a reliable detection even with the occurrence of noise.

Peeling Decoder
[1]
${\tt Input}:$ Number of sampling groups $C$, number of observations P, offsets $\{\mathbf{d}_{c;t}\}_{CP}$, positive integers $l$ and $b$, and the unit time length $t_0$
${\tt Initialize}:$ $\mathbf{T}_{c}[j]\gets$ array for all bins with element to be all 1
${\tt Initialize}:$ $\mathcal{P}\gets$ empty list of second-order Pauli error rates ($\widehat{\alpha},\widehat{{p}}^{(2)}_{\widehat{\alpha}}$)
Run $\textbf{BC}(\{\mathbf{d}_{c;t}\},l,b,C,P,t_0)$ to construct sampling matrices $\{\mathbf{M}_c\}$ and bins $\{\mathbf{U}_c[j]\}$
$\text{hassingle}\gets$  false
$c\in[C]$ and $\emph{j}\in\mathbb{F}_2^{b}$
$(\widehat{\mathfrak{B}},\widehat{\alpha},\widehat{{p}}^{(2)}_{\widehat{\alpha}}) \leftarrow \textbf{BD}(\mathbf{U}_c[\emph{j}],\mathbf{D}_c,\mathbf{T}_{c}[\emph{j}])$
$\widehat{\mathfrak{B}}= \textrm{single-ton}$
$\text{hassingle}\gets$ true
$\mathcal{P}\gets$ ($\widehat{\alpha},\widehat{{p}}^{(2)}_{\widehat{\alpha}}$)
$c'\in[C]$ and $c'\neq c$
${\tt Locate~bin~index}$ $\emph{j}_{c'} \leftarrow \mathbf{M}_{c'}^T \widehat{\alpha}$
$\mathbf{U}_{c'}[\emph{j}_{c'}] \leftarrow \mathbf{U}_{c'}[\emph{j}_{c'}] - \widehat{p}_{\widehat{\alpha}}(-1)^{\langle\mathbf{D}_{c'},\widehat{\alpha}\rangle_p}$
${\widehat{\mathfrak{B}}}\neq \textrm{single-ton}$
continue to next $\emph{j}$
Return: $\mathcal{P}$

Even though the algorithms are rephrased from [Harper et al., 2021], the model of noise from desired fidelity terms varies in our implementation of this Pauli error rate estimator.
In this work, this estimator plays as a subroutine utilized to process the estimated fidelity from the fidelity estimation.
Therefore, we cannot impose a similar assumption as stated in [Harper et al., 2021] on the estimated fidelity.
On the other hand, the theorem and corresponding proofs can be modified in our implementation to coordinate with the fidelity estimator's form.

According to the fidelity estimator, there exist some instructions, such as the ratio estimation, that are inherently biased.
Generally, we denote the noise of the estimated fidelity $f^{(2)}$ by $w=\Delta+\omega$, where $\Delta=\mathbb{E}(w)$ is a constant that represents the bias of the estimation.
In this sense, the noise $\omega$ is random with zero mean, and we will further assume that $\omega$ is zero-mean Gaussian noise.
This is due to the central limit theorem as there include many samples of fidelity estimation in the first step.
Even though we have little knowledge about this noise, the whole noise $w$ is bounded according to Lemma <ref>.
For convenience, we denote the bound by $\mathcal{B}>0$.
Thus we model this noise as a constant belonging to an interval $\Delta\in[-\mathcal{B},\mathcal{B}]$ along with zero-mean Gaussian noise $\omega\sim \mathcal{N}(0,\mathcal{B}^2)$.

We must consider even further how this form of noise behaves in every bin.
A bin is constructed from up to $B$ fidelity terms, and all these fidelity terms contain independent random noise.
Thus the Gaussian part of bins' noise that comes from the linear combination of different $\omega$ is $W\sim\mathcal{N}(0,\frac{\mathcal{B}^2}{B})$.
As for the constant parts, we can treat them as biases of the underlying error rates.
From the Walsh-Hadamard transform, the size of bias $\Delta^{(p)}_\alpha$ for an error rate $p_\alpha$ is no larger than $\frac{B\mathcal{B}}{4^n}$.
In this case, we divide the general noise from the noisy estimation of Pauli fidelity into the two parts mentioned above.

In our execution of this error rate estimator, the main challenge comes from the noisy inputs.
The major part of dealing with the ubiquitous noise is the robust bin detector.
Based on this, the analysis of the failure or exception for running this estimator will focus on the resilience of the bin detector.
Primarily, we rephrase the tail bounds in [Li et al., 2014] since we will employ this bound heavily in our proof.

Given $\textbf{g},\textbf{k} \in \mathbb{R}^N$ where $\textbf{k}$ is an isotropic Gaussian random variable $\textbf{k}\sim\mathcal{N}(0,\nu^2\mathbbm{1}_N)$, then the following tail bound holds:

\begin{align}
\Pr\left(\frac{1}{N}\|\textbf{g}+\textbf{k}\|^2\geq\tau_1\right)\leq&\mathrm{e}^{-\frac{N}{4}\left(\sqrt{2\tau_1/\nu^2-1}-\sqrt{1+2\theta_0}\right)^2}\label{tailbound1}\\
\Pr\left(\frac{1}{N}\|\textbf{g}+\textbf{k}\|^2\leq\tau_2\right)\leq&\mathrm{e}^{-\frac{N}{4}\frac{\left(1+\theta_0-\tau_2/\nu^2\right)^2}{1+2\theta_0}}\label{tailbound2},
\end{align}

for $\tau_1, \tau_2$, and $\theta_0$ satisfying

\begin{gather*}
\tau_1\geq\nu^2(1+\theta_0),\ \tau_2\leq\nu^2(1+\theta_0),\ \theta_0 = \frac{\|\textbf{g}\|^2}{N\nu^2}.
\end{gather*}

As stated in Algorithm <ref>, the procedure aims to distinguish every input bin set as zero-ton, single-ton, or multi-ton.
Thus the number of type errors is six by enumerating all pairs between the preceding three.
To be more clear, we follow the notation in [Harper et al., 2021], and let $\mathfrak{B}$ be the real type and $\hat{\mathfrak{B}}$ to be the detected type.
Besides that, we also need to consider the case in which a bin is detected to be single-ton correctly while the detector estimates the error rate with either a value with a huge deviation or a wrong index.
Let us define these failures explicitly.

The execution of the bin detector (Algorithm <ref>) is regarded  as  a failure if the output falls into one of the following cases,

* Confusion between single-ton and zero-ton bins

\begin{align*}
f_1 = \Prob{\widehat{\mathfrak{B}}=\texttt{z}\,|\,\mathfrak{B}=\texttt{s}}\ \ f_2=\Prob{\widehat{\mathfrak{B}}=\texttt{s}\,|\,\mathfrak{B}=\texttt{z}}
\end{align*}

* Confusion between single-ton and multi-ton bins

\begin{align*}
f_3=\Prob{\widehat{\mathfrak{B}}=\texttt{m}\,|\,\mathfrak{B}=\texttt{s}}\ \ f_4=\Prob{\widehat{\mathfrak{B}}=\texttt{s}\,|\,\mathfrak{B}=\texttt{m}}
\end{align*}

* Confusion between zero-ton and multi-ton bins

\begin{align*}
f_5=\Prob{\widehat{\mathfrak{B}} = \texttt{m}\, |\, \mathfrak{B} = \texttt{z}}\ \ f_6=\Prob{\widehat{\mathfrak{B}} = \texttt{z}\, |\, \mathfrak{B} = \texttt{m}}
\end{align*}

* Errors on the estimated value or index

\begin{align*}
f_7=\Prob{\widehat{\mathfrak{B}}=\texttt{s},\widehat{\alpha}\not=\alpha\, |\,\mathfrak{B}=\texttt{s}, \alpha}\ \ f_8=\Prob{\widehat{\mathfrak{B}}=\texttt{s},\alpha,|\hat{p}_\alpha - p_\alpha|> \frac{2\mathcal{B}}{\sqrt{B}} \,|\,\mathfrak{B}=\texttt{s}, \alpha, p_\alpha}
\end{align*}

We then show a comprehensive lemma to illustrate the robustness in most executions of this detector.
The lemma states an exponentially vanishing failure probability with regard to some mild conditions that we expect to rule out later by the total probability formula.
In the proof, we will separately analyze the probability of all these eight types of failures.

Suppose Assumption A1 & A2 are satisfied, and suppose $B\geq O\left(\frac{\mathcal{B}^2}{\epsilon_0^2}\right)$.
With the fidelity extracted from Algorithm <ref> and bins constructed from Algorithm <ref>, the detector will deal with noisy bins.
Let $E$ denote the event that an arbitrary bin detection with inputs as those in <Ref> get failed as defined in Definition <ref>.
Let $D$ be the event that every bin contains at most $P_1$ nonzero terms
Let $V$ be the event that all the prior bin detection executes successfully.
Let $H$ be the event that the peeling graph is cycle-free.
Let $K$ denote that all the peelings do not rule out the randomness and that the bias parts remain limited dependence.

\begin{gather}
\Prob{E|D,V,H,K} \leq  \mathrm{e}^{-O(n)}.
\end{gather}

During the proof we denote $\frac{\mathcal{B}}{\sqrt{B}}$ by $\nu$.
The original noise in each bin can be divided as zero-mean Gaussian noise $\mathcal{N}(0,\nu^2)$ and a bias with variance less than $\frac{B^2\nu^2}{N}$ given the event $K$ as illustrated in Lemma <ref>.
From Assumption A1 & A2 and the constraint that $B\geq O\left(\frac{\mathcal{B}^2}{\epsilon_0^2}\right)$, we choose $B$ such that $\epsilon_0\geq \frac{4\mathcal{B}}{\sqrt{B}}$.
In the remainder of this proof, we will keep this constraint on $\epsilon_0$.

We will follow the Definition <ref> to check  all the possibility of variants of detection failures.
For the first one, we consider the case that a bin detection identify a single-ton bin as a zero-ton bin.
The failure rate is defined as follows,

\begin{align}
f_1 = \Prob{\widehat{\mathfrak{B}}=\texttt{z}\,|\,\mathfrak{B}=\texttt{s}}=\Prob{\frac{1}{P_1}\|\mathbf{U}_{\texttt{s}}\|^2\leq T(1+\gamma_1)\nu^2}.
\end{align}

According to the lemma 7 in [Harper et al., 2021], the maintained array $\mathbf{T}$ keeps track of the variance information of Gaussian noise in each bin.
Thus, the above failure rate can be bounded by Lemma <ref>, Lemma <ref>, and Markov inequality.

\begin{align}
f_1=\Prob{\frac{1}{P_1}\|\mathbf{U}_{\texttt{s}}\|^2\leq T(1+\gamma_1)\nu^2}\leq& \Prob{\frac{1}{P_1}\|\mathbf{U}_{\texttt{s}}\|^2\leq T(1+\gamma_1)\nu^2\,|\,\|\Delta \mathbf{W}\|_\infty\leq\frac{2\mathcal{B}}{B}}+\Prob{\|\Delta \mathbf{W}\|_\infty\geq\frac{2\mathcal{B}}{B}}\notag\\
\leq& \mathrm{e}^{-\frac{P_1}{4}\frac{((\epsilon_0-2\nu/\sqrt{B})^2/T\nu^2-\gamma_1)^2}{1+2(\epsilon_0-2\nu/\sqrt{B})^2/T\nu^2}}+\frac{P_1B^3}{N}.
\end{align}

Since the number $N=4^n$, $B$ is a polynomial parameter and $P_1= O(n)$, the failure rate is exponentially low.

Then, we consider the case that we recognize a zero-ton bin as a single-ton bin.

\begin{align}
f_2=\Prob{\widehat{\mathfrak{B}}=\texttt{s}\,|\,\mathfrak{B}=\texttt{z}}=\Prob{\frac{1}{P_1}\|\Delta \mathbf{W}+\mathbf{W}\|^2\geq T(1+\gamma_1)\nu^2}.
\end{align}

Similarly, we can use the formula of the total probability to bound this failure rate.

\begin{align}
f_2\leq& \Prob{\frac{1}{P_1}\|\Delta \mathbf{W}+\mathbf{W}\|^2\geq T(1+\gamma_1)\nu^2\,|\,\|\Delta \mathbf{W}\|_\infty\leq\frac{2\mathcal{B}}{B}}+\Prob{ \|\Delta \mathbf{W}\|_\infty\geq\frac{2\mathcal{B}}{B}}\notag\\
\leq& \mathrm{e}^{-\frac{P_1}{4}\left(\sqrt{1+2\gamma_1}-\sqrt{1+\frac{8}{BT}}\right)^2}+\frac{P_1B^3}{N}.
\end{align}

As for the failure rate $f_3$, it measures the probability of identifying a single-ton bin as a multi-ton bin.
For clearness, we denote the bin contains error rate $p_m$ actually, and the estimated rate is $\hat{p}_{\hat{\alpha}}$.
And the sign vector added to peel the error rate is represented by $\mathbf{s}_{\widehat{\alpha}}$

\begin{align}
f_3=\Prob{\widehat{\mathfrak{B}}=\texttt{m}\,|\,\mathfrak{B}=\texttt{s}}=\Prob{\frac{1}{P_1}\left\|\mathbf{U}_{\texttt{s}}-\widehat{p}_{\widehat{\alpha}}\mathbf{s}_{\widehat{\alpha}}\right\|^2 \geq  T(1+\gamma_2)\nu^2}.
\end{align}

Since here involves multiple error rates, we only consider the case that the estimated error rate satisfies that $\alpha=\widehat{\alpha}$ and $\abs{p_\alpha-\hat{p}_{\widehat{\alpha}}}\leq \nu$.
By Lemma <ref>, we have

\begin{align}
f_3\leq f_7+f_8+\mathrm{e}^{{-\frac{P_1}{4}\left(\sqrt{1+2\gamma_2}-\sqrt{1+\frac{2(\sqrt{B}+2)^2}{BT}}\right)^2}}+\frac{P_1B^3}{N}.
\end{align}

As for the next failure, we consider the case that the bin contains several nonzero terms while the detector identifies it as a single-ton bin.

\begin{align}
f_4=\Prob{\widehat{\mathfrak{B}}=\texttt{s}\,|\,\mathfrak{B}=\texttt{m}}=\Prob{\frac{1}{P_1}\left\|\mathbf{U}_\texttt{m}-\widehat{p}_{\widehat{\alpha}}\mathbf{s}_{\widehat{\alpha}}\right\|^2 \leq T(1+\gamma_2)\nu^2}.
\end{align}

We then decompose the bin vector by the Gaussian noise part $\mathbf{W}$, the bias part $\Delta \mathbf{W}$, the error rates $\mathbf{\beta}$, and the sign matrix $\mathbf{s}$.
The failure rate can be calculated as follows,

\begin{align}
f_4\leq \Prob{\frac{1}{P_1}\left\|\mathbf{U}_\texttt{m}-\widehat{p}_{\widehat{\alpha}}\mathbf{s}_{\widehat{\alpha}}\right\|^2 \leq T(1+\gamma_2)\nu^2\,|\,\frac{1}{P_1} \|\mathbf{s}\mathbf{\beta}+\Delta \mathbf{W}\|^2\geq 2T\gamma_2\nu^2}+\Prob{\frac{1}{P_1}\|\mathbf{s}\mathbf{\beta}+\Delta \mathbf{W}\|^2\leq 2T\gamma_2\nu^2}.
\end{align}

The first term in the right-hand side can be easily bounded by Lemma <ref> since $T$ tracks the variance size of every element in $\mathbf{W}$,

\begin{align}
\Prob{\frac{1}{P_1}\left\|\mathbf{U}_\texttt{m}-\widehat{p}_{\widehat{\alpha}}\mathbf{s}_{\widehat{\alpha}}\right\|^2 \leq T(1+\gamma_2)\nu^2\,|\,\frac{1}{P_1} \|\mathbf{s}\mathbf{\beta}+\Delta \mathbf{W}\|^2\geq 2T\gamma_2\nu^2}\leq \mathrm{e}^{-\frac{P_1}{4}\frac{\gamma_2^2}{1+4\gamma_2}}.
\end{align}

As for the second term, we also need to decompose the remaining two terms,

\begin{align}
\Prob{\frac{1}{P_1}\|\mathbf{s}\mathbf{\beta}+\Delta \mathbf{W}\|^2\leq 2T\gamma_2\nu^2}\leq  \Prob{\frac{1}{P_1}\|\mathbf{s}\mathbf{\beta}+\Delta \mathbf{W}\|^2\leq 2T\gamma_2\nu^2\,|\,\frac{1}{P_1}\|\mathbf{s}\beta\|^2\geq3T\gamma_2\nu^2}+\Prob{\frac{1}{P_1}\|\mathbf{s}\beta\|^2\leq3T\gamma_2\nu^2}.
\end{align}

Note that $\Delta W$ is a random variable with the  mean to be zero  from the summation of multiple Bernoulli variables.
Thus it is an approximate Gaussian variable by the central limit theorem.
The variance falls into the interval $\left[\frac{B}{N}\mathcal{B}^2,\frac{2B}{N}\mathcal{B}^2\right]$.
By Lemma <ref>, the first term can be bounded by

\begin{align}
\Prob{\frac{1}{P_1}\|\mathbf{s}\mathbf{\beta}+\Delta \mathbf{W}\|^2\leq 2T\gamma_2\nu^2\,|\,\frac{1}{P_1}\|\mathbf{s}\beta\|^2\geq3T\gamma_2\nu^2}\leq \mathrm{e}^{-\frac{P_1}{4}\frac{(1+\gamma_2N/2B^2)^2}{1+3\gamma_2N/B^2}}\ll1.
\end{align}

The second term can follow the analysis in Appendix E of [Li et al., 2014], and we get

\begin{align}
\Prob{\frac{1}{P_1}\|\mathbf{s}\beta\|^2\leq3T\gamma_2\nu^2}\leq N\mathrm{e}^{-\frac{P_1\epsilon}{4}\left(1-\frac{3T\gamma_2\nu^2}{\epsilon_0^2}\right)}.
\end{align}

Even $N=4^n$, we can choose a proper $P_1$ and guarantee that this is still an exponential vanishing failure rate.
Therefore, we have

\begin{align}
f_4\leq \mathrm{e}^{-O(n)}.
\end{align}

The following two failure types are easy to handle as they can be reduced to the previous rates.
For the failure that the procedure identifies a zero-ton bin as a multi-ton bin, this bin must pass the zero-ton verification first.
Thus, we have

\begin{align}
f_5=\Prob{\widehat{\mathfrak{B}} = \texttt{m}\, |\, \mathfrak{B} = \texttt{z}}\leq f_2=\Prob{\widehat{\mathfrak{B}}=\texttt{s}\,|\,\mathfrak{B}=\texttt{z}}
\end{align}

If the detector recognizes a multi-ton as a zero-ton bin, the probability of this failure is very similar to $f_4$, and we can claim that

\begin{align}
f_6=\Prob{\widehat{\mathfrak{B}} = \texttt{z}\, |\, \mathfrak{B} = \texttt{m}}\leq \mathrm{e}^{-O(n)}.
\end{align}

Then we need to consider the remaining index and value errors.
As the sign for every single-ton bin is generated from the inner product of offsets and the index as follows,

\begin{align}
\sgn{U_{c;t}[j]}=\sgn{(-1)^{\langle d_{c;t},\alpha\rangle_p}p_\alpha+\Delta W+W}=\langle d_{c;t},\alpha\rangle_p\oplus Z_t,\ \forall\,t\in\{P_1+1,\cdots,P_1+P_2\},
\end{align}

where sgn is a function that returns the sign of input and the Bernoulli variable $Z$ denotes whether the bias and noise flips the sign of this bin.
By deliberately choosing these $P_2$ offsets, we can use these bit strings to construct a generating matrix $\mathbf{G}$ of some error correcting code,

\begin{align}
\sgn{\mathbf{U}_c[j]}=\langle\mathbf{G}_c,\alpha\rangle\oplus\mathbf{Z}.
\end{align}

From the error-correcting codes' perspective, there are $2n$ logical bits, and the code words have length $P_2$.
Thus, the code rate is $R=\frac{2n}{P_2}\leq1$.
Besides, we denote the code protects error with Hamming weight up to $\beta P_2$ with a positive parameter $\beta\in[0,1]$.
To decide whether a code with a given distance can protect the index information, we need to figure out the property of the Bernoulli variables $\mathbf{Z}$.
By definition, the variable $Z=1$ if and only if $\abs{\Delta W+W}\geq p_\alpha$ and they contain opposite signs.
Therefore, we have

\begin{align}
\Prob{Z=1}\leq \frac{1}{2}\Prob{\abs{\Delta W+W}\geq p_\alpha}\leq\Prob{W\geq p_\alpha-\abs{\Delta W}\,|\,\abs{\Delta W}\leq \frac{2\mathcal{B}}{B}}+\Prob{\abs{\Delta W}\geq \frac{2\mathcal{B}}{B}}\notag.
\end{align}

Due to the fact that $W$ is Gaussian noise with zero mean and variance $T\nu^2$, there exists a tail bound for Gaussian random variables.

\begin{align}
\Prob{W\geq p_\alpha-\abs{\Delta W}\,|\,\abs{\Delta W}\leq \frac{2\mathcal{B}}{B}}\leq \mathrm{e}^{-\frac{\left(\epsilon_0-\frac{2\nu}{\sqrt{B}}\right)^2}{2T\nu^2}}\leq\mathrm{e}^{-\frac{\left(4-\frac{2}{\sqrt{B}}\right)^2}{2T}}.
\end{align}

The last inequality comes from the assumption about the smallest  Pauli error rate.
Combining this tail bound and the bound of the bias, the success probability of the Bernoulli variable is

\begin{align}
\Prob{Z=1}\leq\mathrm{e}^{-\frac{\left(4-\frac{2}{\sqrt{B}}\right)^2}{2T}}+\frac{B^3}{N}\ll\frac{1}{2}.
\end{align}

We denote this success rate by $\mathbb{P}$, and the index failure rate is therefore bounded by Hoeffding's inequality,

\begin{align}
f_7=\Prob{\widehat{\mathfrak{B}}=\texttt{s},\widehat{\alpha}\not=\alpha\, |\,\mathfrak{B}=\texttt{s}, \alpha}\leq \mathrm{e}^{-2P_2(\beta-\mathbb{P})}.
\end{align}

The value estimation comes from the single-ton search step in Algorithm <ref>.

\begin{align}
f_8=&\Prob{\widehat{\mathfrak{B}}=\texttt{s},\alpha,|\hat{p}_\alpha - p_\alpha|> 2\nu \,|\,\mathfrak{B}=\texttt{s}, \alpha, p_\alpha}=\Prob{\abs{\frac{\mathbf{s}_{\alpha}^T\mathbf{U}_{\texttt{s},\alpha}}{P_1}-p_\alpha}\geq 2\nu}\notag\\
\leq&\Prob{\abs{\frac{\mathbf{s}_{\alpha}^T}{P_1}(\Delta \mathbf{W}+\mathbf{W})}\geq 2\nu\,|\,\abs{ \frac{\mathbf{s}_{\alpha}^T\Delta\mathbf{W}}{P_1}}\leq \nu}+\Prob{\abs{ \frac{\mathbf{s}_{\alpha}^T\Delta\mathbf{W}}{P_1}}\geq \nu}\notag\\
\leq& \mathrm{e}^{-\frac{P_1}{2T}}+\frac{4B^2}{P_1N}.
\end{align}

Till now, we have completed the analysis of all types of failures that would occur in the detection of bins.
In the implementation of our algorithm, we choose $P_1 P_2\sim O(n)$ and $N=4^n$ is an exponentially increasing parameter.
Therefore, by the union bound, we have that

\begin{align}
\Prob{E|D,V,H,K} \leq \sum_{i=1}^8f_i\leq  \mathrm{e}^{-O(n)}.
\end{align}

Suppose Assumptions A1 & A2 are satisfied.
Suppose that every bin contains at most $\frac{n}{2}$ nonzero error rates and that every previous bin detection ran successfully before. In an arbitrary peeled bin set $U_{c_0}[j_0]$, the bias in every bin keeps bounded dependence, and the variance after the peeling will be bounded as follows,

\begin{align}
\text{Var}(\Delta W_{c_0}[j_0])\leq\frac{4B}{N}\mathcal{B}^2.
\end{align}

This above-mentioned observation succeeds with probability at least $1-\mathrm{e}^{-O(n)}$.

We prove this lemma inductively.
In the beginning, we first focus on the bias from an arbitrary bin before peeling.
According to Lemma <ref>, the Pauli fidelity terms are carrying bounded noise $w\in[-\mathcal{B},\mathcal{B}]$, and we can regard the bias of fidelity as still falling into this bound naturally.
The procedure constructs every bin by a specific linear combination of multiple fidelity terms, and Lemma <ref> tells the form of the resulting bins.
Let us consider the bias $\Delta W$ in a bin $U_{c;t}[j]$, which is in the form of the summation of constants with random signs,

\begin{align}
\Delta W_{c;t}[j]=\sum_{\alpha:\mathbf{M}^T_c\alpha=j}(-1)^{\langle d_{c;t},\alpha\rangle_p}\Delta_\alpha^{(p)},
\end{align}

where these $\Delta^{(p)}$ are the effective biases of every error rate by Walsh-Hadamard tranform.
Since all the bias terms are pair-wise independent, we can calculate the variance,

\begin{align}
\text{Var}(\Delta W_{c;t}[j])\leq \frac{B}{N}\mathcal{B}^2.
\end{align}

Then we need to analyze the effect of peeling processing on the biases of bins.
The peeling decoder algorithm points out that the single-ton bin will be averaged over random offsets to peel the corresponding error rate term in other bins.
The bias in a single-ton bin is propagated to the target bin during every peeling.

For brevity, we denote the $l$ single-ton bins and the target bin by $\{U_{c_1}[j_1],\cdots,U_{c_l}[j_l]\}$ and $U_{c_0}[j_0]$, respectively.
Averaging over those single-ton bins with $P_1$ randomly chosen offsets, and the target bin after the peeling is

\begin{align}\label{eq:peeling}
U_{c_0;t_0}[j_0]'= U_{c_0;t_0}[j_0]+\frac{1}{P_1}\sum_{k=1}^l\sum_{t=0}^{P_1}(-1)^{\langle d_{c_k;t}+d_{c_0;t_0},\alpha_k\rangle_p}U_{c_k;t}[j_k].
\end{align}

When the previous detection works perfectly, the right-hand side of Eq. (<ref>) becomes a single-ton bin, and all the other terms are either Gaussian noise or biases.
By the induction assumption, we can find all the involving bins have biases with variance less than  $\frac{4B}{N}\mathcal{B}$.
Note that the bias terms span all possible Pauli indices in a bin, and the average over different random signs is not completely independent.
The bias linked to the nonzero error rate will be propagated faithfully.
Thus we have

\begin{align}
\text{Var}(\Delta W_{c_0;t_0}[j_0]')\leq\frac{B}{N}\mathcal{B}^2+\frac{3lB^2}{N^2}\mathcal{B}^2+\frac{4lB}{P_1N}\mathcal{B}^2\leq \frac{4B}{N}\mathcal{B}^2,
\end{align}

where $2l\leq n\leq P_1\sim O(n)$.
The second term comes from the covariance of the desired nonzero terms' biases, and the third term is the independent part of summation.

In the above analysis, we have presumed two things.
For the first, we assume there is only one bias term from the peeling that performs dependently with the original bias terms.
This actually contains two aspects that both the peeling bin and the peeled bin have untouched biases with the corresponding Pauli index.
These rely on the constraint that the accumulation of all previous signs that come from peelings would not cause any collapses to the signs of nonzero terms in a bin.
Then we need to bound the probability that the accumulation of random signs in multiple peelings rules out the randomness.
This is possible because the propagating terms carry one more random sign during every peeling process.
Random signs must come from the inner product of the corresponding offset and a nonzero bit string in the kernel space of the subsampling matrix.
To check whether these random signs will annihilate themselves, we can construct a directed group with nodes representing bins and edges representing peeling processes.
Therefore, only the peeling that lays upstream in the corresponding path will affect the target bin's bias part.
Note there are $BC$ bin sets, and each bin set will be affected by at most $BC$ bin sets.
By the fact that the sum of all previous bitstrings must be zero to eliminate the random sign or equal to the original nonzero terms' indices, we have the failure rate $\Prob{K^c}$ for the above variance bound to be

\begin{align}
\Prob{K^c}\leq\frac{nB^3C^2}{2N}.
\end{align}

Since $C$ is a constant and $N$ is exponentially large, the failure rate is exponentially vanishing.

§.§ Theoretic bounds

The abovementioned lemmas help to remove obstacles over the engagement of the fidelity estimator and the error rate transformer.
Consequently, we can develop a comprehensive proposition to demonstrate the validity of combining these two subroutines.
This is based on the following assumptions.

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)$.

According to Lemma <ref>, each round of the fidelity estimator results in fidelity terms with errors in the size at most $\mathcal{B}=\frac{\sigma\epsilon^2(f^{(1)}_ir^{(0)}_i+r^{(0,1)}_i)}{(f^{(0)}_i-\epsilon r^{(0)}_i)t_0^2}$, where $\sigma$ is a constant related to the choice of fitting times.
Hence, we reduce the proposition to a claim that with bounded-noise fidelity, the transformer can estimate all the $s$ error rates correctly with their support indices.
This claim can be verified by a series of lemmas.

We first consider the failure probability and error bounds of the transformer.
Since the peeling decoder is the main routine for this transformer, the failure probability can be tracked as follows,

\begin{align}
\mathbb{P}_F
\leq \Prob{\text{Peeling decoder fails}\big|{E}_{\rm bin}^c}
+\Prob{E_{\rm bin}|D,H}+\Prob{D^c}+\Prob{H^c} \label{error_rate}.
\end{align}

The event ${E}_{\rm bin}$ denotes that there exist some failing bin detections.
On the contrary, ${E}_{\rm bin}^c$ denotes that no bin detection error occurred in the entire execution of <Ref>.
$D$ and $H$ denote that the maximum number of nonzero terms is at most $\frac{n}{2}$ and that the peeling routines are all cycle-free, respectively.
According to the Proposition 4 in [Li et al., 2014], we know the first term is vanishing with $s$ as $O\left(\frac{1}{s}\right)$.
The third term is exponentially vanishing according to Lemma 8 in [Harper et al., 2021], and the fourth term is at most $O\left(\frac{\log^{\log\log s}s}{s}\right)$ as stated in Lemma 6 of [Li et al., 2014].
And the remaining is the second term.

To bound the second term, we have the following equation,

\begin{align}
\Prob{E_{\text{bin}}\big| D,H}=1-\Prob{E_{\text{bin}}^c\big| D,H}=1-\Prob{\bigcap_{i=1}^M E_i^c\big| D,H},
\end{align}

where $M$ denotes the number of bin detection, and $E_i^c$ denotes the event that the $i$-th detection succeeds.
Thus, the question is derived to bound the probability of no anomaly detection.
By the definition of conditional probability, we have

\begin{align}
\Prob{\bigcap_{i=1}^M E_i^c\big| D,H}=\Prob{E_M^c\bigg|\bigcap_{i=1}^{M-1}E_i^c,D,H}\Prob{\bigcap_{i=1}^{M-1}E_i^c\bigg|D,H}
\end{align}

By Lemma <ref> and the fact that $B\geq O\left(\frac{\mathcal{B}^2}{\epsilon^2_0}\right)$, every term above is exactly what we bounded.
Therefore, all these terms are close to 1 with an exponentially decaying gap

\begin{gather}
\Prob{E^c_k\big|\bigcap_{i=1}^{k-1}E^c_i,D,H}\geq 1-\mathrm{e}^{-O(P_1)}.
\end{gather}

And we have that

\begin{align}
\Prob{E_{\text{bin}}\big| D,H}\leq 1-(1-\mathrm{e}^{-O(n)})^M\leq BCs\mathrm{e}^{-O(n)}\leq \mathrm{e}^{-O(n)}.
\end{align}

Thus, we get the stated failure probability.

We then consider the query number of the fidelity estimator, which represents the number of experiments needed to estimate the Pauli information of a Hamiltonian.
As illustrated by Proposition 4 in [Li et al., 2014], the peeling decoder succeeds with the probability of at least $1-O\left(\frac{1}{s}\right)$ adopting $C=O(1)$ subsampling groups and $B=\Theta(s)$ bin sets in every group.
Therefore, to prepare these bins, the number of queried eigenvalues is $BPC=O(Ps)$.
Note in every bin set, so we choose $P=P_1+P_2=O(n)$.
Therefore, we get the illustrated query number.

As for the error bound, since the event that the procedure estimates the error rates with noise larger than $O\left(\frac{\epsilon^2}{t_0^2\sqrt{s}}\right)$ is classified as value error, the above analysis of the failure probability demonstrates that all error rates fall into the stated interval in successful execution.
Therefore, the absolute values of decomposition parameters can be estimated by the square root of these Pauli error rates, and we get the stated precision.

§.§ With Prior Knowledge

From an experimental view of system's Hamiltonian learning, we will always hold some prior knowledge about the system's underlying dynamics.
For example, we can find out some significant interactions of the system that are gained from the experimental details or the designed properties of the device.
Formally speaking, given the prior knowledge, we can expect that the unknown Hamiltonian contains the corresponding terms on the Pauli basis with unknown parameters, which we refer to as the structure information.
Indeed, it is not always true that this prior information covers all the possible structure, so we cannot fully rely on the prior structure to estimate the Hamiltonian.
In this subsection, we would like to discuss how the partial prior knowledge about the structure would boost this Hamiltonian learning method, especially in the Pauli error rate estimation.

During the bin detection, the structure information can serve as a lower bound for the nonzero-term detection.
For example, given that there exists one known term in the detected bin $U$, we shall accept the detection only if it is found to be a single-ton or multi-ton bin.
Otherwise, the algorithm can output an "exception", and we shall adaptively adjust the execution parameter like $\gamma_1$ in Algorithm <ref>.
And the same happens when we have two or more terms in one bin.
As for the peeling step, since this procedure would invoke the bin detector as a subroutine, the prior information would not directly improve this step.

As discussed above, the prior information about the possible structure can help to improve the accuracy of the overall estimation.
However, since our method is designed without prior information, the benefits we get from the structure knowledge remains in a heuristic way and can be mainly observed in numerical results.

§ SIGN ESTIMATION

We have briefly summarized the sign estimation procedure in Sec <ref>, including the motivation and basic ideas.
In this section, we will illustrate the protocol to estimate signs of all decomposed parameters.
As this sign estimation plays the third step in the whole Hamiltonian estimation procedure, this estimation employs some information gained from previous Pauli estimation.

In the above section, we regard the Hamiltonian channel as a processing matrix form in Eq. (<ref>), or, more specifically, as a Pauli channel.
That is due to the fact that our target information is faithfully stored in diagonal terms.
On the contrary, we estimate the sign information mainly from off-diagonal terms in this section, so we need to care for the whole channel.
With the expansion in Eq. (<ref>), we consider the SPAM as a local Pauli eigenstate $\rho_\gamma$ and a Pauli operator $P_\beta$ with the measurement outcome as follows,

\begin{align}\label{eq:originlineq}
\Tr(\mathcal{H}_t(\rho_\gamma)P_\beta)=\Tr(\rho_\gamma P_\beta)+it\sum_{\alpha\in{\sf P}^n}s_\alpha\Tr(\rho_\gamma P_\alpha P_\beta-P_\alpha\rho_\gamma P_\beta)+o(t^2),
\end{align}

where we want to solve all decomposition parameters $\bf s$.

In order to solve all the parameters to get the sign information, we construct a family of linear equations to extract all the interesting parameters.
Firstly, we want to rule out the effect of the constant term in Eq (<ref>).
We use variant $t$ and fit the measurement results on these $t$.
Therefore, we get the first-order version of the Hamiltonian evolution.

\begin{align}\label{eq:lineq}
\Tr(\mathcal{H}^{(1)}(\rho_\gamma)P_\beta)=i\sum_{\alpha\in{\sf P}^n}s_\alpha \Tr(\rho_\gamma [P_\alpha, P_\beta]),
\end{align}

where $[P_\alpha, P_\beta]$ is the commutator notation.
Moreover, since we know the support of nonzero ${\bf s}^\star$ exactly from the last section, it helps to reduce the number of unknown variables and equally the number of equations needed.
Sequentially, we choose multiple state and measurement settings to construct the necessary equations.
The collection of SPAM settings is denoted by $\{(\rho_1,M_1),\cdots,(\rho_m,M_m)\}$.
Collect these first-order measurement results, and this vector is denoted by $\mathcal{M}$.
For the coefficient matrix $\Phi$, we define $\Phi_{j,\alpha}\coloneqq i\Tr(\rho_{\gamma_j}[P_\alpha,M_{\beta_j}])$.
As we pick $m$ pairs states and measurements, and $\Phi$ only covers $s=\abs{{\bf s}^\star}$ nonzero parameters, the matrix is in the shape of $m\times s$.
Therefore, we construct the following equations based on Eq. (<ref>)

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

The challenge arises when the observations $\mathcal{M}$ are noisy.
This is possible as we have some noisy SPAMs, and the linear regression brings systematic errors due to the higher-order terms.
Consider the case that the protocol can only estimate $\mathcal{M}$ with the additive noise $\omega$, so the equation becomes

\begin{gather}\label{eq:noisyeq}
\hat{\mathcal{M}}=\Phi\cdot{\bf s}^\star+\omega.
\end{gather}

Indeed solving the equation directly is not applicable due to the unknown noise. Instead, we adopt the method from compressed sensing [Candes and Tao, 2005, Candes et al., 2006, Baraniuk et al., 2008] and consider the problem of finding a feasible vector $x$ satisfying the following constraints

\begin{align*}
\begin{split}
\min_x \|x\|_1&\\
||\Phi x- \hat{\mathcal{M}}||_2&\leq \epsilon
\end{split}.
\end{align*}

This problem is also illustrated in the main text as Eq. (<ref>).
Notably, this problem has a nonempty feasible set for every $\epsilon\geq\|\omega\|_2$ since the ideal value of ${\bf s}^\star$ is a feasible solution.
We always choose $\epsilon=\sqrt{m}\left(\frac{\sigma\tau}{t_1}+o(t_1)\right)$, which is the upper bound of $\|\omega\|_2$.
We denote this problem as sign optimization.

We will show that the solution $x$ is a well-qualified estimation of the ideal parameters given the measurements and states are randomly chosen from the Pauli group and local eigenstates thereof.
In the following, we first give a definition of the restricted isometry property introduced by [Candes and Tao, 2005, Baraniuk et al., 2008].
Then we show how to prove this property for the randomly sampled process matrices in our specific case and how to bound the noise effects.

We say that a matrix $\Phi$ satisfies the Restricted Isometry Property (RIP) of order $s$ if there exists a
$\delta_s\in(0,1)$ such that

\begin{gather}
(1-\delta_s)\|x_T\|_2^2\leq\|\Phi x\|_2^2\leq(1+\delta_s)\|x_T\|_2^2
\end{gather}

holds for all $x$ on all possible support $T$ such that $\#(T)\leq s$.

Suppose we know the exact support of the vector ${\bf s}^\star$ on the Pauli group.

Thus, with the number of equations $m\geq \Omega(s)$, the solution $x^\star$ from Algorithm <ref> is close to the ideal parameter ${\bf s}^\star$ with probability $1-\mathrm{e}^{-O(m-s)}$,

\begin{gather}\label{eq:spambound}
\|x^\star-{\bf s}^\star\|_2\leq C\cdot \epsilon,
\end{gather}

where $C$ is a constant related to the choice of the sampling set of SPAM.

We first figure out the size of noise $\omega$ that comes from SPAM errors and linear regressions in (<ref>).
According to Assumption <ref>, every measurement carries noise at most $\tau$.
Therefore, we have

\begin{gather}
\Tr(\mathcal{H}_t(\tilde{\rho}_\gamma)\tilde{P}_\beta)=\Tr(\rho_\gamma P_\beta)+it\sum_{\alpha\in{\sf P}^n}s_\alpha\Tr(\rho_\gamma P_\alpha P_\beta-P_\alpha\rho_\gamma P_\beta)+o(t^2)+\tau,
\end{gather}

which satisfies for all $t$.
Considering the ordinary least square, we have

\begin{gather}
\abs{\Tr(\mathcal{H}^{(1)}_t(\tilde{\rho}_\gamma)\tilde{P}_\beta)-\Tr(\mathcal{H}^{(1)}_t(\rho_\gamma)P_\beta)}\leq\frac{\sum_{i=1}^K\abs{i\cdot t_1-\Bar{t}}}{\sum_{i=1}^K(i\cdot t_1-\Bar{t})^2}\cdot\tau+o(t),
\end{gather}

where $K$ denotes the number of evolving times used to extract the first-order term and $\Bar{t}$ is the averaged time length.
As we choose the $t$ from $t_1$ up to $Kt_1$ arithmetically, the foregoing bound is derived to $\frac{\sigma\tau}{t_1}+o(t_1)$ for some regression constant $\sigma$, which is very similar to Eq. (<ref>).
Therefore, the norm of the total noise term $\omega$ can be bounded by

\begin{gather}
\|\omega\|_2\leq\sqrt{m}\left(\frac{\sigma\tau}{t_1}+o(t_1)\right).
\end{gather}

According to Algorithm <ref>, the chosen $\epsilon$ makes sure that $\epsilon\geq\|\omega\|_2$, which guarantees the existence of the solution of the optimization problem.

We then prove the approximate RIP for the constructed $\Phi$ matrix.
The proof is inspired by [Shabani et al., 2011, Baraniuk et al., 2008].
In the noisy version, we randomly choose the local eigenstates of Pauli operators and Pauli measurements, so $\Phi$ is a $m\times s$ random matrix.
By the fact that all these states and measurements are in a finite dimension system, the matrix $\Phi$ is always bounded.
Consider an arbitrary $s$-sparse vector $x$, and we have the following upper and lower bounds,

\begin{gather}
\frac{w_l}{m}\|x\|^2\leq x^T\phi_i\phi_i^Tx\leq\frac{w_u}{m}\|x\|^2,\ \ \forall\,x\in\mathbb{R}^s\\
l\|x\|^2\leq \mathbb{E}_{\Phi}\|\Phi x\|^2\leq u\|x\|^2,\ \ \forall\,x\in\mathbb{R}^s,
\end{gather}

where we use $w_l,w_u,l,u$ to show the upper and lower bounds, and we denote the $i$th row of $\Phi$ by $\phi_i$.
Implement Hoeffding's inequalities on the random matrix $\Phi$, we have

\begin{align}
\Prob{\|\Phi x\|^2\geq (1+\delta_s)\|x\|^2}\leq&\mathrm{e}^{-2m(1+\delta_s-u)^2/(w_u-w_l)^2}\\
\Prob{\|\Phi x\|^2\leq (1-\delta_s)\|x\|^2}\leq&\mathrm{e}^{-2m(1-\delta_s-l)^2/(w_u-w_l)^2},
\end{align}

for some positive $\delta_s$.
Combining these two, we get,

\begin{align}\label{eq:rip}
\Prob{\abs{\|\Phi x\|^2-\|x\|^2}\geq \delta_s\|x\|^2}\leq 2\mathrm{e}^{-2m(\delta_s+\mu_0)/(w_u-w_l)^2},
\end{align}

where $\mu_0$ denotes $\min(1-u,l-1)$.

Now we need to generalize this inequality to every vector in a specific $s$-sparse support.
In our case, that is the parameter support we learned from stage 1.
According to Lemma 5.1 in [Baraniuk et al., 2008],we can use the counting and union bound to prove with probability $1-\mathrm{e}^{-2m(\delta_s+\mu_0)/(w_u-w_l)^2+s\log(12/\delta_s)}$, every $x$ on the support gained from stage 1 satisfies

\begin{gather}
(1-\delta_s)\|x\|_2^2\leq\|\Phi x\|_2^2\leq(1+\delta_s)\|x\|_2^2.
\end{gather}

Note that we can decide the sparse support, so it is not necessary to further count all possible supports.
We need to choose $m>s$ to get a vanishing failure probability, and this is consistent with the intuition that the number of linear equations must be larger than that of unknown variables.

Suppose we get a solution $x^\star$ in the feasible set of the problem in Eq. (<ref>).
Therefore, according to Theorem 1.3 in [Candès, 2008] and the definition of RIP matrix in [Candes and Tao, 2005], the result of Eq. (<ref>) satisfies the following distance bound toward the real parameter ${\bf s}^\star$

\begin{align}
\|x^\star-{\bf s}^\star\|_2\leq \frac{2\sqrt{1+\delta_s}}{1-(1+\sqrt{2})\delta_s}\cdot\epsilon. \notag
\end{align}

Note the constant $\delta_s$ is artificially picked.
Therefore, we can always choose a proper $\delta_s<\sqrt{2}-1$.

To actually rule out the effects of noise on the sign estimation, we need to make sure the additive noise $\omega$ is bounded by some bars that depend on the lower bound of the nonzero decomposition parameters.
This can be achieved by calculating the propagation of the noise and the chosen number of equations $m$.
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)}$.

As proved in Lemma <ref>, the process of solving the linear equations will not enlarge the errors.
The solved decomposition parameters $x^\star$ satisfy that following bound with a probability of at least $1-\mathrm{e}^{-O(m-s)}$

\begin{gather*}
\|x^\star-{\bf s}^\star\|_2\leq C\cdot\sqrt{m}\left(\frac{\sigma\tau}{t_1}+o(t_1)\right),
\end{gather*}

where $C$ comes from the constant in Lemma <ref>.
According to Assumption A1 & A2 and <ref>, $m$ determines the upper bound of norm distance, which is supposed to be smaller than $O(\sqrt{\epsilon_0})$.
Therefore, this error indicates that there is no sign flip between the estimated parameters $x^\star$ and the actual decomposition vector ${\bf s}^\star$.

Similar to Remark <ref>, we need to choose some proper $t_1$ to attune the two types of noise.
Thus, in our model, we would choose $t_1\sim\sqrt{\sigma\tau}$.

§ MAIN THEOREM

Till now, we have introduced our protocol that combines the fidelity estimator, Pauli error rates estimator, and the sign estimator.
The validity of this protocol, therefore, relies on those guarantees of subroutines.
In this section, we introduce the main theorem that demonstrates the correctness and efficiency of our protocol working on the general $n$-qubit quantum systems.

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 $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.

According to Proposition <ref>, Algorithm  <ref> and <ref> with these parameters can estimate the absolute values of decomposition parameters up to the accuracy of $\|\abs{\hat{s}^\star}-\abs{s^\star}\|_\infty\leq O\left(\frac{\epsilon}{t_0\sqrt[4]{s}}\right)$.
Moreover, the sign estimator will extract perfect sign information from Proposition <ref>.
Therefore, the illustrated accuracy can be achieved.
The failure probability comes from the union bound of these two stages, which are $\delta+\left(\frac{1}{s}\right)$ and $\mathrm{e}^{-O(m-s)}$, respectively.
As $m=O(s)$, we can claim the stated failure probability.

As for the measurement complexity, it is gained from simple counting.
In the first stage, we run the cascading circuit with altogether $l\cdot\abs{\kappa}$ sequences.
The set $\kappa$ includes exponentially increasing sequence lengths from 1 to $O\left(\frac{1}{\Delta}\right)$ where $\Delta$ refers to the smallest fidelity residual of the detected channel.
In the second stage, the protocol needs to prepare $m$ equations, which are equal to $m$ measurements.
Therefore, the complexity $M_Q$ satisfies

\begin{gather*}
M_Q=l\cdot \abs{\kappa}+m=\tilde{O}\left(\frac{sn}{\epsilon^4}\right).
\end{gather*}

The classical complexity is also combined from the foregoing two stages.
In the first stage, the post-processing is mainly contributed by the peeling algorithm, of which the complexity is stated as $O(sn^2)$ in Theorem 1 of [Harper et al., 2021].
Also, to extract the desired fidelity terms, the procedure needs $O(sn)$ fittings.
In the second stage, the optimization problem illustrated in Eq. (<ref>) is of size $m\times s$.
Since we have confirmed the existence of the solution, the complexity is $\text{poly}(s)$.
Therefore, we have

\begin{gather}
\end{gather}

§ IMPLEMENTATION AND NUMERICAL RESULTS

This section serves as the further analysis and numerical detection of some concerns that are raised in the main text.

§.§ Threshold Behaviors

In this section, we perform the exhaustive numerical simulation to show how we choose the bin parameter $b$ in each case of the simulations introduced in <Ref>.
Even though the procedure runs with an inputted parameter, $b$, this parameter is supposed to rely on the structure of the target Hamiltonian to guarantee the successful execution.
Hence, we have to employ a trial process to determine the proper $b$ for each Hamiltonian.
More clearly, we increase the parameter $b$ and run the procedure on the target Hamiltonian to find the proper choice of $b_0$ such that the estimation with $b\geq b_0$ is steadily close to the ideal values.
In Figure <ref>, we exhibit the simulation results of running the procedure exhaustively to detect the behaviors of different $b$ on different systems.

This figure includes the thorough results for the trial processes for successfully executing the learning procedure.
In figure (a), we show the execution of random TFIM Hamiltonian with multiple system sizes.
The metric for errors is chosen as the relative 1-norm distance.
Note in the case of $n=1$, all the choices of $b$ perform a very similar estimation distribution, which is due to the only one nonzero term in this case. In the case of $n=2$, the reconstruction becomes constantly close to the ideal values when $b\geq4$, which implies the proper choice, $b=4$. For $n=3$, the steady reconstruction shows up when the procedure chooses $b\geq4$.
As for cases of $n=4$, $n=5$, and $n=7$, the steady estimations are recognized when $b\geq 5$. These trials solve the concerns of choosing proper $b$ raised in the TFIM simulation in Figure <ref>.
In figure (b), we focus on the $b$ choosing problem of the molecular Hamiltonian learning.
In the case of learning $H_2$ Hamiltonian, the reconstruction errors decay fast by increasing $b$ up to 5.
After $b\geq 6$, the estimation performs well, which indicates that $b=6$ is a proper choice.
In the case of $H_3$ and $LiH$ molecules, the proper choices are $b=6$ and $b=7$, respectively.

In Figure <ref>, we display the trial results for the random TFIM Hamiltonian with different system sizes.
In the case of $n=1$, all these $b$ perform well due to the fact that there is only one nonzero parameter in the Hamiltonian.
In the case of the $n=2$ system, these results are more typical, which is very similar to Figure <ref>.
According to the box plot, the procedure performs much better when $b\geq2$, and the performance keeps steady after $b$ grows up to 4.
Therefore, the procedure will choose $b=4$ as a proper parameter to execute on the 2-qubit system.
In all other cases, this process works very similarly.
Based on the understanding of the requirement of the parameter $b$ and the observations on Figure <ref> and Figure <ref>, parameters for random TFIM Hamiltonian learning are just the same as those used in <Ref>.

In addition to the random TFIM systems, we also check the threshold behaviors for molecular systems.
The numerical results of the trial processes are shown in Figure <ref>(b).
A key property of molecular Hamiltonian is that they contain much more nonzero terms than the Ising models.
Hence, the threshold behaviors are generally asking for larger $b$.
In the first plot, the reconstruction errors of $H_2$ molecules decay rapidly from $b=2$ to $b=5$.
After the procedure increases $b$ up to 6, the reconstruction becomes steadily close to the ideal decomposition parameters.
Based on these results, the procedure regards $b=6$ as a proper choice of the bin parameter, and we exhibit the corresponding distribution of $b=6$ estimations in Figure <ref>.
Similarly, we run the trial simulation for $H_3$ and $LiH$ molecules, respectively.
By witnessing the reconstruction with steady errors, the procedure can then determine the proper parameter in each case.

§.§ Higher-Order Fitting

During the simulation for the numerical results, we found there exist some lower bounds of the errors after we detect the distribution of our reconstructions for systems with comparable large sizes.
This means there are some systematic errors in some processes of the procedure.
In this section, we will show these errors are mainly contributed by the higher-order terms during the fitting of the second-order fidelity.

Firstly, we need to consider the definite expression of the fidelity terms of a Hamiltonian channel.
According to Eq. (<ref>), we need to expand this series to even higher-order terms.

\begin{align}
\mathcal{H}_t(\rho)
=&\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]\\
&+it^3\sum_{\alpha,\beta,\gamma\in{\sf P}^n}s_\alpha s_\beta s_\gamma\left[\frac{1}{6}(\rho P_\alpha P_\beta P_\gamma-P_\alpha P_\beta P_\gamma\rho)+\frac{1}{2}(P_\alpha\rho P_\beta P_\gamma-P_\alpha P_\beta \rho P_\gamma)\right]\notag\\
&+t^4\sum_{\alpha,\beta,\gamma,\delta\in{\sf P}^n}s_\alpha s_\beta s_\gamma s_\delta\left[\frac{1}{24}(\rho P_\alpha P_\beta P_\gamma P_\delta+P_\alpha P_\beta P_\gamma P_\delta\rho)-\frac{1}{6}(P_\alpha \rho P_\beta P_\gamma P_\delta+P_\alpha P_\beta P_\gamma \rho P_\delta)+\frac{1}{4}P_\alpha P_\beta\rho P_\gamma P_\delta\right]\notag\\
\end{align}

Note that these Pauli indices are all equivalent, and we can use this property to determine the effects on diagonal terms (Pauli terms).
For the first-order terms, they simply cause no diagonal effects. For the second-order terms, the diagonal effects from them are depicted by Eq. (<ref>), which is the major source of learning the Hamiltonian in our protocol.
As for the third-order terms, there are also zero effects for them on the diagonal terms.
The first terms obviously contribute no effects on the Pauli terms since they annihilate themselves on the Pauli terms.
The second terms have no effects either due to the rotational symmetry of the indices.
The fourth-order terms contribute significantly to the diagonal terms.
Apparently, so far, the odd-order terms contain no effects of the diagonal terms, and this can be easily extended to all the odd-order terms.
This guides us to further improve our fitting process by considering other even-order terms.

We also show the numerical results to verify this idea explicitly.

This figure depicts the effects of different fitting settings on the estimation of Hamiltonians.
To rule out the statistical fluctuations in the detection, we run different fitting processes on 50 random 4-qubit TFIM Hamiltonians for 10 rounds each.
For $m=2$, it is the ordinary fitting used for the vast estimation exhibited in the main text.
For $m=4$, we consider the effects of fourth-order terms besides the second-order terms in the fitting process. This modification asks for a little more time steps for regression while bringing much better estimation.
As for $m=6$, we consider the further sixth-order terms, and this modification requires even more time steps to implement six-order regression, which increases the numbers of oracle used.

In Figure <ref>, we exhibit the reconstruction cases with different fitting processes and fixed fidelity noise from finite shot of circuit.
In this figure, $m$ denotes the highest terms that we consider in the fitting.
As can be observed, the estimation is better when considering the fourth-order terms than the ordinary fitting with a mild sacrifice of the measurement complexity due to more time step needed for regression.
This proves the previous idea about systematic errors.
On the other hand, when the procedure keeps fetching higher-order terms, both the errors and complexity grow due to the effects of the circuit noise.
Even though the systematic errors are suppressed by considering higher-order terms, the circuit noise would be enlarged by the higher-order regression.
Going from $m=4$ to $m=6$, the circuit noise is now the dominant contributor of the overall errors.
Therefore, we propose an effective solution to significantly alleviate the systematic errors in the fitting process, while we still need to suppress the circuit noise by further increase the shot numbers in order to balance the systematic errors and circuit noise of fidelity.