User-specified random sampling of quantum channels and its applications

TL;DR

This paper introduces a novel method for sampling quantum channels using Hamiltonian Monte Carlo, leveraging a new parameterization of the channel space, with applications demonstrated in three scenarios.

Contribution

It develops an exact parameterization of quantum channels enabling efficient Hamiltonian Monte Carlo sampling, adapting classical state sampling techniques to quantum channels.

Findings

Successful implementation of the sampling method

Application to three quantum channel sampling tasks

High-quality samples generated from user-specified distributions

Abstract

Random samples of quantum channels have many applications in quantum information processing tasks. Due to the Choi--Jamio\l{}kowski isomorphism, there is a well-known correspondence between channels and states, and one can imagine adapting \emph{state} sampling methods to sample quantum channels. Here, we discuss such an adaptation, using the Hamiltonian Monte Carlo method, a well-known classical method capable of producing high quality samples from arbitrary, user-specified distributions. Its implementation requires an exact parameterization of the space of quantum channels, with no superfluous parameters and no constraints. We construct such a parameterization, and demonstrate its use in three common channel sampling applications.

Figures

Tables

Table 1. Table 1: Simulated data for the first example in Sec.  IV.1 . Each entry n k ( i ) superscript subscript 𝑛 𝑘 𝑖 n_{k}^{(i)} represents the number of clicks in detector Π k subscript Π 𝑘 \Pi_{k} when the input state ρ ( i ) superscript 𝜌 𝑖 \rho^{(i)} is sent.

${\rho }^{\left(i\right)}$

${\Pi }_k$ 1 2 3 4 1 7 5 6 6 2 3 13 4 4 3 2 7 8 7 4 3 7 8 6

Table 2. Table 2: Simulated data for the second example in Sec.  IV.1 . Each entry n k ( i ) superscript subscript 𝑛 𝑘 𝑖 n_{k}^{(i)} represents the number of clicks in detector Π k subscript Π 𝑘 \Pi_{k} when the input state ρ ( i ) superscript 𝜌 𝑖 \rho^{(i)} is sent.

${\rho }^{\left(i\right)}$

${\Pi }_k$ 1 2 3 4 5 6 7 8 9 1 8 1 4 1 8 2 1 2 0 2 3 10 3 1 1 2 3 3 1 3 1 1 9 1 2 2 3 5 3 4 8 3 3 4 0 1 4 3 1 5 2 3 3 1 10 2 4 2 0 6 3 4 4 2 0 8 2 2 2 7 3 2 4 2 0 2 9 4 1 8 2 0 4 0 0 2 6 8 5 9 3 3 3 1 5 0 1 3 8

Table 3. Table 3: Simulated data for the examples in Sec.  IV.2 . Each entry n k ( i ) superscript subscript 𝑛 𝑘 𝑖 n_{k}^{(i)} represents the number of clicks in detector Π k subscript Π 𝑘 \Pi_{k} when the input state ρ ( i ) superscript 𝜌 𝑖 \rho^{(i)} is sent.

${\rho }^{\left(i\right)}$

${\Pi }_k$ 1 2 3 4 1 9 4 4 7 2 6 6 3 9 3 3 5 10 6 4 8 4 5 7

Table 4. Table 4: Comparison of results based on AIC, BIC, and RBR criteria, with different number of measured copies N 𝑁 N . Candidate models are nested channel families (see main text). Each row below collects the counts for each family of true channels; each column collects the counts for the model that a criterion selects as the best fit for the data.

		true	# cases where the best-fit model is
		family	Dephasing	Pauli	SUnital	Unital	General
		Dephasing	947	43	10	0	0
		Pauli	583	408	9	0	0
		SUnital	629	319	52	0	0
		Unital	562	405	31	2	0
		General	596	372	30	2	0
		Dephasing	983	17	0	0	0
		Pauli	721	279	0	0	0
		SUnital	795	200	5	0	0
		Unital	741	256	3	0	0
		General	762	235	3	0	0
		Dephasing	712	62	93	56	77
		Pauli	224	271	173	140	192
		SUnital	239	151	315	117	178
		Unital	191	164	215	234	196
$N=20$	RBR     BIC   AIC	General	166	156	180	164	334
		Dephasing	935	42	21	2	0
		Pauli	305	644	45	6	0
		SUnital	343	429	213	13	2
		Unital	257	534	148	57	4
		General	278	531	115	33	43
		Dephasing	995	4	1	0	0
		Pauli	539	461	0	0	0
		SUnital	646	328	26	0	0
		Unital	601	393	6	0	0
		General	586	404	10	0	0
		Dephasing	823	52	68	22	35
		Pauli	158	404	139	129	170
		SUnital	155	197	376	140	132
		Unital	91	185	222	343	159
$N=50$	RBR     BIC   AIC	General	75	178	139	176	432
		Dephasing	938	41	18	1	2
		Pauli	173	733	72	16	6
		SUnital	141	368	455	30	6
		Unital	87	409	264	215	25
		General	78	442	147	93	240
		Dephasing	998	2	0	0	0
		Pauli	367	631	2	0	0
		SUnital	427	471	102	0	0
		Unital	357	593	42	8	0
		General	363	609	26	2	0
		Dephasing	905	43	33	11	8
		Pauli	129	578	105	94	94
		SUnital	76	210	506	126	82
		Unital	55	196	223	394	132
$N=100$	RBR     BIC   AIC	General	34	132	118	190	526

Equations

$$\ast (|\psi ⟩)\equiv |\overline{\psi }⟩\in \mathcal{H},\text{for }|\psi ⟩\in \mathcal{H},$$

$$⟨\overline{\psi }|\overline{\varphi }⟩=⟨\varphi |\psi ⟩,\forall |\psi ⟩,|\varphi ⟩\in \mathcal{H},$$

$$*\Bigl{(}\sum_{i}c_{i}|\psi_{i}\rangle\Bigr{)}=\sum_{i}c_{i}^{*}|\overline{\psi}_{i}\rangle,$$

$$*\Bigl{(}\sum_{ij}c_{ij}|\psi_{i}\rangle\langle\phi_{j}|\Bigr{)}\equiv\sum_{ij}c_{ij}^{*}|\overline{\psi}_{i}\rangle\langle\overline{\phi}_{j}|.$$

$$\text{vec}(|\psi ⟩⟨\varphi |)\equiv \ast (|\varphi ⟩)\otimes |\psi ⟩=|\overline{\varphi }⟩\otimes |\psi ⟩=|\overline{\varphi }\psi ⟩,$$

$$\text{vec}(ABC)=(C^{\mathcal{T}}\otimes A)\text{vec}(B).$$

$${\rho }_{\mathcal{E}}$$

$$=(\text{mathbbm}1\otimes \mathcal{E})\left(|\text{mathbbm}1⟩⟩⟨⟨\text{mathbbm}1|\right),$$

$$\mathcal{E}(X)={\text{tr}}_1\{{\rho }_{\mathcal{E}}(X^{\mathcal{T}}\otimes \text{mathbbm}1)\}.$$

$${\text{tr}}_2({\rho }_{\mathcal{E}})=\text{mathbbm}1,$$

$$p_{k}^{(i)}=\textrm{tr}\{\Pi_{k}^{(i)}\mathcal{E}(\rho^{(i)})\}=\textrm{tr}\Bigl{\{}\rho_{\mathcal{E}}\Lambda_{k}^{(i)}\Bigr{\}},$$

$$L(D|{\rho }_{\mathcal{E}})=\prod _{i}L(D^{(i)}|{\rho }_{\mathcal{E}})=\prod _i\left[\prod _k(p_k^{(i)}{)}^{n_k^{(i)}}\right],$$

$$A=\left(\begin{array}{cccc}|& |& & |\\ {\phi }_{11}& {\phi }_{12}& \dots & {\phi }_{dd}\\ |& |& & |\end{array}\right),$$

$${\phi }_i\equiv \left(\begin{array}{c}{\phi }_{i1}\\ {\phi }_{i2}\\ ⋮\\ {\phi }_{id}\end{array}\right),\text{for }i=1,2,\dots ,d,$$

$${\phi }_i^{†}{\phi }_j={\delta }_{ij}\text{for }i,j=1,2,\dots ,d.$$

$$\Phi \equiv \left(\begin{array}{cccc}|& |& & |\\ {\phi }_1& {\phi }_2& \dots & {\phi }_d\\ |& |& & |\end{array}\right).$$

$$\Psi \equiv P\Phi =\left(\begin{array}{cccc}|& |& & |\\ {\psi }_1& {\psi }_2& \dots & {\psi }_d\\ |& |& & |\end{array}\right)$$

