Hybrid Quantum-Classical Algorithms and Quantum Error Mitigation

The variational quantum eigensolver (VQE) is a hybrid quantum-classical algorithm for computing the ground state and the ground state energy of a Hamiltonian H of interest. In the seminal work by Peruzzo et al.,⁹⁾ the VQE algorithm was theoretically proposed and experimentally demonstrated for finding the ground state energy of the HeH⁺ molecule using a two-qubit photonic quantum processor. We note that the conventional approach for finding eigenstates of a Hamiltonian with a universal quantum computer is by adiabatic state preparation and quantum phase estimation (QPE).⁵¹⁾ We refer to McArdle et al.⁴⁸⁾ and Cao et al.⁴⁹⁾ for the review.

The VQE algorithm relies on the Rayleigh–Ritz variational principle, i.e., for any parameterised quantum state \(|\varphi(\vec{\theta})\rangle\), we have \begin{equation} \min_{\vec{\theta}}\langle\varphi(\vec{\theta})|H |\varphi(\vec{\theta})\rangle \geq E_{G}, \end{equation} (2) where \(E_{G}\) is the ground state energy of the Hamiltonian H and the minimisation is over all parameters \(\vec{\theta}\).¹²⁾ As we will explain shortly, we can efficiently calculate \(\langle\varphi(\vec{\theta})|H |\varphi(\vec{\theta})\rangle\) with quantum processors. Therefore, by optimising parameters \(\vec{\theta}\) via a classical computer, considering \(E(\vec{\theta}) =\langle\varphi(\vec{\theta})|H |\varphi(\vec{\theta})\rangle\) as a to-be-minimised cost function, we can approximate the ground state energy and the ground state.

Now, we explain how to measure \(E(\vec{\theta}) =\langle\varphi(\vec{\theta})|H |\varphi(\vec{\theta})\rangle\) for a Hamiltonian H. As Pauli operators and products of them \(\{I, X, Y, Z \}^{\otimes N_{q}}\) form the complete basis for operators, any Hamiltonian can be expanded as \begin{equation} \begin{split} H &= \sum_{\alpha} f_{\alpha} P_{\alpha},\\ P_{\alpha} &\in \{I, X, Y, Z \}^{\otimes N_{q}}, \end{split} \end{equation} (3) where \(f_{\alpha}\) are real coefficients and \(\{X, Y, Z \}\) are Pauli operators. While an arbitrary \(N_{q}\)-qubit operator may have exponential terms in the expansion, Hamiltonians in reality are generally sparse so that the expansion only has the number of terms polynomial to \(N_{q}\). For example, the Hamiltonian of the 1D Ising model is \begin{equation} H = h\sum_{i} Z_{i} Z_{i+1} + \lambda \sum_{i} X_{i}, \end{equation} (4) where the number of terms is linear to the number of qubits. This also holds true for other systems. For example, to simulate the electronic structure of molecules, we consider the fermionic Hamiltonian \begin{equation} H_{f} = \sum_{ij} t_{ij} a_{i}^{\dagger} a_{j} +\sum_{ijkl} u_{ijkl} a^{\dagger}_{i} a_{k}^{\dagger} a_{l} a_{j}, \end{equation} (5) where \(a_{j}^{\dagger}\) (\(a_{j}\)) denotes the creation (annihilation) operator for a fermion in the j-th orbital, and \(t_{ij}\) and \(u_{ijkl}\) are the one and two electron interactions, which can be efficiently calculated by integrating the basis set wave functions.^52,53) The Hamiltonian consists of a polynomial number of terms consisting of products of fermionic operators. They can be further mapped to qubit operators via encoding methods such as Jordan–Wigner, parity, and Bravi–Kitaev encodings.^51,54,55) For example, the Jordan–Wigner transformation is defined as \begin{equation} \begin{split} a_{i}^{\dagger} &\rightarrow I^{\otimes i-1} \otimes \sigma^{-} \otimes Z^{\otimes N-i},\\ a_{i} &\rightarrow I^{\otimes i-1} \otimes \sigma^{+} \otimes Z^{\otimes N-i}, \end{split} \end{equation} (6) where \(\sigma^{\pm} = (X\mp i Y)/2\). We can therefore map the fermion Hamiltonian to the form of Eq. (3) with a polynomial number of Pauli operators.

By assuming a Pauli decomposition of the Hamiltonian \(H =\sum_{\alpha} f_{\alpha} P_{\alpha}\), the cost function \(E(\vec{\theta})\) becomes \begin{equation} E(\vec{\theta}) = \sum_{\alpha} f_{\alpha} \langle\varphi(\vec{\theta})|P_{\alpha}|\varphi(\vec{\theta})\rangle, \end{equation} (7) where each term \(\langle\varphi(\vec{\theta})|P_{\alpha}|\varphi(\vec{\theta})\rangle\) is evaluated by calculating the expectation value of each Pauli operator \(P_{\alpha}\). Note that the measurement of \(P_{\alpha}\) can be fully parallelised by using many quantum processors.

After we have obtained \(E(\vec{\theta})\), we update the parameters by using a classical computer to minimise the cost function. As an example, the gradient descent method updates the parameter settings as \begin{equation} \vec{\theta}^{(n+1)} = \vec{\theta}^{(n)}-a \nabla E(\vec{\theta}^{(n)}), \end{equation} (8) where \(\vec{\theta}^{(n)}\) and \(\vec{\theta}^{(n+1)}\) denote the parameters at the n-step and the \(n+1\)-step, respectively, \(a>0\) is a parameter determining the step size, and \(\nabla E(\vec{\theta}^{(n)})\) is the gradient of the cost function at \(\vec{\theta}^{(n)}\). The gradient descent method deterministically decreases the cost function to a local minimum. The concept of the VQE is illustrated in Fig. 2. In practice, it is important to opt for a fast and accurate optimisation method properly to reach the global minimum or feasible solution in a reasonable time. Also, the optimisation method should be robust to physical noise and shot noise in the quantum hardware. In addition to gradient based methods, another type of optimisation method is via direct search of the cost function. While the gradient may be more sensitive to physical noise — it typically vanishes exponentially in the number of qubits,⁵⁶⁾ direct search is believed to be more robust to physical noise, which may necessitate less repetitions.⁵⁷⁾ We refer to McArdle et al.⁴⁸⁾ for the review of classical optimisation algorithms.

»View larger version

Figure 2. Schematic of the variational quantum eigensolver. The expectation values of the Pauli operators \(P_{\alpha}\) are measured for the ansatz state, and the expectation value of the Hamiltonian is measured as a cost function on the classical computer. The cost function is sent to a classical optimiser to update the parameters.

Whether the VQE algorithm works also depends on the choice of the ansatz. To have an efficient quantum simulation algorithm, we need to use a suitable ansatz for the problem. If the ansatz state cannot express the solution, e.g., when the solution is a highly entangled state but the ansatz can only generate low entangled states, it cannot find the correct solution. In literature, several different types of ansätze are proposed for different purposes. For example, the unitary coupled cluster ansatz is known as a suitable physically inspired ansatz for electronic structure problems in chemistry.^9,58,59) However, the unitary coupled cluster ansatz generally necessitates a complicated form of quantum circuit with gates applied on multiple number of qubits, where each multi-qubit gate could be decomposed as a sequence of two-qubit gates, and the number of multi-qubit gates is quadratic to the number of qubits and the number of electrons. Since the unitary coupled cluster ansatz involves many general two-qubit gates, it may be hard to implement on noisy quantum devices with short coherence time and limited connectivity. This problem might be circumvented by leveraging so-called “hardware efficient ansätze”. This family of ansätze might be more experimentally feasible,^9,10,13) because they are constructed based on realisable demands on connectivity and gate operations that correspond to real quantum devices. However, a hardware efficient ansatz does not reflect the details of the simulated quantum system, and it has been shown that exponentially vanishing gradients (so-called barren plateaus) are liable to occur for randomly initialised parameters.⁶⁰⁾ There are several other methods proposed to circumvent the vanishing gradient problem.^61–64) We refer to McArdle et al.⁴⁸⁾ for a more detailed discussion for ansatz construction.

The disadvantage of the VQE method is that the correctness of the solution relies on the heuristic choice of the ansatz and the optimisation may be caught by a local instead of global minimum. Furthermore, the total number of measurements is \(O(\epsilon^{-2})\) for reaching a precision ϵ due to shot noise,^9,12) which is quadratically worse than the conventional QPE algorithm that uses universal quantum computers.

The VQE algorithm has been experimentally demonstrated by several groups.^{9,10,65–70)} To date, the hydrogen chain was simulated with a 12-qubit superconducting system.⁷¹⁾

Now we introduce the basic algorithms for variational quantum simulation (AQS), in particular, for simulating real¹⁴⁾ and imaginary¹⁸⁾ time evolution. The real time evolution of a quantum system can be described via the Schrödinger equation as \begin{equation} \frac{d |\psi(t)\rangle}{dt} = - i H |\psi(t)\rangle, \end{equation} (9) where H is the Hamiltonian and \(|\psi(t)\rangle\) is the time-dependent state. The conventional approach for simulating the evolution is to realize the time evolution \(e^{-iH t}\) as a unitary circuit and the state at time t is obtained by applying the unitary to the state.^6,72–74) The circuit depth generally increases polynomially with respect to the evolution time t. Instead, variational quantum simulation algorithms assume that the quantum state \(|\psi(t)\rangle\) is represented by an ansatz quantum circuit, \(|\varphi (\vec{\theta}(t))\rangle = U(\vec{\theta}(t))|\varphi_{\textit{ref}}\rangle\), and the time evolution of the Schrödinger equation of the state \(|\psi(t)\rangle\) is mapped to the evolution of the parameters \(\vec{\theta}(t)\).

Different variational principles can be used to have different evolution equations of the parameters. The three most conventional variational principles are — The Dirac and Frenkel variational principle,^75,76) McLachlan's variational principle,⁷⁷⁾ and the time-dependent variational principle.^78,79) The Dirac and Frenkel variational principle is not suitable for variational quantum simulation because the equation of the parameters may involve complex solutions, which contradicts with the requirement that parameters are real. Although the other two variational principles both give real solutions, it is shown that the time-dependent variational principle could be more unstable and it cannot be applied for evolution of density matrices and imaginary time evolution. On the contrary, McLachlan's principle generally produces stable solutions and it is also applicable to all the other scenarios beyond real time simulation. We refer a detailed study of the three variational principles to Yuan et al.¹⁹⁾

Now, we focus on McLachlan's variational principle⁷⁷⁾ and show how to derive the evolution of the parameters that effectively simulates the time evolution. McLachlan's principle requires one to minimise the distance between the ideal evolution and the evolution of the parametrised trial state as \begin{equation} \delta \|(\partial/\partial t + iH)|\varphi (\vec{\theta}(t))\rangle \| = 0, \end{equation} (10) where \(\|{|}\varphi\rangle \| =\langle\varphi|\varphi\rangle\) and δ is the variation over the derivative of the parameters \(\dot{\theta}_{j}\). With real parameters \(\vec{\theta}\), the solution gives the evolution of the parameters \begin{equation} \sum_{j} M_{k,j}\dot{\theta}_{j} = V_{k}, \end{equation} (11) with coefficients \begin{equation} \begin{split} M_{k,j} &= \mathop{\text{Re}}\nolimits \left(\frac{\partial \langle\varphi(\vec{\theta}(t))|}{\partial \theta_{k}}\frac{\partial |\varphi(\vec{\theta}(t))\rangle}{\partial \theta_{j}}\right),\\ V_{k} &= -\mathop{\text{Im}}\nolimits\left(\langle\varphi(\vec{\theta}(t))|H \frac{\partial |\varphi(\vec{\theta}(t))\rangle}{\partial \theta_{k}} \right), \end{split} \end{equation} (12) where \(\mathop{\text{Re}}\nolimits(\cdot)\) and \(\mathop{\text{Im}}\nolimits(\cdot)\) are the real and imaginary parts, respectively. We refer to Appendix A for the derivation.

Next, we describe the variational imaginary time simulation algorithm.¹⁸⁾ The normalised Wick-rotated Schrödinger equation⁸⁰⁾ is obtained by replacing t in Eq. (9) with = \(i\tau\), \begin{equation} \frac{d |\psi(\tau)\rangle}{d\tau} = - (H -\langle H\rangle)|\psi(\tau)\rangle, \end{equation} (13) where \(\langle H\rangle =\langle \psi(\tau)|H|\psi(\tau)\rangle\) is included for preserving the norm of the state \(|\psi(\tau)\rangle\). Notably, the imaginary time evolution can be leveraged for preparing a Gibbs state and for discovering the ground state of quantum systems.^18,19) Following the same procedure for real time evolution, we first make use of McLachlan's principle, \begin{equation} \delta \|(\partial/\partial \tau +H-\langle H\rangle)|\varphi(\vec{\theta}(\tau))\rangle \| = 0, \end{equation} (14) which then gives the evolution of the parameters: \begin{equation} \sum_{j} M_{k,j} \dot{\theta}_{j} = C_{k}, \end{equation} (15) with M defined in Eq. (12) and C defined by \begin{equation} C_{k} = -{\mathop{\text{Re}}\nolimits}\left(\langle\varphi(\vec{\theta}(\tau))|H \frac{\partial |\varphi(\vec{\theta}(\tau))\rangle}{\partial \theta_{k}}\right) = - \frac{1}{2} \frac{\partial E(\vec{\theta})}{\partial \theta_{k}}, \end{equation} (16) with \(E(\vec{\theta}) = (\langle\varphi(\vec{\theta}(\tau))|H |\varphi(\vec{\theta}(\tau))\rangle)\). Note that the C vector is related to the gradient of the energy, implying that variational imaginary time simulation can be regarded as a generalisation of the gradient descent method. It has been numerically found that the variational imaginary time simulation algorithm might be less sensitive to local minima in contrast to simple gradient descent methods.¹⁸⁾ In addition, when imaginary time evolution does reach a minimum that is not the ground state, it tends to be an excited eigenstate of the Hamiltonian, which can be thus exploited for finding general eigen-spectra.³²⁾ It has recently been observed that an equivalent formulation of the variational imaginary time algorithm can be obtained by exploiting the quantum natural gradient approach.^81–83) Notice that since the M matrix has to be measured, variational imaginary time simulation needs more measurements than the conventional gradient descent method. However, for an increasing system size and simulation time, the number of measurements required for M matrix is asymptotically negligible.⁸⁴⁾

