Quantum gradient descent algorithms for nonequilibrium steady states and linear algebraic systems

The dynamics of density matrix, $\rho\in\mathbb{C}^{N\times N}$, $N=2^{n}$, is governed by a quantum master equation in the Lindblad form lindbald1976on ,

where $\mathcal{L}$ is the Liouvillian superoperator such that

$t$ is the time and $H$ is the Hamiltonian operator. The term $[H,\rho]$ describes the unitary part of the dynamics, while $\mathcal{D}\rho$ gives the non-unitary time evolution governed by dissipative superoperators,

Here each $L_{k}$ is the $k$th jump operator associated with a dissipative channel occurring at rate $\mu_{k}$. Let us consider the special but quite common case in which the Hamiltonian $H$ and the jump operator $L_{k}$ is usually expressed as the sum of tensor products lloyd1996universal ; mahdian2020hybrid ,

Under the Choi-Jamiolkowski isomorphism havel2003robust ; yoshioka2019constructing ; kamakari2021digital , $\rho=\sum_{i}p_{i}|\Psi_{i}\rangle\langle\Psi_{i}|\rightarrow|\rho\rangle=\sum_{i}p_{i}|\Psi_{i}\rangle\otimes|\Psi_{i}\rangle$, Eq. (3) is written as,

with the Liouville operator

where $L_{k}^{*}$ denotes the complex conjugate of $L_{k}$ and $I_{N}$ is the $N\times N$ identity matrix. In this representation, the time evolution is formally solved as $|\rho(t)\rangle=e^{\mathcal{H}t}|\rho(0)\rangle$ with the initial state $|\rho(0)\rangle$.

The NESS $\rho_{\textrm{ss}}$, $\dot{\rho}_{\textrm{ss}}=\mathcal{L}\rho_{\textrm{ss}}=0$ or $|\rho_{\textrm{ss}}\rangle=\lim_{t\rightarrow\infty}|\rho(t)\rangle$ schirmer2010stabilizing ; minganti2018spectral , corresponds to a pure state $|\rho_{\textrm{ss}}\rangle$ satisfying $\mathcal{H}|\rho_{\textrm{ss}}\rangle=0$. Since the operator (III.1) is not hermitian, we consider instead the hermitian operator $\mathcal{H}^{{\dagger}}\mathcal{H}\succeq 0$. It is straightforward to see that $|\rho_{\textrm{ss}}\rangle$ also be the ground state of $\mathcal{H}^{{\dagger}}\mathcal{H}$. Thus, the objective function is defined as the expectation value of $\mathcal{H}^{{\dagger}}\mathcal{H}$,

We here use QGD to minimize $f(|\rho\rangle)$. The QGD iterative process can be regraded as,

Therefore, according to Eqs. (6) and (III.1), $D$ can be expressed as a sum of $M=\mathcal{O}(poly(n))$ local terms, $D=\sum_{m=0}^{M-1}d_{m}\hat{D}[m]\in\mathbb{C}^{N^{2}\times N^{2}}$, where the unitary $\hat{D}[m]$ consists of some local Pauli operators. After normalization, we obtain the iterative equation such that

where the normalized constant,

can be estimated from the expectation value of $\hat{D}^{{\dagger}}[m_{1}]\hat{D}[m_{2}]$. The Algorithm 1 below and the Fig. (1.a) outline the preparation of NESS with QGD.

Our QGD algorithm requires two quantum registers, the register 1 and register 2. The input is a tensor product state $|\Psi_{0}\rangle=|0^{\widetilde{m}}\rangle_{1}|\rho^{(0)}\rangle_{2}$, $\widetilde{m}=\log M$. The state $|\rho^{(0)}\rangle$ is an initial guess state picked from easily prepared states such as $|\rho^{(0)}\rangle_{2}=|0^{2n}\rangle_{2}$.

In step 1, we employ the amplitude encoding method to prepare a $\widetilde{m}$-qubit input state $\sum_{m=0}^{M-1}\frac{d_{m}}{\sqrt{\mathcal{N}_{D}}}|m\rangle$, where $\mathcal{N}_{D}=\sum_{m=0}^{M-1}d_{m}^{2}$. The preparation takes $\mathcal{O}(poly(\widetilde{m}))$ steps if $d_{m}$ and $\mathcal{N}_{D}$ can be efficiently computed by classical algorithms Grover2002Creating ; Soklakov2006Efficient . The preparation process can be carried out by applying the following unitary gate on state $|0^{\widetilde{m}}\rangle$,

