Quantum Process Tomography of Unitary Maps from Time-Delayed Measurements

We now describe how to reconstruct a single-qubit Hamiltonian from a time series of measurements. The Bloch sphere picture Nielsen and Chuang (2010) provides a geometric perspective and insight into single-qubit quantum states and operations. The Bloch vector $\vec{r}\in\mathbb{R}^{3}$ associated with $\psi\in\mathbb{C}^{2}$ is defined via the relation

For pure states as considered here, $\vec{r}$ is a unit vector. We can parametrize any single-qubit Hamiltonian (up to multiples of the identity map, which are irrelevant here) as

with $\vec{h}\in\mathbb{R}^{3}$. The corresponding time evolution operator is the rotation

when representing $\vec{h}=\omega\vec{v}/2$ by a unit vector $\vec{v}\in\mathbb{R}^{3}$ and $\omega\in\mathbb{R}$. On the Bloch sphere, this operator corresponds to a classical rotation about $\vec{v}$, as illustrated in Fig. 3a and described by Rodrigues’ rotation formula ($\theta=\omega t$):

$\mathsf{U}_{\theta,\vec{v}}\in\mathrm{SO}(3)$ is a rotation matrix (parametrized by $\theta$ and $\vec{v}$) applied to $\vec{r}$. In the above context of Takens embedding with $U=U(\Delta t)$, we may equivalently work with $\mathsf{U}=\mathsf{U}_{\omega\Delta t,\vec{v}}$ and matrix powers thereof.

The time trajectory effected by the rotation applied to $\vec{r}$ results in a circle embedded within the Bloch sphere. Knowing the circle would allow to determine $\vec{v}$, and the dynamics on the circle to determine $\omega$. By construction we also know one point on the circle already, namely $\vec{r}$ (the initial condition), but we only have access to the expectation value of an observable $M$ for $t>0$. In the following, we denote the “measurement direction” by $\vec{m}\in\mathbb{R}^{3}$, i.e., $M$ is parametrized as $M=\vec{m}\cdot\vec{\sigma}$. Algorithm 1 facilitates a recovery of $\omega$ and $\vec{v}$ using the projections of the time trajectory on $\vec{m}$. The main idea is to first reconstruct $\omega$ based on the time dependence, and then the unit vector $\vec{v}$.

As caveat, the solution to this problem is not unique, due to the four possible signs of the coefficients $\alpha_{2}$ and $\alpha_{3}$ in the algorithm. Fig. 3a visualizes how two rotations can generate the same projection onto the $z$-axis, with $\vec{h}^{\prime}$ resulting from $(\alpha_{2},\alpha_{3})\to-(\alpha_{2},\alpha_{3})$. The two equations for $\alpha_{2}$ and $\alpha_{3}$ are shown in Fig. 3b. (In the example, $\alpha_{1}=-0.2051$, and hence the radius of the norm constraint circle is close to $1$.) In order to decide between these distinct solutions, one could perform one further measurement in a different basis, or may use a-priori knowledge about the Hamiltonian.

We remark that the nonlinear optimization in the first step of the algorithm might get trapped in a local minimum, which could be resolved by restarting the optimization. On the other hand, when neglecting uncertainties associated with the measurements, a correct solution is indicated by a zero residual both in steps 1 and 3. We assume that a range of realistic frequencies $\omega$ is known beforehand, such that the Nyquist condition (given the non-uniform sampling points $t_{q}$) holds Marvasti (2001).

As straightforward approach for determining a unitary time step matrix $U$ which matches the measurement averages, one could start from the natural parametrization $U=\operatorname{e}^{-iH}$ (with the time step already absorbed into $H$), and then optimize the Hermitian matrix $H$ directly. However, this method encounters the difficulty of a complicated optimization landscape with many local minima, in particular when using gradient descent-based approaches.

Here we discuss an alternative numerical approach tailored towards generic (dense) few-qubit Hamiltonians. As discussed in Sect. IV.1, for single qubits a transformation $\psi^{\prime}=U\psi$ by a unitary matrix $U\in\mathrm{SU}(2)$ can equivalently be described by a spatial rotation in three dimensions on the Bloch sphere: $\vec{r}^{\prime}=\mathsf{U}r$, with $\mathsf{U}\in\mathrm{SO}(3)$ given by Rodrigues’ formula (10). An equivalent mapping from $U$ to $\mathsf{U}$ is via