Given the current parameters \(\vec{\theta}\), we now show how to efficiently measure the M, V, and C terms with quantum circuits. Suppose \(|\varphi(\vec{\theta}(t))\rangle = U_{N}(\theta_{N})\cdots U_{k}(\theta_{k})\cdots U_{1}(\theta_{1})|\varphi_{\textit{ref}}\rangle\) with the derivative of each unitary \(U_{k}(\theta_{k})\) expressed as \begin{equation} \frac{\partial U_{k}(\theta_{k})}{\partial \theta_{k}} = \sum_{i} g_{k,i} U_{k} \sigma_{k,i}, \end{equation} (17) where \(\sigma_{k,i}\) are unitary operators and \(g_{k,i}\) are complex coefficients. For instance, assuming \(U_{k} (\theta_{k}) = e^{-i\theta_{k} X}\), we have \(\partial U_{k} (\theta_{k})/\partial\theta_{k} = -i U_{k} X\) and hence \(g_{k,i} = -i\) and \(\sigma_{k,i} = X\). Thus, the derivative of the ansatz state is \begin{equation} \frac{\partial |\varphi(\vec{\theta}(t))\rangle}{\partial \theta_{k}} = \sum_{i} g_{k,i} U_{k,i} |\varphi_{\textit{ref}}\rangle, \end{equation} (18) where we defined \begin{equation} U_{k,i} = U_{N} U_{N-1} \cdots U_{k+1} U_{k} \sigma_{k,i} \cdots U_{2} U_{1}. \end{equation} (19) Now each \(M_{k,j}\) can be written as \begin{equation} M_{k,j} = \sum_{i,q}\mathop{\text{Re}}\nolimits (g^{*}_{k,i}g_{j,q} \langle\varphi_{\textit{ref}}| U^{\dagger}_{k,i} U_{j,q} |\varphi_{\textit{ref}}\rangle). \end{equation} (20) Supposing the Hamiltonian is decomposed as \(H =\sum_{\alpha} f_{\alpha} \sigma_{\alpha}\), with real coefficients \(f_{\alpha}\) and Pauli operators \(\sigma_{\alpha}\), we have \(V_{k}\) and \(C_{k}\) as \begin{equation} \begin{split} V_{k} &= - \sum_{i,\alpha} \mathop{\text{Re}}\nolimits (ig^{*}_{k,i} f_{\alpha} \langle\varphi_{\textit{ref}}| U^{\dagger}_{k,i}\sigma_{\alpha}U |\varphi_{\textit{ref}}\rangle),\\ C_{k} &= - \sum_{i,\alpha} \mathop{\text{Re}}\nolimits (g^{*}_{k,i} f_{\alpha} \langle\varphi_{\textit{ref}}| U^{\dagger}_{k,i}\sigma_{\alpha}U |\varphi_{\textit{ref}}\rangle). \end{split} \end{equation} (21) Note that each term constituting M, V, and C can be written as a sum of terms as \begin{equation} a \mathop{\text{Re}}\nolimits (e^{i\theta} \langle\varphi_{\textit{ref}}| \mathcal{V} |\varphi_{\textit{ref}}\rangle), \end{equation} (22) where \(a,\theta \in\mathbb{R}\) depend on the coefficients \(g_{k,i}\) and \(f_{j}\), and \(\mathcal{V}\) is either \(U^{\dagger}_{k,i}U_{j,q}\) or \(U^{\dagger}_{k,i}\sigma_{\alpha}U\). Then we can calculate M, C, and V with the quantum circuits shown in Fig. 3. Refer to Appendix B for a detailed explanation about construction of the quantum circuit. Notice that comprehensive analysis on the sampling cost for M, V, and C has been given by van Straaten and Koczor.⁸⁴⁾ We summarise the variational algorithm for simulating real (imaginary) time evolution.

»View larger version

Figure 3. Quantum circuits for computing (a) \(\mathop{\text{Re}}\nolimits(e^{i\theta}\langle\varphi_{\textit{ref}}|U_{k,i}^{\dagger} U_{j,q}|\varphi_{\textit{ref}}\rangle)\) and (b) \(\mathop{\text{Re}}\nolimits(e^{i\theta}\langle\varphi_{\textit{ref}}| U_{k,i}^{\dagger}\sigma_{j} U|\varphi_{\textit{ref}}\rangle)\).^14,18,19)

The schematic figure is also shown in Fig. 4. We can also simulate time-dependent Hamiltonian evolution by using the time-dependent Hamiltonian at each step. The accuracy of the simulation can be computed at each step by the distance between the evolution of the ansatz state and that of the ideal evolution.¹⁹⁾ In the case of the real time simulation, we have \begin{equation} \begin{split} \Delta &= \|(\partial/\partial t + iH)|\varphi (\vec{\theta}(t))\rangle \|^{2}\\ &= \sum_{k,j} M_{k,j} \dot{\theta}_{k} \dot{\theta}_{j} -2 \sum_{k} V_{k} \dot{\theta}_{k} +\langle H^{2}\rangle, \end{split} \end{equation} (23) which is a function of M, V, and \(\langle H^{2}\rangle\). By minimising the simulation error Δ, several strategies have been proposed to adaptively construct the optimal ansatz circuit for variational quantum simulation.^85,86) Similar arguments also hold for imaginary time evolution. Note that a variational real time simulation algorithm was demonstrated in an experiment using 4 superconducting qubits.⁸⁷⁾ In this experiment, adiabatic quantum computing was simulated, which was used for discovering eigenstates of an Ising Hamiltonian.

»View larger version

Figure 4. Schematic of the variational quantum simulation algorithm. The elements of M matrix and V (C) vectors are measured for the ansatz state. The results are sent to the classical computer to solve \(M\dot{\vec{\theta}} = V\) (C) to update the parameters, which will be fed to quantum processors.

The quantum approximate optimisation algorithm (QAOA)¹³⁾ was initially proposed for solving classical optimisation problems. The algorithm works by mapping the classical problem to a Hamiltonian \(H_{P}\) so that the ground state corresponds to the solution. Since we try to solve a classical optimisation problem, and the Hamiltonian \(H_{P}\) is diagonal in the computational basis as \(H_{P} =\sum_{\alpha} f_{\alpha} P_{\alpha}\) where \(P_{\alpha}\in \{I, Z \}^{\otimes N_{q}}\).

As an example, we consider the Boolean satisfiability problem, which aims to find solutions of Boolean variables so that all given clauses in a propositional formula are true. The j-th Boolean variable denoted as \(x_{j}\) takes a value either 1 or 0, each of which corresponds to true and false. The Boolean satisfiability problem consists of \(x_{i}\lor x_{j}\), \(x_{i}\land x_{j}\), and \(\bar{x}\) operations. \(x_{i}\lor x_{j}\) becomes 1 when either of \(x_{i}\) or \(x_{j}\) is 1, \(x_{i}\land x_{j}\) is a product of \(x_{i}\) and \(x_{j}\), and \(\bar{x}\) operations flip the value x. An example of Boolean satisfiability problem is \((x_{1}\lor x_{2})\land (\bar{x}_{1}\lor x_{2})\land (\bar{x}_{1}\lor \bar{x}_{2})\) and the solution that all the clauses are true (with value 1) is \(x_{1} = 0\) and \(x_{2} = 1\).

The general form of this problem (with clauses involving three or more booleans) is generally NP-hard and has a wide range of applications in computer science and cryptography.⁸⁸⁾ To map the Boolean satisfiability problem to a Hamiltonian, we first get the Hamiltonian with its ground state being the solution of each clause. For example, the Hamiltonian for the first clause \((x_{1}\lor x_{2})\) is \begin{equation} H_{1} = \frac{1}{4} (I-Z_{1})(I-Z_{2}). \end{equation} (24) After having the Hamiltonian \(H_{k}\) for each clause, the Hamiltonian for all clauses is expressed as \(H_{P} =\sum_{k} H_{k}\).

Since the solution corresponds to the ground state of \(H_{P}\), we can try to solve the problem by using a quantum algorithm for searching for the ground state. The variational approach for solving such problems is first studied by Farhi et al.¹³⁾ who introduced the ansatz \begin{equation} U(\vec{\theta}_{1}, \vec{\theta}_{2}) = \prod_{k = 1}^{D}e^{-i \theta_{1}^{(k)} H_{X}} e^{-i \theta_{2}^{(k)} H_{P}}, \end{equation} (25) with \(H_{X} =\sum_{j = 1}^{N_{q}} X_{j}\), D being the number of repetitions of the ansatz quantum circuit, \(\vec{\theta}_{1} = (\theta_{1}^{(1)},\theta_{1}^{(2)},\ldots)\), and \(\vec{\theta}_{2} = (\theta_{2}^{(1)},\theta_{2}^{(2)},\ldots)\). With initial states \(|{-},-,\ldots\rangle\) and \(|{-}\rangle ={\frac{1}{\sqrt{2}}}(|0\rangle-|1\rangle)\), we obtain the ansatz state \(|\varphi(\vec{\theta}_{1},\vec{\theta}_{2})\rangle = U(\vec{\theta}_{1},\vec{\theta}_{2})|{-},-,\ldots\rangle\). Directly optimising the parameters might be challenging for a fixed Hamiltonian, and Farhi et al.¹³⁾ suggested to consider gradually changing a Hamiltonian in each step of optimisation as \begin{equation} H(t) = \left(1-\frac{t}{T} \right) H_{X} +\frac{t}{T} H_{P}, \end{equation} (26) with \(H(0) = H_{X}\) and \(H(T) = H_{P}\). Since \(|{-},-,\ldots\rangle\) is the ground state of \(H(0) = H_{X}\) and solution is the ground state of \(H(T) = H_{P}\), the optimisation emulates the process of adiabatic state preparation with \(D\rightarrow\infty\). With sufficiently large D, and adaptive optimisations of the parameters for each optimisation step t, the algorithm may still be able to find the ground state solution.

A recent thorough study of the performance of the QAOA on MaxCut problems can be found in Zhou et al.⁸⁹⁾ Utilising nonadiabatic mechanisms, a heuristic strategy was proposed to learn the parameters exponentially faster than the conventional approach. Meanwhile, the QAOA was implemented with 40 trapped-ion qubits.⁹⁰⁾

Now we introduce the application of variational quantum algorithms in machine learning. In general, machine learning provides a universal approach to learn the pattern of the given data and to predict or reproduce new data. For example, in supervised learning, the given data are described by \(\{(\vec{x}_{1}, y_{1}), (\vec{x}_{2}, y_{2}),\ldots, (\vec{x}_{N_{D}}, y_{N_{D}})\}\) with \(N_{D}\) being the number of the data. Then we try to construct the model \(y = f(\vec{x})\) so that it predicts the output \(y_{\textit{new}} = f(\vec{x}_{\textit{new}})\) for any new data \(\vec{x}_{\textit{new}}\). Supposing the model \(y = f(\vec{x},\vec{\theta})\) has parameters \(\vec{\theta}\), the training process is to tune the parameters to minimise a cost function, such as \begin{equation} C_{\textit{ML}}(\vec{\theta}) = \sum_{k} |y_{k}- f(\vec{x}_{k}, \vec{\theta}) |^{2}. \end{equation} (27) When the cost function \(C_{\textit{ML}}(\vec{\theta})\) is minimised to a small value, it indicates that the given data are well-modelled. In practice, there are different ways to choose the model. For example, a linear regression model is, \begin{equation} f(\vec{x}, \vec{\theta}) = \vec{w}\cdot \vec{x}+b, \end{equation} (28) with real parameters \(\vec{w}\) and b, and \(\vec{\theta} = (\vec{w}, b)\). In deep leaning, i.e., a multi-layer neural network, the model is described as \begin{equation} \begin{split} f(\vec{x}, \vec{\theta}) &= \mathbf{g}_{N_{d}}\circ \mathbf{u}_{N_{d}}\cdots \circ \mathbf{u}_{2} \circ \mathbf{g}_{1} \circ \mathbf{u}_{1} (\vec{x}),\\ \mathbf{u}_{k} (\vec{x}) &= W_{k} \vec{x}+ \vec{b}_{k}, \end{split} \end{equation} (29) where \(N_{d}\) is the depth of the neural network, \(\mathbf{g}_{k}\) is a nonlinear activation function and \(W_{k}\) and \(\vec{b}_{k}\) are a parametrised matrix and vector, respectively. To circumvent overfitting, we can add the norm of the parameters to the cost function to restrict the degree of freedom of the parameters, which is called regularisation.

Whether machine learning works or not is highly dependent on the choice of the model. Since quantum states could efficiently represent multipartite correlations that admit no efficient classical representation, quantum machine learning protocols are proposed involving quantum neural networks that consist of variational quantum circuits. Consequently, we can leverage a quantum neural network to dramatically enhance the representability of the model.^20,21) In addition, as the ansatz circuit is a unitary operator, the norm of the quantum state is necessarily unity, and this constraint may lead to a natural regularisation of the parameters to avoid overfitting.²⁰⁾ The schematic figures for classical and quantum neural networks are shown in Fig. 5.