where the elements $\{w_{0,1},\cdots,w_{M-1,M-1}\}$ are arbitrary so long as $W$ is unitary.

In step 2, we implement the $\widetilde{m}$-qubit-controlled unitary,

on state $|\Psi_{1}\rangle$ to generate the entangled state

The non-unitary gradient operator $D$ is realized by an extended unitary gate on a larger Hilbert space.

Without loss of generality, $\hat{D}[m]$ can be expressed as $\hat{D}[m]=I\otimes\hat{P}_{1}\otimes\cdots\otimes\hat{P}_{J}\otimes I$, $\quad J<2n$, where the Pauli operator $\hat{P}_{j}\in\{X,Y,Z\}$ acts on $j$th qubit. The multi-qubit-controlled unitary gate $\wedge_{\widetilde{m}}\hat{D}$ can be written as

Namely, $\wedge_{\widetilde{m}}\hat{D}$ can be decomposed into multi-qubit-controlled Pauli operators $\wedge_{\widetilde{m}}U$, $U\in\{\hat{P}_{1},\cdots,\hat{P}_{J}\}$. Specifically, $\wedge_{\widetilde{m}}\hat{D}$ can be implemented by $\mathcal{O}(MJ)$ multiple-qubit-controlled $\textrm{SU}(2)$ unitary operators.

Last step, we apply $\widetilde{m}$ Hadamard gates on register 1 and the system state now becomes a separable one,

Suppose that the steps 1-3 are repeated $S$ times. The final state of the system is given by

When the register 1 is in state $|0^{\widetilde{m}}\rangle$, the state $|\rho^{(S)}\rangle$ could be a better approximation of the NESS $|\rho_{\textrm{ss}}\rangle$ if the $|\rho^{(S)}\rangle$ satisfies $|f(|\rho^{(S)}\rangle)|\leq\varepsilon$. Otherwise, repeat the overall procedures again until the convergence condition is satisfied. The success probability of obtaining $|0^{\widetilde{m}}\rangle$ is given by

which decreases exponentially with respect to the number of iterations. Performing quantum amplitude amplification brassard2000quantum , we have that the number of measurements is at most

where $L=\lfloor\pi/4\theta\rfloor$ and $\sin^{2}\theta=P_{\textrm{suc}}$. The final success probability is at least $\max(P_{\textrm{suc}},1-P_{\textrm{suc}})$.

In addition to focusing on the NESS of open quantum systems, it is also important to provide an efficient scheme to compute expectation values of physical observables. Our algorithm yields a quantum state $|\rho_{\textrm{ss}}\rangle$ related to the density matrix $\rho_{\textrm{ss}}$. Note that in our approach the state is normalized, i.e., $\langle\rho|\rho\rangle=1$. In general, the condition on the trace $\textrm{Tr}[\rho]=\langle I_{N}|\rho\rangle=1$ is not automatically fulfilled. Thus, the expectation value of an arbitrary observable $\mathcal{M}$ must be estimated as,

with respect to the density matrix $\rho_{\textrm{ss}}$ for a given observable $\mathcal{M}$. Here, the super-operator $\hat{\mathcal{M}}=\mathcal{M}\otimes I_{N}$ and the maximally entangled state

which is prepared by applying the unitary $U_{I_{N}}$ on an initial state $|0^{2n}\rangle$. The quantum circuit of the unitary $U_{I_{N}}$ consists of the Hadamard gate $H$ and CNOT gate, see Fig. (2.b.) In particular, by setting $\mathcal{M}=I_{N}$, the trace of NESS $\rho_{\textrm{ss}}$ can be calculated by $\textrm{Tr}[\rho_{\textrm{ss}}]=\langle I_{N}|\rho_{\textrm{ss}}\rangle$.

We propose two strategies to estimate the expectation value (19) from the output state of the algorithm. While the swap test calculates the square of the inner product buhrman2001quantum ; fanizza2020beyond , the magnitude and sign of the inner product are also required in our specific cases.

