Variational Quantum Process Tomography of Non-Unitaries

Quantum process tomography is a fundamental and critical benchmarking and certification tool that is capable of fully characterizing an unknown quantum process. Standard quantum process tomography suffers from an exponentially scaling number of measurements and complicated data post-processing due to the curse of dimensionality. On the other hand, non-unitary operators are more realistic cases. In this work, we put forward a variational quantum process tomography method based on the supervised quantum machine learning framework. It approximates the unknown non-unitary quantum process utilizing a relatively shallow depth parametric quantum circuit and fewer input states. Numerically, we verified our method by reconstructing the non-unitary quantum mappings up to eight qubits in two cases: the weighted sum of the randomly generated quantum circuits and the imaginary time evolution of the Heisenberg XXZ spin chain Hamiltonian. Results show that those quantum processes could be reconstructed with high fidelities (>99%) and shallow depth parametric quantum circuits (d≤8), while the number of input states required is at least two orders of magnitude less than the demands of the standard quantum process tomography. Our work shows the potential of the variational quantum process tomography method in characterizing non-unitary operators.


Introduction
Quantum process tomography is a fundamental and indispensable technique in quantum information processing [1]. Recently, it has been increasingly crucial in benchmarking and verifying the performance of a quantum device and its dynamics when the system goes larger.
However, standard quantum process tomography (SQPT) [2,3] is cursed by the exponentially exploding Hilbert space dimension used to represent the quantum process, which needs to prepare an informationally complete set of input states and perform the standard quantum state tomography on the corresponding output quantum states [3][4][5]. Specifically, it requires 4 n input states and 4 2n quantum measurements for an n-qubit quantum process. Such exponential overhead severely constrains the problem size on which SQPT can be feasibly conducted. Currently, SQPT has only been experimentally implemented up to three qubits [6][7][8][9][10][11]. Many alternative proposals are put forward to address these problems, including compressed sensing tomography [12,13] and ansatz-based tomography [14,15], but at the cost of assuming specific structures of the unknown quantum process. While other benchmarking protocols only partially characterize a quantum system, including randomized benchmarking [16,17], direct fidelity estimation [18], and so on. Therefore, quantum process tomography still plays an indispensable role since it is a comprehensive characterization of quantum systems with full information and no prior assumptions. On the other hand, realistic quantum processes are often non-unitary operators. So characterizing a non-unitary quantum process is a practical problem.
Recently, quantum machine learning combined with parametric quantum circuit (PQC) has been making considerable progress [19][20][21][22][23][24]. Reference [25] showed that a gradientbased quantum machine learning algorithm could efficiently extract the information of certain quantum states encoded in a PQC, after which the unknown quantum state can be reconstructed classically with high fidelity using the optimal parameters of the PQC. Further in Ref. [26], PQC with relatively shallow depth demonstrated the capability of characterizing unitary mappings with fewer input states.
In this work, we propose a variational quantum process tomography (VQPT) method of non-unitaries based on the supervised quantum machine learning framework, extending the work in Ref. [26]. As shown in Figure 1a, we utilize a d-depth PQC C( θ) to approximate the unknown quantum process denoted by O, where θ is the parameter list encoding the information to be optimized in this PQC. As shown in Ref. [26], when O is a unitary, we randomly prepare a training set of N quantum states |ψ j to train the PQC, each of which is separately fed into the unknown quantum process O and the PQC C( θ). Regarding non-unitary cases, we put forward two methods in this paper-(1) superposing unitary PQCs to approximate the non-unitary; (2) transforming the non-unitary problem into a unitary one-to tackle the problem. Then, similarly, as long as each pair of output quantum states, O |ψ j and C( θ) |ψ j are identical, and that N is large enough, the unknown quantum process should be approximated by C( θ) no matter O is unitary or not. As a result, all the information of O is stored in the parameters θ, and we can systematically reconstruct O from those parameters using a classical computer. We numerically verified our approach by reconstructing two typical non-unitary processes: the weighted sum of the randomly generated quantum circuits and the imaginary time evolution of the Heisenberg XXZ spin chain Hamiltonian, respectively. Numerical results show that we could reconstruct a non-unitary quantum process up to 8 qubits with an average gate fidelity higher than 99%, and the number of required input quantum states is smaller than that required by SQPT by at least 2 orders of magnitude.  Figure 1. (a) The general framework of our VQPT method. By training the PQC with the quantum states in the training set Ψ and validation set Φ and optimizing the parameters in C( θ) based on a gradient-descent approach, the PQC C( θ) gradually approximates the physical quantum process O.
(b) The structure of the PQC. It begins and ends with a single-qubit layer. Each two-qubit layer is counted as a depth, and the circuit contains d depths of operations and ends with a single-qubit layer.
Compared to other tomography approaches, our method has several advantages. First, combined with variational optimization, it avoids exponential measurements and complicated quantum state tomography subroutines in SQPT. It only requires measuring a single qubit for each configuration, hence is less prone to errors with relatively shallow circuit and simple measurement [27,28]. Second, with a small number of input states, our approach greatly reduces the exponential state preparation overheads and reconstructs O with high fidelities. And last, numerical simulations demonstrate the capability of our method on non-unitary quantum processes up to eight-qubit cases, further proving its generalities and potentials. This paper is organized as follows. In Seciton 2, we introduce the scheme of our VQPT for non-unitary processes. In Secitons 3 and 4, we give two cases for solving the non-unitary quantum process tomography. We verified our method with numerical simulations of non-unitary quantum process tomography for the weighted sum of the randomly generated quantum circuit in Seciton 3 and the imaginary time evolution of a Heisenberg XXZ spin chain Hamiltonian in Seciton 4. We conclude in Section 5.