»View larger version

Figure 5. Comparison between (a) classical neural networks and (b) quantum neural networks used for supervised learning. The figure (b) is the quantum circuit proposed in quantum circuit leaning.²⁰⁾

Note that there are quantum machine learning algorithms^91–94) that are based on universal quantum computing without variational quantum circuits. While these algorithms are proven to have exponential speedups over classical algorithms under certain conditions, they generally require a deep quantum circuit and are not suitable for NISQ devices.

In this section, we illustrate five examples of quantum machine learning algorithms. The first two algorithms are introduced for leaning classical data for solving a classical problem; the latter three are introduced for learning quantum data.

3.2.2 Data-driven quantum circuit learning

Data-driven quantum circuit learning (DDQCL) implements a generative model by using a variational quantum circuit.²¹⁾ We briefly summarise the classical generative model. Suppose the data are described with \(\{\vec{x}_{1},\vec{x}_{2},\ldots,\vec{x}_{N_{D}} \}\), where \(N_{D}\) is the number of data and we assume \(\vec{x}\) has a binary representation, i.e., \(\vec{x} = \{x_{1},x_{2},\ldots,x_{N_{E}}\}\) with \(x_{j}\in \{-1, 1\}\) and \(N_{E}\) being the length of the binary string. We can assume the data is generated from an unknown probability distribution \(p_{D}(\vec{x})\), and the task of a generative model is to construct a model \(p_{M}(\vec{x}|\vec{\theta})\) to approximate the probability distribution \(p_{D}(\vec{x})\). The joint probability that the data is generated from \(p_{M}(\vec{x}|\vec{\theta})\) is \begin{equation} \mathcal{L}(\vec{\theta}) = \prod_{n = 1}^{N_{D}} p_{M}(\vec{x}_{n}| \vec{\theta}), \end{equation} (31) which is the so-called likelihood function. By maximising \(\mathcal{L}(\vec{\theta})\), we can obtain the optimised model probability distribution for the data. Alternatively, we can adopt the cost function \begin{equation} C_{\textit{DD}}(\vec{\theta}) = - \frac{1}{N_{D}} \sum_{n = 1}^{N_{D}} \log[\max(\varepsilon, p_{M}(\vec{x}_{n}| \vec{\theta}))], \end{equation} (32) where a small value \(\varepsilon>0\) is introduced for avoiding singularities of the cost function. For a classical generative model, an example model \(p_{M}(\vec{x}|\vec{\theta})\) could be generated via a Boltzmann machine on a classical computer.

For DDQCL, the model \(p_{M}(\vec{x}|\vec{\theta})\) is obtained from a quantum computer, for example, by the projection probability \begin{equation} p_{Q}(\vec{x}| \vec{\theta}) = |\langle \vec{x} | \varphi(\vec{\theta})\rangle|^{2}, \end{equation} (33) where \(|\varphi(\vec{\theta})\rangle\) is an ansatz state generated on the variational quantum circuit, \(|\vec{x}\rangle = |x_{1}, x_{2},\ldots, x_{N_{E}}\rangle\) is the computational basis, and the probability is defined according to the Born Rule. The quantum generative model is also referred to as the Born machine.⁹⁵⁾ Since quantum states can efficiently represent complex multipartite correlations and the model can be efficiently sampled by measuring the prepared state, DDQCL may be able to represent generative models that are classically challenging. In experiment, DDQCL has been applied to successfully learn Greenberger–Horne–Zeilinger states (GHZ) states and coherent thermal states by using a 4-qubit trapped ion device.²¹⁾ The schematic figure for sampling \(p_{Q}(\vec{x}|\vec{\theta})\) is shown in Fig. 6.

»View larger version

Figure 6. A quantum circuit for sampling \(p_{Q}(\vec{x}|\vec{\theta})\).

3.2.3 Quantum generative adversarial networks

Quantum generative adversarial networks (QuGANs)^22,96) are a quantum analogue of generative adversarial networks (GANs). Conventional GANs are composed of three parts — true data, generator, and discriminator, as shown in Fig. 7(a). The generator competes with the discriminator, where the former tries to produce fake data and the latter tries to determine whether the input data are true or fake. By optimising both the generator and the discriminator, the generator can learn the distribution of the true data until the discriminator cannot tell the difference between the true data and the fake data from the generator. GANs are generalised to quantum computing by replacing each part with a quantum system.²²⁾ By using a quantum circuit as the generator, we can represent the \(N_{q}\)-dimensional Hilbert space with \(\log N_{q}\) qubits and compute sparse and low-rank matrices with \(O(\text{poly}(\log N_{q}))\) steps. Suppose the true data are described with an ensemble of quantum states as a density matrix \(\rho_{\textit{true}}\), and the discriminator implements a quantum measurement. The schematic figure for quantum GANs is shown in Fig. 7(b).

»View larger version

Figure 7. (Color online) Schematic diagrams for (a) classical GANs and (b) quantum GANs. In classical GANs, the generator and the discriminator consist of classical neural networks, while quantum neural networks are used in quantum GANS.

Suppose the fake density matrix from the generator is produced from a parametrised quantum circuit as \(\rho_{G}(\vec{\theta}_{G})\) with parameters \(\vec{\theta}_{G}\), which aims to learn the true density matrix \(\rho_{\textit{true}}\). By using ancillary qubits, we can express density matrices with pure states and we refer to Dallaire-Demers and Killoran⁹⁶⁾ for detailed ansatz constructions. The discriminator implements a parametrised positive-operator valued measure \(\{P^{t}(\vec{\theta}_{D}), P^{f}(\vec{\theta}_{D})\}\), with parameters \(\vec{\theta}_{D}\) and \(P^{t}(\vec{\theta}_{D})+P^{f}(\vec{\theta}_{D}) = I\), to distinguish between \(\rho_{\textit{true}}\) and \(\rho_{G}(\vec{\theta}_{G})\). Assuming that \(\rho_{\textit{true}}\) and \(\rho_{G}(\vec{\theta}_{G})\) are sent to the discriminator randomly with equal probability, the probability that the discriminator fails is \begin{align} P_{\textit{fail}} &= \frac{1}{2}(\mathop{\text{Tr}}\nolimits[\rho_{G}(\vec{\theta}_{G}) P^{t}(\vec{\theta}_{D})]+\mathop{\text{Tr}}\nolimits[\rho_{\textit{true}} P^{f}(\vec{\theta}_{D})]),\\ &= \frac{1}{2}(C_{\textit{GA}}(\vec{\theta}_{D}, \vec{\theta}_{G})+1) \end{align} (34) with \(C_{\textit{GA}}(\vec{\theta}_{D},\vec{\theta}_{G}) =\mathop{\text{Tr}}\nolimits[(\rho_{G}(\vec{\theta}_{G}) -\rho_{\textit{true}}) P^{t}(\vec{\theta}_{D})]\) being proportional to the trace distance between \(\rho_{G}(\vec{\theta}_{G})\) and \(\rho_{\textit{true}}\) when optimised over \(P^{t}(\vec{\theta}_{D})\). With fixed parameters \(\vec{\theta}_{G}\) for the generator, we optimise the discriminator \(\vec{\theta}_{D}\) to minimise the failure probability \(P_{\textit{fail}}\) or \(C_{\textit{GA}}(\vec{\theta}_{D},\vec{\theta}_{G})\). With optimised discriminator, we fix \(\vec{\theta}_{D}\) and in turn optimise the generator \(\vec{\theta}_{G}\) to maximise the failure probability. By repeating this process, the parameters \(\vec{\theta}_{D}\) and \(\vec{\theta}_{G}\) arrives at an equilibrium with \(\rho_{G}(\vec{\theta}_{G}) =\rho_{\textit{true}}\), indicating a successful learning of the data.

3.2.4 Quantum autoencoder for quantum data compression

The classical autoencoder is used for compressing classical data as shown in Fig. 8(a). For an input set of training data \(\{\vec{x}_{1},\vec{x}_{2},\ldots,\vec{x}_{N_{D}} \}\), the autoencoder \(\mathcal{E}\) first encodes it to \(\{\vec{x}_{1}',\vec{x}_{2}',\ldots,\vec{x}_{N_{D}}' \}\) with each \(\vec{x}_{i}'\) being a smaller vector than \(\vec{x}_{i}\). A decoder \(\mathcal{D}\) is then applied to transform the data to \(\{\vec{x}_{1}'',\vec{x}_{2}'',\ldots,\vec{x}_{N_{D}}''\}\), with each \(\vec{x}_{i}''\) having the same size as \(\vec{x}_{i}\). The task of an autoencoder is to compress the input data into smaller size, with the requirement that the input data can be recovered from the compressed one with \(\sum_{k} \|\vec{x}_{k}''-\vec{x}_{k} \|^{2}\leq\varepsilon\) where ε is a desired accuracy.

»View larger version

Figure 8. Schematic figures for (a) classical autoencoder and (b) quantum autoencoder. The encoding operation \(\mathcal{E}\) compress the data, and the decoding operation \(\mathcal{D}\) decodes the data to the original dimension.

Quantum autoencoder implements a similar task for compressing quantum states.^23,97) Here, we consider the scheme proposed in Romero et al.²³⁾ Consider the case with input states of \(n+m\) qubits being compressed to states of n qubits. Suppose the input states are an ensemble \(\{p_{k}, |\phi_{k}\rangle_{\textit{AB}} \}\) with probabilities \(p_{k}\) and unknown states \(|\phi_{k}\rangle_{\textit{AB}}\). Here the subsystems A and B consists of n and m qubits, respectively. With a compression circuit \(U(\vec{\theta})\) and a post-selection measurement (for example in the computational basis) on the m-qubits of system B, as shown in Fig. 8(b), each input state \(|\phi_{k}\rangle_{\textit{AB}}\) of \(n+m\) qubits is mapped to a state \(\rho_{k}^{\textit{comp}}\) of n qubits. The inverse process then decodes the compressed state and map each \(\rho_{k}^{\textit{comp}}\) to an output state \(\rho_{k}^{\textit{out}}\). The average fidelity between the inputs and outputs is \begin{equation} C_{\textit{AE}}^{(1)}(\vec{\theta}) = \sum_{k} p_{k} F(|\phi_{k}\rangle, \rho^{\textit{out}}_{k}(\vec{\theta})), \end{equation} (35) which serves as a cost function and can be computed via the SWAP test or the destructive SWAP test circuit. Refer to Appendix C for details. When \(C_{\textit{AE}}^{(1)}(\vec{\theta})\) is maximised to a desired accuracy, it indicates a successful compression of the input state ensemble. Meanwhile, the successful compression and decoding, i.e., \(C_{\textit{AE}}^{(1)} = 1\), are achieved if and only if \begin{equation} U(\vec{\theta}) |\phi_{k}\rangle_{\textit{AB}} = |\phi_{k}^{\textit{comp}}\rangle_{A} \otimes |\bar{0}\rangle_{B},\ \forall k, \end{equation} (36) where \(|\phi_{k}^{\textit{comp}}\rangle_{A}\) corresponds to the compressed state, and \(|\bar{0}\rangle_{B}\) is some fixed reference state. Thus we can also adopt a simpler cost function \begin{equation} C_{\textit{AE}}^{(2)}(\vec{\theta}) = \sum_{k} p_{k} \mathop{\text{Tr}}\nolimits[U(\vec{\theta}) |\phi_{k}\rangle\langle\phi_{k}|_{\textit{AB}} U(\vec{\theta})^{\dagger} I_{A} \otimes |\bar{0}\rangle \langle\bar{0}|_{B}], \end{equation} (37) which can be computed as the probability projected to \(I_{A}\otimes |\bar{0}\rangle\langle\bar{0}|_{B}\).

Variational quantum algorithms can be used for solving matrix-vector multiplication and solving linear systems of equations. The task of matrix-vector multiplication is to obtain \(|v_{\mathcal{M}}\rangle =\mathcal{M} |v_{0}\rangle/\|\mathcal{M} |v_{0}\rangle\|\), where \(\mathcal{M}\) is a sparse matrix, \(|v_{0}\rangle\) is a given state vector, and \(\|{|}\psi\rangle \| =\sqrt{\langle\psi|\psi\rangle}\). Meanwhile, linear systems of equation is to solve a linear equation \(\mathcal{M} |v_{\mathcal{M}^{-1}}\rangle = |v_{0}\rangle\) to have \(|v_{\mathcal{M}^{-1}}\rangle =\mathcal{M}^{-1} |v_{0}\rangle\). There exist several variational quantum algorithms for implementing these tasks.^{30,31,104,105)} Herein, we illustrate the methods introduced in Bravo-Prieto et al.³⁰⁾ and Xu et al.³¹⁾

We first consider matrix-vector multiplication proposed by Xu et al.,³¹⁾ where the solution \(|v_{\mathcal{M}}\rangle\) corresponds to the ground state of the Hamiltonian, \begin{equation} H_{\mathcal{M}} = I - \frac{\mathcal{M}|v_{0}\rangle\langle v_{0}|\mathcal{M}^{\dagger}}{\| \mathcal{M} |v_{0}\rangle\|^{2}}. \end{equation} (42) When \(\mathcal{M}\) is a sparse matrix and the circuit for preparing \(|v_{0}\rangle\) is known, the expectation value \(E_{\mathcal{M}} =\langle\varphi(\vec{\theta})| H_{\mathcal{M}} |\varphi(\vec{\theta})\rangle\) can be efficiently evaluated by using the Hadamard test or the SWAP test circuit. Thus by minimising the expectation value \(E_{\mathcal{M}}\), we can find the ground state \(|v_{\mathcal{M}}\rangle\). Because the ground state has a unique ground state energy 0, we can further know whether the discovered state is the ground state or not, unlike the conventional VQE where the ground state energy is generally unknown. The matrix-vector multiplication algorithm can be applied for implementing Hamiltonian simulation. Suppose we want to simulate a time evolution operator \(U =\exp(- i H t)\) via Trotterisation \(U\approx\prod\exp(- i H\delta t)\). By setting \(M = 1- i H\delta t\), we can approximate each \(\exp(- i H\delta t)\) and hence the whole evolution as \(U = M^{N_{S}}+O(t^{2}/N_{S})\), where \(N_{S} = t/\delta t\). Therefore, by subsequently applying matrix-vector multiplication, we can simulate the time evolution of quantum systems.

