# Practical limitations of quantum data propagation on noisy quantum
processors

Gaurav Saxena<EMAIL_ADDRESS>LG Electronics Toronto AI Lab, Toronto,
Ontario M5V 1M3, Canada  Ahmed Shalabi  LG Electronics Toronto AI Lab,
Toronto, Ontario M5V 1M3, Canada  Thi Ha Kyaw  LG Electronics Toronto AI Lab,
Toronto, Ontario M5V 1M3, Canada

###### Abstract

The variational quantum imaginary time evolution algorithm is efficient in
finding the ground state of a quantum Hamiltonian. This algorithm involves
solving a system of linear equations in a classical computer and the solution
is then used to propagate a quantum wavefunction. Here, we show that owing to
the noisy nature of current quantum processors, such a quantum algorithm or
the family of quantum algorithms will require single- and two-qubit gates with
very low error probability to produce reliable results. Failure to meet such
condition will result in erroneous quantum data propagation even for a
relatively small quantum circuit ansatz. Specifically, we provide the upper
bounds on how the relative error in variational parameters’ propagation scales
with the probability of noise in quantum hardware. We also present an exact
expression of how the relative error in variational parameter propagation
scales with the probability of partially depolarizing noise.

##  I Introduction

Quantum computers promise a quantum advantage in numerous applications such as
factoring [1], solving linear systems [2], quantum simulations [3], and
quantum chemistry/material discovery [4]. The aforementioned quantum
applications sometimes require millions of physical qubits and gates [5],
which current quantum hardware does not possess [6] and hence hinders the
direct experience of any proposed practical quantum speedup [7]. In the
current noisy intermediate-scale quantum (NISQ) era [8], many interesting
practical applications are evaluated through hybrid quantum-classical
algorithms otherwise known as NISQ algorithms [9, 10, 11]. These algorithms
are designed for short-depth quantum circuits with limited numbers of qubits
and gates. One such important class of algorithms are variational quantum
algorithms (VQAs) [12, 13, 14]. VQAs use short-depth parameterized quantum
circuits to perform some quantum computation and then leverage classical
computing power for optimization in a feedback fashion. It is, however, an
open question whether NISQ computers will yield a quantum advantage for
practically meaningful problems.

One significant family of VQA algorithms is the quantum imaginary time
evolution (QITE) [15, 16, 17]. QITE finds its uses in various applications
including efficiently finding the ground state of a Hamiltonian, diffusion
quantum Monte Carlo, etc [16, 15, 17, 18, 19, 20, 21, 22, 23]. QITE works by
propagating quantum data to evolve a quantum state. In general, imaginary time
evolution or a Wick rotation [24] is a simple but powerful mathematical
technique that involves mapping to the time coordinate $\tau=it$. While
simulating imaginary time evolution can be done classically, the resources
required scale exponentially with the size of the Hilbert space in question
[25]. Performing _quantum_ imaginary time evolution involves non-unitary
evolution and thus, implementing it using digital quantum computers is hard as
quantum processors only allow for unitary operations via quantum gates. To
address this limitation, a hybrid quantum-classical variational approach was
proposed to implement the quantum imaginary time evolution [16, 15, 17].
Imaginary time evolution of the Wick-rotated Schrödinger equation is
implemented by encoding an initial quantum state as a parameterized quantum
state $|\varphi(\tau)\rangle\approx|\phi(\theta(\tau))\rangle$ where $\theta$
is a real valued parameter vector
$\theta(\tau)=\left(\theta_{1}(\tau),\theta_{2}(\tau),\ldots,\theta_{N}(\tau)\right)$.
To simulate non-unitary dynamics, an approximate formulation of imaginary time
evolution derived from a generalized McLachlan’s variational principle is used
to evolve the parameterized variational circuit [21]. The final step is to
solve a linear equation of the form, (see Fig. 1(a) and section II)

$M\dot{\mathbf{\theta}}=Y\,.$ (1)

The solution $\dot{\theta}$ then dictates the evolution of variational
parameters of the ansatz circuit for the next time step:
$\theta(\tau+\Delta\tau)\approx\theta(\tau)+\dot{\theta}\times\Delta\tau$. The
matrix $M$ and the vector $Y$ depend on the ansatz circuit and the given
Hamiltonian. Similarly, for the real-time evolution, we have the same equation
where $Y$ is replaced by a different vector $V$ and the equation becomes
$M\dot{\mathbf{\theta}}=V$ (see Section II for details).

Figure 1: (a) A simple illustration to show that owing to noisy nature of the
current quantum processors, the outcomes from the QPU are misestimated and
thus cannot be relied upon without any error mitigation. (b) Plot of the upper
bound on the relative error in the rate of change of the parameters with
respect to the probability with which the depolarizing noise occurs for
different ansatzes with same condition number but different number of
parameters $N$, when $N=6$ and $N=14$ can be seen as simulating H2 and LiH
molecules in the minimum basis set using two layers of the hardware-efficient
ansatz [26], respectively.

Even though the family of QITE algorithms is well-studied in the literature
and various improved algorithms of the same kind have been developed/proposed
recently [27, 28, 29], none of the works have assessed [30, 31, 16, 32] its
performance on noisy quantum processors. Given the current noisy nature of
quantum processors, it becomes crucial to understand how well the quantum
imaginary (real) time evolution algorithm performs in absence of any error
correction or mitigation technique. In particular, we are interested in the
impact of hardware errors incurred during each iteration on the exact rate of
change of parameters, $\dot{\theta}$ used to estimate the next set of
variational parameters $\theta(\tau+\Delta\tau)$.

In this Article, we provide the quantitative answer to the above question by
providing upper bounds on the error incurred in exactly solving
$M\dot{\theta}=Y$ (see Fig. 1(b)). This paper is arranged as follows. First,
in Section II, we provide a brief review to the quantum imaginary time
evolution for completeness. Then in Section III, we give the details of our
error model. Section IV presents our main results followed by Discussion in
Section V.

##  II The Quantum Imaginary Time Evolution Algorithm for mixed states

In this section, we will present an overview of the QITE algorithm for mixed
states. For a more comprehensive overview, we refer the reader to Ref. [21]
and citations therein. Before we begin, a quick aside on notation. $R_{k}$
denotes the $k$-th rotation gate (which is a single or two-qubit gate) and
$\mathcal{R}_{k}$ denotes the $k$-th quantum gate acting on the whole system.
$\mathcal{R}_{k}(\tau)$ on an operator $\tau$ is to be understood as the
action of $R_{k}$ gate on $\tau$ by conjugation, i.e., $R_{k}\tau
R_{k}^{\dagger}$. $R_{k}\tau$ is to be understood as the multiplication of the
operators.

Under imaginary time evolution where the imaginary time $\tau:=it$, an initial
mixed state $\rho(0)$ evolves according to

$\displaystyle\rho(\tau)=\frac{e^{-H\tau}\rho(0)e^{-H\tau}}{\operatorname{Tr}\left[e^{-2H\tau}\rho(0)\right]}.$
(2)

The above time evolution can be rewritten as

$\displaystyle\frac{d\rho(\tau)}{d\tau}=-(\\{H,\rho(\tau)\\}-2\langle
H\rangle\rho(\tau)).$ (3)

Encoding the state $\rho$ in a parameterized circuit gives us the following
time evolution equation for a parameterized mixed state $\rho(\theta(\tau))$

$\displaystyle\sum_{i}\frac{\partial\rho(\vec{\theta})}{\partial\theta_{k}}\dot{\theta}_{k}=-(\\{H,\rho(\vec{\theta}(\tau))\\}-2\langle
H\rangle\rho(\vec{\theta}(\tau))).$ (4)

The evolution of the $\theta$ parameters can be derived by calculating the
variation $\delta\rho(\theta)$ and using the MacLachlan’s variational
principle

$\displaystyle\delta\|d\rho/d\tau+(\\{H,\rho\\}-2\langle H\rangle\rho)\|=0,$
(5)

which yields the following matrix equation for the time evolution of
$\vec{\theta}$

$\displaystyle\sum_{q}M_{k,q}\dot{\theta}_{q}=Y_{k}.$ (6)

The matrix elements of $M$ are defined as

$\displaystyle
M_{k,q}:=\operatorname{Tr}\left[\left(\frac{\partial\rho(\vec{\theta}(t))}{\partial\theta_{k}}\right)^{\dagger}\frac{\partial\rho(\vec{\theta}(t))}{\partial\theta_{q}}\right],$
(7)

and for $Y$

$\displaystyle
Y_{k}:=-\operatorname{Tr}\left[\left(\frac{\partial\rho(\vec{\theta})}{\partial\theta_{k}}\right)^{\dagger}(\\{H,\rho(\vec{\theta}(t))\\}-2\langle
H\rangle\rho(\vec{\theta}(t)))\right]$ (8)

It is important to note that when $\vec{\theta}$ is real, $M$ and $Y$ are also
real and $Y$ will simplify to

$\displaystyle
Y_{k}=-\operatorname{Tr}\left[\frac{\partial\rho(\vec{\theta}(t))}{\partial\theta_{k}}\\{H,\rho(\vec{\theta}(t))\\}\right]$
(9)

To implement the imaginary time evolution for mixed states, we need to measure
the coefficients $M_{k,q}$ and $Y_{k}$. We utilize the general circuit shown
in Fig. 2 to measure $\frac{\partial\rho}{\partial\theta_{k}}$. The generator
of the circuit can be expressed as

$\displaystyle\mathcal{L}(\rho)=\sum_{k}g_{k}S_{k}\rho T_{k}^{\dagger}$ (10)

where $S_{k}$ and $T_{k}$ are unitary operators and $g_{k}$ are coefficients.
Using the above we can write

$\displaystyle\frac{\partial\mathcal{R}_{k}(\rho)}{\partial\lambda_{k}}=\sum_{i}r_{k,i}S_{k,i}\rho
T_{k,i}^{\dagger}.$ (11)

From there we can rewrite the derivatives
$\frac{\partial\rho}{\partial\theta_{k}}$ to obtain the following expressions
for $M$ and $Y$ [21]111Note two minor corrections in Fig. 3 of Ref. [21].
First one is that $U_{k}:=S_{k}^{\dagger}T_{k}$, and the second is that in the
second system, $T_{q}$ will come after $U_{q}$ and not $S_{q}$ (as given in
Fig. 3 of [21]) as the definition of $M_{k,q}$ involves taking the conjugate
of one term.:

$\displaystyle\begin{split}M_{k,q}&=\sum_{i,j}\frac{1}{2}\\{r^{*}_{k,i}r_{q,j}{\rm
Tr}\left[\left(\mathcal{R}_{k,i}\rho_{0}\right)^{\dagger}\left(\mathcal{R}_{q,j}\rho_{0}\right)\right]+\rm{h.c.},\\\
Y_{k}&=-\sum_{i}\frac{1}{2}\\{r_{k,i}{\rm
Tr}\left[(\mathcal{R}_{k,i}\rho_{0})(\rho H)\right]\\}\end{split}$ (12)

where $\rho$ is the final state given by

$\displaystyle\rho=\mathcal{R}_{N}\circ\mathcal{R}_{N-1}\circ\cdots\circ\mathcal{R}_{1}(\rho_{0})$
(13)

and

$\displaystyle\mathcal{R}_{k,i}\rho_{0}=\mathcal{R}_{N}\cdots\mathcal{R}_{k+1}[S_{k,i}\left(\mathcal{R}_{k-1}\cdots\mathcal{R}_{2}\mathcal{R}_{1}(\rho_{0})\right)T_{k,i}^{\dagger}]\,.$
(14)

Further, to simulate the real time evolution the same setup can be used [21].
The main difference in the real and imaginary time evolution is in the matrix
equation to be solved where $Y$ is replaced by $V$ for the real evolution,
i.e., the linear equation for matrix inversion becomes $M\dot{\theta}=V$ where
V is given by

$\displaystyle
V_{k}=\sum_{i,j}\frac{1}{2}\left\\{r_{k,i}g_{j}\operatorname{Tr}\left[\left(\mathcal{R}_{k,i}\rho_{0}\right)\left(S_{j}\rho
T_{j}^{\dagger}\right)\right]+\text{ h.c. }\right\\}.$ (15)

##  III Error Models

We consider the errors of the following probabilistic form:

$\displaystyle\mathcal{E}(\rho):=p\mathcal{N}(\rho)+(1-p)\rho$ (16)

where $\mathcal{N}$ is some quantum channel, representing the hardware error
occurring with probability $p$ in the quantum circuit. We assume that the
error $\mathcal{E}$ occurs after every gate implementation and is a global
Markovian error (see Fig. 2). Considering such general Markovian error allows
us to estimate the error scaling in the worst-case scenario. Such error models
are commonly studied in literature, however most studies include local errors
on all qubits and not global [34, 35, 32]. By considering global errors, we
also account for cases where noise model could be correlated. Note that in
this article, we do not consider temporal correlations, and thus have avoided
all non-Markovian errors in this analysis. We further impose a physical
condition that the channel $\mathcal{N}\in{\rm CPTP}(A_{0}\to A_{1})$ is such
that $\|\mathcal{N}(X)\|_{2}\leq\|X\|_{2}\;\forall\,X\in\mathfrak{L}(A_{0})$,
where $\|\cdot\|_{2}$ denotes the Frobenius norm 222We denote $\|\cdot\|_{2}$
and $\|\cdot\|$ interchangeably to indicate the Frobenius norm defined as
$\|A\|_{2}=\sqrt{{\rm Tr}(A^{*}A)}$, $\mathfrak{L}(A_{0})$ denotes the set of
all linear operators in the Hilbert space $A_{0}$, and CPTP$(A_{0}\to
A_{1})$333We denote a system and its corresponding Hilbert space using an
uppercase letter with a numerical subscript, like $A_{0}$ or $A_{1}$. denotes
the set of completely positive and trace preserving maps that map operators in
the Hilbert space $A_{0}$ to operators in the Hilbert space $A_{1}$. With this
condition, we also leave the algorithmic analyses due to the presence of
replacement channels along with non-Markovian noise [38, 39, 40, 41] as an
open problem.

Figure 2: Error setting in our studies. The above circuit evaluates the
elements of the matrix $M$ (see Sec. II for the form of $M_{k,q}$ element)
with global error $\mathcal{E}$ occurring after every gate where
$\mathcal{E}(\rho):=p\mathcal{N}(\rho)+(1-p)\rho$. The above circuit can be
decomposed into three parts. The first part/line is a single qubit (in the
state $|0\rangle+e^{i\theta}|1\rangle$) that controls the circuits below. The
second and third parts/lines each have the same set of gates and the same
number of multiple qubits, and represent the same ansatz circuits. The
difference in the second and third lines at some sites (here, $k$-th and
$q$-th site) arise because of the differentiation with respect to the
respective parameter at those sites which is required for computing the matrix
element. For representing the differentiation with respect to the respective
parameter, the gate at that site is replaced by an appropriate controlled-
unitary followed by appropriate unitary at that site. In addition, we consider
ansatz circuits with only single and two-qubit gates. We are assuming that the
initial state preparation, the first control qubit, and the measurements are
error-free (fuzzy-free in the diagram). For each application of
controlled-$U_{k}$($U_{q}$) followed by $S_{k}$($T_{q}$) as the $k$($q$)-th
gate in the circuit, the error is applied after the combined gates (red dashed
box).

Notice that the set of all channels $\mathcal{N}$ that obey
$\|\mathcal{N}(X)\|_{2}\leq\|X\|_{2}\;\forall\,X\in\mathfrak{L}(A_{0})$
contain the set of all probabilistic unitary channels, i.e., channels of the
form

$\displaystyle\mathcal{N}(\rho)=\sum_{i}p_{i}\mathcal{U}_{i}(\rho)=\sum_{i}p_{i}U_{i}\rho
U_{i}^{\dagger}\,$ (17)

where $p_{i}\geq 0$ and $\sum_{i}p_{i}=1$. By doing so, we are accounting for
all the probabilistic unitary errors that contain the set of random Pauli
errors. The special cases of which are the (partially and completely)
dephasing and depolarizing channels, which are considered very relevant
experimentally [42, 43].

##  IV Error Bounds

By assuming the initial state preparation and measurements to be error-free
(see Fig. 2) and by considering gate errors as described above in Section III,
one can obtain upper bounds on the error $\epsilon$ (defined as the norm of
the difference between the ideal and erroneous $\dot{\theta}$ in the exact
solution of $M\dot{\theta}=Y$) using fundamental inequalities of perturbation
theory (see Eq. (A42) of the appendix). However, the bound obtained this way
is very loose and blows out very quickly, thereby giving any sensible bound
only for very small probabilities of noise.

(a)

(b)

Figure 3: Plots for Theorem 1 results. In the above plots, we fixed the
condition number of $M$ and $\|M\|$ to be $66.7$ and $0.9977$, respectively,
and varied the number of parameters. These values were obtained by using the
five-parameter ansatz considered for simulating the Hamiltonian of the
Hydrogen molecule and by generating random angles (see Appendix C for
details). Figure (a) is the plot of the constraint in Eq. (19). In this graph,
we plotted the allowed maximum probability with respect to fixed $\delta$
which dictates when $\Delta M\Delta\dot{\theta}$ can be considered negligible.
For each given $N$, the graph stops after a certain $\delta$; it shows that
beyond that $\delta$ the right hand side in Eq. (19) will be negative, thus
leading to invalid bound. To plot Fig. (b), we choose $\delta=0.04$ for which
we calculated the allowed probabilities to be $0.032,0.036,0.052,0.068,0.129$
for $N=5,6,10,12,14$, respectively. Figure (b) is the plot of the relative
error in $\dot{\theta}$ as a function of probability of noise. The plot stops
at the maximum allowed probability as governed from the constraint.

Tight upper bounds can be found by imposing stricter conditions on solving the
linear equations or by knowing more details of the noise. We consider
instances of both these cases and find tight upper bounds.

When we solve the perturbed linear equation $M\dot{\theta}=Y$, we are
essentially solving the equation $(M+\Delta
M)(\dot{\theta}+\Delta\dot{\theta})=Y+\Delta Y$. Now for a fixed very small
$\delta<<1$, if each element of $\Delta M$ and $\Delta\dot{\theta}$ is less
than $\delta$, i.e., $(\Delta M)_{i,j}\leq\delta$ and
$(\dot{\theta})_{j}\leq\delta$, then $(\Delta M\Delta\dot{\theta})_{i}\sim
O(\delta^{2})$ and can safely be ignored as its contribution to $\Delta Y$ can
be considered negligible. By slightly loosening this restriction, we can ask
for $\|\Delta M\|\leq N^{2}\delta$ and $\|\Delta\dot{\theta}\|\leq N\delta$ to
hold for the contribution of $\Delta M\Delta\dot{\theta}$ to be considered
negligible. Considering these constraints, we get the following tight upper
bound for the general noise model.

######  Theorem 1.

For probabilistic errors $\mathcal{E}$ of the form:
$\mathcal{E}(\rho)=p\mathcal{N}(\rho)+(1-p)\rho$, where the channel
$\mathcal{N}$ is such that
$\|\mathcal{N}(X)\|_{2}\leq\|X\|_{2}\;\forall\,X\in\mathfrak{L}(\mathcal{H})$
where $\mathfrak{L}$ denotes the set of linear operators in a Hilbert space
$\mathcal{H}$, we show that the upper bound in the error
$\epsilon:=\|\Delta\dot{\theta}\|$ is given by

$\displaystyle\frac{\epsilon}{\|\dot{\theta}\|}\leq{\rm
cond}(M)\left(\left(1-(1-p)^{2N}\right)\left(1+\frac{N}{2\|M\|}\right)+(1-(1-p)^{N})\left(1+\frac{\sqrt{N}}{\sqrt{2\|Y\|}}\right)\right)$
(18)

for the probability of noise $p$ in the following range

$\displaystyle 0\leq p\leq
1-\left(1-\frac{2N^{2}\delta}{N+2\|M\|}\right)^{\frac{1}{2N}}$ (19)

where $\delta$ in the term on the right is small enough such that the
expression is positive.

The proof of the above theorem has been provided in Appendix A. Note that
beyond the probability range provided in Eq. (19), the contribution of $\Delta
M\Delta\dot{\theta}$ to $\Delta Y$ cannot be considered negligible. We also
numerically demonstrate the constraint of Theorem 1 in Fig. 3(a), where we
show how the allowed probability (for which the main result of Theorem 1 is
valid) changes with varying $\delta$ and number of parameters, $N$. By fixing
$\delta=0.04$, we show in Fig. 3(b) how the relative error scales with the
probability of noise in quantum hardware. Note that for smaller $N$, we can
have larger $\delta$ which will lead to larger allowed probabilities.

Now let us consider the evolution of quantum state in QITE under a specific
noise model, the partially depolarizing noise, which is a major source of
error in NISQ computers [43]. For the depolarizing noise, we use its symmetry
to find the relation between the relative error, the probability of noise, and
the number of parameters.

######  Theorem 2.

For partially depolarizing channels, the upper bound in the error
$\epsilon:=\|\Delta\dot{\theta}\|$ is given by

$\displaystyle\frac{\epsilon}{\|\dot{\theta}\|}=\frac{1-(1-p)^{N}}{(1-p)^{N}}.$
(20)

The proof of this theorem can be found in Appendix B. There are several key
features to note here. First, we now have an exact expression on how the error
in the solution of $M\dot{\theta}=Y$ will scale as opposed to Theorem 1 where
we only had an upper bound. The above exact expression is because (in Appendix
B) we were able to derive exactly the form of $M^{err}$ as a function of $M$
(see Eq. (B60)). Moreover, the relative error under depolarizing noise only
depends on the probability of depolarizing noise and the number of parameters
$N$ and is independent of the vector $Y$ or the matrix $M$, i.e. it is
independent of the Hamiltonian being simulated or the structure of the ansatz
circuit.

It is clear from the expression that as the probability of noise increases,
the error also increases. This leads us to conclude that QITE can perform well
only when the probability of noise is very low, or in other words, under
depolarizing noise, QITE in fact requires very high fidelity quantum gates
without any error mitigation deployed. Moreover, analyzing the behaviour of
error scaling under the depolarizing noise, we can conjecture that as the
probability increases (see Fig. 1(b)), the error scaling for general noise
model must also scale up rapidly with increasing probability of noise. Lastly,
given the exact expression of relative error, we can perform a simple post-
processing (mitigation) to get the accurate $\dot{\theta}$.

From the two theorems above, we see that as the noise probability after each
gate increases, the relative error $\epsilon/\|\dot{\theta}\|$ grows very
rapidly or diverges. This means that as the noise probability increases, the
new parameters to be updated that are required for the next run of the QITE
will be further away from the actual parameters and so the error in estimating
the ground state of a given Hamiltonian will increase. Also note that if we
consider circuits where multiple unitary gates are applied onto different
qubit registers simultaneously, and the probabilistic error is applied after
each such layer, then the expressions of relative errors are similar to the
results in Theorem 1 and 2 except that the number of parameters, $N$, in the
equations will be replaced by the depth, $d$, of the circuit.

In addition, we perform numerical analyses of our results by considering the
paradigmatic Hydrogen molecule H2 in the minimal STO-3G basis. The reduced H2
qubit Hamiltonian in the Bravyi-Kitaev representation is given by [16]

$\displaystyle
H=g_{0}I+g_{1}Z_{0}+g_{2}Z_{1}+g_{3}Z_{0}Z_{1}+g_{4}Y_{0}Y_{1}+g_{5}X_{0}X_{1}\,$
(21)

where $g_{0}=0.2252$, $g_{1}=0.3435$, $g_{2}=-0.4347$, $g_{3}=0.5716$,
$g_{4}=0.0910$, and $g_{5}=0.0910$. In our analysis, we make use of two
ansatzes to understand our error bounds. One is a simple ansatz and the other
ansatz is the universal ansatz [12] as shown in Figs. 4(a) and 4(b) of the
Appendix, respectively. Notice that the simple ansatz has 5 parameters and the
universal ansatz has 9 parameters. Using these ansatzes, we then plot the
error vs probability figures throughout. For example, Fig. 1(b) was obtained
by using the 5-parameter simple ansatz circuit undergoing depolarizing noise
to get the reference line $(N=5,$ black solid line$)$ with random angles to
generate one matrix $M$ with a reference condition number, cond$(M)$. While
fixing the same cond$(M)$ and $\|M\|$, we can now start to vary the number of
parameters $N$. By doing so, we effectively change the structure of the ansatz
circuit. Alternatively, when $N<5$, we can set some $\theta=0$ and vary the
other parameters to get the same condition number while using the same ansatz
circuit. The same strategy is applied for the general probabilistic errors
(Fig. 3(b) of the main text and Fig. 6 of Appendix C), and for the
depolarizing error (Fig. 1(b) of the main text or Fig. 5 of Appendix C). From
all of the presented numerical results, we find that the error in updating the
variational parameters increases drastically with an increase in error
probability in the physical hardware. Since depolarizing noise commonly occurs
in superconducting qubits architecture, we expect that such behaviour could
become a bottleneck in performing QITE in high-depth quantum circuits without
appropriate post processing or error mitigation.

Finally, it is important to note that one does not have to solve,
$M\dot{\theta}=Y$ via matrix inversion or any method in particular. One can
use any numerical method to solve the linear system of equations. Our goal is
to understand how the relative error in $\dot{\theta}$ scales with the
probability of noise in quantum hardware. In doing so, we found the upper
bound of relative error and showed that it is directly proportional to the
condition number of $M$ (see Eq.(17)). The condition number tells us how
sensitive the solution of a linear solver is given a slight perturbation in
the matrix $M$ or vector $Y$. Now, concerning the computation of condition
number of $M$, computing the norm of the inverse/pseudo inverse of $M$ might
be difficult, specially for large matrices. However, it is possible to
approximate the condition number of a matrix without computing the norm of its
inverse [44], thus resolving the issue of computing the inverse/condition
number of a matrix.

##  V Discussion

In this work we presented upper bounds on the error incurred during a single
iteration of quantum imaginary time evolution in the presence of global
Markovian noise in the quantum circuit. We assumed that the error can occur
after every gate in the ansatz circuit with a probability $p$. Furthermore, we
derived an exact expression for the error incurred specifically due to
partially depolarizing channels. We note that these error analyses are also
applicable to real-time quantum data evolution. This is because the
variational problem for real time evolution is similar since the matrix M in
the linear problem of solving for $\dot{\theta}$ is exactly the same. i.e.,
the results from the two theorems would not differ except to a different
dependence on the norm of $V$ (equation (15)). Our results primarily show that
to implement the variational quantum data evolution both real- and imaginary-
time evolution, we would require gates with very low error probability due to
the sensitivity of matrix inversion to external environmental perturbations.
We also show that scalability is a problem with this approach because of how
the error scales with the number of parameters in the circuit and the
probability of gate errors. In fact, when we consider multiple unitary gates
applying onto different qubit registries, $N$, the notation we used for the
number of variational parameters, can be taken as the quantum circuit depth.
Hence, the results in Fig. 1 still hold with $N$ then being the depth of the
circuit. Furthermore, our results can be used to compare different ansatzes in
order to choose a well-suited ansatz for the given problem. In other words, if
the amount of error that can be tolerated is fixed, one can find the maximum
probability of allowed noise in the system and thus, relevant error mitigation
techniques can be employed in order to restrict error probability going beyond
a certain threshold.

This work is important considering the rising interest in extracting quantum
advantage from NISQ algorithms and processors and which makes it critical to
do an in-depth error analysis in order to truly understand whether the NISQ
algorithms really stand a chance in showing quantum advantage. We found that
the error analyses performed on QITE algorithm, which as we stated previously
is one of the most efficient NISQ algorithms to find the ground state of a
Hamiltonian, were either too basic, or assumed very strict conditions, or were
based on simple numerical simulations [27, 28, 29, 30, 16, 32]. For instance,
in Ref. [30], the error analysis presented was done by considering two
different evolutions of the state and how those two evolutions diverged. No
hardware noise was considered in getting the theoretical results. In Ref.
[16], no specific hardware noise or even general error model was considered.
By assuming $p=10^{-4}$ error rate per gate, the authors numerically concluded
that “the imaginary time algorithm can perform significantly better than the
gradient descent”. However, our analysis is much more general. To show this
let us take the same numbers as considered in Ref. [16], and consider that the
probability of depolarizing noise after every gate is $10^{-4}$. Our results
(shown in Fig.1(b)) suggest that the error incurred would be very small and
so, with this noise probability, we can infer that the imaginary time
evolution algorithm will perform well - which is in agreement with the results
of Ref. [16].

Finally, it is yet to be seen how various quantum error mitigation techniques
[45, 46] would help to suppress error in variational parameters propagation
from solving a linear system of equations as we have presented here. We
conjecture that the standard zero noise extrapolation method [47, 48, 49, 50]
would fail due to the need of multiple quantum circuit layers for its
extrapolation technique to work. However, we expect probabilistic error
cancellation (PEC) [47, 51, 52, 53] might help to reduce the growth of errors
in the evolution of the circuit parameters. Exact details on the amount of
error suppression with respect to the overhead needed for PEC are left for a
future work.

## Code Availability

Code and data generated to construct the plots will be made available upon
request.

## Acknowledgements

We thank Kevin Ferreira (Senior Director of LG Electronics Toronto AI Lab) for
his constant support throughout this work.

## References

  * Shor [1994] P. Shor, Algorithms for quantum computation: discrete logarithms and factoring, in  _Proceedings 35th Annual Symposium on Foundations of Computer Science_ (Santa Fe, NM, USA, 1994) pp. 124–134.
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum Algorithm for Linear Systems of Equations, Phys. Rev. Lett. 103, 150502 (2009).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, _et al._ , Quantum Chemistry in the Age of Quantum Computing, Chem. Rev. 119, 10856 (2019).
  * Lee  _et al._ [2021] J. Lee, D. W. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and R. Babbush, Even More Efficient Quantum Computations of Chemistry Through Tensor Hypercontraction, PRX Quantum 2, 030305 (2021).
  * Lin and Tong [2022] L. Lin and Y. Tong, Heisenberg-Limited Ground-State Energy Estimation for Early Fault-Tolerant Quantum Computers, PRX Quantum 3, 010318 (2022).
  * Hoefler  _et al._ [2023] T. Hoefler, T. Häner, and M. Troyer, Disentangling hype from practicality: On realistically achieving quantum advantage, Commun. ACM 66, 82–87 (2023).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, _et al._ , Noisy intermediate-scale quantum algorithms, Rev. Mod. Phys. 94, 015004 (2022).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, _et al._ , The Variational Quantum Eigensolver: A review of methods and best practices, Phys. Rep. 986, 1 (2022).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Variational quantum algorithms, Nat. Rev. Phys. 3, 625 (2021).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5, 1 (2014).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New J. Phys. 18, 023023 (2016).
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Phys. Rev. A 92, 042303 (2015).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nat. Phys. 16, 205 (2020).
  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, Variational ansatz-based quantum simulation of imaginary time evolution, npj Quantum Inf. 5, 1 (2019).
  * Nishi  _et al._ [2021] H. Nishi, T. Kosugi, and Y.-i. Matsushita, Implementation of quantum imaginary-time evolution method on NISQ devices by introducing nonlocal approximation, npj Quantum Inf. 7, 1 (2021).
  * Foulkes  _et al._ [2001] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Quantum monte carlo simulations of solids, Rev. Mod. Phys. 73, 33 (2001).
  * Yeter-Aydeniz  _et al._ [2020] K. Yeter-Aydeniz, R. C. Pooser, and G. Siopsis, Practical quantum computation of chemical and nuclear energy levels using quantum imaginary time evolution and Lanczos algorithms, npj Quantum Inf. 6, 1 (2020).
  * Tsuchimochi  _et al._ [2023] T. Tsuchimochi, Y. Ryo, S. L. Ten-no, and K. Sasasako, Improved algorithms of quantum imaginary time evolution for ground and excited states of molecular systems, Journal of Chemical Theory and Computation 19, 503 (2023).
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, Theory of variational quantum simulation, Quantum 3, 191 (2019), 1812.08767v4 .
  * Endo  _et al._ [2020] S. Endo, J. Sun, Y. Li, S. C. Benjamin, and X. Yuan, Variational quantum simulation of general processes, Phys. Rev. Lett. 125, 010501 (2020).
  * Gomes  _et al._ [2020] N. Gomes, F. Zhang, N. F. Berthusen, C.-Z. Wang, K.-M. Ho, P. P. Orth, and Y. Yao, Efficient Step-Merged Quantum Imaginary Time Evolution Algorithm for Quantum Chemistry, J. Chem. Theory Comput. 16, 6256 (2020).
  * Wick [1954] G. C. Wick, Properties of bethe-salpeter wave functions, Phys. Rev. 96, 1124 (1954).
  * Lee  _et al._ [2023] S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, M. Kastoryano, R. Babbush, J. Preskill, D. R. Reichman, E. T. Campbell, E. F. Valeev, L. Lin, and G. K.-L. Chan, Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry, Nature Communications 14, 1952 (2023).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017).
  * Kyaw  _et al._ [2023] T. H. Kyaw, M. B. Soley, B. Allen, P. Bergold, C. Sun, V. S. Batista, and A. Aspuru-Guzik, Boosting quantum amplitude exponentially in variational quantum algorithms, Quantum Sci. Technol. 9, 01LT01 (2023).
  * Zhang  _et al._ [2023] Z.-J. Zhang, J. Sun, X. Yuan, and M.-H. Yung, Low-Depth Hamiltonian Simulation by an Adaptive Product Formula, Phys. Rev. Lett. 130, 040601 (2023).
  * Benedetti  _et al._ [2021] M. Benedetti, M. Fiorentini, and M. Lubasch, Hardware-efficient variational quantum algorithms for time evolution, Phys. Rev. Res. 3, 033083 (2021).
  * Zoufal  _et al._ [2023] C. Zoufal, D. Sutter, and S. Woerner, Error bounds for variational quantum time evolution, Phys. Rev. Appl. 20, 044059 (2023).
  * Stilck França and García-Patrón [2021] D. Stilck França and R. García-Patrón, Limitations of optimization algorithms on noisy quantum devices, Nat. Phys. 17, 1221 (2021).
  * Wang  _et al._ [2021] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Noise-induced barren plateaus in variational quantum algorithms, Nat. Commun. 12, 1 (2021).
  * Note [1] Note two minor corrections in Fig. 3 of Ref. [21]. First one is that $U_{k}:=S_{k}^{\dagger}T_{k}$, and the second is that in the second system, $T_{q}$ will come after $U_{q}$ and not $S_{q}$ (as given in Fig. 3 of [21]) as the definition of $M_{k,q}$ involves taking the conjugate of one term.
  * Xue  _et al._ [2021] C. Xue, Z.-Y. Chen, Y.-C. Wu, and G.-P. Guo, Effects of Quantum Noise on Quantum Approximate Optimization Algorithm, Chin. Phys. Lett. 38, 030302 (2021).
  * Marshall  _et al._ [2020] J. Marshall, F. Wudarski, S. Hadfield, and T. Hogg, Characterizing local noise in QAOA circuits, IOP SciNotes 1, 025208 (2020).
  * Note [2] We denote $\|\cdot\|_{2}$ and $\|\cdot\|$ interchangeably to indicate the Frobenius norm defined as $\|A\|_{2}=\sqrt{{\rm Tr}(A^{*}A)}$.
  * Note [3] We denote a system and its corresponding Hilbert space using an uppercase letter with a numerical subscript, like $A_{0}$ or $A_{1}$.
  * Bastidas  _et al._ [2018] V. M. Bastidas, T. H. Kyaw, J. Tangpanitanon, G. Romero, L.-C. Kwek, and D. G. Angelakis, Floquet stroboscopic divisibility in non-Markovian dynamics, New J. Phys. 20, 093004 (2018).
  * Kyaw  _et al._ [2017] T. H. Kyaw, S. Allende, L.-C. Kwek, and G. Romero, Parity-preserving light-matter system mediates effective two-body interactions, Quantum Sci. Technol. 2, 025007 (2017).
  * Gour [2019] G. Gour, Comparison of Quantum Channels by Superchannels, IEEE Trans. Inf. Theory 65, 5880 (2019).
  * Regula and Takagi [2021] B. Regula and R. Takagi, Fundamental limitations on distillation of quantum channel resources, Nat. Commun. 12, 1 (2021).
  * Mi  _et al._ [2022] X. Mi, M. Sonner, M. Y. Niu, K. W. Lee, B. Foxen, R. Acharya, I. Aleiner, T. I. Andersen, F. Arute, K. Arya, _et al._ , Noise-resilient edge modes on a chain of superconducting qubits, Science 378, 785 (2022).
  * Georgopoulos  _et al._ [2021] K. Georgopoulos, C. Emary, and P. Zuliani, Modeling and simulating the noisy behavior of near-term quantum computers, Phys. Rev. A 104, 062432 (2021).
  * von Petersdorff [2018] T. von Petersdorff, Errors for linear systems, https://www.math.umd.edu/~petersd/466/linsysterrn.pdf/ (2018).
  * Cai  _et al._ [2023] Z. Cai, R. Babbush, S. C. Benjamin, S. Endo, W. J. Huggins, Y. Li, J. R. McClean, and T. E. O’Brien, Quantum error mitigation, Rev. Mod. Phys. 95, 045005 (2023).
  * Quek  _et al._ [2022] Y. Quek, D. S. França, S. Khatri, J. J. Meyer, and J. Eisert, Exponentially tighter bounds on limitations of quantum error mitigation, arXiv 10.48550/arXiv.2210.11505 (2022), 2210.11505 .
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error Mitigation for Short-Depth Quantum Circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient Variational Quantum Simulator Incorporating Active Error Minimization, Phys. Rev. X 7, 021050 (2017).
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Error mitigation extends the computational reach of a noisy quantum processor, Nature 567, 491 (2019).
  * Giurgica-Tiron  _et al._ [2020] T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, Digital zero noise extrapolation for quantum error mitigation, in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (IEEE, Denver, CO, USA, 2020) pp. 306–316.
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, Practical Quantum Error Mitigation for Near-Future Applications, Phys. Rev. X 8, 031027 (2018).
  * Song  _et al._ [2019] C. Song, J. Cui, H. Wang, J. Hao, H. Feng, and Y. Li, Quantum computation with universal error mitigation on a superconducting quantum processor, Sci. Adv. 5, eaaw5686 (2019).
  * Zhang  _et al._ [2020] S. Zhang, Y. Lu, K. Zhang, W. Chen, Y. Li, J.-N. Zhang, and K. Kim, Error-mitigated quantum gates exceeding physical fidelities in a trapped-ion system, Nat. Commun. 11, 1 (2020).
  * Higham [2002] N. J. Higham,  _Accuracy and Stability of Numerical Algorithms_, 2nd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2002).
  * Sim  _et al._ [2019] S. Sim, P. D. Johnson, and A. Aspuru-Guzik, Expressibility and Entangling Capability of Parameterized Quantum Circuits for Hybrid Quantum-Classical Algorithms, Adv. Quantum Technol. 2, 1900070 (2019).