Variational Algorithm for Approximating Non-Unitary Quantum Processes with Parametric Quantum Circuits
There are three main components in variational quantum algorithms: ansatzes, cost function, and optimization methods. Figure 2 gives an illustrative framework of variational quantum algorithms (VQA) [29]. It is a quantum-classical hybrid architecture. In this section, we describe the variational algorithm for non-unitary quantum processes from the above three aspects.

Ansatz of the Parametric Quantum Circuit
Ansatz is an essential aspect of a VQA. It encodes the parameters θ, and can be further trained to minimize the cost function. The specific structure of an ansatz varies from each other depending on the task. Typical ansatzes include hardware-efficient ansatz [30], variational Hamiltonian ansatz [31], and so on.
Here in our work, the design of our PQC is shown in Figure 1b, where interlaced layers of single-qubit gates and two-qubits CNOT gates are used. Such an organization facilitates the quick generation of entanglement between qubits, thus making it possible to approximate complicated quantum processes. Specifically, each two-qubit layer is counted as a depth, varying between odd and even depth. Each single-qubit layer contains three rotational gates (R z , R y , and R z ) on each qubit, where R y and R z are defined as The sequences R z , R y , and R z ensure that arbitrary single qubit rotations can be realized with appropriate parameters. Our PQC ends with a single-qubit layer. As a result, the total number of parameters for such a circuit with n qubits and d depths is 3n(d + 1).
It is noted that in practice, the design of the PQC should consider the underlying quantum hardware, especially the choice and pattern of the two-qubit gates. Hence the ansatz is not fixed, and our method is a generalized framework with various circuit ansatzes.

Cost Function on the Training Set
In our method, we build a cost function on the training set, which evaluates the distance between C( θ) and the target O. Concretely, we first randomly generate a training set of N quantum states, denoted as Ψ = {|ψ 1 , |ψ 2 , . . ., |ψ N }. Here, a random quantum state is generated by applying an R y gate with random parameters onto each qubit and controlled phase (CZ) gates with random control and target qubits. Each state |ψ j is fed into the unknown quantum process O and the PQC C( θ), with the output quantum states denoted as |ψ ideal j = O |ψ j and |ψ out j = C( θ) |ψ j . Then, we compute the Euclidean distance between |ψ ideal j and |ψ out j , which is The inner product on the right-hand side of Equation (2) can be efficiently computed with a quantum computer using a generalized SWAP-test algorithm [26]. Then, the cost function f ( θ) is defined as the summation of the distance obtained over all input states, which is namely, f ( θ) is the mean square error between the two sets of output quantum states.