Here we also consider the algorithms for solving linear equations independently proposed by Xu et al.³¹⁾ and Bravo-Prieto et al.¹⁰¹⁾ The solution is mapped to the ground state of the Hamiltonian \begin{equation} H_{\mathcal{M}^{-1}} = \mathcal{M}^{\dagger} (I- |v_{0}\rangle\langle v_{0}|) \mathcal{M}, \end{equation} (43) with the smallest eigenvalue 0 as well. This Hamiltonian was firstly proposed by Subaşı et al.,¹⁰⁶⁾ who applied adiabatic algorithms to find the ground state with a universal quantum computer. With the variational method, solving the ground state is similar to the one for matrix-vector multiplication. Furthermore, Xu et al.³¹⁾ used Hamiltonian morphing optimisation for avoiding local minima, and Bravo-Prieto et al.¹⁰¹⁾ used a local cost function for circumventing a barren plateau issue and run a simulation with up to 50 qubits. Refer to Appendix D for these optimisation methods.

Next we show how to find excited states and excited energy spectra of a Hamiltonian. Calculating excited energy spectra is important for studying many-body quantum physics problems. For example, it can be used to study chemical reaction dynamics, important for creating new drugs and new methodologies for mass production of beneficial materials.^53,107) In addition, evaluation of excited states enables us to calculate the photodissociation rates and absorption bands, which are essential for designing solar cells and investigating their dynamics.^108,109) There exist several VQAs for evaluating excited states and excited energy spectra.^{32,33,36,110–117)} These algorithms can be used as a subroutine for other applications, e.g., the Green's function,^35,118) non-adiabatic coupling, and Berry's phase¹¹⁹⁾ and simulating real time evolution¹²⁰⁾ etc. In this section, we review three VQAs — the overlap-based method, the subspace expansion method, and the contraction VQE method — for calculating excited state energy, and the application in calculating the Green's function.

Circuit recompilation aims to approximate a given quantum circuit with ones that are compatible with a practical experiment hardware. Concerning noisy gates or restricted set of realisable gates, the compiler runs to reduce the circuit noise or the implementation cost, as shown in Fig. 9. For example, an arbitrary two-qubit unitary is generally not directly supported on a practical quantum hardware and we need to compile the unitary into a sequence of realisable gates. While a naive decomposition of the unitary may induce too many unnecessary gates for hardware with specific topological structures, finding efficient and simple decomposition of quantum circuits is vital for near-term quantum computing.

»View larger version

Figure 9. Schematic figure for circuit recompilation algorithms. A variational quantum circuit with hardware constraints tries to approximate the target quantum circuit.

Denote the target unitary as \(U_{T}\) and the variational circuit as \(U(\vec{\theta})\), which consists of gates compatible with the hardware. Circuit recompilation is to tune the parameters \(\vec{\theta}\) so that \(U_{T}\approx U(\vec{\theta})\). We can mathematically define a metric of the distance between \(U_{T}\) and \(U(\vec{\theta})\) via the average gate infidelity (AGI)^127,128) \begin{equation} C_{I}(\vec{\theta}) = 1-\int d \varphi |\langle \varphi|U_{T}^{\dagger} U(\vec{\theta})|\varphi\rangle|^{2} \end{equation} (63) where pure states φ are randomly chosen according to the Haar measure. By construction, AGI indicates how different two unitary gates are on average for randomly sampled pure states, and it vanishes when the circuit is perfectly re-compiled. The AGI method has been applied for compiling high fidelity CNOT gate with cross-resonance gate suffering from crosstalk and single qubit operations. Furthermore, a high-fidelity four-qubit syndrome extraction circuit was recompiled to be achievable with simultaneous cross resonance drives under crosstalk.¹²⁸⁾

Another related cost function¹²⁷⁾ is \begin{equation} C_{\textit{HS}}(\vec{\theta}) = 1-\frac{1}{d}|{\mathop{\text{Tr}}\nolimits}(U_{T}^{\dagger} U(\vec{\theta}))|^{2}, \end{equation} (64) which can be efficiently calculated with the “Hilbert–Schmidt test” circuit¹²⁷⁾ by using two copies of states. In particular, we have \begin{equation} \frac{1}{d^{2}}|\mathop{\text{Tr}}\nolimits[U_{T}^{\dagger} U(\vec{\theta})]|^{2} = |\langle \Psi^{+}| U_{T} \otimes U^{*}(\vec{\theta})|\Psi^{+}\rangle|^{2}, \end{equation} (65) where \(|\Psi^{+}\rangle = d^{-1}\sum_{i} |i\rangle\otimes |i\rangle\) is the maximally entangled state with d being the dimension of the system. Then, \(|\mathop{\text{Tr}}\nolimits[U_{T}^{\dagger} U(\vec{\theta})]|^{2}/d^{2}\) can be evaluated by applying \(U_{T}\otimes U^{*}(\vec{\theta})\) to \(|\Psi^{+}\rangle\) and measure the probability to observe \(|\Psi^{+}\rangle\). Note that \(C_{\textit{HS}}(\vec{\theta})\) and AGI is related,^129,130) \begin{equation} C_{\textit{HS}}(\vec{\theta}) = \frac{d+1}{d} C_{I}(\vec{\theta}). \end{equation} (66) Thus \(C_{\textit{HS}}(\vec{\theta})\) can be used as an alternative cost function for circuit recompilation. Note that \(C_{I}(\vec{\theta})\) and \(C_{\textit{HS}}(\vec{\theta})\) generally exponentially vanish for an increasing system size. This problem is solved by using the weighted average of this global cost function and the corresponding local cost function (see Appendix D for an explanation of global and local costs). Then simple 9-qubit unitaries were successfully recompiled for noiseless simulator and for the Rigetti and IBM experimental processor.¹²⁷⁾

Meanwhile, circuit recompilation can be implemented for specific input states, which may reduce the optimisation complexity.^127,131) The goal is to find the recompiled circuit \(U(\vec{\theta})\) so that \(U_{T} |\psi_{\textit{in}}\rangle\approx U(\vec{\theta}) |\psi_{\textit{in}}\rangle\) holds for the state \(|\psi_{\textit{in}}\rangle\) and the target unitary \(U_{T}\). When \(|\psi_{\textit{in}}\rangle\) is a product state, we can easily define the Hamiltonian \(H_{\textit{in}}\) whose ground state is \(|\psi_{\textit{in}}\rangle\). For example, for \(|\psi_{\textit{in}}\rangle = |00\cdots 0\rangle\), we can set \(H_{\textit{in}} = -\sum_{j} Z_{j}\), where \(Z_{j}\) is the Pauli Z operator acting on the jth qubit. Then, by minimising the expectation value of \(H_{\textit{in}}\) for the trial state \(|\varphi(\vec{\theta})\rangle = U^{\dagger}(\vec{\theta}) U_{T} |\psi_{\textit{in}}\rangle\), we can obtain \(U(\vec{\theta})\) such that \(U^{\dagger}(\vec{\theta})U_{T} |\psi_{\textit{in}}\rangle\approx |\psi_{\textit{in}}\rangle\) and hence \(U_{T} |\psi_{\textit{in}}\rangle\approx U(\vec{\theta}) |\psi_{\textit{in}}\rangle\). We can leverage the conventional VQE method or variational imaginary time simulation for the optimisation and the fidelity of the recompiled state is lower-bounded as \begin{equation} F_{R} \geq 1- \frac{\delta}{E_{1}-E_{0}}, \end{equation} (67) where \(\delta =\langle \varphi(\vec{\theta})|H_{\textit{in}}|\varphi(\vec{\theta})\rangle-E_{0}\) is the deviation of the cost function from the ground state energy, and \(E_{0}\) and \(E_{1}\) are the ground state and first excited state energy. Note that \(E_{0} = -N_{q}\) and \(E_{1} = 2-N_{q}\) when the unitary acts on \(N_{q}\) qubits. This algorithm successfully recompiled a unitary operator involving 7 qubit system into another quantum circuit with different topology and dramatically reduced the number of two-qubit gates to 72 from 144, while the number of single-qubit gates increased to 77 from 42.¹³¹⁾

Quantum metrology aims to discover the optimal setup for probing a parameter with the minimal statistical shot noise.^132–134) The basic setup of quantum metrology is as follows — firstly, we prepare the initial probe state \(|\psi_{p}\rangle\) and evolve the state under the Hamiltonian \(H(\omega)\) with ω being the target parameter. After time t, the probe state can be described as \(|\psi_{p}(\omega, t)\rangle\), which is measured and analysed to extract the information of ω. Typically, ω can be a magnetic field, for example, with Hamiltonian \(H(\omega) =\omega \sum_{j}\sigma_{z}^{(j)}\) with \(\sigma_{z}^{(j)}\) being the Pauli Z operator, and the measurement is a Ramsey type measurement. When separable states are used as probes, the statistical error of the parameter behaves as \(\delta\omega \propto 1/\sqrt{N_{q}}\) with \(N_{q}\) being the number of qubits. This scaling is called the standard quantum limit (SQL).¹³⁴⁾ Notably the scaling can be improved by using entangled states, such as GHZ states, symmetric Dicke states, and squeezed states. In the absence of noise or for specific types of noise in the evolution under the Hamiltonian, the optimal strategy has been revealed. For example, with no environmental noise, the optimal probe state has been proved to be the GHZ state which achieves the Heisenberg limit \(\delta\omega \propto 1/N_{q}\).^133–135) However, in the presence of general types of noise, analytical arguments about the optimal strategy are usually very hard.

Variational-state quantum metrology is to use a quantum computer to find the optimal quantum state for quantum metrology with noisy hardware. Different proposals have been studied with either general variational quantum circuits¹³⁶⁾ or specific experimental setups,¹³⁷⁾ i.e., optical tweezer arrays of neutral atoms.^138,139) In general, suppose the initial probe state is created on a quantum device, as \(|\varphi_{p} (\vec{\theta})\rangle\). By evolving the state under the Hamiltonian \(H(\omega)\) followed by a proper measurement, we can define a cost function to reflect the metrology performance. In particular, we can either use the quantum Fisher information (QFI)¹³⁶⁾ or the spin squeezing parameter.^137,140) Here, we focus on the QFI, which characterises the minimum uncertainty of the estimated parameter as \begin{equation} \delta \omega \geq \frac{1}{\sqrt{N_{s} F_{Q}[\rho_{P}(\omega, t)]}}, \end{equation} (68) where \(F_{Q}[\rho_{P}(\omega, t)]\) is the QFI for the noisy output state \(\rho_{P}(\omega, t)\) after the evolution time t and \(N_{s}\) is the number of samples. Denoting \(T = N_{s}t\) to be the effective total time of all \(N_{s}\) samples, we can define a cost function as \begin{equation} C_{\textit{ME}}(\vec{\theta},t) = \frac{T}{t} F_{Q}[\rho_{P} (\omega, t, \vec{\theta})]. \end{equation} (69) For fixed T, we aim to optimise \(\vec{\theta}\) and t to maximise \(C_{\textit{ME}}(\vec{\theta},t)\). We show the schematic of this method in Fig. 10. Although QFI of mixed states cannot be directly evaluated in an efficient way, classical Fisher information (CFI) defined for a fixed measurement basis lower bounds QFI and it can be computed efficiently. We note that CFI is equivalent to QFI when the measurement basis is optimal, so QFI can in principle be obtained by optimising the measurement.

»View larger version

Figure 10. Schematic figure for variational quantum-state metrology. After the probe state is prepared by the variational quantum circuit, it evolves under the Hamiltonian \(H(\omega)\) with environmental noise, and is measured to evaluate the cost function.

With variational-state quantum metrology, a highly asymmetric state has been discovered for a 9-qubit system that outperforms previous results.¹³⁶⁾ Interestingly, even though the Hamiltonian and the noise model is symmetric under the permutation of qubits, there is a symmetry breaking for the optimal solution. Note that unlike conventional analytical approaches employed in quantum metrology, we do not have to know the noise model of the quantum device to obtain the optimised state. Recently, this algorithm was generalised to multi-parameter estimation.¹⁴¹⁾

Quantum error correction (QEC) makes use of a larger number of physical qubits to encode logical qubits to protect them against physical errors. In general, conventional formulation of QEC does not take into account the experimental implementation of the code. However, in the NISQ era, for example, the set of possible operations and the qubit topology are restricted for each physical hardware. Therefore, hardware-friendly implementation of QEC, tailored to actual experiments, is crucial for near-term quantum computers.^142–144) Here, we illustrate two examples^143,144) for realising QEC on NISQ computers.

3.7.2 Variational quantum error corrector (QVECTOR)