##  Appendix A Proof of Theorem 1 (Upper bound for the general noise model)

Theorem 1 statement. For probabilistic errors $\mathcal{E}$ of the form:
$\mathcal{E}(\rho)=p\mathcal{N}(\rho)+(1-p)\rho$, where the channel
$\mathcal{N}$ is such that
$\|\mathcal{N}(X)\|_{2}\leq\|X\|_{2}\;\forall\,X\in\mathfrak{L}(\mathcal{H})$
where $\mathfrak{L}$ denotes the set of linear operators in a Hilbert space
$\mathcal{H}$, we show that the upper bound in the error
$\epsilon:=\|\Delta\dot{\theta}\|$ is given by

$\displaystyle\frac{\epsilon}{\|\dot{\theta}\|}\leq{\rm
cond}(M)\left(\left(1-(1-p)^{2N}\right)\left(1+\frac{N}{2\|M\|}\right)+(1-(1-p)^{N})\left(1+\frac{\sqrt{N}}{\sqrt{2\|Y\|}}\right)\right)$
(A1)

for the probability of noise $p$ in the following range

$\displaystyle 0\leq p\leq
1-\left(1-\frac{2N^{2}\delta}{N+2\|M\|}\right)^{\frac{1}{2N}}.$ (A2)

###### Proof.