Strategy 1. Our first strategy calculates the real and imaginary parts separately based on the modified swap test zhao2019quantum or Hadamard test aharonov2008a .

Start with an initial state $\frac{1}{\sqrt{2}}(|0\rangle+\zeta|1\rangle)$ with $\zeta=1$ or $\zeta=i$. We prepare the state $|\phi_{0}\rangle=\frac{1}{\sqrt{2}}(|0\rangle|I_{N}\rangle+\zeta|1\rangle|\rho_{\textrm{ss}}\rangle)$ by applying the control unitary $\wedge_{1}U$ (see Proposition 1) on the state $(1/\sqrt{2})[|0\rangle+\zeta|1\rangle]|0^{2n}\rangle$. After applying the Hadamard gate on the first qubit, we obtain the state

Now we measure the qubit system in the basis of Pauli operator $Z$. The measurement output is a random variable $\in\mathbb{Z}=\{r_{\pm}=\pm 1\}$, where the output $r_{+}$ $(r_{-})$ corresponds to the measurement output state $|0\rangle$ $(|1\rangle)$. The expectation value of $\mathbb{Z}$ is given by

If $\zeta=1$, we have $\mathbb{E}(\mathbb{Z})=\textrm{Re}[\langle I_{N}|\rho_{\textrm{ss}}\rangle]$. If $\zeta=i$, then $\mathbb{E}(\mathbb{Z})=\textrm{Im}[\langle I|\rho_{\textrm{ss}}\rangle]$. Let $P_{+}$ be the probability of measurement outcome $r_{+}$. To estimate the probability $P_{+}$, we perform independent trials of the Bernoulli test. Then sample $R$ times and record the number $R_{0}$ of observed $|0\rangle$. The frequency $R_{0}/R$ then gives an estimator for $P_{+}$ up to a sampling error $\tilde{\varepsilon}$. From Hoeffding’s inequality hoeffding1963 , we have

We obtain the number of measurements, $R=\mathcal{O}(\widetilde{\varepsilon}^{-2}\log\delta^{-1})$. Therefor we have the following conclusion.

Proposition 1. Let $\wedge_{1}U=|0\rangle\langle 0|\otimes U_{I_{N}}+|1\rangle\langle 1|\otimes U_{\rho_{\textrm{ss}}}$ be the controlled unitary gate with $U_{I_{N}}|0^{2n}\rangle=|I_{N}\rangle$ and $U_{\rho_{\textrm{ss}}}|0^{2n}\rangle=|\rho_{\textrm{ss}}\rangle$. There exists a quantum algorithm that calculates $\textrm{Re}[\langle I_{N}|\rho_{\textrm{ss}}\rangle]$ and $\textrm{Im}[\langle I_{N}|\rho_{\textrm{ss}}\rangle]$ to accuracy $\widetilde{\varepsilon}$ with success probability at least $1-\delta$ by using $\mathcal{O}(\widetilde{\varepsilon}^{-2}\log\delta^{-1})$ queries of $\wedge_{1}U$.

Next we estimate the trace average $\textrm{Tr}[\mathcal{M}\rho]=\langle I_{N}|\hat{\mathcal{M}}|\rho\rangle$.

Proposition 2. Denote the hermitian operator $\mathbb{M}=X\otimes\hat{\mathcal{M}}$ and the controlled unitary gate $\wedge_{1}U=|0\rangle\langle 0|\otimes U_{I_{N}}+|1\rangle\langle 1|\otimes U_{\rho_{\textrm{ss}}}$ with $U_{I_{N}}|0^{2n}\rangle=|I_{N}\rangle$ and $U_{\rho_{\textrm{ss}}}|0^{2n}\rangle=|\rho_{\textrm{ss}}\rangle$. The success probability to calculate $\langle I_{N}|\hat{\mathcal{M}}|\rho_{\textrm{ss}}\rangle$ to accuracy $\widetilde{\varepsilon}$ via measuring the expectation value of $\mathbb{M}$ on state $|\phi_{0}\rangle=\frac{1}{\sqrt{2}}(|0\rangle|I_{N}\rangle+\zeta|1\rangle|\rho_{\textrm{ss}}\rangle)$ is at least $1-\delta$ by using $\mathcal{O}(\widetilde{\varepsilon}^{-2}\log\delta^{-1})$ queries of $\wedge_{1}U$.