Gradient-Based Optimization
The cost function f ( θ) is a hybrid quantum-classical function, and its gradient can be estimated based on the chain rule, where it contains functions to be evaluated with a quantum computer and functions to be evaluated on a classical computer: S in our case means the generalized SWAP-test circuit with parametric quantum gates, and ∂S ( θ)/∂ θ j can be computed using the parameter-shift rule [32], where θ j denotes the j-th parameter in the parameter list θ, and θ ± j = θ j ± π 2 . Hence, the gradient of the cost function f ( θ) can be computed following Ref. [33], which proposed a method to embed Equation (4) into a classical automatic differentiation framework, such that the gradient of the hybrid quantum-classical cost function can be automatically computed using a hybrid quantum-classical computer. The gradient can then be fed into a gradient-based optimizer to minimize the cost function f ( θ).

Evaluation Criteria
After the training, we evaluate experimental results between O and C( θ) using the average gate fidelity [34]. The average gate fidelity F avg between the rebuilt C( θ) and the actual O is given by In addition, we borrow the idea from supervised machine learning and design a validation set to test the generalization ability of the training outcomes. Since during training, we do not know whether our PQC is expressive enough to represent O, and whether the number of input states is enough to ensure convergence to O. Moreover, in practice, we may also have the problem of overfitting such that the optimal C( θ) is very distinct from O, but the cost function f ( θ) has already converged to 0. To overcome these problems, we generate another set of N input states denoted as Φ = {|φ 1 , |φ 2 , . . ., |φ N }, i.e., the validation set, which is independent of the training set. After the training process, we feed each |φ j into the unknown quantum process and the resulting optimal PQC, obtaining two outputs |φ ideal j = O |φ j and |φ out j = C( θ) |φ j . Then, we compute the quantum fidelity between |φ ideal j and |φ out j efficiently through the SWAP-test [35] on a quantum computer and summarize all the instances of the validation set, which is defined as the accuracy It is noted that accuracy is a faithful tool. If accuracy is close to one, the PQC C( θ) we obtained can be well generalized to the new input validation states. The values of accuracy and F avg are indeed strongly correlated in our numerical simulation. Therefore, we can pick out the simulation with a larger accuracy value as a more faithful reconstruction of O. Moreover, the accuracy is an efficient evaluation criterion. Since the direct characterization of the distance between O and C( θ) as required in Equation (5) scales exponentially with the number of qubits n, it is possible to determine whether the training is successful or not based on the accuracy, without resorting to the complete characterization of F avg .

Case Study I: Weighted Sum of the Randomly Generated Quantum Circuits
For an n-qubit randomly generated quantum circuit (RQC) [36], we first apply Hadamard gates to initialize the state to a symmetric superposition. Then, the circuit is organized by depth, including CZ gates alternating between odd and even configurations to entangle neighboring qubits and randomly chosen single-qubit gate (T, √ X or √ Y). Finally, Hadamard gates are applied to each qubit. It is noted that such randomly generated quantum circuits are hard for efficient simulation on a classical computer [36,37].
Here, we take the weighted sum of RQC as the target non-unitary operator O where the parameter p denotes the weighted probability that the operator is an identity matrix.

Methods
For the non-unitary operator in Equation (7), we correspondingly utilize the weighted sum of the unitary PQCs to approximate the target non-unitary quantum process where p i are the weighting parameters and C i ( θ i ) are the corresponding individual unitary PQCs as shown in Figure 3. Here in the known structure case, we assume that C ( θ, p) = pI + (1 − p)C( θ).  When dealing with unitary mappings, the cost function in Equation (3) merely involves the third term with variational parameters since |ψ ideal j and |ψ out j are both norm-1 state vectors. In terms of such non-unitary cases, these two output states cannot be kept at norm one, depending on the circuit operations C . Therefore, the cost function involves both circuit parameters θ i and weighting parameters p i : Hence, the gradient calculation and consequent optimization are related to the collaborative optimization of both parameters θ i and p i . In Section 3.2, we give concrete examples of the non-unitary quantum process and the specific form of cost function and its gradient.