We know that

$\displaystyle\begin{split}M_{k,q}^{\rm err}&=\sum_{i,j}r^{*}_{k,i}r_{q,j}{\rm
Tr}\left[\left(\mathcal{R}_{k,i}\rho_{0}\right)^{\dagger}\left(\mathcal{R}_{q,j}\rho_{0}\right)\right],\;{\rm
and}\\\ Y_{k}^{\rm err}&=-\sum_{i}\\{r_{k,i}{\rm
Tr}\left[(\mathcal{R}_{k,i}\rho_{0})(\tilde{\rho}H+H\tilde{\rho})\right]\\}\end{split}$
(A3)

where $\tilde{\rho}$ denotes the erroneous $\rho$ that we get from the circuit
subjected to errors. Note that when $\theta$ is real, then $M_{k,q}$ and
$Y_{k}$ are also real.

For single-qubit gates, we have the following forms of $r$ coefficients and
the $S$ and $T$ matrices. Suppose we are differentiating at the $k$-th gate
with respect to the $k$-th parameter, then

$\displaystyle r_{k,1}=-\frac{i}{2}\,,$
$\displaystyle\;S_{k,1}=P_{k}R_{k}\,,\;T_{k,1}=R_{k}$ (A4) $\displaystyle
r_{k,2}=\frac{i}{2}\,,$ $\displaystyle\;S_{k,2}=R_{k}\,,\;T_{k,2}=P_{k}R_{k}$
(A5)

Note, in a multi-qubit circuit, $P_{k}$ is a compact notation of the form
$I\otimes\cdots\otimes I\otimes P\otimes I\otimes\cdots\otimes I$ where $P$ is
the Pauli matrix that occurs in the $k$-th rotation gate. If there is a two-
qubit gate at the $k$-th position, then the $r$ coefficients and the $S$ and
$T$ matrices would have the following form