Variational quantum error corrector (QVECTOR) is for discovering a device-tailored quantum error correcting code.¹⁴⁴⁾ Different from the compiler for preparing a target logical state,¹⁴³⁾ QVECTOR aims to discover the optimal encode circuit that preserves the quantum state under noise. As shown in Fig. 11, the circuit \(U_{S}\) is used for preparing the to-be-encoded k-qubit state \(|\varphi\rangle\), \(V(\vec{\theta}_{1})\) and \(V^{\dagger}(\vec{\theta}_{1})\) are the noisy encoding and decoding circuits on \(n\geq k\) qubits, and \(W(\vec{\theta}_{2})\) is noisy gates for state recovery, which operates on \(n+r\) qubits. This quantum circuit corresponds to creating the \([n, k]\) quantum codes with r additional qubits. The parameters \(\vec{\theta}_{1}\) and \(\vec{\theta}_{2}\) are optimised to maximise the average code fidelity, \begin{equation} C_{\textit{QV}}(\vec{\theta}_{1}, \vec{\theta}_{2}) = \int_{\psi \in \mathcal{C}} \langle \psi| \mathcal{E}_{R} (\vec{\theta}_{1}, \vec{\theta}_{2})(|\psi\rangle \langle \psi|) |\psi\rangle\,d\psi, \end{equation} (72) where \(\mathcal{E}_{R} (\vec{\theta}_{1},\vec{\theta}_{2})\) is the encoding, recovery, and decoding operations \(V(\vec{\theta}_{1})\), \(W(\vec{\theta}_{2})\), and \(V^{\dagger}(\vec{\theta}_{1})\), and the integration is calculated following the Haar distribution with \(|\psi\rangle = U_{S} |0\rangle^{\otimes k}\). In practice, we can set \(U_{S}\) to be the unitary 2-design for efficiently calculating the average fidelity.¹⁴⁵⁾ With numerical simulation, QVECTOR learned the three qubit code¹⁴⁶⁾ under the phase damping noise, resulting in a six times longer \(T_{2}\) than conventional methods. In the presence of amplitude and phase damping noise, QVECTOR learned the quantum code outperforming the five qubit stabilizer code.

»View larger version

Figure 11. Schematic figure of the variational circuit for QVECTOR.

Dissipative-system variational quantum eigensolver (dVQE) is for obtaining the non-equilibrium steady state (NESS) in an open quantum system.¹⁴⁷⁾ In practice, quantum systems inevitably interact with their environments, and quantum states decohere owing to the noise. Therefore, the ability to simulate open quantum systems is indispensable for studying practical quantum phenomena. Notably, investigating NESS of open quantum systems is very important, e.g., in revealing the transport mechanism in nano-scale devices such as single-atom junctions.¹⁴⁸⁾

The time independent Markovian open quantum dynamics can be described by the Lindblad master equation, \begin{equation} \frac{d}{dt}\rho = - i [H, \rho]+ \mathcal{L}(\rho), \end{equation} (73) where ρ is the system state, H is the Hamiltonian, \(\mathcal{L}(\rho) =\sum_{k} (2L_{k}\rho L_{k}^{\dagger} -L_{k}^{\dagger} L_{k}\rho-\rho L_{k}^{\dagger} L_{k})\) and \(L_{k}\) is a Lindblad operator. Denoting \(\mathcal{L}'(\rho) = -i [H,\rho]+\mathcal{L}(\rho)\), the non-equilibrium steady state corresponds to the state with \(\mathcal{L}'(\rho) = 0\). To find the steady state, we first vectorise the state \(\rho =\sum_{n m}\rho_{n m} |n\rangle\langle m|\) to \begin{equation} |\rho\rangle\!\rangle = \frac{1}{\mathcal{N}}\sum_{n m} \rho_{nm} |n\rangle |m\rangle, \end{equation} (74) where \(\mathcal{N}\) is a normalisation factor. With vectorised \(\mathcal{L}'\), the steady state satisfies \(\mathcal{L}' |\rho_{\textit{ness}}\rangle\!\rangle = 0\) and hence we have \begin{equation} \langle\!\langle \rho_{\textit{ness}}| \mathcal{L}'^{\dagger} \mathcal{L}' |\rho_{\textit{ness}} \rangle\!\rangle = 0, \end{equation} (75) where \begin{equation} \begin{split} \mathcal{L}' &= (-i(H\otimes I-I\otimes H^{T})+ \mathcal{D}_{v}),\\ \mathcal{D}_{v} &= \sum_{k} \left(L_{k} \otimes L_{k}^{*} -\frac{1}{2} L_{k}^{\dagger} L_{k} \otimes I -\frac{1}{2} I \otimes L_{k}^{T} L_{k}^{*} \right). \end{split} \end{equation} (76) Here \(A^{*}\) denotes the complex conjugate of an operator A.

As \(|\rho(\vec{\theta})\rangle\!\rangle\) is a vector, we can express this vector representation of the quantum state through using a pure trial state on a variational quantum circuit by imposing proper constraints so that the corresponding density matrix \(\rho(\vec{\theta})\) will be physical, i.e., positive semi-definiteness and hermiteness have to be satisfied. The trace of the state needs to be unity but this condition will be considered when the expectation value of a given observable is measured. We can prepare the vectorised physical state on a variational quantum circuit as \begin{equation} \begin{split} |\rho(\vec{\theta}) \rangle\!\rangle &= U(\vec{\theta}_{1}) \otimes U^{*}(\vec{\theta}_{1}) |\rho_{d}(\vec{\theta}_{2}) \rangle\!\rangle,\\ |\rho_{d}(\vec{\theta}_{2}) \rangle\!\rangle &= \frac{1}{\mathcal{N}} \sum_{j} \alpha_{j}(\vec{\theta}_{2}) |j\rangle_{s} \otimes |j\rangle_{a}, \end{split} \end{equation} (77) where \(\vec{\theta}\equiv \{\vec{\theta}_{1},\vec{\theta}_{2} \}\), \(\alpha_{j}(\vec{\theta_{2}})>0\), and \(\sum_{j}\alpha_{j}(\vec{\theta_{2}}) = 1\), \(|j\rangle\) is in the computational basis, and the corresponding density matrix is \(\rho(\vec{\theta}) =\sum_{j}\alpha_{j} U(\vec{\theta}_{1}) |j\rangle_{s}\langle j|_{s} U^{\dagger} (\vec{\theta}_{1})\). Here, \(\mathcal{N} =\sqrt{\sum_{j}\alpha_{j}^{2}}\) ensures the normalisation. The subscripts s and a denote system and ancilla, respectively. We refer to Yoshioka et al.¹⁴⁷⁾ for the detailed ansatz construction. By preparing a trial state \(|\rho(\vec{\theta})\rangle\!\rangle\), and minimising the cost function \(C_{\textit{NE}}(\vec{\theta}) =\langle\!\langle\rho(\vec{\theta}|\mathcal{L}'^{\dagger}\mathcal{L}' |\rho (\vec{\theta})\rangle\!\rangle\), we can obtain an approximation of NESS. After the optimal parameters \(\vec{\theta}_{1}^{(op)}\) and \(\vec{\theta}_{2}^{(op)}\) are found, we can measure the expectation value of any given observable O for \(\rho(\vec{\theta})\) by randomly generating the state \(U(\vec{\theta}_{1}^{(op)})|j\rangle\) with probability \(\alpha_{j}(\vec{\theta}_{2}^{(op)})\), and averaging the measurement outcomes. This method was demonstrated for 8-qubit dissipative Ising model with a 16-qubit classical simulation.

There are other applications, such as VQAs for nonlinear problems,¹⁴⁹⁾ fidelity estimation,¹⁵⁰⁾ factoring,¹⁵¹⁾ singular value decomposition,^101,152) quantum foundations,¹⁵³⁾ circuit QED simulation,¹⁵⁴⁾ and Gibbs state preparation.^152,155,156)

We first show how to generalise the simulation algorithm for real and imaginary time evolution from pure states to mixed states.¹⁹⁾ The main idea is again to consider a parametrised representation of mixed states and map the dynamics to the evolution of the parameters. As we are considering mixed states, only McLachlan's variational principle applies, which leads to the evolution of parameters with information determined by the density matrix.¹⁹⁾ Although conventional quantum computers operate on pure states, we can also represent mixed states with their purifications by using ancilla qubits.

In this section, we review the variational quantum simulation algorithms for general processes, including the generalised time evolution, its application in solving linear algebra tasks, and simulating open system dynamics.³⁴⁾

The variational quantum simulation algorithm for imaginary time evolution can be applied for preparing a Gibbs state.¹⁹⁾ Starting at the maximally mixed state \(I_{d}/d\) with dimension d, imaginary time evolution of the state with Hamiltonian H and time τ prepares the state \(e^{-2H\tau}\), which corresponds to the Gibbs state at temperature \(T = 1/2\tau\). However, because the identity operator is invariant under any unitary transformation, we cannot directly input the maximally mixed state \(I_{d}/d\) to an unitary ansatz circuit to generate the Gibbs state. One resolution is to consider the purification of the Gibbs state. At time \(\tau = 0\), we input the maximally entangled state \(|\psi_{\text{max}}\rangle ={\frac{1}{\sqrt{d}}}\sum_{i = 1}^{d} |i\rangle_{s} |i\rangle_{a}\), where s and a denote the target system and the ancilla system for purification. It is easy to verify that \(\mathop{\text{Tr}}\nolimits_{a}[|\psi_{\text{max}}\rangle\langle\psi_{\text{max}}|] = I_{s}/d\) with \(\mathop{\text{Tr}}\nolimits_{a}[\cdot]\) denoting partial trace of ancilla. To realise imaginary time evolution on the system s, i.e., \(e^{-H_{s}\tau}|\psi_{\text{max}}\rangle/\|{e}^{-H_{s}\tau}|\psi_{\text{max}}\rangle\|\), we apply a joint unitary ansatz circuit on \(|\psi_{\text{max}}\rangle\) as shown in Fig. 12. Then we variatonally simulate imaginary time evolution by the evolution of the parameters of the joint unitary ansatz. The application of Gibbs state preparation in quantum Boltzmann machines was recently studied by Zoufal et al.¹⁶⁰⁾

»View larger version

Figure 12. The ansatz quantum circuit for obtaining the Gibbs state for \(N_{q} = 2\).