for $\alpha,\beta=1,2,3$, where we have used the relation (7) together with the orthogonality relation of the Pauli matrices: $\frac{1}{2}\operatorname{tr}[\sigma^{\alpha}\sigma^{\beta}]=\delta_{\alpha\beta}$. For our purposes, the Bloch picture has the advantage that measurement averages depend linearly on the Bloch vector, and hence on $\mathsf{U}$, e.g., $\braket{\psi^{\prime}|\sigma^{z}|\psi^{\prime}}=\vec{e}_{z}\cdot\vec{r}^{\prime}=\vec{e}_{z}\cdot(\mathsf{U}r)$. In practice, we first optimize the entries of $\mathsf{U}$ to reproduce the measurement data (under the constraint that $\mathsf{U}$ is an orthogonal matrix), and afterwards find $U$ related to $\mathsf{U}$ via Eq. (11).

Generalization to a larger number of qubits is feasible via tensor products of Pauli and identity matrices (in other words, Pauli strings). The analogue of (11) for $n$ qubits, with $U\in\mathrm{SU}(2^{n})$, is $\mathsf{U}\in\mathrm{SO}(4^{n}-1)$ with entries

where $\alpha$ and $\beta$ are now index tuples from the set $\{0,1,2,3\}^{\otimes n}\backslash\{(0,\dots,0)\}$, using the convention that $\sigma^{0}$ is the $2\times 2$ identity matrix. We exclude the tuple $(0,\dots,0)$ here since it conveys no additional information: the identity matrix is always mapped to itself by unitary conjugation.

Specifying an element of $\mathrm{SO}(4^{n}-1)$ involves more free parameters than for a unitary matrix from $\mathrm{SU}(2^{n})$ in case $n\geq 2$, thus the optimization might find an orthogonal $\mathsf{U}$ which reproduces the measurement data, but does not originate from a $U\in\mathrm{SU}(2^{n})$ via (12). We circumvent this representability issue as follows, focusing on the case of two qubits here. Every $U\in\mathrm{SU}(4)$ can be decomposed as Kraus and Cirac (2001); Zhang et al. (2003)

with the “entanglement” gate

and single-qubit unitaries $u^{\mathrm{a}},u^{\mathrm{b}},u^{\mathrm{c}},u^{\mathrm{d}}\in\mathrm{SU}(2)$. The single qubit gates can be handled as before, i.e., represented by orthogonal rotation matrices $\mathsf{u}^{\mathrm{a}},\mathsf{u}^{\mathrm{b}},\mathsf{u}^{\mathrm{c}},\mathsf{u}^{\mathrm{d}}\in\mathrm{SO}(3)$. We find the Bloch representation of $R$ via (12) (cf. Gamel (2016); Huang and Mendl (2021)); the parameters $\theta_{1},\theta_{2},\theta_{3}$ appear in the matrix entries solely as $\cos(\theta_{j})$ and $\sin(\theta_{j})$ for $j=1,2,3$. In summary, we express the decomposition (13) in terms of to-be found orthogonal matrices in the Bloch picture.

After switching to the described Bloch representation, we “relax” the condition that an involved (real) matrix $\mathsf{U}$ is orthogonal, by admitting any real matrix, but adding the term

to the overall cost function, where $\lVert{\cdot}\rVert_{\mathsf{F}}$ denotes the Frobenius norm. As advantage, we bypass a parametrization of $\mathsf{U}$ to enforce strict orthogonality and avoid local minima in the optimization.

By choosing $\gamma=2$ in Eq. (4), the time steps $t_{q}=\gamma^{q}$ are integers, and we can generate corresponding powers of $U$ by defining $\mathsf{U}_{0}=\mathsf{U}$, $\mathsf{U}_{q}=\mathsf{U}_{q-1}^{2}$ for $q=1,\dots,2d$, such that

In our case, the $\mathsf{U}_{q}$’s are separate matrices, which are set in relation via an additional cost function term

For the case of two-qubit Hamiltonians, we additionally substitute two-dimensional vectors $(c_{j},s_{j})$ for $(\cos(\theta_{j}),\sin(\theta_{j}))$, again to avoid local minima. The condition $c_{j}^{2}+s_{j}^{2}=1$ translates to another penalty term in the overall cost function:

An intriguing possibility of the time-delay measurements is the identification of the overall Hamiltonian based on measurements restricted to a subsystem. This scenario will actually require the preparation of more than one exactly known initial state. As minimal example, consider a two-qubit system, the choice between two initial states, being able to select a measurement basis via the gate $C$, and time-delayed measurements on one of the qubits while the other qubit is inaccessible, as illustrated in Fig. 4.

