Quantum algorithm for calculation of transition amplitudes in hybrid tensor networks

Quantum computers are expected to be capable of executing classically intractable calculations [1, 2, 3, 4, 5, 6, 7, 8]. It has been reported that quantum computers can outperform classical computers in some tasks [9, 10, 11]. However, quantum resource limitations are obstacles to the practical application of quantum computers. Current quantum computers are so-called noisy intermediate-scale quantum (NISQ) devices [12], and we can control only tens to hundreds of noisy qubits on them [13, 14, 15, 4, 16, 17, 18]. Their hardware limitation makes it difficult to apply quantum computers to practical tasks that require large numbers of qubits or deep quantum circuits [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31].

The hybrid tensor network approach has recently been proposed to overcome the limitations of NISQ devices  [19]. The approach enables the treatment of a larger number of quantum states than the number of qubits of actual quantum devices by representing quantum states as a combination of conventional classical tensors and quantum tensors. A quantum tensor has upper and lower indices, for example, $\psi_{i_{1},i_{2}\dots i_{K}}^{j_{1},j_{2}\dots j_{L}}$, which represents $K$-qubit systems indexed by an $L$-bit string. In other words, the quantum state is defined in terms of $L$ classical bits $(j_{1},j_{2},\dots,j_{L})$ as

where $\psi_{i_{1},i_{2},\dots,i_{K}}^{j_{1},j_{2},\dots,j_{L}}\in\mathbb{C}$, $\sum_{i_{1},i_{2},\dots,i_{K}}|\psi_{i_{1},i_{2},\dots,i_{K}}^{j_{1},j_{2},\dots,j_{L}}|^{2}=1$, and $\ket{i_{1}i_{2}\dots i_{K}}$ is a computational basis of the $K$-qubit Hilbert space. One of the most vital forms of the hybrid tensor network is the hybrid tree-tensor network, where a network of quantum and classical tensors composes a tree graph. While contraction of general hybrid tensor networks can be costly, hybrid tree-tensor networks can be contracted efficiently to obtain expectation values [19]. In this paper, we mainly discuss a two-layer hybrid tree-tensor network with only quantum tensors, which is called a quantum-quantum tree-tensor network, since it captures the essential properties of hybrid tree-tensor networks and the use of classical tensors restricts the range of representation to classically tractable quantum states such as matrix product states [32]. A quantum-quantum tree-tensor network that expresses $k$ subsystems of $n$-qubit states is represented as

where $\ket{\psi_{HT}}$ is the unnormalized state, $\psi_{i_{1},i_{2},\dots,i_{k}}=\braket{i_{1}i_{2}\dots i_{k}}{\psi}$ is the probability amplitude of a $k$-qubit state $\ket{\psi}$, and $\ket{\varphi^{i_{m}}}(m=1,2,\dots,k)$ is an $n$-qubit state of the $m$-th subsystem indexed by a classical bit $i_{m}$. Figure 1 shows the network diagram of $\ket{\psi_{HT}}$. While the state $\ket{\psi_{HT}}$ represents a $kn$-qubit state, we can efficiently calculate the expectation value of an observable $O=\bigotimes_{m=1}^{k}O_{m}$, where $O_{m}$ operates on $\ket{\varphi^{i_{m}}}$ via proper tensor contractions, using a quantum computer with $O(\max(k,n))$ qubits. The contraction for evaluating the expectation value $\braket{O}=\bra{\psi_{HT}}O\ket{\psi_{HT}}$ (without normalization) can be implemented as follows. First, we measure $O_{m}$ with multiple states and construct operators $M_{m}^{i^{\prime}_{m}i_{m}}=\bra{\varphi^{i^{\prime}_{m}}}O_{m}\ket{\varphi^{i_{m}}}$. Note that each $M_{m}$ is a $2\times 2$ Hermitian matrix when $O_{m}$ is Hermitian. Then, since $M_{m}$ is Hermitian, we can measure the expectation value of $\bigotimes_{m=1}^{k}M_{m}$ for the state $\ket{\psi}$, which is equal to $\braket{O}$.