$\displaystyle\begin{split}r_{k,1}=\frac{i}{2}\,,&\;S_{k,1}=|0\rangle\langle
0|\otimes I\,,\;T_{k,1}=|1\rangle\langle 1|\otimes P_{k}R_{k}\\\
r_{k,2}=-\frac{i}{2}\,,&\;S_{k,2}=|1\rangle\langle 1|\otimes
P_{k}R_{k}\,,\;T_{k,2}=|0\rangle\langle 0|\otimes I\\\
r_{k,3}=-\frac{i}{2}\,,&\;S_{k,3}=|1\rangle\langle 1|\otimes
P_{k}R_{k}\,,\;T_{k,3}=|1\rangle\langle 1|\otimes R_{k}\\\
r_{k,4}=\frac{i}{2}\,,&\;S_{k,4}=|1\rangle\langle 1|\otimes
R_{k}\,,\;T_{k,4}=|1\rangle\langle 1|\otimes R_{k}P_{k}\end{split}$ (A6)

or equivalently,

$\displaystyle\begin{split}r_{k,1}=\frac{i}{2}\,,\;&\;S_{k,1}=|0\rangle\langle
0|\otimes I+|1\rangle\langle 1|\otimes R_{k}\,,\;T_{k,1}=|1\rangle\langle
1|\otimes P_{k}R_{k}\\\ r_{k,2}=-\frac{i}{2}\,,\;&\;S_{k,2}=|1\rangle\langle
1|\otimes P_{k}R_{k}\,,\;T_{k,2}=|0\rangle\langle 0|\otimes I+|1\rangle\langle
1|\otimes R_{k}\end{split}$ (A7)

When the gates at the $k$-th and the $q$-th are both single-qubit or two-
qubit, then,

$\displaystyle\begin{split}r^{*}_{k,1}r_{q,1}&=r^{*}_{k,2}r_{q,2}=1/4\,\;\textrm{\,
and,}\\\ r^{*}_{k,2}r_{q,1}&=r^{*}_{k,1}r_{q,2}=-1/4\,.\end{split}$ (A8)

If one of the two gates is a single-qubit gate and the other one is a two
qubit gate, then the signs gets reversed in the above equation.

First let us only consider single-qubit gates at the $k$-th and $q$-th sites
in the circuit. Then,

$\displaystyle\begin{split}{\rm
Tr}[(\mathcal{R}_{k,i}\rho_{0})^{\dagger}(\mathcal{R}_{q,j}\rho_{0})]={\rm
Tr}\Big{[}&\Big{(}p^{N-k+1}\mathcal{D}\circ\mathcal{R}_{N}\circ\mathcal{D}\circ\cdots\mathcal{D}\mathcal{R}_{k+1}\mathcal{D}(T_{k,i}\tilde{\rho}_{k-1}S^{\dagger}_{k,i})\\\
&\;+p^{N-k}(1-p)({\rm
all\;}N-k{\rm\;combinations\;of\;}\mathcal{D}{\rm\;after\;gates})(T_{k,i}\tilde{\rho}_{k-1}S^{\dagger}_{k,i})\\\
&\;+p^{N-k-1}(1-p)^{2}({\rm
all\;}N-k-1{\rm\;combinations\;of\;}\mathcal{D}{\rm\;after\;gates})(T_{k,i}\tilde{\rho}_{k-1}S^{\dagger}_{k,i})\\\
&\;+\cdots+(1-p)^{N-k+1}(\mathcal{R}_{N}\cdots\mathcal{R}_{k+1}(T_{k,i}\tilde{\rho}_{k-1}S^{\dagger}_{k,i}))\Big{)}\\\
&\Big{(}p^{N-q+1}\mathcal{D}\circ\mathcal{R}_{N}\circ\mathcal{D}\circ\cdots\mathcal{D}\mathcal{R}_{q+1}\mathcal{D}(S_{q,j}\tilde{\rho}_{q-1}T^{\dagger}_{q,j})\\\
&\;+p^{N-q}(1-p)({\rm
all\;}N-q{\rm\;combinations\;of\;}\mathcal{D}{\rm\;after\;gates})(S_{q,j}\tilde{\rho}_{q-1}T^{\dagger}_{q,j})\\\
&\;+p^{N-q-1}(1-p)^{2}({\rm
all\;}N-q-1{\rm\;combinations\;of\;}\mathcal{D}{\rm\;after\;gates})(S_{q,j}\tilde{\rho}_{k-1}T^{\dagger}_{q,j})\\\
&\;+\cdots+(1-p)^{N-q+1}(\mathcal{R}_{N}\cdots\mathcal{R}_{q+1}(S_{q,j}\tilde{\rho}_{q-1}T^{\dagger}_{q,j}))\Big{)}\Big{]}\end{split}$
(A9)

where $\tilde{\rho}_{k-1}$ has the following form (expression of
$\tilde{\rho}_{q-1}$ will follow similarly)

$\displaystyle\begin{split}\tilde{\rho}_{k-1}&=p^{k-1}\mathcal{D}\circ\mathcal{R}_{k-1}\circ\mathcal{D}\circ\mathcal{R}_{k-2}\circ\mathcal{D}\circ\cdots\circ\mathcal{R}_{2}\circ\mathcal{D}\mathcal{R}_{1}(\rho_{0})\\\
&\;+p^{k-2}(1-p)(({\rm
sum\;of\;all\;possible\;}k-2{\rm\;combinations\;of\;}\mathcal{D}{\rm\;after\;gates})(\rho_{0}))\\\
&\;+p^{k-3}(1-p)^{2}(({\rm
sum\;of\;all\;possible\;}k-3{\rm\;combinations\;of\;}\mathcal{D}{\rm\;after\;gates})(\rho_{0}))\\\
&\;+\cdots+(1-p)^{k-1}\mathcal{R}_{k-1}\cdots\mathcal{R}_{1}(\rho_{0})\,.\end{split}$
(A10)

We know that

$\displaystyle M_{k,q}^{\rm err}=\frac{1}{4}\left({\rm
Tr}[(\mathcal{R}_{k,1}\rho_{0})^{\dagger}(\mathcal{R}_{q,1}\rho_{0})]-{\rm
Tr}[(\mathcal{R}_{k,1}\rho_{0})^{\dagger}(\mathcal{R}_{q,2}\rho_{0})]-{\rm
Tr}[(\mathcal{R}_{k,2}\rho_{0})^{\dagger}(\mathcal{R}_{q,1}\rho_{0})]+{\rm
Tr}[(\mathcal{R}_{k,2}\rho_{0})^{\dagger}(\mathcal{R}_{q,2}\rho_{0})]\right)$
(A11)

which can be expressed as

$\displaystyle\begin{split}M_{k,q}^{\rm
err}=\frac{1}{4}\Big{(}&\sum_{l}p_{l}{\rm
Tr}[\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})]\\\
&\;+(1-p)^{2N}{\rm
Tr}\left[\left(\mathcal{R}_{N}\cdots\mathcal{R}_{k+1}(\rho_{k}P_{k}-P_{k}\rho_{k})\right)\left(\mathcal{R}_{N}\cdots\mathcal{R}_{q+1}(P_{q}\rho_{q}-\rho_{q}P_{q})\right)\right]\Big{)}\end{split}$
(A12)

where $p_{l}$ are probabilities which are of the form $p^{a}(1-p)^{b}$ such
that $\sum_{l}p_{l}=1-(1-p)^{2N}$. The first term in the RHS of the Eq. (A12)
is a general expression of writing the terms where $\mathcal{D}$ occurs
atleast once either in $\mathcal{N}_{l}$, $\mathcal{M}_{l}$,
$\tilde{\rho}_{{}_{k,l}}$, and $\tilde{\rho}_{{}_{q,l}}$. That is,
$\mathcal{N}_{l}$ is of the form
$\mathcal{D}^{a_{N}}\circ\mathcal{R}_{N}\circ\mathcal{D}^{a_{N-1}}\circ\mathcal{R}_{N-1}\circ\cdots\circ\mathcal{D}^{a_{k+1}}\circ\mathcal{R}_{k+1}$
where $a_{i}={0,1}$, and $\tilde{\rho}_{{}_{k,l}}$ has the form
$\mathcal{D}^{b_{k-1}}\circ\mathcal{R}_{k-1}\circ\mathcal{D}^{b_{k-1}}\circ\mathcal{R}_{k-2}\circ\cdots\circ\mathcal{D}^{b_{1}}\mathcal{R}_{1}(\rho_{0})$
where $b_{i}={0,1}$ and the probabilities have been taken care of in $p_{l}$,
and similarly for $\mathcal{M}_{l}$ and $\tilde{\rho}_{{}_{q,l}}$. Using the
Cauchy-Schwarz inequality we can write

$\displaystyle\begin{split}{\rm
Tr}[\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}&-P_{k}\tilde{\rho}_{{}_{k,l}})\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})]\\\
&\leq\sqrt{{\rm
Tr}[\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\mathcal{N}_{l}(P_{k}\tilde{\rho}_{{}_{k,l}}-\tilde{\rho}_{{}_{k,l}}P_{k})]}\sqrt{{\rm
Tr}[\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})\mathcal{M}_{l}(\tilde{\rho}_{{}_{q,l}}P_{q}-P_{q}\tilde{\rho}_{{}_{q,l}})]}\\\
&=\|\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\|_{2}\|\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})\|_{2}\end{split}$
(A13)

The first inequality follows since in the LHS we are taking the product of two
anti-Hermitian matrices and if we multiply both the matrices by $i$ and take a
minus sign out of trace, we essentially get LHS to be equal to be less than or
equal to the trace of two Hermitian matrices and then when we apply the
Cauchy-Schwarz inequality, we get the above bound. Furthermore, it also holds
that

$\displaystyle{\rm
Tr}[\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\mathcal{M}_{l}(\tilde{\rho}_{{}_{q,l}}P_{q}-P_{q}\tilde{\rho}_{{}_{q,l}})]\leq\|\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\|_{2}\|\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})\|_{2}$
(A14)

and therefore,

$\displaystyle|{\rm
Tr}[\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})]|\leq\|\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\|_{2}\|\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})\|_{2}$
(A15)

Next, note that for any operator, say $X$, we have

$\displaystyle\|\mathcal{N}_{l}(X)\|_{2}\leq\|X\|_{2}=\|\mathcal{R}_{N}\circ\mathcal{R}_{N-1}\circ\cdots\circ\mathcal{R}_{k+1}(X)\|_{2}$
(A16)

Now, let us denote
$\mathcal{R}_{N}\circ\mathcal{R}_{N-1}\circ\cdots\circ\mathcal{R}_{k+1}$ as
$\mathcal{R}$, then have the following