For the numerical experiment shown in Fig. 8 below, we use $X$-, $Y$- and $Z$-basis measurements on the top qubit and minimize the mean squared error between the observed measurement averages and the prediction based on a general Ansatz for the (traceless) Hamiltonian, i.e., $15$ real parameters.

Finally, we consider quantum systems defined on a lattice, and assume an Ising-type Hamiltonian as in Eq. (1) (instead of a generic Hamiltonian) to avoid an exponential growth of the number of parameters.

For the purpose of reconstructing the Hamiltonian parameters, the Schrödinger time evolution can be well approximated via a second order Strang splitting method McLachlan and Quispel (2002), i.e., the Time-Evolving Block Decimation (TEBD) algorithm Vidal (2004). Interpreting the resulting layout of single- and two-site unitary operators as forming a quantum circuit, the reconstruction task amounts to a variational circuit optimization to reproduce the reference measurement averages. Specifically, when considering the interaction and local field parts of the Hamiltonian

by themselves, the individual terms in each of the sums pairwise commute. Fig. 5 illustrates the corresponding quantum circuit for a one-dimensional lattice; the constructing works for higher-dimensional lattices as well.

With the algorithm described in Sec. IV.1 we are able to reconstruct a single-qubit Hamiltonian up to numerical precision. Fig. 6 shows the results obtained for the $\omega$ and least squares optimizations. Following the Takens framework, we use $2d+1$ time points, where $d=3$ is the number of Hamiltonian parameters to be reconstructed. Specifically, we have used the time points $t_{q}=0.3\times(1.3)^{q}$ for $q=0,\dots,6$ here. With this, we are able to find $\omega$, $\alpha_{1}$ and $\kappa$ up to an error of $\sim 10^{-15}$. Thus $7$ measurements, plus a further measurement in a different basis to distinguish between the possible solutions due to symmetry, are sufficient for the reconstruction.

In this subsection we consider a generic single- or two-qubit Hamiltonian $H$. Specifically, we draw the coefficients of $H$ with respect to the Pauli basis from the standard normal (Gaussian) distribution.

We first focus on a single-qubit $2\times 2$ Hamiltonian. Fig. 7a shows the result of directly fitting the vector $\vec{h}\in\mathbb{R}^{3}$ (see Eq. (8)), with the KL divergence between the predicted and actual measurement probabilities as cost function. This approach is compared to the “relaxation” procedure (described in Sect. IV.2) in Fig. 7b, using $\gamma=2$ and $\Delta t=0.1$. The manifold for the Takens embedding has dimension $d=3$ for a single qubit, and hence $2d+1=7$ time points should be used. However, we face the difficulty that the Hamiltonian is not uniquely specified solely by Pauli-$Z$ measurements according to the discussion in Sect. IV.1. For this reason, we include Pauli-$X$ measurement data as well, and reduce the number of time points to $4$ (to compensate for this additional source of data). For both versions we use the Adam optimizer Kingma and Ba (2015). The direct fitting method might get trapped in a local minimum and hence not be able to find the ground-truth vector $\vec{h}$, which leads to a large variation around the median. This issue is ameliorated by the relaxation method.

In Fig. 7c visualizes the loss function and relative errors when applying the relaxation procedure to reconstruct the unitary time step matrix and corresponding Hamiltonian of a two-qubit system ($d=15$ real parameters). Specifically, for parameter optimization we express Eq. (13) using the Bloch picture, i.e., in terms of orthogonal rotation matrices for the single-qubit unitaries, and likewise using the Bloch picture analogue of the entanglement gate (14), parametrized by two-dimensional vectors $(c_{j},s_{j})$ for $(\cos(\theta_{j}),\sin(\theta_{j}))$. The condition $c_{j}^{2}+s_{j}^{2}=1$ translates to another penalty term in the overall cost function, see Eq. (18). We set $\Delta t=0.05$, $\gamma=2$ and use $6$ time points for this experiment. The reference measurement data at each time point are the expectation values of the nine observables $\{I\otimes\sigma^{\alpha},\sigma^{\alpha}\otimes I,\sigma^{\alpha}\otimes\sigma^{\alpha}\}_{\alpha\in\{x,y,z\}}$. We use $54$ measurement data points in total (instead of $2d+1=31$) since we found that the additional information improves the reconstruction. As last step (after the main optimization), we find the Hamiltonian entries giving rise to the computed $\mathsf{U}$ via the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.

According to Fig. 7c, the median relative errors of $\mathsf{U}$ and the Hamiltonian are below $10^{-3}$, but there are still cases where the method cannot find the reference solution at all, even tough the loss value is small. An explanation could be that $\mathsf{U}$ and $H$ are not uniquely determined. We leave a clarification of this point for future work.