The approach allows for simulations beyond the scale of the quantum hardware. For example, the energy and spin-spin correlation functions of electrons can be calculated with this approach. However, the approach has a serious problem preventing it from being used in a range of applications: it can only calculate the expectation value of observables. In other words, there is a problem in the approach when the contracted operator $M_{m}$ is non-Hermitian. Hereinafter, we denote the Hermitian and non-Hermitian contracted operators as $M_{m}$ and $N_{m}$, respectively. The reason for the problem is that the number of terms to be measured increases exponentially with $k$ when $\bigotimes_{m=1}^{k}N_{m}$ is calculated naively. Specifically, non-Hermitian operators are decomposed into sums of Pauli operators, as in $\bigotimes_{m=1}^{k}N_{m}=\bigotimes_{m=1}^{k}(h_{Im}I_{m}+h_{Xm}X_{m}+h_{Ym}Y_{m}+h_{Zm}Z_{m})$, where $I_{m}$ is the identity operator, $X_{m}$, $Y_{m}$, and $Z_{m}$ are Pauli operators which act on the $m$-th qubit, and $h_{\alpha}(\alpha\in\{I_{m},X_{m},Y_{m},Z_{m}\})$ are the corresponding coefficients, with up to $4^{k}$ terms appearing. One example where the problem occurs is in the calculation of the transition amplitude related to Green’s functions and photon emissions/absorptions [33, 34]. Overlap of two quantum states, which is exploited in subroutines in many algorithms [35, 36, 37, 38], is a special case of the transition amplitude. Thus, the difficulty of computing the expectation value of a non-Hermitian operator limits the applicability of the hybrid tensor network approach.

In this study, we propose a method for calculating transition amplitudes with the hybrid tensor network approach. The main point of the method is the treatment of tensor products of non-Hermitian operators. In a naive calculation, an exponential number of terms in $\bigotimes_{m=1}^{k}N_{m}$ will appear. We propose two ways to avoid the problem. One is a Monte-Carlo contraction method and the other is a singular value decomposition (SVD) contraction method, and the second method is the main proposal in this paper. Although the first method can avoid having to measure all the terms whose number increases exponentially with $k$, the second method is exponentially more efficient than the first one in terms of the sampling cost.

In the following, we give an overview of quantum-quantum tree-tensor networks [19] for obtaining the expectation values of observables in Sec. II. Then, the method of calculating transition amplitudes and overlaps is described in Sec. III; the contraction of quantum tensors in subsystems is discussed in Sec. III.1 and the contraction of non-Hermitian matrices in Sec. III.2. We discuss the future applications of the method in Sec. IV.

We present an overview of the hybrid tensor network simulation on the state described by Eq. (2) introduced in Ref. [19]. Letting the observable $O$ be $O=\bigotimes_{m=1}^{k}O_{m}$ and $O_{m}=\bigotimes_{r=1}^{n}O_{mr}$ ($r=1,2,\dots,n$), the expectation value of the observable including the normalization constant can be described as

and

where $A$ is a normalization constant, $\vec{i}=(i_{1},i_{2},\dots,i_{k})$ with $i_{m}$ taking either $0$ or $1$, $M_{m}^{i_{m}^{\prime}i_{m}}=\bra{\varphi^{i_{m}^{\prime}}}O_{m}\ket{\varphi^{i_{m}}}$, and $M_{Am}^{i_{m}^{\prime}i_{m}}=\braket{\varphi^{i_{m}^{\prime}}}{\varphi^{i_{m}}}$. $M_{m}^{i_{m}^{\prime}i_{m}}$ ($M_{Am}^{i_{m}^{\prime}i_{m}}$) is an element of the $2\times 2$ matrix $M_{m}$ ($M_{Am}$). When $O_{m}$ is assumed to be an observable, i.e., a Hermitian operator, $M_{m}$ also becomes a Hermitian operator. $M_{Am}$ is a Hermitian operator since $M_{Am}$ is a special case of $O_{m}=I_{m}$ in $M_{m}$, where $I_{m}$ is the identity operator.