$\displaystyle\|\mathcal{R}(X)\|_{2}$
$\displaystyle=\|\mathcal{R}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\|_{2}$
(A17)
$\displaystyle=\|\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}}\|_{2}$
(A18)
$\displaystyle=\|P_{k}\tilde{\rho}_{{}_{k,l}}P_{k}-\tilde{\rho}_{{}_{k,l}}\|_{2}$
(A19) $\displaystyle\leq\sqrt{2}$ (A20)

Therefore, for single-qubit gates (at $k$-th and $q$-th sites) we have

$\displaystyle|{\rm
Tr}[\mathcal{N}_{l}(\tilde{\rho}_{{}_{k,l}}P_{k}-P_{k}\tilde{\rho}_{{}_{k,l}})\mathcal{M}_{l}(P_{q}\tilde{\rho}_{{}_{q,l}}-\tilde{\rho}_{{}_{q,l}}P_{q})]|\leq
2\,.$ (A21)

Now, let us consider that the gates at the $k$-th site and $q$-site are two-
qubit gates. Let $W_{k}=|0\rangle\langle 0|\otimes I+|1\rangle\langle
1|\otimes R_{k}$ and $W_{q}=|0\rangle\langle 0|\otimes I+|1\rangle\langle
1|\otimes R_{q}$ be two unitary matrices. Then $M^{\rm err}_{k,q}$ can be
expressed as follows

$\displaystyle\begin{split}M^{\rm err}_{k,q}&=\frac{1}{4}\left({\rm
Tr}[(\mathcal{R}_{k,1}\rho_{0}-\mathcal{R}_{k,2}\rho_{0})^{\dagger}(\mathcal{R}_{q,1}\rho_{0}-\mathcal{R}_{q,2}\rho_{0})]\right)\\\
&=\frac{1}{4}\Big{(}\sum_{l}p_{l}{\rm
Tr}\left[\mathcal{N}_{l}\left(T_{k,1}\tilde{\rho}_{k-1}S^{\dagger}_{k,1}-T_{k,2}\tilde{\rho}_{k-1}S^{\dagger}_{k,2}\right)\mathcal{M}_{l}\left(S_{q,1}\tilde{\rho}_{q-1}T^{\dagger}_{q,1}-S_{q,2}\tilde{\rho}_{q-1}T^{\dagger}_{q,2}\right)\right]+\\\
&(1-p)^{2N}{\rm
Tr}\left[\left(\mathcal{R}_{N}\cdots\mathcal{R}_{k+1}\left(T_{k,1}{\rho}_{k-1}S^{\dagger}_{k,1}-T_{k,2}{\rho}_{k-1}S^{\dagger}_{k,2}\right)\right)\left(\mathcal{R}_{N}\cdots\mathcal{R}_{q+1}\left(S_{q,1}{\rho}_{q-1}T^{\dagger}_{q,1}-S_{q,2}{\rho}_{q-1}T^{\dagger}_{q,2}\right)\right)\right]\Big{)}\\\
&=\frac{1}{4}\left(\sum_{l}p_{l}{\rm
Tr}\left[\mathcal{N}_{l}\left(T_{k,1}\tilde{\rho}_{k-1}S^{\dagger}_{k,1}-T_{k,2}\tilde{\rho}_{k-1}S^{\dagger}_{k,2}\right)\mathcal{M}_{l}\left(S_{q,1}\tilde{\rho}_{q-1}T^{\dagger}_{q,1}-S_{q,2}\tilde{\rho}_{q-1}T^{\dagger}_{q,2}\right)\right]\right)+(1-p)^{2N}M_{k,q}\end{split}$
(A22)

Also, for two-qubits gates we get

$\displaystyle\begin{split}S_{k,1}\tilde{\rho}_{k-1}T^{\dagger}_{k,1}-S_{k,2}\tilde{\rho}_{k-1}T^{\dagger}_{k,2}&=(|0\rangle\langle
0|\otimes I+|1\rangle\langle 1|\otimes
R_{k})\tilde{\rho}_{k-1}(|1\rangle\langle 1|\otimes
R_{k}^{\dagger}P_{k})-(|1\rangle\langle 1|\otimes
P_{k}R_{k})\tilde{\rho}_{k-1}(|0\rangle\langle 0|\otimes I+|1\rangle\langle
1|\otimes R_{k}^{\dagger})\\\ &=\tilde{\rho}_{k}(|1\rangle\langle 1|\otimes
P_{k})-(|1\rangle\langle 1|\otimes P_{k})\tilde{\rho}_{k}\end{split}$ (A23)

Similarly, we can obtain the expression for the gate at $q$-th site as

$\displaystyle
S_{q,1}\tilde{\rho}_{q-1}T^{\dagger}_{q,1}-S_{q,2}\tilde{\rho}_{q-1}T^{\dagger}_{q,2}=\tilde{\rho}_{q}(|1\rangle\langle
1|\otimes P_{q})-(|1\rangle\langle 1|\otimes P_{q})\tilde{\rho}_{q}$ (A24)

Therefore, we can write

$\displaystyle\begin{split}M_{k,q}^{\rm
err}=\frac{1}{4}\left(\sum_{l}p_{l}{\rm
Tr}\big{[}\mathcal{N}_{l}\left(\left(|1\rangle\langle 1|\otimes
P_{k}\right)\tilde{\rho}_{k}-\tilde{\rho}_{k}\left(|1\rangle\langle 1|\otimes
P_{k}\right)\right)\mathcal{M}_{l}\left(\tilde{\rho}_{q}\left(|1\rangle\langle
1|\otimes P_{q}\right)-\left(|1\rangle\langle 1|\otimes
P_{q}\right)\tilde{\rho}_{q}\right)\big{]}\right)+(1-p)^{2N}M_{k,q}\end{split}$
(A25)

Let us now simplify the above expression by finding an upper bound to the
first terms in RHS

$\displaystyle\begin{split}|{\rm
Tr}\big{[}\mathcal{N}_{l}\big{(}\left(|1\rangle\langle 1|\otimes
P_{k}\right)&\tilde{\rho}_{k}-\tilde{\rho}_{k}\left(|1\rangle\langle 1|\otimes
P_{k}\right)\big{)}\mathcal{M}_{l}\big{(}\tilde{\rho}_{q}\left(|1\rangle\langle
1|\otimes P_{q}\right)-\left(|1\rangle\langle 1|\otimes
P_{q}\right)\tilde{\rho}_{q}\big{)}\big{]}|\\\
&\leq\|\mathcal{N}_{l}\big{(}\left(|1\rangle\langle 1|\otimes
P_{k}\right)\tilde{\rho}_{k}-\tilde{\rho}_{k}\left(|1\rangle\langle 1|\otimes
P_{k}\right)\big{)}\|_{2}\|\mathcal{M}_{l}\big{(}\tilde{\rho}_{q}\left(|1\rangle\langle
1|\otimes P_{q}\right)-\left(|1\rangle\langle 1|\otimes
P_{q}\right)\tilde{\rho}_{q}\big{)}\|_{2}\end{split}$ (A26)

Moreover, we have

$\displaystyle\begin{split}\|\mathcal{N}_{l}\big{(}\left(|1\rangle\langle
1|\otimes P_{k}\right)\tilde{\rho}_{k}-\tilde{\rho}_{k}\left(|1\rangle\langle
1|\otimes P_{k}\right)\big{)}\|_{2}&\leq\|\left(|1\rangle\langle 1|\otimes
P_{k}\right)\tilde{\rho}_{k}-\tilde{\rho}_{k}\left(|1\rangle\langle 1|\otimes
P_{k}\right)\|_{2}\\\ &=\sqrt{{\rm Tr}\left[\left(\left(|1\rangle\langle
1|\otimes P_{k}\right)\tilde{\rho}_{k}-\tilde{\rho}_{k}\left(|1\rangle\langle
1|\otimes P_{k}\right)\right)\left(\tilde{\rho}_{k}\left(|1\rangle\langle
1|\otimes P_{k}\right)-\left(|1\rangle\langle 1|\otimes
P_{k}\right)\tilde{\rho}_{k}\right)\right]}\\\ &=\sqrt{2}\sqrt{{\rm
Tr}[\tilde{\rho}^{2}(|1\rangle\langle 1|\otimes I)]-{\rm
Tr}[\tilde{\rho}(|1\rangle\langle 1|\otimes
P_{k})\tilde{\rho}(|1\rangle\langle 1|\otimes P_{k})]}\\\
&\leq\sqrt{2}\end{split}$ (A27)

Therefore, for two-qubit gates (at $k$-th and $q$-th sites) we have

$\displaystyle|{\rm Tr}\big{[}\mathcal{N}_{l}\big{(}\left(|1\rangle\langle
1|\otimes P_{k}\right)$
$\displaystyle\tilde{\rho}_{k}-\tilde{\rho}_{k}\left(|1\rangle\langle
1|\otimes
P_{k}\right)\big{)}\mathcal{M}_{l}\big{(}\tilde{\rho}_{q}\left(|1\rangle\langle
1|\otimes P_{q}\right)-\left(|1\rangle\langle 1|\otimes
P_{q}\right)\tilde{\rho}_{q}\big{)}\big{]}|\leq 2$ (A28)

Thus, we can bound $M^{\rm err}_{k,q}$ as

$\displaystyle{-(1-(1-p)^{2N}))\frac{1}{2}+(1-p)^{2N}M_{k,q}\leq}M^{\rm
err}_{k,q}\leq(1-(1-p)^{2N}))\frac{1}{2}+(1-p)^{2N}M_{k,q}$ (A29)

which implies

$\displaystyle|M^{\rm err}_{k,q}-(1-p)^{2N}M_{k,q}|\leq\frac{1-(1-p)^{2N}}{2}$
(A30)

Therefore, we can say that for the matrix $M^{\rm err}-(1-p)^{2N}M$, we have

$\displaystyle\|M^{\rm err}-(1-p)^{2N}M\|_{2}\leq(1-(1-p)^{2N})\frac{N}{2}$
(A31)

Using this, we can get an upper bound on $\|\Delta M\|$ as given below

$\displaystyle\|M^{\rm err}-M\|$ $\displaystyle=\|M^{\rm
err}-(1-p)^{2N}M+(1-p)^{2N}M-M\|$ (A32) $\displaystyle\leq\|M^{\rm
err}-(1-p)^{2N}M\|+(1-(1-p)^{2N})\|M\|$ (A33)
$\displaystyle\leq(1-(1-p)^{2N})\left(\frac{N}{2}+\|M\|\right)$ (A34)

Hence,