• Fixed weighted sum of PQCs
A specific six-qubit randomly generated circuit is organized as shown in Figure 4a. Here, we first considered the weighted summation as a fixed and known structure to validate the capability of our VQPT method, where the weighting parameter is set at p = 0.1. Hence, we utilized the ensemble of identity process I and the PQC C( θ) to approximate the non-unitary operator Numerically, we conducted the scalability tests up to eight-qubit cases, as shown in Figure 4b (Data details are listed in Table A1 in Appendix A). Results show a perfect reconstruction of the weighted summation operator with F avg > 99.99%. Meanwhile, compared to the number of input states needed in SQPT (dashed line), our method has at least two orders of magnitude fewer demands on the input states (rectangle marks). In Figure 4c, we give the variational optimization details of (4, 5, 8), (5,6,12) and (6,7,35) configurations, where we denote (n, d, N) as an n-qubit, d-depth, and N-input PQC configuration in numerical simulation. The cost function f ( θ) and the fidelity F avg both converge to the optimums.  shows the maximum F avg reached during repetitive trials, and the logarithmic right-axis represents the corresponding number of input states N needed, while the red rectangular denotes the input state in our PQC, which is at least two orders of magnitude fewer than that in SQPT. It is noted that the data presented is valid for both fixed and variational weighting parameter conditions, since we utilized the same PQC configuration and achieved the same results. (c) VQPT on the weighted sum of the randomly generated quantum circuit on n-qubit, d-depth, and N-input configurations (denoted as (n, d, N)). The weighting probability p = 0.1 is a fixed and known parameter. (d) VQPT on the weighted sum of the randomly generated quantum circuit on (6, 7, 35) cases with unknown weighting parameters p. We initialized the parameter p init = 0.1, and we presented two trials with different target values p opt = 0.05 and p opt = 0.20 (red dotted lines). Results demonstrated convergence to the optimal values, and the maximum F avg is over 99%.

• Variational weighted sum of PQCs
Further, we studied the case where the weighting parameter p is unknown, As shown in Section 3.1, when conducting tomography on a non-unitary operator with an unknown weighting parameter, we could expand the specific cost function form with both circuit parameters θ and weighting parameter p. The first item ψ ideal j |ψ ideal j in Equation (9) is a fixed value, and the last two items involve these parameters and can be further expanded with Equation (11): Based on the expanded cost function items in Equations (12) and (13), we can calculate the corresponding gradient and conduct consequent optimization, where both θ and p are variational targets. We give two examples of (6, 7, 35) with different weighting parameters p opt = 0.05 and p opt = 0.20, respectively, where the initial value is set as p init = 0.1. Numerical results in Figure 4d show that our method is capable of optimizing the variational circuit parameters θ (solid lines) together with the weighting parameter p (triangle and rectangular marks) based on the flexible cost function f ( θ, p). Moreover, we conducted the scalability tests on unknown weighting parameter p cases, which achieved F avg > 99.99% on all conditions as listed in Table A1 in Appendix A.
Here, we give a two-qubit case with the target matrix and the reconstructed matrix shown in Figure 5.

Case Study II: Imaginary Time Evolution of Heisenberg XXZ Spin Chain
We take the Hamiltonian of the Heisenberg XXZ spin chain [38] in a magnetic field as our example, which is written as, Here n is the number of spins (qubits), J is the tunneling strength, ∆ is the interaction strength, and h is the magnetization strength. The imaginary time evolutionary operator with time τ [39] is denoted as In the simulations, we fixed h = 0.1 and J = 1.

Methods
Mathematically, we could transform the target non-unitary processes into unitary ones and afterward utilize the method in Ref. [26] to conduct consequent unitary learning, as shown in Figure 6. Here, we mainly introduce the following two methods.

• Unitary dilation
Every contraction operator on a small Hilbert space has a unitary dilation in an extended Hilbert space, which is guaranteed by the Sz.-Nagy dilation theorem [40].

Theorem 1.
A contraction operator T applied up to N times on a smaller space H can be equivalent to a unitary U T applied up to N times on a larger space K, with where P H is the projection operator into space H.
In a particular case, the Sz.-Nagy theorem has a minimal dilation when the extended space K has the smallest dimension. Here, we give an example of a minimal unitary dilation with N = 1: where D T = (I − T † T) 1 2 . Therefore, we can utilize Equation (17) to construct the extended unitary of the original non-unitary operator by adding one ancilla qubit.

• Unitary decomposition
Reference [41] shows that any quantum operator can be exactly decomposed as a linear combination of at most four unitary operators.
Specifically, any operator O can be decomposed into a Hermitian and an anti-Hermitian component: such that O = S + A. Further, the Hermitian and anti-Hermitian components can be expressed as the sum of two unitary operators u i : as long as the expansion parameter approaches zero. With this decomposition, we can write the action of O as a sum of at most four unitary operators u i , regardless of whether O itself is unitary or not.