The procedure for constructing $M_{m}$ and $M_{Am}$ depends on how the indices $i_{m}$ of the wave function $\ket{\varphi^{i_{m}}}$ are mapped. We suppose there are two cases of the mapping; one case is where the indices $i_{m}$ are mapped to unitary gates, i.e., $\ket{\varphi^{i_{m}}}=U_{Cm}^{i_{m}}\ket{0}^{\otimes n}$; the other case is where the indices $i_{m}$ are mapped to the initial wave function as $\ket{\varphi^{i_{m}}}=U_{Cm}\ket{{i_{m}}}\ket{0}^{\otimes n-1}$, where $U_{Cm}^{i_{m}}$ and $U_{Cm}$ are unitary operators with polynomial depth in the $m$-th subsystem. The second case can be regarded as a special example of the first case, since $\ket{\varphi^{i_{m}}}=U_{Cm}\ket{{i_{m}}}\ket{0}^{\otimes n-1}=U_{Cm}(X_{1})^{i_{m}}\ket{0}^{\otimes n}$ and we can think of $U_{Cm}(X_{1})^{i_{m}}$ as $U_{Cm}^{i_{m}}$, where $X_{1}$ is a Pauli $X$ operator which acts on the first qubit. Note that the first method needs a Hadamard test circuit while $M_{m}$ and $M_{Am}$ can be efficiently constructed via direct measurements in the second case, as will be described later.

First, we consider the construction of $M_{m}$ in the case of $\ket{\varphi^{i_{m}}}=U_{Cm}^{i_{m}}\ket{0}^{\otimes n}$. Since the procedure for measuring the diagonal elements is relatively straightforward, we will focus on the measurement of the non-diagonal elements. Figure 2(a) shows a quantum circuit to obtain the matrix element $M_{m}^{i_{m}^{\prime}i_{m}}$. The procedure for constructing $M_{m}$ is as follows. First, we prepare initial states. We use a Hadamard test to prepare $\ket{\varphi^{i_{m}^{\prime}}}=U_{Cm}^{i_{m}^{\prime}}\ket{0}^{\otimes n}$ and $\ket{\varphi^{i_{m}}}=U_{Cm}^{i_{m}}\ket{0}^{\otimes n}$. The ancilla qubit is initialized to $\frac{\ket{0}+e^{i\alpha}\ket{1}}{\sqrt{2}}$, where $\alpha$ is the phase. We set $\alpha=0$ ($\alpha=\frac{\pi}{2}$) to obtain the real (imaginary) part of $M_{m}^{i_{m}^{\prime}i_{m}}$. Since $M_{m}$ is a Hermitian matrix, only measurements of $\mathrm{Re}(M_{m}^{01})$, and $\mathrm{Im}(M_{m}^{01})$ are required. Then, we measure on a computational basis. We use the fact that $O_{mr}$ is a Hermitian operator and has a spectral decomposition: $O_{mr}=U_{mr}^{{\dagger}}D_{mr}U_{mr}$, where $U_{mr}$ is a unitary matrix and $D_{mr}$ is a diagonal matrix. Also, we assign elements of $D_{mr}$ to the measurement results. More concretely, denoting $D_{mr}=\mathrm{diag}[\lambda^{(mr)}_{j_{r}=0},\lambda^{(mr)}_{j_{r}=1}]$, the measured value is computed as $\prod_{r=1}^{n}\lambda^{(mr)}_{j_{r}}$ corresponding to the measured outcome $\vec{j}=(j_{1},j_{2},\dots,j_{n})$.

We explain how to construct $M_{Am}$. Figure 2(b) shows the circuit to construct $M_{Am}$. In the case of $\ket{\varphi^{i_{m}}}=U_{Cm}^{i_{m}}\ket{0}^{\otimes n}$, since the $\ket{\varphi^{i_{m}}}$ are non-orthogonal to each other, we have $A\neq 1$. Therefore, $M_{Am}$ must be calculated. The circuit in Fig. 2(b) is a modified version of that in Fig. 2(a), except that system measurements are not necessary; we can construct $M_{Am}$ by using the same construction procedure as for $M_{m}$.