$\displaystyle\frac{\|\Delta
M\|}{\|M\|}\leq\left(1-(1-p)^{2N}\right)\left(\frac{N}{2\|M\|}+1\right)\,.$
(A35)

Similarly, we can find the upper bound on $Y^{\rm err}_{k,q}$ to be

$\displaystyle Y^{\rm
err}_{k,q}\leq(1-(1-p)^{N}))\frac{1}{\sqrt{2}}+(1-p)^{N}Y_{k,q}$ (A36)

and from this we get

$\displaystyle|Y^{\rm
err}_{k,q}-(1-p)^{N}Y_{k,q}|\leq(1-(1-p)^{N})\frac{1}{\sqrt{2}}$ (A37)

Following steps similar to the upper bound in Eq. (A35), we get the upper
bound on $\|\Delta Y\|$ to be

$\displaystyle\frac{\|\Delta
Y\|}{\|Y\|}\leq(1-(1-p)^{N})\left(1+\frac{\sqrt{N}}{\sqrt{2\|Y\|}}\right)$
(A38)

Finally, to derive upper bounds on the error in $\dot{\theta}$ we utilize the
formalism of studying perturbations to matrix linear equations to provide
bounds on the solution of our matrix equation $M\dot{\theta}=Y$ where both $M$
and $Y$ are perturbed. The ($p$-)condition number (denoted ${\rm cond}_{p}$)
of a matrix $Z$ is a defined as

$\displaystyle{\rm cond}_{p}(Z):=\|Z\|_{p}\|Z^{-1}\|_{p}$ (A39)

where $\|\cdot\|_{p}$ represents the matrix $p$-norm. If $Z$ is not
invertible, then ${\rm cond}_{p}(Z)=\infty$. In our work, we use $\|\cdot\|$
to mean $\|\cdot\|_{2}$ defined as $\|Z\|_{2}:=\sqrt{{\rm Tr}(Z^{*}Z)}$ (also
known as the Frobenius norm), and ${\rm cond}(Z)$ to mean ${\rm cond}_{2}(Z)$,
unless otherwise stated. If we consider a perturbation in both $M$ and $Y$ of
the form

$\displaystyle(M+\Delta M)(\Delta\dot{\theta}+\dot{\theta})=Y+\Delta Y$ (A40)

then the solution to the perturbed system obeys the following inequality [54]

$\displaystyle\frac{\|\Delta\dot{\theta}\|}{\|\dot{\theta}\|}\leq\frac{{\rm
cond}(M)}{1-{\rm cond}(M)\frac{\|\Delta M\|}{\|M\|}}\left(\frac{\|\Delta
M\|}{\|M\|}+\frac{\|\Delta Y\|}{\|Y\|}\right)\,.$ (A41)

Using the above inequality, we calculate the upper bounds under general noise
model and show the impact of generalized unitary errors in any given circuit
amounts to perturbing $M$ and $Y$.

By plugging in Eqs. (A35) and (A38) into equation (A41), we can bound the
error, $\epsilon:=\|\Delta\dot{\theta}\|$ as

$\displaystyle\epsilon\leq\frac{\|\dot{\theta}\|{\rm cond}(M)}{1-{\rm
cond}(M)\left(\left(1-(1-p)^{2N}\right)\left(1+\frac{N}{2\|M\|}\right)\right)}\left(\left(1-(1-p)^{2N}\right)\left(1+\frac{N}{2\|M\|}\right)+(1-(1-p)^{N})\left(1+\frac{\sqrt{N}}{\sqrt{2\|Y\|}}\right)\right)$
(A42)

which is valid under the constraint that

$\displaystyle 1\leq{\rm cond}(M)\leq\frac{2\|M\|}{(1-(1-p)^{2N})(N+2\|M\|)},$
(A43)

for a given probability of noise where cond$(\cdot)$ represents the condition
number of the matrix and $N$ is the total number of parameters in the circuit.
Equivalently, the above constraint can be expressed as a constraint on
probability provided the condition number of M is known as follows

$\displaystyle 0\leq p\leq 1-\left(1-\frac{2\|M\|}{{\rm
cond}(M)(N+2\|M\|)}\right)^{\frac{1}{2N}},$ (A44)

The above constraints arise due to the fact that the right hand side of Eq.
(A42) needs to be positive. The constraint on condition number can be
understood as follows: for a given probability of noise, the upper bound on
the error can be found for condition numbers between one and the bound in Eq.
(A43). Similarly, we can interpret the constraint on the allowed probability
range. Further analyzing the constraint on the condition number implies that
as the probability of noise decreases, the upper bound can be found for
matrices with larger condition numbers, and therefore, as the noise
probability tends to zero, the upper bound can be found for almost all the
cases (as the bound on the condition number tends to infinity). Also note that
the blow-out in the error as probability of noise is increasing is coming from
the fundamental result of perturbation theory and not from noise in QITE
algorithm itself. Therefore, it is crucial to find tighter bounds that are
valid for larger probabilities of noise.

To find tighter bounds, we need to impose certain physical constraints or
assumptions. The above general result is derived without any constraint. Now,
if we assume that the contribution of $\Delta M\Delta\dot{\theta}$ to $\Delta
Y$ in the perturbed equation $(M+\Delta
M)(\dot{\theta}+\Delta\dot{\theta})=Y+\Delta Y$ can be considered negligible
when $\|\Delta M\|\leq N^{2}\delta$ and $\|\Delta\dot{\theta}\|\leq N\delta$
for some fixed very small $\delta<<1$, then we can find bounds tighter than
the one presented above.

So, for the above constraint to work, we need the following to hold:

$\displaystyle\|\Delta
M\|\leq\|M\|(1-(1-p)^{2N})\left(1+\frac{N}{2\|M\|}\right)\leq N^{2}\delta$
(A45)

where the first inequality is coming form Eq. (A35) and the second inequality
is imposed so that the contribution of $\Delta M\Delta\dot{\theta}$ can be
considered negligible. By this assumption we then have

$\displaystyle M\Delta\dot{\theta}+\Delta M\dot{\theta}=\Delta Y$ (A46)

from which we get

$\displaystyle\|\Delta\dot{\theta}\|\leq\|M^{-1}\|(\|\Delta
Y\|+\|M\|\|\dot{\theta}\|)\,.$ (A47)

From the above equation and by using $\|Y\|\leq\|M\|\|\dot{\theta}\|$, we get

$\displaystyle\frac{\|\Delta\dot{\theta}\|}{\|\dot{\theta}\|}$
$\displaystyle\leq{\rm cond}(M)\left(\frac{\|\Delta M\|}{\|M\|}+\frac{\|\Delta
Y\|}{\|Y\|}\right)$ (A48) $\displaystyle\leq{\rm
cond}(M)\left(\left(1-(1-p)^{2N}\right)\left(1+\frac{N}{2\|M\|}\right)+(1-(1-p)^{N})\left(1+\frac{\sqrt{N}}{\sqrt{2\|Y\|}}\right)\right)$
(A49)

which is valid when the probability of noise is in the following range:

$\displaystyle 0\leq p\leq
1-\left(1-\frac{2N^{2}\delta}{N+2\|M\|}\right)^{\frac{1}{2N}}.$ (A50)

The above constraint on M is straightforward to get from Eq. (A45). ∎

##  Appendix B Proof of Theorem 2 (Error under depolarizing noise)

Theorem 2 statement. For partially depolarizing channels, the upper bound in
the error $\epsilon:=\|\Delta\dot{\theta}\|$ is given by

$\displaystyle\frac{\epsilon}{\|\dot{\theta}\|}=\frac{1-(1-p)^{N}}{(1-p)^{N}}.$
(B51)

###### Proof.

The action of the depolarizing channel is given by

$\displaystyle\mathcal{E}(X)=p{\rm Tr}[X]\frac{I}{d}+(1-p)X\,.$ (B52)

Assuming that the depolarizing noise affects the circuit after the application
of each gate, we get the following expression for ${\rm
Tr}\left[\left(\mathcal{R}_{k,i}\rho_{0}\right)^{\dagger}\left(\mathcal{R}_{q,j}\rho_{0}\right)\right]$

$\displaystyle\begin{split}{\rm
Tr}\left[\left(\mathcal{R}_{k,i}\rho_{0}\right)^{\dagger}\left(\mathcal{R}_{q,j}\rho_{0}\right)\right]=\;&\left(p(1-p)^{k-1}(a_{0}+a_{1}p+\ldots+a_{N-k}p^{N-k}){\rm
Tr}[P_{k}\rho_{k-1}]\right)\Bigg{(}\frac{1}{d}p(1-p)^{q-1}(b_{0}+b_{1}p+\ldots+b_{N-q}p^{N-q})\times\\\
&\qquad\qquad\quad{\rm
Tr}\left[P_{q}\rho_{q-1}\right]\Bigg{)}+\frac{(1-p)^{N}}{d}\\\
+&\left(1-p)^{N-k+1}(1-(1-p)^{k-1}\right)\Bigg{(}\frac{(1-p)^{N-q+1}(1-(1-p)^{q-1})}{d^{2}}{\rm
Tr}[\mathcal{R}_{N}\ldots\mathcal{R}_{k+1}(P_{k})\mathcal{R}_{N}\ldots\mathcal{R}_{q+1}(P_{q})]\\\
&\qquad\qquad\quad+\frac{(1-p)^{N}}{d}{\rm
Tr}\left[\mathcal{R}_{N}\circ\cdots\circ\mathcal{R}_{k+1}\left(P_{k}\right)\mathcal{R}_{N}\circ\cdots\circ\mathcal{R}_{q+1}\left(T_{q,j}\rho_{q-1}S^{\dagger}_{q,j}\right)\right]\Bigg{)}\\\
+&\left(1-p\right)^{N}\Bigg{(}\frac{p}{d}(1-p)^{q-1}(b_{0}+b_{1}p+\ldots+b_{N-q}p^{N-q}){\rm
Tr}[P_{q}\rho_{q-1}]{\rm Tr}[P_{k}\rho_{k-1}]\\\
&\qquad\qquad\quad+\frac{(1-p)^{N-q+1}(1-(1-p)^{q-1})}{d}{\rm
Tr}\left[R_{N}\circ\cdots\circ
R_{k+1}\left(S_{k,i}\rho_{k-1}T^{\dagger}_{k,i}\right)R_{N}\circ\cdots\circ
R_{q+1}\left(P_{q}\right)\right]\\\ &\qquad\qquad\quad+(1-p)^{N}{\rm
Tr}\left[R_{N}\circ\cdots\circ
R_{k+1}\left(S_{k,i}\rho_{k-1}T^{\dagger}_{k,i}\right)R_{N}\circ\cdots\circ
R_{q+1}\left(T_{q,j}\rho_{q-1}S^{\dagger}_{q,j}\right)\right]\Bigg{)}\;.\end{split}$
(B53)