$$\displaystyle\psi_{dk}\equiv{\left\{\begin{array}[]{ll}\textrm{e}^{\textrm{i}\phi_{k}}(\cos\theta_{k-1})S_{k}&\textrm{for }k=1,\ldots,K_{d}\\ 0&\textrm{for }k=K_{d}+1,\ldots,d^{3}\end{array}\right.},$$

$$\sum _{k=1}^m|{\psi }_{dk}{|}^2=S_m^2,\text{for any }m=1,2,\dots ,K_d,$$

$$v_{nk}$$

$$\displaystyle=\frac{\mbox{\footnotesize$1$}}{\mbox{\footnotesize$S_{n+1}$}}{\left\{\begin{array}[]{ll}\psi_{dk}\frac{\mbox{\footnotesize$\cos\theta_{n}$}}{\mbox{\footnotesize$\sin\theta_{n}$}}&\textrm{for }k=1,\ldots,n\\[4.30554pt] \psi_{d(n+1)}\frac{\mbox{\footnotesize$-\sin\theta_{n}$}}{\mbox{\footnotesize$\cos\theta_{n}$}}&\textrm{for }k=n+1\\[4.30554pt] 0&\textrm{for }k=n+2,\ldots,d^{3}\end{array}\right.}.$$

$$S_{n+1}{v}_n^{†}{\psi }_d$$

$$=\cos {\theta }_n\sin {\theta }_n{S}_{n+1}^2-\sin {\theta }_n\cos {\theta }_n{S}_{n+1}^2=0.$$

$$\left(\begin{array}{cccc}|& |& & |\\ {\psi }_1& {\psi }_2& ...& {\psi }_{d-1}\\ |& |& & |\end{array}\right)\equiv V\tilde{\Psi },$$

$${\text{tr}}_1({\rho }_{\mathcal{E}})=\text{mathbbm}1.$$

$$\rho =\frac{1}{2}(\text{mathbbm}1+\mathbf{\sigma }\cdot \mathbf{s}+\mathbf{t}\cdot \mathbf{\tau }+\mathbf{\sigma }\cdot \mathbf{C}\cdot \mathbf{\tau }),$$

$${\rho }_{\mathcal{E}}=\frac{1}{2}(\text{mathbbm}1+\mathbf{\sigma }\cdot \mathbf{C}\cdot \mathbf{\tau }).$$

$${\mathbf{v}}_1=(-1,-1,-1),$$

$${\mathbf{v}}_2=(-1,1,1),$$

$${\mathbf{v}}_3=(1,-1,1),$$

$${\mathbf{v}}_4=(1,1,-1),$$

$$(c_1,c_2,c_3)={\alpha }_1{\mathbf{v}}_1+{\alpha }_2{\mathbf{v}}_2+{\alpha }_3{\mathbf{v}}_3+{\alpha }_4{\mathbf{v}}_4,$$

$${\alpha }_1={\cos }^2{\theta }_1$$

Full text

Keywords:

quantum tomography, quantum channels, quantum parameter estimation, Choi-Jamiołkowski isomorphism, error regions, plausible region, random sampling, Monte Carlo methods

I Introduction

Quantum channels, or completely positive (CP) and trace-preserving (TP) maps, are a central concept in describing the dynamics of quantum systems. They form the basic models for imperfect quantum operations used for quantum information processing (QIP). Random—according to some specified distribution—samples of quantum channels are needed in many QIP tasks, including the evaluation of the distributional average of channel-related quantities, the computation of error bars for quantum process tomography, the exploration of typical properties of quantum channels, the numerical optimization of functions of channels over a complicated landscape, and others.

Sampling from specific distributions over the quantum state space is a well-studied problem, with many different approaches, including the Monte Carlo (MC) technique for arbitrary distributions, and other methods for sampling from specific distributions Shang+1:15 ; Seah+1:15 ; Zyczkowski98 ; Zyczkowski99 ; Zyczkowski01 ; Blume10 ; Huszar12 ; YS19 . Due to the Choi–Jamiołkowski isomorphism Choi ; Jamiolkowski , which gives a correspondence between CP channels and states, these state sampling methods can be adapted to sample quantum channels. Indeed, in the recent work by Thinh et al. Thinh18 , a Metropolis–Hasting Markov chain (MHMC) MC approach was used to sample channels from arbitrary distributions, by sampling the purification of the Choi–Jamiołkowski state corresponding to the channel. References Geometry and Bruzda09 discuss a procedure for generating samples of quantum channels with a specific distribution on the channel space by making use of the channel-state correspondence. An alternative way of generating the same distribution of channels is to couple the input state and an ancilla initially in an arbitrary pure state by a Haar-random unitary operator and then taking the partial trace over the ancilla. See Ref. Bruzda10 for more discussion on the properties of the distribution of channels generated by these two procedures.

As a general method for sampling from arbitrary distributions, MC methods stand out in their wide-ranging applicability and efficiency. The MHMC variety of MC methods used, for example, in Ref. Thinh18 , however, suffer from strong correlations between sample points, and one requires large samples for reliable answers not biased by these correlations. This was observed, for instance, in the MHMC state sampling algorithm of Ref. Shang+1:15 . A significant improvement in the quality of the samples was seen when we switched to the Hamiltonian Monte Carlo (HMC) approach Seah+1:15 , reaffirming the advantage of HMC over MHMC MC also observed in other settings Neal96 ; Hajian07 ; Porter14 ; Neal11 ; Duane87 .

The HMC method requires the availability of a parameterization of the domain space with exactly the right number of parameters, with no superfluous parameters and no constraints. The parameterization of the channel/state space used in Ref. Thinh18 , which has superfluous parameters, cannot be used for HMC. The exact parameterization of states used in the HMC algorithm in Ref. Seah+1:15 gives, through the Choi–Jamiołkowski isomorphism, a parameterization of the set of all CP, but not necessarily TP, maps. The TP property has to be imposed as an explicit constraint, thus rendering the parameterization unsuitable in a HMC algorithm for sampling CPTP channels.

In this work, we construct one exact parameterization of the space of CPTP maps, with no superfluous parameters, and no constraints. This can then be used in a HMC procedure for sampling from arbitrary, user-specified, distributions over the channel space. To illustrate the usefulness of our parameterization and the HMC algorithm, we apply our methods to three quantum sampling problems. Our examples are focused on problems in quantum process tomography, reflecting the interests of the authors; our parameterization and the HMC method, however, are just as useful for sampling problems in other areas of QIP. As an aside, our construction exactly parameterizes the space of all bipartite mixed quantum states with the completely mixed state for one of the parties.

Here is the brief outline of our paper. We first review the Choi–Jamiołkowski isomorphism in Sec. II. Section III explains our main contribution: the exact parameterization of the space of CPTP channels. In Sec. IV, we illustrate the use of our parameterization in a HMC sampling algorithm through three examples from quantum process tomography: (a) the construction of error regions in process estimation; (b) marginal likelihood for estimating specific properties of the channel; (c) model selection among candidate channel families. The reader is referred to Ref. Seah+1:15 or Appendix A for an introduction to the HMC algorithm used here. We conclude in Sec. V.

II The channel-state duality

There are many ways of writing the CPTP map of a quantum channel. Given our desire to make the connection with the sampling of quantum states, we make use of the channel-state duality and describe the quantum channel by a state via the Choi–Jamiołkowski isomorphism. Here, we remind the reader of this isomorphism, and, in the process, define the notation used throughout the article.

We begin with the $d$-dimensional Hilbert space $\mathcal{H}$ describing the state vectors (pure states) of the system. We define a map $*:\mathcal{H}\rightarrow\mathcal{H}$,

[TABLE]

such that

[TABLE]

and $*$ is “$*$-linear”, i.e.,

[TABLE]

where $c_{i}^{*}$ is the complex conjugate of $c_{i}$. Note that Eq. (2) specifies the $*$ map only up to a unitary transformation of no consequence. One specific realisation of the $*$ map, and what we use in our numerical examples below, is to first pick a basis $\{|i\rangle\}$ on $\mathcal{H}$, define $|\overline{i}\rangle\equiv|i\rangle$, and then extend the action of $*$ to arbitrary vectors using the $*$-linearity property. See also Sec. 3.1 in Ref. MUB10 for qubit examples of the $*$ map.

We extend the action of the $*$ map to adjoint vectors, $*(\langle\psi|)=\langle\overline{\psi}|=(|\overline{\psi}\rangle)^{\dagger}=[*(|\psi\rangle)]^{\dagger}$, and further to the set of operators on $\mathcal{H}$, denoted as $\mathcal{B}(\mathcal{H})$,

[TABLE]

We write $*(X)\equiv\overline{X}$, for any $X\in\mathcal{B}(\mathcal{H})$. Note that $\overline{X^{\dagger}}=(\overline{X})^{\dagger}$, and we denote $X^{\mathcal{T}}\equiv\overline{X}^{\dagger}$, a basis-independent transpose operation. If $X$ is non-negative, then so is $X^{\mathcal{T}}$.

Using the $*$ map, we define the vectorization map, a linear map from operators to vectors in a vector space $\mathcal{V}$, $\textrm{vec}:\mathcal{B}(\mathcal{H})\rightarrow\mathcal{V}$,

[TABLE]

for any $|\psi\rangle,|\phi\rangle\in\mathcal{H}$ and extended to all operators by linearity. We write, for any $X\in\mathcal{B}(H)$, $\textrm{vec}(X)\equiv|X\rangle\!\rangle\in\mathcal{V}$. Note the useful identity,

[TABLE]

Also, if $\{|i\rangle\}$ is an orthonormal basis for $\mathcal{H}$, then so is $\{|\overline{i}\rangle\}$. Consequently, the vectorized identity operator, $|\mathbbm{1}\rangle\!\rangle=\textrm{vec}(\mathbbm{1})=\sum_{i=1}^{d}|\overline{i}i\rangle$, can be regarded as a bipartite maximally entangled (unnormalized) state on $\mathcal{H}\otimes\mathcal{H}$.

Now, we are ready to state the channel-state duality. Consider a CP map, $\mathcal{E}:\mathcal{B}(\mathcal{H})\rightarrow\mathcal{B}(\mathcal{H})$, acting as $\mathcal{E}(\cdot)=\sum_{a}E_{a}(\cdot)E_{a}^{\dagger}$ for a (nonunique) set of Kraus operators $\{E_{a}\}$. We define

[TABLE]

where we have used the identity in Eq. (6); the $\mathbbm{1}$ in $\mathbbm{1}\otimes\mathcal{E}$ denotes the identity map. Thus defined, $\rho_{\mathcal{E}}$ is a nonnegative operator on $\mathcal{V}$; it can also be regarded as an unnormalized state (density operator) on the bipartite Hilbert space $\mathcal{H}\otimes\mathcal{H}\equiv\mathcal{H}_{1}\otimes\mathcal{H}_{2}$, labelling the two subsystems by $1$ and $2$. In the latter picture, one regards $|\mathbbm{1}\rangle\!\rangle\langle\!\langle\mathbbm{1}|$ as the density operator for a maximally entangled state on $\mathcal{H}\otimes\mathcal{H}$, and $\rho_{\mathcal{E}}$ is the density operator that results from the action of the map $\mathbbm{1}\otimes\mathcal{E}$ on it.

That $\rho_{\mathcal{E}}$ is invariant under a change of Kraus representation for the $\mathcal{E}$ is manifest in the last line of Eq. (II). We can turn the logic around: Any bipartite state on $\mathcal{H}\otimes\mathcal{H}$ possesses a spectral decomposition into eigenvectors, and the identification of those eigenvectors, with their corresponding (square root of the) eigenvalues, as vectorized Kraus operators immediately gives an associated CP map on $\mathcal{B}(\mathcal{H})$. Equation (II) hence states a duality between CP maps $\mathcal{E}$ and states $\rho_{\mathcal{E}}\geq 0$. $\rho_{\mathcal{E}}$ is sometimes called the “Choi state” of the CP map $\mathcal{E}$. Observe that

[TABLE]

We are primarily interested in CP maps that are also TP. In this case, the state $\rho_{\mathcal{E}}$ dual to the CP and TP channel satisfies the partial trace condition,

[TABLE]

i.e., $\mathcal{E}$ is CPTP if and only if $\rho_{\mathcal{E}}\geq 0$ and $\textrm{tr}_{2}(\rho_{\mathcal{E}})=\mathbbm{1}$. A simple count verifies that we have just the right number of parameters: A CP $\mathcal{E}$ is represented by $d^{4}$ real parameters—a positivity-preserving map that specifies how a $d^{2}$-element basis of operators on $\mathcal{H}$ is mapped back to itself—and this is the same number of real parameters needed to specify an unnormalized non-negative $\rho_{\mathcal{E}}$; the TP condition removes $d^{2}$ parameters, leaving $d^{2}(d^{2}-1)$ real parameters for a CPTP map, i.e., a quantum channel. Note that the set of $\rho_{\mathcal{E}}$s corresponding to quantum channels form a convex set of states, each with trace $d$. We denote the convex set of all $\rho_{\mathcal{E}}$ that satisfy Eq. (9) by $\mathscr{S}_{\mathrm{TP}}$, and refer to $\rho_{\mathcal{E}}\in\mathscr{S}_{\mathrm{TP}}$ as a TP state.

This duality between quantum channels and states enables us to sample quantum channels with algorithms for sampling quantum states (see the next section). Furthermore, the problem of process tomography—the estimation of the full description of a quantum channel acting on a quantum system—can be re-cast as that of state tomography. As the applications of our channel sampling algorithm discussed below are related to estimating quantum channels, we use the remainder of this section to recall this connection between state and process tomography, stemming from the channel-state duality MLEreview .

Quantum process tomography seeks to discover the full description of some unknown quantum channel $\mathcal{E}$, through $N$ uses of the channel. Standard strategies involve choosing a set of input states $\{\rho^{(i)}\}$, sending $N^{(i)}$ copies of state $\rho^{(i)}$ through the channel $\mathcal{E}$, and then measuring the output state using a POVM $\Pi^{(i)}\equiv\{\Pi_{k}^{(i)}\}$. For each $i$, the tomographic outcome probabilities come from the Born rule,

[TABLE]

where $\Lambda_{k}^{(i)}\equiv(\rho^{(i)})^{\mathcal{T}}\otimes\Pi_{k}^{(i)}$. Written in this manner, the expression for $p_{k}^{(i)}$ reminds one of the situation of state tomography of $\rho_{\mathcal{E}}$, where the set $\{\Lambda_{k}^{(i)}\}$ forms a pseudo-POVM in that $\Lambda_{k}^{(i)}\geq 0\ \forall k,i$, and $\sum_{k}\Lambda_{k}^{(i)}=\rho^{(i)}\otimes\mathbbm{1}$ for any $i$. Note that $\sum_{k}p_{k}^{(i)}=1$, as guaranteed by the TP condition in Eq. (9) together with the normalization $\textrm{tr}(\rho^{(i)})=1$.

The likelihood function for the data $D=\{D^{(i)}=(n_{1}^{(i)},n_{2}^{(i)},\ldots)\}$—$n_{k}^{(i)}$ denotes the number of clicks in detector $\Pi_{k}^{(i)}$ when $\rho^{(i)}$ is sent, and ${\sum_{k}n_{k}^{(i)}=N^{(i)}}$—collected is

[TABLE]

where we omit the combinatorial factors that are needed for proper normalization but are not important here. Disregarding quantum constraints, the likelihood is maximized, over all $\{p_{k}^{(i)}\}$, by setting {p_{k}^{(i)}=\frac{\mbox{\footnotesizen_{k}^{(i)}}}{\mbox{\footnotesizeN^{(i)}}}}; with quantum constraints, a constrained maximization of $L(D|\rho_{\mathcal{E}})$ over all permissible probabilities—those $p^{(i)}_{k}$s that could have come from a nonnegative $\rho_{\mathcal{E}}$ and which satisfy ${\sum_{k}p_{k}^{(i)}=1\ \forall i}$—yields what is known as the maximum-likelihood estimator (MLE) for $\rho_{\mathcal{E}}$ MLEreview .

III Parameterizing channels

III.1 Arbitrary channels

To obtain a sample of quantum channels according to some specified distribution, we generate Choi states $\rho_{\mathcal{E}}$ with the HMC algorithm. The HMC method demands a parameterization of the state space (in this case the space of $\rho_{\mathcal{E}}$) with no superfluous parameters and no external constraints. In Ref. Seah+1:15 , the ability to sample quantum states with the HMC algorithm was demonstrated using a parameterization of the full quantum state space. Because of the TP condition, sampling of quantum channels demands a parameterization of, not the full quantum state space as in Ref. Seah+1:15 , but only of the set $\mathscr{S}_{\mathrm{TP}}$ of TP states. Here, as our central result, we explain how to accomplish this.

We first choose a product basis $\{|\overline{i}j\rangle\}_{i,j=1}^{d}$ on $\mathcal{H}\otimes\mathcal{H}$ and represent $\rho_{\mathcal{E}}$ as a $d^{2}\times d^{2}$ matrix—also denoted as $\rho_{\mathcal{E}}$, to simplify notation—with complex entries. Positivity of $\rho_{\mathcal{E}}$ means that we can write $\rho_{\mathcal{E}}=A^{\dagger}A$, where $A$ is a $d^{2}\times d^{2}$ upper triangular complex matrix with real entries in the last column. The $d^{2}$ columns of $A$ are labelled using a double index,

[TABLE]

so that $\rho_{\mathcal{E}}=\sum_{ijkl=1}^{d}\varphi_{ij}^{\dagger}\varphi_{kl}|\overline{i}j\rangle\langle\overline{k}l|$, as the abstract, basis-independent object. Stacking the columns of $A$ to form columns with $d^{3}$ entries,

[TABLE]

permits writing the TP condition in Eq. (9), that is ${\textrm{tr}_{2}(\rho_{\mathcal{E}})=\sum_{ij}{\left(\sum_{k}\varphi_{ik}^{\dagger}\varphi_{jk}\right)}|\overline{i}\rangle\langle\overline{j}|=\mathbbm{1}=\sum_{ij}\delta_{ij}|\overline{i}\rangle\langle\overline{j}|}$, as an orthonormality condition on the $\varphi_{i}$s,

[TABLE]

Hence, to sample quantum channels, we simply need to find a parameterization for the orthonormal set $\{\varphi_{i}\}_{i=1}^{d}$.

Let us count the number of parameters needed. Since $A$ is upper triangular, $\varphi_{ik}$ has $(ik)+1$ generically nonzero entries, where $(ik)\equiv(i-1)d+(k-1)$ is a $d$-nary number. Each $\varphi_{i}$ thus has $K_{i}\equiv\sum_{k}[(ik)+1]=id^{2}-\frac{1}{2}d(d-1)$ nonzero entries. These nonzero entries are all complex, except for the $d^{2}$ of them in $\varphi_{dd}$, which are real. The orthonormality conditions on the $\varphi_{i}$s remove $d^{2}$ real parameters. Altogether then, the $\varphi_{i}$s are described by $2\sum_{i}K_{i}-d^{2}-d^{2}=d^{2}(d^{2}-1)$ real parameters, exactly the number needed to describe a quantum channel.

To specify an appropriate parametrization of the $\varphi_{i}$ set, it is convenient to reshuffle the rows of $\varphi_{i}$ so that all the identically-zero entries of each $\varphi_{i}$ are collected together. We first define the matrix

[TABLE]

Observe that the orthonormality conditions on the $\varphi_{i}$s translate into the requirement that $\Phi^{\dagger}\Phi=\mathbbm{1}$. Let $P$ be a $d^{3}\times d^{3}$ permutation matrix such that

[TABLE]

has columns $\psi_{i}$s, each of which is a reshuffled $\varphi_{i}$ with all identically zero entries located below the generically nonzero ones, i.e., the $k$th entry of $\psi_{i}$, which we denote as $\psi_{ik}$, is generally nonzero for $k=1,\ldots,K_{i}$, and zero for $k=K_{i}+1,\ldots,d^{3}$. Such a $P$ matrix exists because $A$ is upper triangular. Requiring $\Phi^{\dagger}\Phi=\mathbbm{1}$ is equivalent to demanding $\Psi^{\dagger}\Psi=\Phi^{\dagger}P^{-1}P\Phi=\mathbbm{1}$.

We are now ready to state the parameterization for the $\psi_{i}$s, thereby giving a parameterization for $\mathscr{S}_{\mathrm{TP}}$. We begin with $\psi_{d}$, parameterizing it with spherical coordinates so that it is normalized,

[TABLE]

where $\theta_{0}\equiv 0$ fixed, and the $S_{k}$s are recursively defined as $S_{k}=(\sin\theta_{k})S_{k+1}$, with $S_{K_{d}}=1$. Here, the $\phi_{k}$s for the $\psi_{dk}$s that come from the real entries of $\varphi_{dd}$ are understood to be set to zero (which ones they are, depends on the choice of $P$). $\psi_{d}$ is hence parameterized by real parameters $\theta_{1},\ldots,\theta_{K_{d}-1}$, and $K_{d}-d^{2}$ $\phi_{k}$ (real) parameters, giving $2K_{d}-1-d^{2}$ real parameters in all. Note the identity,

[TABLE]

so that the norm-square of $\psi_{d}$ is simply $\psi_{d}^{\dagger}\psi_{d}=\sum_{k=1}^{K_{d}}|\psi_{dk}|^{2}=S_{K_{d}}^{2}=1$, i.e., $\psi_{d}$ has length 1.

Next, let $v_{n}$, for $n=1,\ldots,K_{d-1}-1$, be the $d^{3}$-long column vector with the $k$th entry defined as

[TABLE]

Observe that $v_{n}$ is orthogonal to $\psi_{d}$, for every $n$, since

[TABLE]

One can check, in a similar manner, that the $v_{n}$ column vectors form an orthonormal set.

The span of $\{v_{n}\}_{n=1}^{K_{d-1}-1}$ lies in the orthogonal subspace of $\psi_{d}$. $\psi_{1},\psi_{2},\ldots\psi_{d-1}$ are to be orthogonal to $\psi_{d}$, so we can set them to be in the linear span of $\{v_{n}\}$. Note the both $\psi_{d-1}$ and $v_{K_{d-1}-1}$ have the same number ($=K_{d-1}$) of nonzero entries, the largest among the $\psi_{i}$s ($i=1,\ldots,{d-1}$) and $v_{n}$s. Specifically, we define

[TABLE]

where $V$ is the (non-square) matrix with columns $v_{1},v_{2},\ldots,v_{K_{d-1}-1}$. $\widetilde{\Psi}$ is defined such that its columns are the coefficients of the $\psi_{i}$s when expressed as a linear combination of the $v_{n}$s, i.e., $\psi_{i}=V\widetilde{\psi}_{i}=\sum_{n}\widetilde{\psi}_{in}v_{n}$, where $\widetilde{\psi}_{i}$ is the $i$th column of $\widetilde{\Psi}$, and $\widetilde{\psi}_{in}$ are its entries. Note that $V$ is a $d^{3}\times(K_{d-1}-1)$ matrix with the last $d^{3}-K_{d-1}$ rows completely zero, while $\widetilde{\Psi}$ is a $(K_{d-1}-1)\times(d-1)$ matrix.

Observe that the orthonormality of the $\psi_{i}$s, for $i=1,\ldots,{d-1}$ is equivalent to the orthonormality of the columns of $\widetilde{\Psi}$, i.e., $\widetilde{\Psi}^{\dagger}\widetilde{\Psi}=\mathbbm{1}$. This is then the same problem as before, for $\Psi$, with now one fewer column. We hence repeat the procedure above, parameterizing $\widetilde{\psi}_{d-1}$ using a new set of spherical coordinates ($\theta$s and $\phi$s; note that none of the $\phi$s are set to zero as the $\psi_{i\neq d}$s are generally complex), defining new $v$ vectors orthogonal to it, getting a new $\widetilde{\Psi}$, and so forth. We do this recursively until all $\psi_{i}$s are parameterized.

Let us check that the recursive procedure yields the right number of parameters for the full set of orthonormal $\psi_{i}$s. As mentioned earlier, in the first round, $\psi_{d}$ (and the $V$ there) is parameterized by $2K_{d}-1-d^{2}$ parameters, that subtraction of $d^{2}$ coming from the $d^{2}$ zero $\phi_{k}$s done for $\psi_{d}$ only. In the next round, $\psi_{d-1}$ is parameterized by an additional (on top of the ones that go into $V$) $2(K_{d-1}-1)-1$ real parameters; in yet the next round, $\psi_{d-2}$ is parameterized by an additional $2(K_{d-2}-2)-1$ real parameters; and so forth. Altogether then, we have $-d^{2}+\sum_{i=0}^{d-1}{\left[2(K_{d-i}-i)-1\right]}=d^{4}-d^{2}$ real parameters, exactly the right number needed for parameterizing $d$-dimensional quantum channels.

To illustrate how one applies the above parameterization, the case of qutrit channels is discussed in Appendix B. In the following sections, we make use of our parameterization in a HMC algorithm to sample quantum channels according to specified distributions, and demonstrate the usefulness of these samples in different applications. Before we get to that, however, let us mention a parameterization designed specifically for unital qubit channels, useful for one of our examples below.

III.2 Unital qubit channels

A useful class of quantum channels is the set of unital channels, those that preserve the identity operator, $\mathcal{E}(\mathbbm{1})=\mathbbm{1}$. The unitality condition can be stated in terms of the Choi state as the requirement

[TABLE]

A unital quantum channel thus has $\rho_{\mathcal{E}}$ such that $\textrm{tr}_{i}(\rho_{\mathcal{E}})=\mathbbm{1}$ for $i=1,2$, stating both the TP and unitality conditions. This is generally a difficult pair of conditions to impose, for a parameterization of unital channels with exactly the right number of parameters, as needed for HMC.

For unital qubit channels, however, this can be done in a straightforward manner, as we describe here MagicBasis . The Choi state of a qubit channel is a two-qubit state. Any two-qubit state (normalized to trace 2) can be written as

[TABLE]

where $\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of Pauli operators for the first qubit and $\bm{\tau}=(\tau_{x},\tau_{y},\tau_{z})$ is the vector of Pauli operators for the second qubit. (Here, the word “vector” is used in the physicist’s sense of a three-dimensional spatial vector.) $\bm{s}$ and $\bm{t}$ are the Bloch vectors for qubits 1 and 2, respectively; $\bm{C}$ is a dyadic, representable by a $3\times 3$ matrix of real numbers corresponding to the coefficients of $\sigma_{i}\tau_{j}$, for $i,j=x,y,z$. The TP condition requires $\bm{s}=0$; the unitality condition demands $\bm{t}=0$. The Choi state of a unital qubit channel thus takes the form

[TABLE]

Up to local unitary transformation, the dyadic $\bm{C}$ can always be chosen to be diagonal $\bm{C}_{\text{diag}}$. For $\rho_{\mathcal{E}}$ to be positive semi-definite, the three diagonal entries of $\bm{C}_{\text{diag}}$ must lie within a tetrahedron with the vertices

[TABLE]

where each vertex corresponds to one of four pairwise orthogonal maximally entangled two-qubit states. We parameterize the three entries of $\bm{C}_{\text{diag}}$ by the convex combination of the four vertices

[TABLE]

where

[TABLE]

Generally, the dyadic $\bm{C}$ can be written as

[TABLE]

where $\bm{R_{1}}$ and $\bm{R_{2}}$ are the rotation matrices representing the local unitary transformations (equivalently, spatial rotations in the Bloch-ball picture) of qubits 1 and 2, respectively. $\bm{R_{1}}$ and $\bm{R_{2}}$ can each be parameterized by three rotation angles. Altogether, we have a parameterization of the set of all unital qubit channels, specified by nine angle parameters.

IV Applications

The HMC algorithm is a method for generating random samples from any target distribution by making use of pseudo-Hamiltonian dynamics in a mock phase space. Upon identifying the parameters in the parameterizations given in Sec. III—which satisfy the requirements of not having superfluous parameters and no constraints—as the position variables in the mock phase space, we can employ the HMC algorithm. We give a brief review of the HMC algorithm in Appendix A. A more detailed discussion can be found in Ref. Seah+1:15 . In this section, we demonstrate the use of random samples of channels in three applications related to process tomography. That the examples are related to tomography simply reflects the authors’ original motivation and source of interest in the matter of channel sampling. The channel parameterization invented here and the resulting ability to sample according to a user-specified distribution using a HMC algorithm are applicable beyond tomography tasks.

IV.1 Error regions for process estimation

Whether one chooses to use the MLE or some other estimator for $\rho_{\mathcal{E}}$, the point estimator will not coincide exactly with the true $\rho_{\mathcal{E}}$ with finite data. It is important then to endow the point estimators with error regions expressing the uncertainty in our knowledge of the identity of the channel. Here, we adopt as error regions the notion of smallest credible regions (SCRs) proposed in Ref. OER13 . SCRs were originally proposed for the estimation of quantum states, whether they are TP states or not, but completely analogous notions can be defined for $\mathscr{S}_{\mathrm{TP}}$. Here, we examine the construction of SCRs for the task of quantum process estimation, as an application of our channel sampling algorithm. We first recall a few key points about SCRs pertinent to our discussion here; the reader is referred to OER13 for further details.

The SCR is the region—a set of states—in $\mathscr{S}_{\mathrm{TP}}$ with the smallest size for a chosen credibility. Size is the prior content of a region in $\mathscr{S}_{\mathrm{TP}}$, i.e., the prior (before any data are taken) probability that the true state is in the region; credibility is the posterior (after incorporating the data) content of that region. The SCRs are bounded-likelihood regions (BLRs), i.e., regions $\mathcal{R}_{\lambda}$ comprising all states with likelihood no smaller than a threshold fraction $\lambda\in[0,1]$ of the maximum likelihood $L_{\max}(D)$,

[TABLE]

with $\mathcal{R}_{0}=\mathscr{S}_{\mathrm{TP}}$. The size $s_{\lambda}$ of the BLR $\mathcal{R}_{\lambda}$ is its prior content, and its credibility $c_{\lambda}$ is its posterior content,

[TABLE]

with ${s_{0}=c_{0}=1}$ when ${\lambda=0}$. The volume element $(\textrm{d}\rho)$ expresses the prior distribution; (\textrm{d}\rho)\frac{\mbox{\footnotesizeL(D|\rho)}}{\mbox{\footnotesizeL(D)}} is the posterior distribution. $L(D)\equiv\int_{\mathcal{R}_{0}}(\textrm{d}\rho)L(D|\rho)$, a normalizing factor, is the likelihood of obtaining the data $D$ for the chosen prior. For tomography problems, it is often natural to state the prior distribution in terms of the POVM-induced probabilities [see Eq. (10)],

[TABLE]

where $w_{0}(p)$ is the prior density, nonzero only for $p\equiv(p_{1}^{(1)},p_{1}^{(2)},\ldots,p_{2}^{(1)},\ldots)$ that corresponds to a $\rho\in\mathscr{S}_{\mathrm{TP}}$, and $(\textrm{d}p)\equiv\textrm{d}p_{1}^{(1)}\textrm{d}p_{1}^{(2)}\cdots~{}$.

To report the error region for an experiment with data $D$, following the scheme of Ref. OER13 , $s_{\lambda}$ and $c_{\lambda}$ are calculated for all values of $\lambda$. The error regions are reported by plotting $s_{\lambda}$ and $c_{\lambda}$ as functions of $\lambda$. For a desired level of credibility, the $\lambda$ value is read off, and the error region is the $\mathcal{R}_{\lambda}$ for that value of $\lambda$. The size and credibility of a BLR [see Eq. (38)] cannot, in general, be computed analytically, due to the complicated integration region. Instead, we make use of MC integration: We generate random samples using HMC according to the prior and posterior distributions; the size and credibility are then the fractions of points contained in the BLR for the two distributions.

A related concept is the plausible region Evans:15 . This is the set of all points in $\mathscr{S}_{\mathrm{TP}}$, for which the data provide evidence in favor of—${L(D|\rho)>L(D)}$. The plausible region is in fact a BLR, with a critical value of $\lambda$,

[TABLE]

Once we have computed the size and credibility curves, we can also identify the plausible region for the data.

As a first example, we look at single-qubit channels. The input states $\rho^{(i)}$ for process tomography are taken to be the tetrahedron states,

[TABLE]

where $\bm{a}_{i}=3^{-1/2}\bm{v}_{i}$ with the vertex vectors of Eq. (III.2). For every $i$, we use the same POVM, the four-outcome tetrahedron measurement, with outcomes

[TABLE]

We simulate data using an amplitude-damping channel described by the Kraus operators

[TABLE]

where $\gamma$, the damping parameter, is set to $0.4$. The matrices above refer to the computational basis. For the matrices appear in the rest of this paper, it should be assumed that they refer to the computational basis as well. 24 copies of each input state $\rho^{(i)}$ are measured (simulated), giving a total of 96 counts over the four input states. The simulated data are reported in Table 1.

For the prior distribution, we choose the conjugate prior,

[TABLE]

where ${\bar{p}=\{\bar{p}_{k}^{(i)}\}}$ corresponds to the Born probabilities [see Eq. (10)] for an amplitude-damping channel with $\gamma=0.5$, expressing our prior belief that that is the actual channel. Figure 1(a) shows the size and credibility curves, obtained from MC integration using 500,000 sample points generated from HMC with the channel parameterization of Sec. III. The critical $\lambda$ value for the plausible region is indicated with a red dashed line, with size value $s=0.2102$ and credibility value $c=0.8586$. The true channel is contained in all BLRs with $\lambda<0.0302$ and $c_{\lambda}>0.5511$, and is thus in the plausible region.

Now, qubit channels are simple to characterize and there are many ways of sampling from the space of qubit channels. It is hence useful to see how our sampling algorithm works for examples beyond the qubit situation, for which proper sampling is more challenging. As a second example, we consider an amplitude-damping qutrit (three-dimensional quantum system) channel with the Kraus operators

[TABLE]

for $\gamma_{1}=0.1$ and $\gamma_{2}=0.5$.

The POVM used is one of the symmetric, informationally complete POVM (SIC-POVM) from the one-parameter family of qutrit SIC-POVMs. It can be described by a set of states $\{\ket{\mu_{i}}\}$; when written in the computational basis, they are given explicitly by

[TABLE]

where $\omega=\textrm{e}^{\textrm{i}{2\pi/3}}$, $\omega^{*}=\omega^{2}$, and $1+\omega+\omega^{2}=0$. The POVM elements are

[TABLE]

The input states are

[TABLE]

For each of the input states, the number of copies measured is 27, giving a total of 243 counts. The simulated data are reported in Table 2.

The prior is the primitive prior, i.e., $w_{0}(p)$ is a constant wherever it is nonzero. Figure 1(b) shows the size and credibility curves, obtained from MC integration with 100,000 sample points using HMC and our channel parameterization. As before, the critical $\lambda$ value for the plausible region is indicated by the vertical dashed line. The size and credibility of the plausible region are $s=0.0032$ and $c=0.9990$ respectively. The true channel is contained in all BLRs with $\lambda<1.1560\times 10^{-6}$, and is thus in the plausible region.

IV.2 Marginal likelihood for channel properties

Often, one is only interested in certain properties of a channel, like the fidelity between the output of the channel and its input, rather than a full channel description in the form of its process matrix. If one could directly measure that one quantity of interest, one expects to accomplish the estimation task with significantly fewer uses of the channel than needed for full tomography. However, a direct measurement of the quantity of interest may be difficult to design and implement, while the process tomography measurement is often standard procedure. Even in the latter case, one should still estimate the quantity of interest directly from the tomography data, rather than first estimating the full process matrix and then computing the quantity of interest from that estimate OEI16 .

The key ingredient in making inferences about a property $F$ of a channel from tomographic data $D$ is the marginal likelihood, obtained by integrating the full likelihood $L(D|p)$ over the irrelevant parameters,

[TABLE]

where $W_{r,D(0)}(F)$ is the integral in the numerator(denominator). $f(p)$ is the function that expresses $F$ in terms of the tomographic probabilities $p$, and $w_{r}(p)$ is the prior density on $p$, which induces a prior density on $F$. \delta\bigl{(}F-f(p)\bigr{)} is the Dirac delta function that enforces ${f(p)=F}$. Once we have the marginal likelihood, we can proceed in an analogous way as in Sec. IV.1 to construct the smallest credible interval (SCI) and the plausible interval for $F$, as well as perform other statistical inference tasks based on the marginal likelihood.

We thus need a general procedure for computing the marginal likelihood $L(D|F)$. In Ref. OEI16 , an iterative algorithm was developed for that purpose, requiring the use of random samples according to specified distributions. The reader is referred to Ref. OEI16 for the full description of the iterative algorithm, and to Appendix C for the details relevant for our examples below. Here, we give only a brief account of the basic ideas. The delta functions in the defining equation (59) are difficult to handle in a numerical evaluation of the integrals. Instead, we evaluate the antiderivatives $P_{r,i}(F)$, with respect to $F$, of $W_{r,i}(F)$,

[TABLE]

with step functions in place of the delta functions. $P_{r,i}$ can be computed by MC integration. The results are closely fitted with several-parameter functions, and then differentiated to give $W_{r,i}$, and hence the marginal likelihood. This procedure works, in principle; in practice, one runs into numerical accuracy problems. If $w_{r}(p)$ has little weight over some range of $F$, a rather generic situation, $P_{r,0}$ will be very flat there, and its derivative cannot be reliably estimated. To overcome this problem, the crux is to note that, because of the delta functions, the marginal likelihood is invariant under the replacement $w_{r}(p)\rightarrow w_{r}(p)g(f(p))$ for any function $g(F)$ positive over the entire range of $F$. We thus have the freedom to choose the $w_{r}(p)$ used to evaluate $L(D|F)$. This freedom of choice is exploited in the iterative procedure described in Ref. OEI16 , where the estimate of $W_{r,0}$ is successively improved by using an $w_{r}(p)$ modified by the previous (possibly inaccurate) estimate of $W_{r,0}$, until the desired convergence level is reached. Each iterative step requires the ability to sample according to the new $w_{r}(p)$; that is where the HMC algorithm, permitting sampling in accordance to a user-specified distribution, comes in.

Below, we carry out the iterative algorithm and compute the marginal likelihood for two common channel properties, average fidelity $F_{\mathrm{avg}}$ and minimum fidelity $F_{\mathrm{min}}$. We make use of the HMC algorithm made possible by our channel parameterization of Sec. III. Both examples are for qubit channels, and use the same (simulated) tomographic data obtained from tetrahedron input states [see Eq. (41)] and the tetrahedron POVM [see Eq. (42)] for the true channel

[TABLE]

a Pauli channel. Here, the $\sigma_{i}$s are the standard Pauli operators, and $(p_{x},p_{y},p_{z})=(0.05,0.15,0.2)$. The data are generated from 96 uses of the channel. The simulated data are reported in Table 3.

We regard the Pauli channel as noise acting on our quantum system. We are interested in the fidelity measures, $F_{\mathrm{avg}}$ and $F_{\mathrm{min}}$, quantifying the effect of this noise channel on our system.

IV.2.1 Average Fidelity

The average fidelity $F_{\mathrm{avg}}$ is defined here as the (squared-)fidelity between the input and output of the channel $\mathcal{E}$, averaged over all input pure states according to the Haar measure. We write $F(\psi,\rho)\equiv\langle\psi|\rho|\psi\rangle$ for the square of the fidelity between a pure state $\psi\equiv|\psi\rangle\langle\psi|$ and an arbitrary state $\rho$. Then, the average fidelity for the channel $\mathcal{E}$ is

[TABLE]

Here, $\textrm{d}U$ is the Haar measure for the space of unitary operators, and $\psi_{0}$ is some fiducial pure state. In arriving at the last line, we have used a standard result of the twirling operation Emerson05 (namely, the expression in the brackets in the second-to-last line), with $q$ given by

[TABLE]

where $O_{i}$s are all the traceless elements of an orthonormal (according to the Hilbert-Schmidt inner product) operator basis, containing an element proportional to the identity operator, for the $d$-dimensional $\mathcal{H}$. In the qubit case, $q$ has the explicit formula,

[TABLE]

where we have chosen the $*$ map such that $|\overline{i}\rangle=|i\rangle$ for $\{|i\rangle\}_{i=0}^{1}$, the $\sigma_{z}$-basis for the qubit (see comment about this choice in the second paragraph of Sec. II).

We use the iterative procedure of Ref. OEI16 to compute the marginal likelihood $L(D|F_{\mathrm{avg}})$, for $F\equiv F_{\mathrm{avg}}$. The final result is shown in Fig. 2; the intermediate steps of the iterative algorithm are described in Appendix D.1. With the marginal likelihood at hand, as an example of its usefulness, we can construct, as in Sec. IV.1, the SCI for our estimate of $F_{\mathrm{avg}}$. Figure 3(a) gives the size and credibility curves, as well as the critical $\lambda$ value for the plausible region. Figure 3(b) shows the SCI for $F_{\mathrm{avg}}$ for different credibility values. The horizontal black line specifies the plausible interval, which includes the true value of ${F_{\mathrm{avg}}=0.7333}$ (indicated with an arrow).

IV.2.2 Minimum fidelity of unital qubit channels

As a second example, also to illustrate the use of the parameterization of the unital qubit channels of Sec. III.2, we look at the minimum, or worst-case, (squared-)fidelity of a unital channel. The minimum fidelity for a channel $\mathcal{E}$ is the fidelity of the output of $\mathcal{E}$ with its (pure) input, minimized over all input states, i.e.,

[TABLE]

In the qubit case, $F_{\mathrm{min}}$ can be written explicitly using the Bloch-ball representation as

[TABLE]

where $\bm{s}$ is the Bloch vector of the input state $\psi$, and $\bm{s}_{\mathcal{E}}$ is that of the output $\mathcal{E}(\psi)$. For a unital qubit channel, $\bm{s}_{\mathcal{E}}$ is the image of a linear map on the Bloch vector: $\bm{s}_{\mathcal{E}}=\bm{Ms}$. The minimum fidelity can thus be written simply as

[TABLE]

where $\mu_{\text{min}}$ is the smallest eigenvalue of $\frac{1}{2}(\bm{M}+\bm{M}^{\dagger})$. This provides the direct connection between the unital qubit channel and $F_{\mathrm{min}}$, and, in particular, allows us to express $F_{\mathrm{min}}$ in terms of the tomographic probabilities associated with a channel $\mathcal{E}$.

Here, we assume the promise that the unknown channel is a unital one; the Pauli channel used to simulated the data is indeed unital. In effect, this unitality assumption restricts the relevant space of Choi states dual to the channels, to a strict subset of $\mathscr{S}_{\mathrm{TP}}$, namely, to those that also satisfy Eq. (30). Any channel sampling is thus done only from this subset. Using the parameterization of Sec. III.2, we employ HMC integration to compute the marginal likelihood $L(D|F_{\mathrm{min}})$. The result is given in Fig. 4; the intermediate steps are provided in Appendix D.2. With this marginal likelihood, one can construct the corresponding SCIs and the plausible region, as well as perform other statistical inferences about the unital qubit channel.

IV.3 Model selection

Often, one may not need the full generality of a CPTP channel to describe the dynamics of a quantum system. Instead, a simpler model with fewer parameters may suffice. Simpler models are computationally easier to work with, are likely more easily motivated from a physical standpoint, and may already describe the tomographic data well. One can phrase this problem as one of model selection in statistics, where the best model, among a few candidate models, is chosen, given the available data. Here, we discuss the quantum problem of model selection for channel families. Our sampling algorithm is used for two purposes here: (1) to evaluate a criterion—based on the notion of relative belief—for the “best” model; (2) to assess and compare the performance of different model selection criteria by testing them on many randomly chosen true channels.

Two criteria for model selection commonly used in classical problems are the Akaike Information Criterion (AIC) AIC and the Bayesian Information Criterion (BIC) BIC . The AIC is based on the quantity (which we denote also as “AIC”),

[TABLE]

where $k$ is the number of parameters in the model and $L_{\max}$ is the maximum value of the likelihood of the model for the data. The best model is the one with the smallest AIC value. The BIC is defined in a similar manner, but uses the value of $N$, the number of copies measured,

[TABLE]

The best model according to this criterion is again the one with the smallest BIC value.

Another approach to model selection is based on the relative belief ratio (RBR) of Ref. Evans:15 . The RBR of a model $M$ is the ratio of its posterior to prior probabilities,

[TABLE]

where

[TABLE]

If the posterior probability for a model $M$ increases after the data, i.e. ${\textrm{RBR}(M|D)>1}$, the data provide evidence in favor of the model; the data provide evidence against the model if ${\textrm{RBR}(M|D)<1}$. It is also useful to have a measure of strength of evidence, since the data might provide evidence in favor of more than one model from our candidate set, and one would like some basis of choosing among those models. The RBR value by itself is not a measure of the strength of evidence (see Ref. Evans:15 for a discussion of various aspects, and also Ref. Evans19 ). We supplement it with the posterior probability

[TABLE]

for the model $M_{0}$ in question, and $M$ ranges over the set of candidate models. If ${\textrm{RBR}(M_{0}|D)>1}$ and $P_{M_{0}}$ is large, then there is strong evidence in favor of $M_{0}$. The best model, according to the RBR criterion of relative belief ratio, is the one with the largest posterior probability $P_{M}$, among all candidate models with ${\textrm{RBR}(M|D)>1}$.

As an example, we consider as candidate models five nested qubit channel families: dephasing channels $\subset$ Pauli channels $\subset$ symmetric unital channels $\subset$ unital channels $\subset$ general CPTP channels; see Fig. 5. The smallest set is the 1-parameter family of dephasing channels,

[TABLE]

and its Choi state is given by

[TABLE]

The set of Pauli channels is a 3-parameter family,

[TABLE]

for $\mathbf{p}\equiv(p_{x},p_{y},p_{z})$, $p_{i}\geq 0$, and $\sum_{i}p_{i}\leq 1$. The Choi state of a Pauli channel is of the form

[TABLE]

where $p_{I}=1-p_{x}-p_{y}-p_{z}$. The 6-parameter family of symmetric unital channels refers to the subset of unital qubit channels such that $\bm{R_{1}}=\bm{R_{2}}$ in Eq. (36). We then have the 9-parameter family of unital qubit channels, and lastly, the 12-parameter set of all CPTP qubit channels.

A natural prior on the model space is one that puts equal weights on each family. This is easily defined by the sampling procedure: the prior sample is constructed by generating 500,000 sample points with the primitive prior for each family. For the dephasing channel, the sample is generated by sampling $p$ uniformly from $[0,1]$. For the Pauli channel, we obtain the sample by generating $(p_{x},p_{y},p_{z})$ uniformly from the 3-simplex. For the symmetric unital, the unital, and the general channels, we make use of HMC and the parameterizations in Sec. III to generate the sample points. Note that in the numerical procedure that generates the samples for, say, the set of Pauli channels, we will never come across a sample point that is exactly a dephasing channel with ${p_{x}=0=p_{y}}$. Thus, even though the channel families are nested sets, one can consider each family to have prior probability of $\frac{1}{5}$. We use this prior to compute the RBR criterion for simulated data of different sizes. For our choice of prior, $P(M)=\frac{1}{5}$ for all models and $P(M|D)$ is calculated by taking the average of $\frac{L(D|\rho)}{L(D)}$ over the sample points for each model. $P_{M_{0}}$ is computed by summing the posterior probabilities $P(M|D)$ of all the models $M$ with the same RBR as model $M_{0}$. Typically, $P_{M_{0}}=P(M_{0}|D)$.

To assess the performance of the three model-selection criteria, for each family of channels, we randomly (according to the primitive prior, as described above) draw 1000 channels. For each channel, we simulate data—with tetrahedron input states and a tetrahedron measurement [see Eqs. (41) and (42)]—for ${N=20}$, $50$, $100$, $1\,000$, $10\,000$, and $100\,000$ copies measured, and evaluate the AIC, BIC, and RBR criteria for that data. Table 4 shows the conclusions when the three criteria are applied to the simulated data. When the number of measured copies is very small, i.e., $N=20$, the results based on AIC and BIC show a strong bias towards simpler (i.e., fewer-parameters) models. In particular, both criteria rarely identify the right model when the true channel comes from the unital or general families. Results based on RBR, however, show significantly more instances where the correct model is identified for the more complex (i.e., more parameters) models. For a moderate number of measured copies, i.e. ${N=1\,000}$, AIC and RBR give equally good results, whereas BIC shows a slight bias towards the simpler models. When the number of measured copies is very large, i.e., ${N=100\,000}$, results based on BIC are most accurate whereas results based on AIC have a slight bias to the more complex models. RBR also performs well in this regime.

Another aspect that we can check easily with our sampling procedure is the bias in the prior. This is particularly important for model selection based on the RBR criterion, to be sure that the probability of drawing a wrong conclusion is low. For example, for data that are typical for a unital channel, if we were to conclude regularly that there is evidence in favor of the general CPTP model and evidence against the unital model, there is bias in favor of the general CPTP model and bias against the unital model. To check for the bias, we draw 1000 random channels from each of the channel families and simulate data based on these true channels. The number of instances where the simulated data provide evidence against each of the four candidate models are calculated. The results are shown in Table 5. As can be seen from the table, there is no significant bias in the prior when $N\geq 100$, and the bias decreases as the number of measured copies increases.

V Conclusions

In this work, we constructed one exact parameterization for the space of CPTP channels. This parameterization has no superfluous parameters, and requires no imposition of any added constraints. These features make it possible to use the parameterization in a HMC algorithm, for producing high-quality—in terms of low correlations—samples of CPTP channels from a user-specified distribution. We demonstrated the usefulness of our parameterization in sampling applications taken from quantum process tomography. The method applies to general quantum channel sampling problems.

While our parameterization serves the purpose, it is, of course, just one of the many parameterizations that could be used in a HMC algorithm for sampling from the quantum channel space. For example, it is conceivable that a useful parameterization of a channel can be given in terms of the marginals and the copula of the respective Choi state copula . This is unexplored territory.

A useful extension of this work will be to discover also an exact parameterization for the case of CPTP and unital channels. As discussed above, this additional requirement of unitality presents difficulties that can be easily overcome only in the qubit situation. The parameterization for the space of CPTP, unital channels beyond the qubit case, remains an open problem. Note that such a parameterization will give also a possibly useful description of the space of all bipartite mixed quantum states with completely mixed states on both the single-party states; our current parameterization gives the larger space of states where only one of the two single-party states is completely mixed.

Acknowledgements.

This work is supported in part by the Ministry of Education, Singapore (through grant number MOE2016-T2-1-130). HKN is also supported by Yale-NUS College (through a start-up grant). The Centre for Quantum Technologies is a Research Centre of Excellence funded by the Ministry of Education and the National Research Foundation of Singapore.

Appendix A Hamiltonian Monte Carlo (HMC)

HMC makes use of pseudo-Hamiltonian dynamics in a mock phase space. The parameters of interest are identified as the position variables $\theta$ and fictitious momentum variables $\vartheta$ are introduced. The Hamiltonian is defined as

[TABLE]

where $w(\theta)$ is the target distribution. Any reasonable target distribution is permitted and, therefore, one can sample in accordance with any $w(\theta)$.

The HMC algorithm generates a set of sample points which follows the target distribution $w(\theta)$. The HMC algorithm is stated as follows Seah+1:15 :

HMC algorithm

Set $j=1$ and choose an arbitrary starting point $\theta^{(1)}$. 2. 2.

Generate $\vartheta^{(j)}$ from a multivariate Gaussian distribution with mean zero and unit variance. 3. 3.

Solve the Hamiltonian equations of motion

[TABLE]

with the initial conditions $(\theta,\vartheta)|_{t=0}=(\theta^{(j)},\vartheta^{(j)})$ to obtain $(\theta^{*},\vartheta^{*})=(\theta,-\vartheta)|_{t=T}$. 4. 4.

Calculate the acceptance ratio

[TABLE] 5. 5.

Draw a random number $b$ uniformly from $[0,1]$. If $b<a$, set $\theta^{(j+1)}=\theta^{*}$; otherwise, set $\theta^{(j+1)}=\theta^{(j)}$. 6. 6.

Set $j=j+1$. If $j$ equals the desired number of samples, escape the loop; otherwise, return to step 2.

Note that the $\vartheta$ distribution in step 2 is proportional to the kinetic-energy factor in $\textrm{e}^{-H(\theta,\vartheta)}$; the HMC algorithm enforces a $\theta$ distribution proportional to the potential-energy factor in $\textrm{e}^{-H(\theta,\vartheta)}$, which is $\textrm{e}^{\log w(\theta)}=w(\theta)$, the target distribution. If the differential equations in (78) can be solved exactly, then the acceptance ratio $a=1$. In practice, the differential equations must be discretized. This is done by the leapfrog method. Due to the discretization error, the acceptance ratio will not be 1 generally. The leapfrog method should be implemented such that the acceptance ratio is around the optimal value of $65\%$ Neal11 .

Appendix B Parameterizing qutrit channels

Here, we report an explicit application of the parameterization of Sec. III, for the case of qutrit channels. We start with the permutation matrix $P$ that reshuffles $\varphi_{i}$s into $\psi_{i}$s, with the identically zero entries located below the generically nonzero ones. A $P$ that can accomplish this is one such that

[TABLE]

After the permutation, we have

[TABLE]

To parameterize the $\psi_{i}$s such that they are orthonormal, we first parameterize $\psi_{3}$, of unit length,

[TABLE]

Recalling that $\varphi_{33}$ [see Eq. (12)] is a $d^{2}$-entry real column, and with the $P$ given above, $\psi_{3,4}$, $\psi_{3,5}$, $\psi_{3,6}$, $\psi_{3,13}$, $\psi_{3,14}$, $\psi_{3,15}$, $\psi_{3,22}$, $\psi_{3,23}$, $\psi_{3,24}$ are real. Thus, $\phi_{4}$, $\phi_{5}$, $\phi_{6}$, $\phi_{13}$, $\phi_{14}$, $\phi_{15}$, $\phi_{22}$, $\phi_{23}$, $\phi_{24}$ are set to zero. Then, we define $\{v_{n}\}_{n=1}^{14}$ which lie in the orthogonal subspace of $\psi_{3}$ as follows,

[TABLE]

To make $\psi_{1}$ and $\psi_{2}$ orthogonal to $\psi_{3}$, we set them to be in the span of $\{v_{n}\}$,

[TABLE]

The orthonormality of $\psi_{1}$ and $\psi_{2}$ is equivalent to the orthonormality of $\widetilde{\psi}_{1}$ and $\widetilde{\psi}_{2}$. We simply need to repeat the previous procedure. We parameterize $\widetilde{\psi}_{2}$ to be of unit length,

[TABLE]

Next, we define $\{u_{n}\}_{n=1}^{4}$, each orthogonal to $\psi_{2}$,

[TABLE]

Finally, to have $\widetilde{\psi}_{1}$ normalized and orthogonal to $\widetilde{\psi}_{2}$, we set

[TABLE]

where

[TABLE]

We check that we have the right number of parameters. The parameters used above are $\theta_{1},...,\theta_{23}$, $\phi_{1},...,\phi_{24}$ (nine of these are set identically to zero), $\widetilde{\theta}_{1},...,\widetilde{\theta}_{13}$, $\widetilde{\phi}_{1},...,\widetilde{\phi}_{14}$, $\bar{\theta}_{1},\bar{\theta}_{2},\bar{\theta}_{3}$, and $\bar{\phi}_{1},...,\bar{\phi}_{4}$, giving a total of $72=3^{2}(3^{2}-1)$ parameters, as needed for specifying qutrit channels.

Appendix C Iterative algorithm for estimating the marginal likelihood

To estimate the marginal likelihood reliably, we follow the procedure in Ref. OEI16 . For the following discussion, we assume

[TABLE]

for the sake of simplicity. First, we note that the integrands in (59) are ill-suited for MC integration due to the presence of the Dirac delta factors. We consider the antiderivatives

[TABLE]

and

[TABLE]

With a sample of $w_{r}(p)$ and \frac{\mbox{\footnotesizew_{r}(p)L(D|p)}}{\mbox{\footnotesizeL(D)}}, we can evaluate the antiderivatives for various values of $F$ and fit them with several-parameters functions. From the fitted functions, we can then calculate the derivatives

[TABLE]

and

[TABLE]

and obtain the marginal likelihood by

[TABLE]

A problem arises when $P_{r,0}(F)$ is very close to a constant over some range of values of $F$. The common situation is that $P_{r,0}(F)$ is very close to zero for a range of values near ${F=0}$ and very close to one for a range of values near ${F=1}$. MC integration is not precise enough to distinguish ${P_{r,0}(F)\gtrsim 0}$ from ${P_{r,0}(F)=0}$ and ${P_{r,0}(F)\lesssim 1}$ from ${P_{r,0}(F)=1}$. As a result, the estimated value of $W_{r,0}(F)$ will be equal to zero over those range of values. We cannot get a reliable estimation of $L(D|F)$ in this situation since $W_{r,0}(F)$ is the denominator in Eq. (94). To overcome this problem, we note that we can do the replacement

[TABLE]

with an arbitrary function ${g(F)>0}$ without changing the value of $L(D|F)$.

The procedure for obtaining a reliable estimation of $L(D|F)$ is as follows:

Sample according to $w_{r}(p)$. Use this sample to calculate $P_{r,0}(F)$. Fit a several-parameters function to $P_{r,0}(F)$ and obtain $W_{r,0}(F)$ by differentiating the fitted function. 2. 2.

Sample according to \widetilde{w}_{r}(p)=\frac{\mbox{\footnotesizew_{r}(p)}}{\mbox{\footnotesizeW_{r,0}\bigl{(}f(p)\bigr{)}}}. Use this sample to calculate

[TABLE]

Fit a several-parameters function to $\widetilde{P}_{r,0}(F)$ and obtain $\widetilde{W}_{r,0}(F)$ by differentiating the fitted function. 3. 3.

Sample according to \frac{\mbox{\footnotesize\widetilde{w}{r}(p)L(D|p)}}{\mbox{\footnotesizeL(D)}}=\frac{\mbox{\footnotesizew{r}(p)L(D|p)}}{\mbox{\footnotesizeW_{r,0}\bigl{(}f(p)\bigr{)}L(D)}}. Use this sample to calculate

[TABLE]

Fit a several-parameters function to $\widetilde{P}_{r,D}(F)$ and obtain $\widetilde{W}_{r,D}(F)$ by differentiating the fitted function. 4. 4.

Obtain the marginal likelihood from

[TABLE]

The reason that we can have a reliable estimation of $L(D|F)$ using $\widetilde{W}_{r,0}(F)$ obtained in step 2 is as follows. Suppose the exact value of $W_{r,0}(F)$ is known, $\widetilde{W}_{r,0}(F)$ will be equal to $1$ and $\widetilde{P}_{r,0}(F)$ will be equal to $F$. If the exact values of $W_{r,0}(F)$ are not known, but we have a good approximation for $W_{r,0}(F)$ from step 1 and use it for the calculation of $\widetilde{P}_{r,0}(F)$ in step 2, the $\widetilde{P}_{r,0}(F)$ that we obtain will still be quite close to $F$ and $\widetilde{W}_{r,0}(F)$ will be nonzero for all range of $F$ values.

In step 1, $P_{r,0}(F)$ can be fitted with a linear combination of regularized incomplete beta functions

[TABLE]

that is

[TABLE]

with the fitting parameters $a_{1},\dots,a_{N-1},$ $b_{1},\dots,b_{N-1},$ $w_{1},\dots,w_{N-1}$. $a_{\text{min}}$ and $b_{\text{min}}$ are fixed by the power laws satisfy by $P_{r,0}(F)$ near $F=0$ and $F=1$,

[TABLE]

and

[TABLE]

In step 2, a truncated Fourier series of the form

[TABLE]

is usually a good fitting function. In step 3, $\widetilde{P}_{r,D}(F)$ can be fitted with a smoothing spline.

Appendix D Intermediate results for the estimation of the marginal likelihood

D.1 Average gate fidelity

The green dots in Fig. 6(a) show the values of $P_{r,0}(F_{\mathrm{avg}})$ obtained by a MC integration with $1\,000\,000$ sample points. The MC integration is not precise enough to distinguish ${P_{r,0}(F_{\mathrm{avg}})\gtrsim 0}$ from ${P_{r,0}(F_{\mathrm{avg}})=0}$ near ${F_{\mathrm{avg}}=\frac{1}{3}}$ and to distinguish ${P_{r,0}(F_{\mathrm{avg}})\lesssim 1}$ from ${P_{r,0}(F_{\mathrm{avg}})=1}$ near ${F_{\mathrm{avg}}=1}$. Therefore, a reliable approximation for W_{r,0}(F_{\mathrm{avg}})=\frac{\mbox{\footnotesize\partial}}{\mbox{\footnotesize\partial F_{\mathrm{avg}}}}P_{r,0}(F_{\mathrm{avg}}) cannot be obtained. To overcome this problem, we follow the procedure stated in Appendix C. First, we fit the green dots with a three-term fitting function of the form of Eq. (C) with ${F=\frac{3}{2}(F_{\mathrm{avg}}-\frac{1}{3})}$, ${a_{\text{min}}=3}$ and ${b_{\text{min}}=\frac{21}{2}}$. The black curve is the fitted curve of $P_{r,0}(F_{\mathrm{avg}})$. The fitting parameters are shown in the inset table. $\widetilde{P}_{r,0}(F_{\mathrm{avg}})$ is obtained from a MC integration with $1\,500\,000$ sample points and shown as the blue dots in Fig. 6(a).

The $\widetilde{P}_{r,0}(F_{\mathrm{avg}})$ is quite close to the straight line ${\frac{3}{2}(F_{\mathrm{avg}}-\frac{1}{3})}$. The $\widetilde{P}_{r,0}(F_{\mathrm{avg}})$ after subtracting the straight line ${\frac{3}{2}(F_{\mathrm{avg}}-\frac{1}{3})}$ is shown as the blue dots in Fig. 6(b). The blue curve shows the fitting curve, a truncated Fourier series whose Fourier amplitudes are reported in Fig. 6(c).

$\widetilde{P}_{r,D}(F_{\mathrm{avg}})$ is evaluated by a MC integration with $1\,500\,000$ sample points and it can be fitted with a smoothing spline. The marginal likelihood shown in Fig. 2 is obtained from the ratio of $\widetilde{W}_{r,D}(F_{\mathrm{avg}})$ and $\widetilde{W}_{r,0}(F_{\mathrm{avg}})$.

D.2 Worst-case fidelity of a unital qubit channel

The green dots in Fig. 7(a) show the values of $P_{r,0}(F_{\mathrm{min}})$ from a MC integration with $1\,000\,000$ sample points. The MC integration is not precise enough to distinguish ${P_{r,0}(F_{\mathrm{min}})\gtrsim 0}$ from ${P_{r,0}(F_{\mathrm{min}})=0}$ near ${F_{\mathrm{min}}=0}$ and to distinguish ${P_{r,0}(F_{\mathrm{min}})\lesssim 1}$ from ${P_{r,0}(F_{\mathrm{min}})=1}$ near ${F_{\mathrm{min}}=1}$. Therefore, a reliable approximation for {W_{r,0}(F_{\mathrm{min}})=\frac{\mbox{\footnotesize\partial}}{\mbox{\footnotesize\partial F_{\mathrm{min}}}}P_{r,0}(F_{\mathrm{min}})} cannot be obtained. To overcome this problem, we follow the procedure stated in Appendix C. First, we fit the green dots with a three-term fitting function of the form in Eq. (C) with ${F=F_{\mathrm{min}}}$, ${a_{\text{min}}=4}$, and ${b_{\text{min}}=\frac{15}{2}}$. The black curve is fitted to the numerical values for $P_{r,0}(F_{\mathrm{min}})$. The fitting parameters are shown in the inset table. $\widetilde{P}_{r,0}(F_{\mathrm{min}})$ is obtained by a MC integration with $1\,500\,000$ sample points and shown as the blue dots in Fig. 7(a).

The values of $\widetilde{P}_{r,0}(F_{\mathrm{min}})$ are quite close to the straight line $F_{\mathrm{min}}$. The corresponding values after subtracting this straight line make up the blue dots in Fig. 7(b). The blue fitting curve is a truncated Fourier series with the Fourier amplitudes of Fig. 7(c).

$\widetilde{P}_{r,D}(F_{\mathrm{min}})$ is evaluated by a MC integration with $1\,500\,000$ sample points and it can be fitted with a smoothing spline. The marginal likelihood shown in Fig. 4 is the ratio of $\widetilde{W}_{r,D}(F_{\mathrm{min}})$ and $\widetilde{W}_{r,0}(F_{\mathrm{min}})$.