• Unitary dilation
Utilizing the unitary dilation method in Equation (17), we can dilate the non-unitary operator into a unitary by adding one ancilla qubit with the minimum dilation. Then, the problem is a trivial unitary one. Firstly, we numerically tested the five-qubit XXZ spin chain evolution with imaginary time τ varies from 0.01 to 0.15. By adding one ancilla qubit, we utilized the (6,7,20) PQC configuration to conduct the tomography. Table 1 shows a satisfying reconstruction of the extended unitary process with high F avg and low f ( θ). Further, we fixed the evolution time τ = 0.01, and conducted numerical simulations from four-qubit to six-qubit XXZ spin chain imaginary time evolution (the corresponding required qubits n ext are extended to five to seven). It can be observed in Table 2 that with a shallow depth PQC and fewer input states, our method achieves satisfying results F avg > 99%. Similarly, in Section 3.2, the number of input states in our PQC N PQC is far fewer than that in SQPT N SQPT .  3 Number of input states in SQPT. 4 Number of parameters in C( θ).

• Unitary decomposition
As shown in Equation (19), we mathematically transformed the non-unitary quantum process tomography into four unitary quantum process tomography utilizing the unitary decomposition method and numerically simulated the cases up to four-qubit. Specifically, we prepared four unitary PQCs C i ( θ i ), i = 1, 2, 3, 4 and superpose the four PQCs as shown in Equation (19) to approximate the non-unitary O XXZ . Therefore, the original non-unitary tomography is transformed into a unitary tomography with fourfold variational parameters. Taking two-qubit imaginary time evolution as an example, we numerically reconstructed the non-unitary quantum process O as the following four unitaries C i ( θ i ) Here, we choose the = 0.05, and by Equations (18) and (19), we can obtain the superposition C( θ) = 10C 1 ( θ 1 ) − 10C 2 ( θ 2 ) + 10iC 3 ( θ 3 ) − 10iC 4 ( θ 4 ) with average gate fidelity 99.90%. More details are presented in Table 2. It is worth noting from Tables 2 and A1 that the accuracy value is a valid leading indicator to pick out the better PQC configurations with higher fidelity F avg values. We give detailed data in Appendix B.

Conclusions
In this work, we propose a non-unitary variational quantum process tomography method based on the quantum machine learning algorithm, which encodes the unknown non-unitary quantum process into a PQC of a certain depth d. A set of randomly generated quantum states are used as the training data to minimize the cost function, and a validation set filters out the instance with the highest F avg as a leading indicator, namely the closest parameter configuration to the unknown quantum process.
We introduce two methods to tackle the non-unitary quantum process tomography. One focuses on the mathematical transformation of the non-unitary problem, while the other attempts to utilize superpositions of PQCs to approximate the non-unitary.
We demonstrated our method by two numerical examples, including the superposition of the random quantum circuits and the imaginary time evolution with the Heisenberg XXZ spin chain from two-qubit to eight-qubit. The results indicate that a faithful reconstruction of O (F avg higher than 99%) can be reached with a relatively low-depth PQC (d ≤ 8), and a relatively small number of training states (at least two orders of magnitude compared to SQPT). Moreover, only a single qubit measurement is required in each configuration instead of measuring in the complete set of computational basis. Our work further proves the potential of the variational quantum process tomography framework on non-unitaries and presents a promising application of using the quantum machine learning algorithm to accelerate the quantum process tomography.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: VQPT Variational quantum process tomography SQPT Standard quantum process tomography PQC Parametric quantum circuit VQA Variational quantum algorithm RQC Randomly-generated quantum circuit

Appendix A
Here, we give the data details of the scalability tests corresponding to the Figure 4a content.

Appendix B
Here, we give detailed data on the analysis of the two evaluation criteria: the average gate fidelity F avg and the accuracy.
As stated in the main text, these two evaluation criteria are strongly correlated. Here, we conduct 100 independent trials and computed the two values. We utilize the Pearson's correlation coefficient r [42] to evaluate such relations in our numerical simulations. The coefficient r ranges from −1 to 1, where r ≥ 0.7 means highly linear positive correlated. It is shown in Figure A1 that the accuracy value on the validation set has a strong correspondence with the final F avg value with the Pearson correlation coefficient r = 0.9588 and r = 0.9286 on two typical non-unitary trials, respectively. Therefore, it is feasible to utilize the accuracy value as a criterion to determine the optimal circuit parameters to reconstruct O with higher F avg .