where $P_{k}(P_{q})$ denotes the Pauli gate at the $k$-th($q$-th) rotation
gate. In the above expression, when we consider $p=0$, we get the ideal (non-
erroneous) scenario. We also noticed that the coefficients
$a_{0},a_{1},\ldots,a_{N-k}$ and $b_{0},b_{1},\ldots,b_{N-q}$ can be found
using Pascal’s triangle. To get these terms, let $n$ be the highest power in
the series. Then, the coefficients come from the $(n+2)$th row of the Pascal’s
triangle by ignoring the first $1$ of that row, and by alternating between the
$+$ and $-$ signs. Suppose, $n=2$, then we have the series as
$a_{0}+a_{1}p+a_{2}p^{2}$ and from Pascal’s triangle we get that $a_{0}=3$,
$a_{1}=-3$, and $a_{2}=1$.

Similarly to find $Y^{\rm err}$, when we include the depolarizing error after
every gate, we get

$\displaystyle\begin{split}{\rm Tr}[(\mathcal{R}_{k,i}\rho_{0})(\rho
H+H\rho)]&=\;\left(p(1-p)^{k-1}(a_{0}+a_{1}p+\ldots+a_{N-k}p^{N-k}){\rm
Tr}[P_{k}\rho_{k-1}]\right)\left(\frac{1}{d}{\rm Tr}[\rho H+H\rho]\right)\\\
&\quad+\;\left(1-p)^{N-k+1}(1-(1-p)^{k-1}\right)\frac{1}{d}{\rm
Tr}[\mathcal{R}_{N}\cdots\mathcal{R}_{k+1}(P_{k})(\rho H+H\rho)]\\\
&\quad+\;(1-p)^{N}{\rm
Tr}[\mathcal{R}_{N}\cdots\mathcal{R}_{k+1}(S_{k,i}\rho_{k-1}T^{\dagger}_{k,i})(\rho
H+H\rho)]\end{split}$ (B54)

Using the equations (B53) and (B54), we can find the elements of $M^{\rm err}$
and $Y^{\rm err}$. For instance, the diagonal elements of $M^{\rm err}$ (when
$k$-th and $q$-th gates are single-qubit gates) can be expressed as

$\displaystyle M^{\rm err}_{k,k}=(1-p)^{2N}\left(\frac{1}{2}-\frac{1}{2}{\rm
Tr}[P_{k}\rho_{k-1}]^{2}\right)\,.$ (B55)

The off-diagonal elements of $M^{\rm err}$, $M_{k,q}^{\rm err}$, where $k<q$,
when $k$-th and $q$-th gates are single-qubit gates, can be expressed as

$\displaystyle\begin{split}M^{err}_{k,q}=\frac{(1-p)^{2n}}{4}\Bigg{(}&{\rm
Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(S_{k,1}\rho_{k-1}T^{\dagger}_{k,1}\right)\left(T_{q,1}\rho_{q-1}S^{\dagger}_{q,1}\right)\right]\\\
-&{\rm Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(S_{k,1}\rho_{k-1}T^{\dagger}_{k,1}\right)\left(T_{q,2}\rho_{q-1}S^{\dagger}_{q,2}\right)\right]\\\
-&{\rm Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(S_{k,2}\rho_{k-1}T^{\dagger}_{k,2}\right)\left(T_{q,1}\rho_{q-1}S^{\dagger}_{q,1}\right)\right]\\\
+&{\rm Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(S_{k,2}\rho_{k-1}T^{\dagger}_{k,2}\right)\left(T_{q,2}\rho_{q-1}S^{\dagger}_{q,2}\right)\right]\Bigg{)}\end{split}$
(B56)

which can be written as

$\displaystyle\begin{split}M^{err}_{k,q}=\frac{(1-p)^{2N}}{4}\Bigg{(}&{\rm
Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(P_{k}\rho_{k}\right)\left(\rho_{q}P_{q}\right)\right]\\\ -&{\rm
Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(P_{k}\rho_{k}\right)\left(P_{q}\rho_{q}\right)\right]\\\ -&{\rm
Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(\rho_{k}P_{k}\right)\left(\rho_{q}P_{q}\right)\right]\\\ +&{\rm
Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}\left(\rho_{k}P_{k}\right)\left(P_{q}\rho_{q}\right)\right]\Bigg{)}\end{split}$
(B57)

Since we are considering that $\rho_{0}$ to be a pure state, the above
expression can be simplified as follows

$\displaystyle M^{err}_{k,q}=\frac{(1-p)^{2N}}{4}\Bigg{(}$ $\displaystyle{\rm
Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}(P_{k})\left(\rho_{q}P_{q}-2\rho_{q}P_{q}\rho_{q}+P_{q}\rho_{q}\right)\right]\Bigg{)}$
(B58)

So, the diagonal and off-diagonal elements can be expressed as

$\displaystyle\begin{split}M^{\rm err}_{k,k}&=\frac{(1-p)^{2N}}{2}\left(1-{\rm
Tr}[P_{k}\rho_{k-1}P_{k}\rho_{k-1}]\right)\\\
M^{err}_{k,q}&=\frac{(1-p)^{2N}}{4}\Bigg{(}{\rm Tr}\left[R_{q}\circ\cdots\circ
R_{k+1}(P_{k})\left(\rho_{q}P_{q}-2\rho_{q}P_{q}\rho_{q}+P_{q}\rho_{q}\right)\right]\Bigg{)}\end{split}$
(B59)

Since the matrix $M^{\rm err}$ is symmetric, we can find the elements
$M_{k,q}^{\rm err}$ when $k>q$. We found that even when the $k$-th and/or
$q$-th gates are two-qubit gates, the same equality holds, i.e.,

$\displaystyle M^{\rm err}_{k,q}=(1-p)^{2N}M_{k,q}\,$ (B60)

which follows owing to the particular form of the depolarizing channels and
due to the alternating signs in the definition of $M$ and $Y$.

Using similar calculations, the elements of $Y^{\rm err}$ can be expressed as

$\displaystyle Y_{k}^{\rm err}=(1-p)^{N}\left(\frac{1}{2}{\rm
Tr}[(\mathcal{R}_{N}\cdots\mathcal{R}_{k+1}(\rho_{k}P_{k}-P_{k}\rho_{k}))\\{\rho,H\\}]\right)\,.$
(B61)

So, we found that under partially depolarizing noise occurring after every
gate, the erroneous $M$ and $Y$ scale by a factor of $(1-p)^{2N}$ and
$(1-p)^{N}$ respectively, i.e.,

$\displaystyle M^{err}$ $\displaystyle=(1-p)^{2N}M$ (B62) $\displaystyle
Y^{err}$ $\displaystyle=(1-p)^{N}Y$ (B63)

With these equalities, given the perturbed equation

$\displaystyle M^{\rm err}\dot{\theta}^{\rm err}=Y^{\rm err}$ (B64)

we can reduce it to

$\displaystyle(1-p)^{N}M\dot{\theta}^{\rm err}=Y\,.$ (B65)

Therefore, we can write

$\displaystyle(1-p)^{N}\dot{\theta}^{\rm err}=\dot{\theta}\,.$ (B66)

From the above expression, we get the relative error to be

$\displaystyle\frac{\|\Delta\dot{\theta}\|}{\|\dot{\theta}\|}=\frac{1-(1-p)^{N}}{(1-p)^{N}}$
(B67)

∎

##  Appendix C Numerical Analyses

We performed the numerical simulation to plot the upper bound of the relative
error (in the rate of change of the parameters) as a function of the
probability of noise. We plotted this error vs probability graph for the
Hamiltonian of the Hydrogen molecule given in the main text. For simulating
this Hamiltonian using a general ansatz, we considered general probabilistic
noise and the depolarizing noise as the noise models.

Then, by choosing two ansatzes one with five parameters (see Fig. 4(a)) and
the other with nine parameters (see Fig. 4(b)), we again obtained the matrices
$M$ and $Y$. The nine parameter ansatz is the universal ansatz taken from Ref.
[55]. Using these ansatz and by generating random angles, we obtained $M$,
$Y$, $\dot{\theta}$. For the case when the noise is partially depolarizing,
the relative error is independent of the circuit and Hamiltonian, and only
depends on $N$ and probability of noise, $p$. The plot between relative error
as a function of probability of depolarizing noise for varying $N$ is given in
Fig. 5. To obtain the plots of relative error with respect to probability of
general hardware noise, we then fixed the condition number and varied the
number of parameters. The condition number that we used to plot Fig. 6(a) was
$66.7239$ and the norm of $M$ was 0.9977, both of which were obtained using
the ansatz in Fig. 4(a) by using the following angles:
$\theta_{1}=1.5249,\,\theta_{2}=2.5142,\,\theta_{3}=0.4457,\,\theta_{4}=1.3250,\,\theta_{5}=2.8769$.
The condition number that we used to plot Fig. 6(b) was $20.22$ which was
obtained using the ansatz in Fig. 4(b) by using the following angles:
$\theta_{1}=1.294,\,\theta_{2}=3.265,\,\theta_{3}=5.187,\,\theta_{4}=3.819,\,\theta_{5}=1.5708,\,\theta_{6}=4.172,\,\theta_{7}=0.618,\,\theta_{8}=3.441,\,\theta_{9}=2.502$.
The Frobenius norm of $M$ and $Y$ was calculated to be $1.689$ and $0.121$,
respectively. In both these figures, $\delta$ was fixed to be $0.04$; we also
fixed the condition number, $\|M\|$ and $\|Y\|$ as obtained above, and varied
$N$. Code/data will be made available upon request. One can see that for a
certain allowed probability of noise, relative error in $\dot{\theta}$ is more
for the case with more number of parameters.

(a)

(b)

Figure 4: (a) Ansatz with five parameters. (b) Ansatz with nine parameters
Figure 5: Plot of relative error (under partially depolarizing noise) with
respect to the probability of depolarizing noise for varying number of
parameters, $N$.

(a)

(b)

Figure 6: Plots of theorem 1 (relative error in the rate of change of
parameters with respect to the probability of general probabilistic noise) for
two different ansatzes. In each of the above plots, we fixed the condition
number of $M$ and varied the number of parameters, $N$. We computed the
condition number and the Frobenius norm of $M$ using the five-parameter and
the nine-parameter ansatz that we considered for simulating the Hamiltonian of
the Hydrogen molecule and by generating random angles. For (a) we fixed the
condition number to be $66.7239$ and norm of M to be 0.9977, and for (b), we
fixed the condition number to be $20.22$ and the norm of M to be 1.689. The
curves are plotted up to their valid probabilities as calculated from Eq.
(19).