Proposition 2 can be proved by noting that

due to the structure of the operator $\mathbb{M}$. Therefore, similar to the proof of Proposition 1, the number of measurements is $\mathcal{O}(\widetilde{\varepsilon}^{-2}\log\delta^{-1})$.

In the last two Propositions we have assumed that the state $|\rho_{\textrm{ss}}\rangle$ can be prepared by applying a unitary $U_{\rho_{\textrm{ss}}}$ on an initial state $|\Phi_{0}\rangle=|0^{2n}\rangle$. However, as shown in Algorithm 1, our state $|\rho_{\textrm{ss}}\rangle$ is prepared via a unitary evolution (denoted as $V$) coupling to a measurement operator $|0^{\widetilde{m}}\rangle\langle 0^{\widetilde{m}}|\otimes I_{N}$. The overall process corresponds to a non-unitary process,

Thus, in order to employ the unitary $\wedge_{1}U$ introduced in the two Propositions, we need to find a unitary $U_{\rho_{\textrm{ss}}}$ to approximate this non-unitary process $\mathcal{E}$. Here, we propose a variational quantum non-unitary process approximation to find the unitary $U_{\rho_{\textrm{ss}}}$. Let $|\Phi_{0}\rangle$ be an initial system state. The evolved (normalized) state under the non-unitary operator $\mathcal{E}$ is given by $|\Phi\rangle=\frac{\mathcal{E}|\Phi_{0}\rangle}{\textrm{Tr}(\mathcal{E}|\Phi_{0}\rangle)}$. The main idea is to train a parameterized quantum circuit $\mathcal{U}(\bm{\alpha})$, $\bm{\alpha}=(\alpha_{1},\alpha_{2},\cdots)$ via minimizing an objective function which quantifies the difference between the states $|\Phi\rangle$ and $\mathcal{U}(\bm{\alpha})|\Phi_{0}\rangle$. $\mathcal{U}(\bm{\alpha})$ consists of single-qubit rotation gates and two-qubit CNOT gates havlivcek2019supervised ; commeau2020variational .

Let us have a detailed analysis on the non-unitary approximation algorithm. Our variational quantum non-unitary process approximation outputs a unitary $U_{\rho_{\textrm{ss}}}$ to approximately construct the non-unitary process $\mathcal{E}$ on a given initial state. A natural objective function is the square of the trace distance between $|\rho_{\textrm{ss}}\rangle$ and $|\rho(\bm{\alpha})\rangle=\mathcal{U}(\bm{\alpha})|\Phi_{0}\rangle$, i.e.,

which can be evaluated by measuring the expectation of the observable $O=I_{N^{2}}-|\Phi_{0}\rangle\langle\Phi_{0}|$ with respect to the state $\mathcal{U}(\bm{\alpha})|\rho_{\textrm{ss}}\rangle$. It is straightforward to verify that the objective function can be expressed as

The objective function is faithful if and only if $|\rho_{\textrm{ss}}\rangle=\mathcal{U}(\bm{\alpha})|\Phi_{0}\rangle$.

The parameters $\bm{\alpha}$ are trained through a gradient-based optimization method. The optimal parameters can be obtained by updating $\bm{\alpha}=\bm{\alpha}-\eta\nabla F(\bm{\alpha})$, where $\eta$ is the learning rate and $\nabla F(\bm{\alpha})=(\partial_{1}F,\partial_{2}F,\cdots)$. The gradient information can be computed on a quantum computer via the finite difference formulae,

where $\Delta\alpha_{i}$ is a small perturbation to $\bm{\alpha}$. The minimization $\bm{\alpha}_{\textrm{min}}=\textrm{arg}\min_{\bm{\alpha}}F(\bm{\alpha})$ gives rise to the optimal unitary $U_{\rho_{\textrm{ss}}}=\mathcal{U}(\bm{\alpha}_{\textrm{min}})$.

Strategy 2. Our second strategy provides an directly estimation of (19) from the output $|\Psi_{\textrm{final}}\rangle$. Moreover, strategy 2 circumvents local measurement and unitary approximations, thus further reduces the computational complexity.