We study the scenario depicted in Fig. 4, namely performing measurements solely on one out of two qubits. The observables are the three Pauli gates here. For the reconstruction to work, we require that two different, precisely characterized initial states can be prepared, which we choose at random for the numerical experiments.

The ground-truth Hamiltonian is similarly constructed at random, by drawing standard normal-distributed coefficients of Pauli strings, which in sum form the Hamiltonian. We use the time step $\Delta t=\frac{1}{5}$ and $\gamma=1.15$ here, and have heuristically found $12$ time-delayed measurements (for each measurement basis) as viable to reliably reconstruct the Hamiltonian. The loss function is the mean squared error between the model prediction for the measurement averages and the ground-truth values. To avoid getting trapped in local minima during the numerical optimization, we use the best out of 10 attempts (in terms of the loss function) with different random starting points.

The specific observables are the three Pauli gates applied to the first qubit. Fig. 8a shows exemplary measurement time trajectories, and Fig. 8b the reconstruction error and loss function (median and $25\%$ quantiles based on $20$ random realizations of the overall setup). One observes that a reliable and precise reconstruction of the Hamiltonian (even up to numerical rounding errors) is possible in principle, when neglecting inaccuracies of the measurement process.

We remark that we have used only a single initial state and less time points for the “relaxation” method (see Fig. 7c), which can explain the seemingly larger errors there.

We first consider the case of a Hamiltonian (1) with uniform parameters (not depending on the lattice site), i.e.,

with $J\in\mathbb{R}$ and $\vec{h}\in\mathbb{R}^{3}$. Thus the task consists of reconstructing $4$ real numbers. The ground-truth values for the following experiment are $J=1$ and $\vec{h}=(0.5,-0.8,1.1)$. $\Lambda$ is a $3\times 4$ lattice with periodic boundary conditions, and the initial state is a single wavefunction with complex random entries (independently normally distributed). We compute the time-evolved quantum state $\psi(t)$ via the KrylovKit Julia package kry (2021), and use the squared entries of $\psi(t)$ as reference Born measurement averages.

For the purpose of reconstructing $J$ and $\vec{h}$, we approximate a time step via a variational quantum circuit as shown in Fig. 5, where the single- and two-qubit gates share their to-be optimized parameters $\tilde{J}$ and $\vec{\tilde{h}}$. For simplicity, we use three uniform time points $\Delta t$, $2\Delta t$, $3\Delta t$ with $\Delta t=0.2$ (instead of time points $t_{q}=\gamma^{q}\Delta t$), such that the quantum circuit for realizing a single time step can be reused.

We quantify the deviation between the exact Born measurement probabilities, $p(t)=\lvert{\psi(t)}\rvert^{2}$, and the probabilities $\tilde{p}(t)$ resulting from the trotterized circuit Ansatz, via the Kullback-Leibler divergence:

Fig. 9a visualizes the optimization progress, i.e., the relative errors of $J$ and $\vec{h}$, and the loss function (21). The darker curves show medians over $100$ trials of random initial states, and the lighter shades $25\%$ quantiles. The dashed horizontal line is the loss function evaluated at the ground truth values of $J$ and $\vec{h}$; it is non-zero due to the Strang splitting approximation of a time step. Interestingly, the optimization arrives at an even smaller loss value starting around training epoch $70$. We interpret this as an artefact of the Strang splitting approximation – note that around this epoch, the relative error of $\vec{h}$ slightly increases again. We have used the RMSProp optimizer Hinton (2014) with a learning rate of $0.005$ here. In summary, the reachable relative error is around $0.02$, and higher accuracy would likely require a smaller time step to reduce the Strang splitting error.

Next, we consider a Hamiltonian (1) with disorder (and external field solely in $x$-direction):

with i.i.d. random coefficients $J_{j,\ell}$ and $h_{j}$. For the numerical experiment, $J_{j,\ell}$ is uniformly distributed in the interval $[0.8,1.2]$, and $h_{j}\sim 0.5\,\mathcal{N}(0,1)$ (normal distribution).

The results are shown in Fig. 9b. As before, we evaluate $100$ realizations of the experiment, with a set of coefficients and random initial state drawn independently for each trial. We take the maximum over the relative errors of $J_{j,\ell}$ for all $j,\ell$ to arrive at the relative error reported in Fig. 9b, and likewise for $h_{j}$. The “reference” loss function is the KL divergence evaluated at the ground truth coefficients. It is non-zero due to Strang splitting errors, and fluctuates due to the different random coefficients at each trial.