Now we discuss the application of the variational quantum simulation algorithms for calculating the Green's function and the spectral function.³⁵⁾ We refer to Sect. 3.4.4 for the review of the Green's function. The definition of the retarded the Green's function at zero temperature is \begin{equation} G_{\alpha \beta}^{(R)}(t) = - i \Theta(t) \langle G| c_{\alpha}(t) c_{\beta}^{\dagger}(0)+c_{\beta}^{\dagger}(0) c_{\alpha}(t) |G\rangle. \end{equation} (92) Here, \(\Theta(t)\) is a Heaviside step function, \(c_{\alpha (\beta)}^{(\dagger)}\) is an annihilation (a creation) fermion operator for the fermionic mode \(\alpha (\beta)\), \(c_{\alpha}(t) = e^{i H t} c_{\alpha} e^{-i H t}\), \(|G\rangle\) is the ground state, and we consider \(t>0\) for simplicity. We can transform fermionic operators to qubit operators^51,54,55) as \begin{equation} c_{\alpha} \rightarrow \sum_{k} f_{k}^{(\alpha)} P_{k},\quad c_{\alpha}^{\dagger} \rightarrow \sum_{k} f_{k}^{(\alpha)*} P_{k}, \end{equation} (93) where \(f_{k}\in \mathbb{C}\) and \(P_{k}\) is a product of Pauli operators. Then the retarded the Green's function can be described as \begin{equation} G_{\alpha \beta}^{(R)}(t) = \sum_{k k'} f_{k}^{(\alpha)} f_{k'}^{(\beta)*} \langle G| e^{i H t} P_{k} e^{-iHt} P_{k'} |G\rangle. \end{equation} (94) To calculate each term \(\langle G| e^{i H t} P_{k} e^{-iHt} P_{k'} |G\rangle\), we first approximate the ground state as \(|\tilde{G}\rangle\) by using either the conventional VQE or variational imaginary time simulation. Then, by using real time variational quantum simulation algorithms, we can approximate \begin{equation} \begin{split} e^{-iHt}|G\rangle &\approx U(\vec{\theta}_{1}(t))|\tilde{G}\rangle,\\ e^{-iHt}P_{k'}|G\rangle &\approx U(\vec{\theta}_{2}(t))P_{k'} |\tilde{G}\rangle. \end{split} \end{equation} (95) Note that, the approximated time evolution operators \(U(\vec{\theta}_{1}(t))\) and \(U(\vec{\theta}_{2}(t))\) are optimised for \(|\tilde{G}\rangle\) and \(P_{k'} |\tilde{G}\rangle\), so they are generally different operators. Finally we have \begin{equation} \langle G| e^{i H t} P_{k} e^{-iHt} P_{k'} |G\rangle \approx \langle\tilde{G}|U(\vec{\theta}_{1}(t))^{\dagger} P_{k} U(\vec{\theta}_{2}(t))P_{k'} |\tilde{G}\rangle, \end{equation} (96) which can be evaluated by the Hadamard test circuit. Note that a mechanism to dramatically reduce the number of control operation from a ancilla qubit was also introduced by Endo et al.³⁵⁾

Furthermore, it is worth noting that, by combining this method with the Gibbs state preparation algorithm in the last section,¹⁹⁾ we can also compute the Green's function for finite temperatures.

There are several other applications of variational quantum simulation algorithms. The first one is for financial problems.¹⁶¹⁾ The partial differential equation for pricing is equivalent to imaginary time evolution, so we can apply variational imaginary time simulation algorithm to simulate it. Secondly, variational imaginary time simulation can be used to simulate non-hermitian transcorrelated Hamiltonian for reducing the resource overhead and improving simulation accuracy for electronic structure calculations.¹⁶²⁾ Finally, it can be used for preparing a Boltzmann distribution with only a single copy of a quantum state.¹⁶³⁾

We first introduce the most simple yet very powerful error mitigation technique based on error extrapolation. The basic idea is to run the circuit with different error rates and use them to extrapolate for the error-free result, which is illustrated in Fig. 14.

»View larger version

Figure 14. (Color) Schematic of quantum error extrapolation. We boost the error rate and obtain other data points. Then we extrapolate the original result and those at boosted error rates to estimate the error-free result. Fitting curves should be chosen depending on the situation. Here, we show the case of two-point (green) and three-point (blue) Richardson extrapolation as examples.

Now, consider quantum error mitigation via the least square fitting.⁴⁰⁾ Although the assumption in literature is good characterisation of error rates and amplification of noise for extrapolation is not used, it may be employed to prepare results corresponding to several error rates. Consider well-characterised different error rates \(\varepsilon_{1},\ldots,\varepsilon_{n'}\), the expansion of Eq. (100) up to the nth order corresponds to a linear equation, \begin{equation} \mathbf{E}\vec{\mathbf{m}} \approx \vec{\mathbf{r}}, \end{equation} (116) where \begin{equation} \mathbf{E}= \begin{pmatrix} 1 &\varepsilon_{1} &\ldots &\varepsilon_{1}^{n}\\ 1 &\varepsilon_{2} &\ldots &\varepsilon_{2}^{n}\\ \vdots &\vdots &\ddots &\vdots\\ 1 &\varepsilon_{n'} &\ldots &\varepsilon_{n'}^{n} \end{pmatrix}, \end{equation} (117) \(\vec{\mathbf{m}} = (\langle M\rangle(0), M_{1},\ldots, M_{n})^{T}\), and \(\vec{\mathbf{r}} = (\langle M\rangle(\varepsilon_{1}),\langle M\rangle(\varepsilon_{2}),\ldots,\langle M\rangle(\varepsilon_{n}))^{T}\). Here the number of error rates \(n'\) can be different from the expansion order n. A solution of the linear equation can be obtained from minimising \begin{equation} \|\mathbf{E}\vec{\mathbf{m}} - \vec{\mathbf{r}}\|^{2}, \end{equation} (118) whose solution corresponds to the Richardson extrapolation method when \(n = n'\). Here we use the vector norm \(\|\vec{\mathbf{v}}\| =\sqrt{\vec{\mathbf{v}}^{T}\cdot\vec{\mathbf{v}}}\).

In practice, there may exist multiple noise parameters from different noisy effects, e.g., \(T_{1}\) and \(T_{2}\) noise. Suppose that there is another noise parameter ν, and we expand the expectation value \(\langle M\rangle(\varepsilon,\nu)\) as \begin{align} \langle M\rangle(\varepsilon, \nu) &= \langle M\rangle(0)+M_{10} \varepsilon +M_{01} \nu\\ &\quad + M_{11} \varepsilon \nu +M_{20} \varepsilon^{2}+M_{02} \nu^{2}+\cdots. \end{align} (119) Suppose we consider the truncation up to the second order, we get the same linear expansion of Eq. (116) with \begin{equation} \mathbf{E}= \begin{pmatrix} 1 &\varepsilon_{1} &\nu_{1} &\varepsilon_{1} \nu_{1} &\varepsilon_{1}^{2} &\nu_{1}^{2}\\ 1 &\varepsilon_{2} &\nu_{2} &\varepsilon_{2} \nu_{2} &\varepsilon_{2}^{2} &\nu_{2}^{2}\\ \vdots &\vdots &\ddots &\vdots &\vdots &\vdots\\ 1 &\varepsilon_{n'} &\nu_{n'} &\varepsilon_{n'} \nu_{n'} &\varepsilon_{n'}^{2} &\nu_{n'}^{2} \end{pmatrix}, \end{equation} (120) and \begin{equation} \vec{\mathbf{m}} = \begin{pmatrix} \langle M\rangle(0)\\ M_{10}\\ M_{01}\\ M_{11}\\ M_{20}\\ M_{02} \end{pmatrix},\quad \vec{\mathbf{r}} = \begin{pmatrix} \langle M\rangle(\varepsilon_{1}, \nu_{1})\\ \langle M\rangle(\varepsilon_{2},\nu_{2})\\ \vdots\\ \langle M\rangle(\varepsilon_{n'},\nu_{n'}) \end{pmatrix}. \end{equation} (121) An estimation of \(\langle M\rangle(0)\) can be similarly obtained by minimising Eq. (118). This method was demonstrated on Rigetti's 8-qubit device.⁴⁰⁾

The quasi-probability method is another error mitigation technique, which was first introduced by Temme et al.¹⁶⁾ for special channels and then generalised to practical Markovian noise by Endo et al.¹⁷⁾ The main idea is for any noise process, its effect can be cancelled out by probabilistically implementing its inverse process. Considering a single noisy gate as \(\mathcal{E}\circ\mathcal{U}\) with ideal gate \(\mathcal{U}\) and noise \(\mathcal{E}\), we may find an operation \(\mathcal{E}^{-1}\) that inverts the noise with \(\mathcal{E}^{-1}\circ\mathcal{E} =\mathcal{I}\). Note that we can mathematically invert the noise \(\mathcal{E}\) for most practical cases, the inverse process is in general unphysical and it cannot be directly realised by applying unitary operations on quantum states. Nevertheless, for any basis set of quantum processes \(\{\mathcal{B}_{i} \}\), we can decompose the inverse channel as \(\mathcal{E}^{-1} =\sum_{i} q_{i}\mathcal{B}_{i}\) and consequently we have \begin{equation} \mathcal{U} = \mathcal{E}^{-1}\circ \mathcal{E} \circ\mathcal{U} = \sum_{i} q_{i} \mathcal{K}_{i}, \end{equation} (122) with \(\mathcal{K}_{i} =\mathcal{B}_{i}\circ\mathcal{E}\circ\mathcal{U}\). Therefore even if we cannot directly implement the inverse process, we can effectively realise it by applying basis operations \(\mathcal{B}_{i}\) to the noisy gate and linearly combining the channels \(\mathcal{K}_{i}\) with real coefficients \(q_{i}\). In particular, supposing we input the gate with state ρ and measure the output with an observable M, the ideal measurement statistic \(\langle M\rangle_{\mathcal{U}(\rho)} =\mathop{\text{Tr}}\nolimits [\mathcal{U}(\rho)M]\) can be obtained as \begin{equation} \langle M\rangle_{\mathcal{U}(\rho)} = \sum_{i} q_{i}\langle M\rangle_{\mathcal{K}_{i}(\rho)}, \end{equation} (123) where each \(\langle M\rangle_{\mathcal{K}_{i}(\rho)} =\mathop{\text{Tr}}\nolimits [\mathcal{K}_{i}(\rho)M] =\mathop{\text{Tr}}\nolimits [\mathcal{B}_{i}\circ\mathcal{E}\circ\mathcal{U}(\rho)M]\) corresponds to the measurement with noisy gates. Instead of measuring all \(\langle M\rangle_{\mathcal{K}_{i}(\rho)}\), we can define \(p_{i} = |q_{i}|/C_{\mathcal{E}}\) with \(C_{\mathcal{E}} =\sum_{i}|q_{i}|\) so that \begin{equation} \langle M\rangle_{\mathcal{U}(\rho)} = C_{\mathcal{E}}\sum_{i} \mathop{\text{sgn}}\nolimits(q_{i})p_{i}\langle M\rangle_{\mathcal{K}_{i}(\rho)}. \end{equation} (124) Therefore, we can probabilistically append the basis operation \(\mathcal{B}_{i}\) to the noisy channel \(\mathcal{E}\circ\mathcal{U}\) with probability \(p_{i}\). Then we measure the output state and multiply the measurement outcome with the parity \(\mathop{\text{sgn}}\nolimits(q_{i})\). The ideal measurement \(\langle M\rangle_{\mathcal{U}(\rho)}\) corresponds to the averaged measurement results multiplied by \(C_{\mathcal{E}}\). Note that, because we multiply \(C_{\mathcal{E}}\) to the averaged result, the variance is amplified with \(\gamma_{Q} = C_{\mathcal{E}}^{2}\), which is regarded as the error mitigation cost.

As an example, we consider error mitigation of the depolarising noise \begin{equation} \mathcal{D}(\rho) = \left(1-\frac{3}{4}p\right)\rho+\frac{p}{4}(X\rho X+Y \rho Y + Z \rho Z), \end{equation} (125) whose inverse process is mathematically given by \begin{equation} \mathcal{D}^{-1}(\rho) = C_{\mathcal{D}}[p_{1} \rho -p_{2} (X\rho X+Y \rho Y + Z \rho Z)], \end{equation} (126) where \(C_{\mathcal{D}} = (p+2)/(2-2p)>1\), \(p_{1} = (4-p)/(2p+4)\), \(p_{2} = p/(2p+4)\), and \(p_{1}+3p_{2} = 1\). Thus, for any ideal unitary \(\mathcal{U}\), it can be written as \begin{align} \mathcal{U} &= \mathcal{D}^{-1} \circ\mathcal{DU}\\ &= C_{\mathcal{D}} [p_{1} \mathcal{I}\circ\mathcal{DU}-p_{2} (\mathcal{X}\circ\mathcal{DU}+\mathcal{Y}\circ\mathcal{DU}+\mathcal{Z}\circ\mathcal{DU})]. \end{align} (127) Here, we denote the I, X, Y, and Z gates as \(\mathcal{I}\), \(\mathcal{X}\), \(\mathcal{Y}\), and \(\mathcal{Z}\), respectively and the noisy gate \(\mathcal{D}\circ \mathcal{U}\) as \(\mathcal{DU}\). To recover the ideal gate \(\mathcal{U}\) from the noisy gate \(\mathcal{DU}\), we append the X, Y, Z gates with probability \(p_{2}\). When the state is measured, the parity corresponding to the generated operation and the constant \(C_{\mathcal{D}}\) is multiplied to the measurement result. By repeating this procedure many times, we can recover the measurement outcome with the ideal gate.

The quasi-probability method can be applied to a general quantum circuit, as illustrated in Fig. 15. Let us assume that the ideal process of the entire quantum circuit is described as \(\prod_{k = 1}^{N_{g}}\mathcal{U}_{k}\), where \(\mathcal{U}_{k}\) corresponds to each gate. Suppose \(\mathcal{U}_{k}\) is decomposed as \begin{equation} \mathcal{U}_{k} = C_{k} \sum_{i_{k}} p_{i_{k}} \mathop{\text{sgn}}\nolimits (q_{i_{k}}) \mathcal{K}_{i_{k}}, \end{equation} (128) and the ideal process can be represented as \begin{align} \prod_{k = 1}^{N_{g}} \mathcal{U}_{k} &= \prod_{k = 1}^{N_{g}}C_{k} \sum_{i_{1} i_{2}\cdots i_{N_{g}}} \prod_{k = 1}^{N_{g}}p_{i_{k}} \prod_{k = 1}^{N_{g}}\mathop{\text{sgn}}\nolimits (q_{i_{k}}) \prod_{k = 1}^{N_{g}} \mathcal{K}_{i_{k}}\\ &= C_{\textit{tot}}\sum_{\vec{i}} p_{\vec{i}} \cdot \mathop{\text{sgn}}\nolimits(q_{\vec{i}})\cdot \mathcal{K}_{\vec{i}}, \end{align} (129) where \(\vec{i} = (i_{1},i_{2},\ldots,i_{N_{g}})\), \(C_{\textit{tot}} =\prod_{k = 1}^{N_{g}}C_{k}\), \(p_{\vec{i}} =\prod_{k = 1}^{N_{g}}p_{i_{k}}\), \(\mathop{\text{sgn}}\nolimits(q_{\vec{i}}) =\prod_{k = 1}^{N_{g}}\mathop{\text{sgn}}\nolimits (q_{i_{k}})\), and \(\mathcal{K}_{\vec{i}} =\prod_{k = 1}^{N_{g}}\mathcal{K}_{i_{k}}\). Therefore, we can similarly realise each noisy process \(\mathcal{K}_{\vec{i}}\) with probability \(p_{\vec{i}}\) and multiply the outcome with \(\mathop{\text{sgn}}\nolimits(q_{\vec{i}})\) to recover the effect of the ideal circuit. We also multiply \(C_{\textit{tot}}\) to the averaged outcome. The error mitigation cost is \(\gamma_{\textit{tot}} = C_{\textit{tot}}^{2}\), which becomes \(\gamma_{\textit{tot}}\approx e^{2 b\varepsilon N_{g}}\) when we assume \(C_{k} = 1+ b\varepsilon\) with positive coefficient b and error probability ε. For stochastic noise, b is generally less than 2.¹⁷⁾ We can see that the cost increases exponentially to \(\varepsilon N_{g}\), i.e., the averaged number of errors in the whole circuit. Therefore, the quasi-probability method is efficient only if \(\varepsilon N_{g} = O(1)\).

»View larger version

Figure 15. (Color online) Schematic of the quasi-probability method. We apply the inverse channel to cancel the noise process \(\mathcal{E}\) with quasi-probability method. The inverse operation specified in the figure corresponds to the case of depolarising channel. For simplicity, we just show the case that two-qubit gate error is a tensor product of depolarising channels. We multiply the product of parities of generated operations \(\mathop{\text{sgn}}\nolimits(q_{\vec{i}})\) to the outcome. We compute the average of product of the parity and the outcome, which is multiplied with the cost \(C_{\textit{tot}}\) to approximate the ideal expectation value of the observable.

In contrast to the extrapolation method, we also need to exactly identify the noise model \(\mathcal{E}\) in quantum circuits to apply the quasi-probability method. However, state preparation and measurement (SPAM) errors in the conventional process tomography will result in errors of the estimated process and hence the effect of error mitigation. Such a problem was addressed in Endo et al.,¹⁷⁾ where gate set tomography (GST) was applied so that the quasi-probability decomposition obtained from GST is free from SPAM errors. The same paper also showed how to use experimentally feasible operations to suppress general Markovian errors. The application of quasi-probability in the presence of non-Markovian errors, such as temporally and spatially correlated errors, was further studied by Huo and Li.¹⁷⁴⁾ Detailed theoretical analysis of the cost using resource theory approach was given by Takagi.¹⁷⁵⁾ The quasi-probability method has been demonstrated using superconducting¹⁷⁶⁾ and trapped ion systems¹⁷⁷⁾ for single- and two-qubit gates.

The quantum subspace expansion (QSE) can not only be used for evaluating excited states, but also for quantum error mitigation in variational quantum eigensolver (VQE).³⁶⁾ Suppose that we already have the approximation of the ground state ρ of the Hamiltonian H via either the VQE or variational imaginary time evolution. The approximation ρ may not be the exact ground state of H due to gate errors or imperfect ansatz and optimisation. To mitigate such errors, we consider a small subset of operators \(S\subset \{I, X, Y, Z \}^{\otimes N_{q}}\) and a subspace \begin{equation} \left\{\rho_{\textit{QEM}} = \sum_{i,j} c_{i} c_{j}^{*} P_{i} \rho P_{j}\,|\, P_{i}\in S,\ c_{i}\in \mathbb{C},\ \mathop{\text{Tr}}\nolimits [\rho_{\textit{QEM}}] = 1 \right\}. \end{equation} (130) The QSE method is to minimise the energy of quantum states from the subspace \begin{equation} \min\nolimits_{\vec{c}} \mathop{\text{Tr}}\nolimits[\rho_{\textit{QEM}} H] \end{equation} (131) over all possible parameters \(\vec{c} = (c_{0}, c_{1},\ldots.)\). Such a minimisation can be reformulated as solving a generalised eigenvalue problem \begin{equation} \tilde{H} \vec{c} = E \tilde{S} \vec{c}, \end{equation} (132) where E is the eigenvalue, \(\tilde{H}_{ji} =\mathop{\text{Tr}}\nolimits[\rho P_{j} H P_{i}]\) and \(\tilde{S}_{ji} =\mathop{\text{Tr}}\nolimits[\rho P_{j} P_{i}]\). Note that both \(\tilde{H}\) and \(\tilde{S}\) can be efficiently measured on a quantum computer, and the problem can be efficiently solved on a classical computer given a small set of S. Once \(\vec{c}\) is determined, we can compute the expectation value of any operator M as \(\sum_{ij}c_{i} c_{j}^{*}\mathop{\text{Tr}}\nolimits[\rho P_{j} M P_{i}]\).

The QSE method is more effective for coherent errors such as over-rotation of quantum gates, and it may not mitigate all local stochastic errors. Specifically, the subspace expansion method is equivalent to optimising states from the subspace \(\{\rho_{\textit{QEM}} =\mathcal{P}\rho\mathcal{P}^{\dagger}/{\mathop{\text{Tr}}\nolimits}[\mathcal{P}\rho\mathcal{P}^{\dagger}]\}\) with \(\mathcal{P} =\sum_{i} c_{i}P_{i}\). Under a singular value decomposition of the operator \(\mathcal{P}\), its effect can be regarded as a combination of state rotation and projection, which may be able to correct coherent errors and stochastic errors, respectively. Small under- and over-rotation of quantum gates could be corrected via the unitary part of \(\mathcal{P}\), which makes QSE effective for coherent errors. For stochastic errors, the projection part of the operator \(\mathcal{P}\) is supposed to project out the errors. For example, considering ρ to be a mixture of the ideal pure state \(|\psi\rangle\langle\psi|\) and the noise σ as \(\rho = (1-\varepsilon)|\psi\rangle\langle\psi|+\varepsilon\sigma\), we can use a projection \(|\psi\rangle\langle\psi|\) to suppress the noise, by finding coefficients to satisfy the condition \(\mathcal{P}\equiv\sum_{i} c_{i}P_{i} = |\psi\rangle\langle\psi|\). However, a general quantum state \(|\psi\rangle\) may be expanded into an exponential number of Pauli terms, such a condition becomes hardly true given a small subset of S and hence the QSE method is less effective for stochastic errors on a general quantum state. Nevertheless, as we shortly discussion in the next section, when the target state preserves certain symmetries either inherently or by encoding via an error correcting code, the QSE method can indeed effectively suppress stochastic errors by efficiently projecting to the subspace with the correct symmetry. It can thus mitigate a large portion of stochastic errors and significantly improve the calculation accuracy when further combined with other QEM techniques.³⁷⁾

In practice, the set S can be chosen to be creation and annihilation operators for fermionic Hamiltonians. The QSE method was experimentally demonstrated for mitigating errors in quantum chemistry calculation of ground and excited state energies of the H₂ molecule.⁶⁵⁾

When the physical system we try to simulate exhibits certain symmetries,^178,179) they can be exploited for QEM using symmetry verification.^37,38) We can design ansätze with unitary components that conserve the particle number and the spin symmetries, such as the unitary coupled cluster ansatz⁶⁹⁾ or certain hardware efficient ansätze.¹⁸⁰⁾ We will focus on symmetries that can be mapped to Pauli operators since they can be efficiently measured. For example, the parity operators for the particle number and for the spins, \(\hat{P}_{N}\) and \(\hat{P}_{N_{\uparrow\downarrow}}\), both have \(\pm 1\) eigenvalues and thus can be mapped to Pauli operators under a suitable encoding scheme.

While the ideal quantum state preserves the symmetries, states prepared by the noisy circuit may break them. Since symmetries are broken only if an error happens, symmetry verification works by discarding the cases that break the symmetries.^37,38) We can implement the measurement of \(\hat{P}_{N}\) and \(\hat{P}_{N_{\uparrow\downarrow}}\) by implementing the parity check circuit with an ancilla qubit interacting with the register qubit as shown in Fig. 16. Since a general Pauli operator is the tensor product of local Pauli operators, the symmetry can also be measured using local Pauli measurement and post-processing in certain cases without using ancillae.¹⁸¹⁾ If the symmetry outcome is not the expected one, the circuit run is discarded. More errors can be detected by considering more general symmetries, such as the one that preserves the particle number. Symmetry verification can also be combined with other QEM methods as will be discussed in Sect. 5.10.

»View larger version

Figure 16. Quantum circuit used for symmetry verification. This quantum circuit is for four register qubits. The ancilla qubit is measured in 0 (1) when the the total particle number is even (odd). Thus, by reading out the ancilla qubit, we can detect errors in the register qubits.

Alternatively, we can effectively implement the parity check circuit via a post-processing approach.³⁸⁾ Here we consider the conventional VQE scenario and the case where \(\hat{P}_{N}\) is measured, and error-free subspace corresponds to \(\hat{P}_{N} |\psi\rangle = m |\psi\rangle\) with \(m\in\{\pm 1\}\). For any noisy state ρ that may have components with both positive and negative eigenvalues of \(\hat{P}_{N}\), projecting the state to the subspace with correct symmetry is \begin{equation} \rho_{m} = \frac{M_{m} \rho M_{m}}{\mathop{\text{Tr}}\nolimits[M_{m} \rho]}, \end{equation} (133) where \(M_{m} =\frac{I+m P_{N}}{2}\). Supposing we measure the Hamiltonian H, which commutes with \(\hat{P}_{N}\), we have \begin{equation} \mathop{\text{Tr}}\nolimits[H \rho_{m}] = \frac{\mathop{\text{Tr}}\nolimits[H \rho]+ m\mathop{\text{Tr}}\nolimits[H \hat{P}_{N} \rho]}{1+m \mathop{\text{Tr}}\nolimits[\hat{P}_{N} \rho]}. \end{equation} (134) Therefore, we can measure \(\mathop{\text{Tr}}\nolimits[H\rho]\), \(\mathop{\text{Tr}}\nolimits[\hat{P}_{N}\rho]\), and \(\mathop{\text{Tr}}\nolimits[H\hat{P}_{N}\rho]\) to effectively obtain the measurement with the post-selected state \(\rho_{m}\). The state of Eq. (133) is consistent with the subspace states in the previous section. When allowing to vary m, optimising m would be equivalent to the QSE procedure. It indeed has been shown that the solution of QSE coincides with the state of Eq. (133) given the symmetry of the target ground state.³⁸⁾

In contrast to the parity-check circuit based approach, the current one does not need the additional ancilla and the parity-check circuit even in the most general cases. However it does need more samples since the error is mitigated via post-processing. For example, suppose the probability that the noisy state is in the correct subspace is \(p =\mathop{\text{Tr}}\nolimits[M_{m}\rho]\). Then to achieve a sampling error ε, the parity-check circuit method requires \(\mathcal{O}(1/(p\varepsilon^{2}))\) samples and the post-processing method requires \(\mathcal{O}(1/(p^{2}\varepsilon^{2}))\) samples. The post-processing approach also applies for general symmetries such as the particle number.¹⁸²⁾ The post-processing approach was experimentally demonstrated by using a two-qubit superconducting processor.¹⁸³⁾

When considering the logical states of a stabiliser quantum error correcting code, it has code symmetries described by the stabilisers. The conventional approach for realising error detection and error correction is based on the parity check circuit with additional ancillary qubits. In McClean et al.,¹⁸⁴⁾ the post-processing approach was extended for realising the error detection and error correction without the parity check circuit. Since an error correcting code has multiple symmetries, we can probabilistically apply each symmetry projection to force the state back into the code space. The post-processing error decoder for the five qubit quantum error correcting code was numerically implemented with an effective threshold of 50% and applied for chemistry simulation.¹⁸⁴⁾

The individual error reduction method makes use of quantum error correction (QEC) on a single qubit and post-processing to mitigate errors.³⁹⁾ Assume that physical noise after each quantum gate is described using the Lindblad master equation \begin{equation} \begin{split} \frac{d \rho}{dt} &= \sum_{k} \mathcal{L}_{k} (\rho)\\ L_{k}(\rho) &= 2L_{k}^{\dagger} \rho L_{k} -L_{k}^{\dagger} L_{k} \rho -\rho L_{k}^{\dagger} L_{k}, \end{split} \end{equation} (135) and the duration of the noise process is τ. Note that \(L_{k}\) operates on the kth qubit. Suppose that the lth qubit is encoded as a logical qubit and the noise on it can be reduced by a factor of \(h_{l}\) by QEC. Then the evolution of the state is described as \begin{equation} \frac{d \rho_{l}}{dt} = \sum_{k \neq l} \mathcal{L}_{k} (\rho)+(1-h_{l})L_{l}(\rho_{l}). \end{equation} (136) Then, the individual error reduction method defines \begin{equation} \tilde{\rho}^{\textit{est}} \equiv \rho^{\textit{noisy}}- \sum_{l} \frac{1}{h_{l}}(\rho^{\textit{noisy}}-\rho^{\textit{noisy}}_{l}), \end{equation} (137) where \(\rho^{\textit{noisy}}\) is the output state from the noisy quantum circuit without QEC described by Eq. (135) and \(\rho_{l}^{\textit{noisy}}\) is the output state with QEC on the lth qubit described by Eq. (136). Then it can be shown that \(\tilde{\rho}^{\textit{est}}\) approximates the ideal quantum state as \begin{equation} \tilde{\rho}^{\textit{est}} = \rho_{\textit{ideal}}+O(\tau^{2}), \end{equation} (138) where \(\rho_{\textit{ideal}}\) is the error-free output state from the quantum circuit. Thus, we can estimate the ideal expectation value \(\langle M\rangle_{\textit{ideal}} =\mathop{\text{Tr}}\nolimits (\rho_{\textit{ideal}} M)\) for an observable M by measuring the expectation value \(\langle M\rangle =\mathop{\text{Tr}}\nolimits[M\rho^{\textit{noisy}}]\) and \(\langle M\rangle_{l} =\mathop{\text{Tr}}\nolimits[M\rho_{l}^{\textit{noisy}}]\), and post-processing them via \begin{equation} \langle \tilde{M}\rangle_{\textit{est}} = \langle M\rangle-\sum_{l} \frac{1}{h_{l}}(\langle M\rangle-\langle M\rangle_{l}). \end{equation} (139) While the individual error reduction method suppresses errors to the first order, higher order effects can be further mitigated by combining it with other error mitigation methods. We also note that this method requires QEC on a single qubit, thus making it harder to implement compared to other QEM methods.

Measurement error mitigation is designed for suppressing errors during the measurement process.^41–43) Suppose that the ideal measurement is described by a set of positive-operator valued measure (POVM) operators \(\{E_{k}\}\) and the probability to obtain the measurement result k is \(p_{k} =\mathop{\text{Tr}}\nolimits[E_{k}\rho]\) for state ρ. Denote the ideal error-free probability distribution as \(\vec{p}_{\textit{ideal}} = (p_{1}, p_{2},\ldots, p_{N_{P}})^{\text{T}}\) with \(N_{P}\) being the number of POVM elements. In the presence of a measurement error, it transforms the probability distribution as \begin{equation} \vec{p}_{\textit{noise}} = \mathrm{N}\cdot \vec{p}_{\textit{ideal}}, \end{equation} (140) with N being the transformation matrix. For example, considering projective measurement of a qubit with probability \(\vec{p} = (p_{0}, p_{1})\), the misalignment measurement error can be described by \begin{equation} \mathrm{N}= \begin{pmatrix} 1-\varepsilon &\eta\\ \varepsilon &1-\eta \end{pmatrix}, \end{equation} (141) with ε and η denoting the error probabilities that the outcome flip from 0 to 1 and vice versa.