Next, we assume the case of $\ket{\varphi^{i_{m}}}=U_{Cm}\ket{{i_{m}}}\ket{0}^{\otimes n-1}$. In this case, we can obtain all the elements of $M_{m}$ only from the results of direct measurements. Figure 2(c) shows a quantum circuit to construct $M_{m}$. The initial states are set to $\ket{a_{m}}\ket{0}^{\otimes n-1}$, where $\ket{a_{m}}$ takes four states as $\ket{0}$, $\ket{1}$, $\ket{+}$ and $\ket{+y}$. $M_{m}$ can be obtained using the corresponding measurement results. Refer to Appendix A for the details of the procedure. Also, since $\ket{\varphi^{i_{m}}}$ is orthogonal, $A=1$ and $M_{Am}$ does not have to be calculated.

Finally, we show the procedure to obtain $\bra{\psi_{HT}}O\ket{\psi_{HT}}$, which can be implemented by contracting $M_{m}$ and $M_{A_{m}}$. Now, denoting $\ket{\psi}=\sum_{\vec{i}}\psi_{\vec{i}}\ket{\vec{i}}=\sum_{i_{1},i_{2},\dots,i_{k}}\psi_{i_{1},i_{2},\dots,i_{k}}\ket{i_{1}i_{2}\dots i_{k}}$, we can rewrite Eq. (3) as

where we have used the fact that $M_{m}$ is a Hermitian operator and has a spectral decomposition, $M_{m}=U_{m}^{{\dagger}}D_{m}U_{m}$. Figure 2(d) shows a quantum circuit to measure $M_{m}$. Henceforth, we will assume $\ket{\psi}=U_{M}\ket{0}^{\otimes n}$, where $U_{M}$ is a unitary operator with polynomial depth. We can compute $A^{2}\braket{O}$ by applying $U_{m}$, measuring in a computational basis, and assigning elements of $D_{m}$ to the measurement results. More concretely, denoting $D_{m}=\mathrm{diag}[\lambda^{(m)}_{i_{m}=0},\lambda^{(m)}_{i_{m}=1}]$, the measured value is computed as $\prod_{m=1}^{k}\lambda^{(m)}_{i_{m}}$ corresponding to the measured outcome $\vec{i}=(i_{1},i_{2},\dots,i_{k})$. In a similar procedure, we can measure $A^{2}$ by contracting $M_{A_{m}}$; hence, we can obtain $\braket{O}$.

We describe the measurement of transition amplitudes with the hybrid tensor network approach. The difference from Sec. II is that we need to contract a non-Hermitian operator $N_{m}$.

We proposed a method to calculate transition amplitudes using a hybrid tensor network. When the hybrid tensor network approach is naively applied to the transition amplitude calculation, the contracted operators become non-Hermitian, and the number of terms to be measured increases exponentially. Our method obtains the expectation value without increasing the number of terms exponentially by using singular value decomposition and a Hadamard test. Our theory can be easily generalized to cases with a mixture of classical and quantum tensors called quantum-classical and classical-quantum tree-tensor networks, and those with deeper tree structures. Moreover, we can easily extend the scenario to the case where the measured operator $O$ is a tensor product of non-Hermitian operators by using the SVD contraction method. This study significantly expands the applicability of the hybrid tensor network.

Future work includes the application of our method to algorithms related to hybrid tensor networks. For example, Deep VQE, [20, 21], which is a large-scale computational algorithm for NISQ devices based on the divide and conquer method, can be treated in the framework of hybrid tensor networks in theory. By applying the proposed method to such algorithms, we can extend the range of applications to various large-scale quantum algorithms.

This work is supported by PRESTO, JST, Grants No. JPMJPR1916 and No. JPMJPR2114; ERATO, JST, Grant No. JPMJER1601; CREST, JST, Grant No. JPMJCR1771; Moonshot R&D, JST, Grant No. JPMJMS2061; MEXT Q-LEAP Grant No. JPMXS0120319794 and JPMXS0118068682. S.E. acknowledges useful discussions with Jinzhao Sun and Xiao Yuan.