Proposition 3. Define the hermitian operator $\hat{\mathbb{M}}=X\otimes I_{M}\otimes\hat{\mathcal{M}}$ and the controlled unitary $\wedge_{1}\widetilde{U}=|0\rangle\langle 0|\otimes V+|1\rangle\langle 1|\otimes I_{M}\otimes U_{I_{N}}$, with $V|0^{\widetilde{m}}\rangle|0^{2n}\rangle=|\Psi_{\textrm{final}}\rangle$ and $U_{I_{N}}|0^{2n}\rangle=|I_{N}\rangle$. Then we can calculate $\langle I_{N}|\hat{\mathcal{M}}|\rho_{\textrm{ss}}\rangle$ to accuracy $\widetilde{\varepsilon}$ via measuring the expectation value of $\hat{\mathbb{M}}$ on the state

The success probability is at least $1-\delta$ by using $\mathcal{O}(\widetilde{\varepsilon}^{-2}\log\delta^{-1})$ queries of $\wedge_{1}\widetilde{U}$.

Concerning the proof of Proposition 3 it is straightforward to check that the expectation of the observable $\hat{\mathbb{M}}$ is given by

In particular, by setting $\hat{\mathcal{M}}=I_{N}$, the measurement returns to the quantity $\textrm{Tr}[\rho_{\textrm{ss}}]=\langle I_{N}|\rho_{\textrm{ss}}\rangle$.

The above strategies do not require the prior knowledge of the spectrum of $\rho$ obtained by quantum phase estimate liang2019quantum and quantum state tomography cramer2010efficient . Thus our approach is experimentally friendly in near-term quantum devices.

In this section we present the resource complexity and error analysis. The computation complexity contains (1) the gate complexity, the number of the basic quantum gates (single and two qubit gates); (2) the qubit complexity, the number of qubits.

Proposition 4. Concerning the computational complexity of the preparation of $|\rho_{\textrm{ss}}\rangle$, the gate complexity of finding the NESS is $\mathcal{O}(nM\widetilde{m}^{2})$, where $M$ is the number of decomposition terms of the non-unitary gradient operator $D$ and $\widetilde{m}=\log M$. The qubit complexity is counted as $\mathcal{O}(\tilde{m}+2n)$ with $n=\log N$.

[Proof]. The dominated source of the computational complexity is the controlled unitary $\mathcal{C}_{\widetilde{m}}(\hat{D})$ in Algorithm 1. The unitary $\mathcal{C}_{\widetilde{m}}(\hat{D})$ can be simulated by a sequence of single qubit controlled gate and Toffoli gates barenco1995elementary ; Long2006Mathematical ; Wen2020One ; xin2017quantum ; wei2018efficient . Specifically, performing Algorithm 1 one time requires $M$ controlled unitary operator $\wedge_{\widetilde{m}}\hat{D}$. Thus the circuit depth of obtaining NESS is $\mathcal{O}(M)$. Under the Pauli decomposition of $D$, the controlled unitary operations $\wedge_{\widetilde{m}}\hat{D}$ can always be rewritten as $\mathcal{O}(n)$ multi-qubit-controlled Pauli operators $\wedge_{\widetilde{m}}U$, $U\in\{X,Y,Z\}$. A general $\wedge_{\widetilde{m}}U$ for any unitary $2\times 2$ matrix $U$ can be approximately simulated by a network (Fig. 2.a) including of unitaries $\wedge_{1}V$, $\wedge_{1}V^{{\dagger}}$ and the Toffoli gate $T$ barenco1995elementary . The number of basic operators scales to $\Theta(\widetilde{m}+1)$. Furthermore, the Toffoli gate $T$ on $\widetilde{m}-1$ qubits can be decomposed into $\mathcal{O}(\widetilde{m}-1)$ single-qubit and CNOT gates. The unitary operators $\wedge_{1}V$ and $\wedge_{1}V^{{\dagger}}$ can be simulated with a cost $\mathcal{O}(1)$. Denote $\mathcal{T}_{\widetilde{m}}$ the gate complexity of decomposing $\wedge_{\widetilde{m}}U$. We have the following recurrence relation, $\mathcal{T}_{\widetilde{m}}=\mathcal{T}_{\widetilde{m}-1}+\mathcal{O}(1)+\mathcal{O}(\widetilde{m}-1)$, which implies that $\mathcal{T}_{\widetilde{m}}=\mathcal{O}(\widetilde{m}^{2})$. Hence, the total gate complexity of simulating $\mathcal{C}_{\widetilde{m}}(\hat{D})$ is roughly