In practice, estimation of the matrix N\(_{\textit{est}}\) can be obtained via quantum detector tomography when there is no detector crosstalk over different qubits,¹⁸⁵⁾ so that N\(_{\textit{est}}\) is a tensor product of transformation matrices. Then we can obtain the ideal error free probability vector via \begin{equation} \vec{p}_{\textit{est}} = \mathrm{N}_{\textit{est}}^{-1}\cdot \vec{p}_{\textit{noise}}. \end{equation} (142) Although, this may produce unphysical probability with negative component due to non-classical noise such as coherent errors. To circumvent this problem, we can adopt \begin{equation} \vec{p}^{*}_{\textit{est}} = \mathop{\mathrm{argmin}}\nolimits_{\vec{p}_{\textit{est}}}(\| \mathrm{N}\cdot \vec{p}_{\textit{est}} -\vec{p}_{\textit{noise}} \|) \end{equation} (143) subject to the constraint that \(\vec{p}_{\textit{est}}\) is normalised and its elements are positive, where \(\|\cdot\|\) denotes Euclidean norm. This method was applied to IBM's five-qubit device to give a significant improvement of results in Maciejewski et al.⁴¹⁾ and Chen et al.⁴²⁾ Recently, measurement error mitigation technique was proposed to further suppress cross-talk errors between qubits during the readout process.¹⁸⁶⁾ This method was demonstrated on the IBM 20-qubit superconducting processor.