As shown in Fig. (1.a), the register 2 contains $2n$ qubits used to store the NESS and the state $|I_{N}\rangle$. The ancillary system requires $\widetilde{m}$ qubits determined by $M$. Thus the total qubit complexity scales to $\mathcal{O}(\widetilde{m}+2n)$. $\Box$

Proposition 5. Concerning the computational complexity of estimating the expectation value $\langle\mathcal{M}\rangle$, the gate and the qubit complexity are $\mathcal{O}(\textrm{poly}(n)+nM\widetilde{m}^{2})$ and $\mathcal{O}(2n)$ ($\mathcal{O}(nM\widetilde{m}^{2})$ and $\mathcal{O}(\widetilde{m}+2n)$) for strategy 1 (strategy 2), respectively.

[Proof]. For both strategies, it is direct to check that the qubit complexity scales to $\mathcal{O}(2n)$ and $\mathcal{O}(\widetilde{m}+2n)$. The gate complexity for strategy 1 is determined by the parameterized quantum circuit $\mathcal{U}(\bm{\alpha}^{*})$ which requires $\mathcal{O}(\textrm{poly}(n))$ basic quantum gates HardwareVQE2017 ; liang2020variational . Therefore, by using the result of Proposition 4, the total gate complexity for strategy 1 is $\mathcal{O}(\textrm{poly}(n)+nM\widetilde{m}^{2})$. For strategy 2, we only require to perform the Algorithm 1 to prepare the state $|\Phi_{\textrm{final}}\rangle$. According to Propositions 3 and 4, we obtain the gate complexity $\mathcal{O}(nM\widetilde{m}^{2})$. $\Box$

Now we determine the upper bound of the error after $S$ time iterations. Suppose that the Algorithm 1 runs accurate enough such that

Let $|\rho_{\textrm{ss}}\rangle$ be the normalized ideal NESS and $|\widetilde{\rho}_{\textrm{ss}}\rangle=|\rho^{S}\rangle$ the actual state.

Proposition 6. Denote $\tau_{1}^{2}=|\langle\phi_{1}|\rho^{(0)}\rangle|^{2}$ the overlap between an initial state $|\rho^{(0)}\rangle$ and the state $|\phi_{1}\rangle$ is the largest eigenstate of operator $D$. The approximation error $\varepsilon$ has an upper bound

scaling with $\tau_{1}^{2}$, iterations $S$, the gap between largest and lowest eigenvalues of $\mathcal{H}^{{\dagger}}\mathcal{H}$, $\kappa=\lambda_{\textrm{max}}-\lambda_{1}$. $\delta_{1}$ and $\delta_{2}$ are the two largest eigenstates of the operator $D$.

[Proof]. Consider the eigenvalue decomposition of the hermitian matrix $D=\sum_{r=1}^{N^{2}}\delta_{r}|\phi_{r}\rangle\langle\phi_{r}|$ with the eigenvalues arranged in order, $\delta_{1}\geq\delta_{2}\geq\cdots\geq\delta_{N^{2}}=\delta_{\textrm{min}}\geq 0$. The eigenvalues of $\mathcal{H}^{{\dagger}}\mathcal{H}$ are then $\lambda_{r}=(1-\delta_{r})/2\gamma$, $\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{N^{2}}=\lambda_{\textrm{max}}$. We express the initial state $|\rho^{(0)}\rangle$ in the basis of the eigenvectors $\{|\phi_{r}\rangle\}$,

with $\tau_{r}$ some coefficients. After applying the non-unitary operator $D$ to $|\rho^{(0)}\rangle$ $S$ times, we have

As $\lim_{S\rightarrow\infty}(\delta_{r}/\delta_{1})^{S}=0$ for all $r$, for larger $S$ we can approximate the ground state of $\mathcal{H}^{{\dagger}}\mathcal{H}$ as

where the normalization constant

The error satisfies