While most error mitigation methods rely on exact or partial information of the noise model, the learning-based QEM method aims to suppress errors via an automatical learning process.^44,45,187) Herein, we illustrate two examples by Czarnik et al.⁴⁵⁾ and Strikis et al.⁴⁴⁾

5.8.1 Quantum error mitigation via Clifford data regression

Assume that we are given test data \(\{x^{\textit{ideal}}_{k}, x^{\textit{noisy}}_{k}\}\) with \(\{x^{\textit{ideal}}_{k}\}\) and \(\{x^{\textit{noisy}}_{k}\}\) being the ideal and noisy expectation values of an observable corresponding to the same target quantum circuits. We aim to learn the relationship between \(\{x^{\textit{ideal}}_{k}\}\) and \(\{x^{\textit{noisy}}_{k} \}\) as \(x^{\textit{ideal}} = g(x^{\textit{noisy}},\vec{\theta})\) via regression, and then employ this relationship to infer the unknown ideal result from any given noisy result.⁴⁴⁾ Here, \(\vec{\theta}\) are free parameters to optimise. Note that the test data \(\{x^{\textit{ideal}}_{k}, x^{\textit{noisy}}_{k}\}\) have to be sampled efficiently. We assume that the ideal error-free results \(\{x_{k}^{\textit{ideal}} \}\) are generated by efficiently simulating Clifford quantum circuits on classical computers, and \(\{x^{\textit{noisy}}_{k} \}\) are sampled from Clifford quantum circuits on real noisy quantum circuits. Given the test data, we minimise the cost function, \begin{equation} \mathrm{C}(\vec{\theta}) = \sum_{k} (g(x^{\textit{noisy}}_{k},\vec{\theta})- x^{\textit{ideal}}_{k})^{2}, \end{equation} (144) to fit the relationship between the ideal and noisy measurement results. We can choose a linear ansatz \(g(x^{\textit{noisy}}_{k},\vec{\theta}) =\theta_{1} x^{\textit{noisy}}_{k}+\theta_{2}\), or a more complex ansatz via a neural network. Once the relationship \(x^{\textit{ideal}} = g(x^{\textit{noisy}},\vec{\theta})\) is found, we input a noisy output \(x^{\textit{noisy}}\) obtained from a general quantum circuits and use the output as an estimation of the ideal result \(x^{\textit{ideal}}\). This method is referred to as Clifford Data Regression⁴⁴⁾ and illustrated in Fig. 17. Clifford Data Regression was demonstrated by using the 16-qubit IBMQ quantum computer, and a 64-qubit classical simulator.⁴⁴⁾ Since the noise effect are naturally incorporated in the learning process, this method may work for both local Markovian noise and correlated noise.

»View larger version

Figure 17. (Color online) Schematic figures for Clifford data regression.⁴⁴⁾ The relationship between noisy results from a quantum computer and ideal results from a classical computer is learned via regression. The training data are sampled from Clifford circuits. Then we use the relationship to predict the error mitigated result for an arbitrary non-Clifford quantum circuit.

5.8.2 Learning-based quasi-probability method

Instead of directly learning the ideal measurement outcome from noisy ones, Strikis et al.⁴⁵⁾ introduced a learning-based method for realising quasi-probability error mitigation without depending on process tomography. Focusing on quantum circuits that consist of arbitrary single-qubit gates and Clifford two-qubit gates, which are sufficient for universal quantum computing. Suppose the dominant error source comes from the two-qubit gates, which are called frame gates. Denote the sequence of single-qubit gates as \(\mathbf{U}= \{U_{1}, U_{2},\ldots \}\), and the ideal and noisy circuits as \(\mathcal{E}^{\textit{ideal}}(\mathbf{U})\) and \(\mathcal{E}^{\textit{noisy}}(\mathbf{U})\), respectively. To recover the ideal process, the quasi-probability method applies single-qubit Pauli operations \(\mathbf{Q}= \{Q^{1}, Q^{2},\ldots \}\) before and after gates of each layer. Denote the noisy circuit with additional operations for error mitigation as \(\mathcal{E}^{\textit{noisy}}(\mathbf{U},\mathbf{Q})\). An example of quantum circuits is shown in Fig. 18. The quasi-probability method is to find coefficients \(q(\mathbf{Q})\) so that \begin{equation} \mathcal{E}^{\textit{est}}(\mathbf{U}) = \sum_{\mathbf{Q}} q(\mathbf{Q}) \mathcal{E}^{\textit{noisy}}(\mathbf{U}, \mathbf{Q}) \end{equation} (145) approximates \(\mathcal{E}^{\textit{ideal}}(\mathbf{U})\). The conventional quasi-probability QEM method relies on the exact error channel tomography to determine the coefficients, which may fail to work for spatially and temporally correlated errors whose gate set tomography is in general hard to efficiently implement.

»View larger version

Figure 18. (Color) An example of quantum circuits for learning-based quasi-probability method.⁴⁵⁾ The gates denoted as Q are for mitigating errors for frame gates \(\mathcal{E}\), state preparation errors \(\mathcal{E}_{\textit{st}}\) and measurement errors \(\mathcal{E}_{\textit{me}}\), and U is for sampling a training set \(\mathbb{T}\). We insert Q and U between every layer of gates.

The learning based approach provides an alternative solution without the tomography of the error channel.⁴⁵⁾ Denoting \(\mathrm{com}^{\textit{ideal}(\textit{est})}(\mathbf{U})\) to be the expectation value of the observable of interest corresponding to \(\mathcal{E}^{\textit{ideal}}(\mathbf{U})\) [\(\mathcal{E}^{\textit{est}}(\mathbf{U})\)], we aim to minimise \begin{equation} \mathrm{C} = \frac{1}{|\mathbb{T}|} \sum_{\mathbf{U}\in \mathbb{T}}|\mathrm{com}^{\textit{ideal}}(\mathbf{U})-\mathrm{com}^{\textit{est}}(\mathbf{U}) |^{2}, \end{equation} (146) over the coefficients \(q(\mathbf{Q})\) with a training set \(\mathbb{T}\). This cost function indicates the distance between the correct expectation values and the error-mitigated expectation values. Similarly to the approach proposed by Czarnik et al.,⁴⁴⁾ the training set \(\mathbb{T}\) is set to be the set of Clifford operations \(\mathbb{C}\). Therefore, the entire quantum circuit becomes a Clifford circuit and we can efficiently calculate the ideal expectation values. It is not hard to see that the error mitigation obtained through this procedure also works for arbitrary non-Clifford single-qubit gates with gate-independent errors, because any operation can be expressed as a linear combination of Clifford operations. Since the two-qubit frame gates are fixed, there is no assumption for the error model of frame gates.

Note that the space \(S_{\mathbf{Q}}\) of all the recovery operations \(\mathbf{Q}\) may exponentially increase with the number of gates of the quantum circuit. There are two ways to overcome this problem via the truncation method and variational optimisation method.⁴⁵⁾ For the truncation method, we convert general errors to Pauli errors via Pauli twirling and truncate the space to a subset whose dimension only increases polynomially with the circuit size. Alternatively, we can use the Monte-Carlo method to evaluate the cost function and variationally search the parameters to minimise the cost function. In practice, we can set the starting point of the optimisation to the one specified with the results of imperfect gate set tomography.

The previous error mitigation methods mostly focus on the scenario of digital quantum simulation with discrete noisy gates. Now we review the stochastic error mitigation method for mitigating errors in continuous processes, which works for analog quantum simulation and digital quantum computing.⁴⁶⁾

Suppose the dynamics of the system of interest is described by the Lindblad master equation as \begin{equation} \frac{d}{dt}\rho = - i [H, \rho]+ \mathcal{L}(\rho), \end{equation} (147) with the ideal evolution Hamiltonian H, noise Lindblad operator \(L_{k}\), and \(\mathcal{L}(\rho) =\sum_{k} (2L_{k}\rho L_{k}^{\dagger} -L_{k}^{\dagger} L_{k}\rho-\rho L_{k}^{\dagger} L_{k})\). Suppose the Lindblad operators are local and their strength is weak. Notice that because the effect of Lindblad operators propagates to the entire system due to time evolution, it is generally hard to suppress such a global effect of physical noise via conventional QEM methods. For simplicity, we also assume that the Hamiltonian is time-independent and we note that the result works for general time-dependent Hamiltonians. Now we express the evolution of the state from t to \(t+\delta t\) as \begin{equation} \rho_{i}(t+\delta t) = \mathcal{E}_{i} (\rho_{i}(t)), \end{equation} (148) with \(i =\mathcal{I},\mathcal{N}\) corresponding to the ideal and noisy process. To emulate the ideal evolution of \(\mathcal{E}_{\mathcal{I}}\) with the noisy evolution \(\mathcal{E}_{\mathcal{N}}\), we can implement the quasi-probability method by applying \begin{equation} \mathcal{E}_{\mathcal{I}} = \mathcal{E}_{\mathcal{Q}} \mathcal{E}_{\mathcal{N}}, \end{equation} (149) with the recovery operation \(\mathcal{E}_{\mathcal{Q}}\) decomposed as \begin{equation} \mathcal{E}_{\mathcal{Q}} = \sum_{i} q_{i} \mathcal{B}_{i} = C_{\mathcal{Q}} \sum_{i} \mathop{\text{sgn}}\nolimits(q_{i}) p_{i} \mathcal{B}_{i}. \end{equation} (150) Here, \(C_{\mathcal{Q}} =\sum_{i} |q_{i}|\), \(p_{i} = |q_{i}|/C_{\mathcal{Q}}\), and \(\mathcal{B}_{i}\) are tensor products of single qubit operators. To simulate the whole ideal evolution of time T, we can continuously apply the recovery \(\mathcal{E}_{\mathcal{Q}}\) with a time interval \(\delta t\) with a cost \(C_{T} = C_{\mathcal{Q}}^{N_{d}}\) where \(N_{d} = T/\delta t\). The schematic figures are shown in Fig. 19.

»View larger version

Figure 19. (Color) Schematic figures for applying quasi-probability operations continuously.⁴⁶⁾ Note that we can obtain the result corresponding to \(\delta t\rightarrow 0\) by using stochastic error mitigation method.

However, continuously applying the recovery operation could be challenging. Note that each recovery operation for a small \(\delta t\) is almost an identity operation. Therefore, instead of continuously applying every weak recovery operations, we can use the Monte Carlo method to stochastically realise strong recovery operations in a similar vein of stochastic Schrödinger equation. The stochastic error mitigation method can be combined with the extrapolation QEM method to further mitigate the model estimation error.⁴⁶⁾

Here we illustrate ways to combine different error mitigation techniques. More specifically we focus on possible combinations of error extrapolation, quasi-probability, and symmetry verification studied by Cai.⁴⁷⁾

5.10.1 Symmetry verification with error extrapolation

Symmetry verification cannot detect errors that stand for transformations within the same symmetry subspace. When acting on the eigenstate of the symmetry, such errors are the errors that commute with the symmetry operator, and we refer to these errors as commuting errors. Similarly we also have anti-commuting errors, and individual occurrence of them would be detectable by symmetry verification. However, when anti-commuting errors occur an even number of times, they will commute with the symmetry and thus cannot be detected. The noisy expectation value of an observable M at the error rate ε after applying symmetry verification will be denoted as \(\langle M\rangle_{s}(\varepsilon)\). As shown in Fig. 20(a), \(\langle M\rangle_{s}(\varepsilon)\) still contains undetectable errors, which can be further mitigated by combining symmetry verification with error extrapolation using \begin{equation} \langle M\rangle_{\textit{est}} = \frac{\alpha \langle M\rangle_{s}(\varepsilon)-\langle M\rangle_{s}(\alpha \varepsilon)}{\alpha -1}. \end{equation} (151) Here we have used linear extrapolation as an example. Its effectiveness has been numerically verified by McArdle et al.³⁷⁾ and its application to Hubbard VQE has been discussed by Cai.¹⁸¹⁾

»View larger version

Figure 20. (Color) Diagrams showing how different error components are suppressed using different combinations of quantum error mitigation techniques: (a) using only symmetry verification, (b) combining quasi-probability and symmetry verification, and (c) combining quasi-probability, symmetry verification and error extrapolation.