Convex optimization of programmable quantum computers

TL;DR

This paper formulates the problem of optimizing program states in programmable quantum gate arrays as a convex optimization task, enabling efficient solutions using semidefinite programming and machine learning techniques.

Contribution

It demonstrates that finding the optimal program state is a convex optimization problem, providing a practical method to improve programmable quantum processors.

Findings

Optimal program states can be computed efficiently using convex optimization.

The approach applies to various quantum processor designs, including teleportation-based and parametric circuits.

This method advances the understanding of quantum processor programmability and simulation accuracy.

Abstract

A fundamental model of quantum computation is the programmable quantum gate array. This is a quantum processor that is fed by a program state that induces a corresponding quantum operation on input states. While being programmable, any finite-dimensional design of this model is known to be nonuniversal, meaning that the processor cannot perfectly simulate an arbitrary quantum channel over the input. Characterizing how close the simulation is and finding the optimal program state have been open questions for the past 20 years. Here, we answer these questions by showing that the search for the optimal program state is a convex optimization problem that can be solved via semidefinite programming and gradient-based methods commonly employed for machine learning. We apply this general result to different types of processors, from a shallow design based on quantum teleportation, to deeper schemes relying on port-based teleportation and parametric quantum circuits.

Figures

Equations

$${\mathcal{E}}_{\pi }(\rho )=\underset{2}{\mathrm{Tr}}\left[Q(\rho \otimes \pi )\right],$$

$$C_{\diamond }(\pi ):={\Vert \mathcal{E}-{\mathcal{E}}_{\pi }\Vert }_{\diamond },$$

$$\text{Find }\tilde{\pi }\text{ such that }{C}_{\diamond }(\tilde{\pi })=\underset{\pi }{\min }{C}_{\diamond }(\pi ).$$

$${\chi }_{{\mathcal{E}}_{\pi }}$$

$$=d^{-1}{\sum }_{ij}|i⟩⟨j|\otimes {\mathrm{Tr}}_2\left[Q(|i⟩⟨j|\otimes \pi )\right],$$

$${\chi }_{\pi }:={\chi }_{{\mathcal{E}}_{\pi }}=\Lambda (\pi ).$$

$$C_{\diamond }(\pi )\le d{C}_1(\pi )\le 2d\sqrt{C_F(\pi )},$$

$$C_1(\pi ):={\Vert {\chi }_{\mathcal{E}}-{\chi }_{\pi }\Vert }_1,$$

$$C_F(\pi )=1-F(\pi {)}^2,$$

$$F(\pi ):={\Vert \sqrt{{\chi }_{\mathcal{E}}}\sqrt{{\chi }_{\pi }}\Vert }_1=\mathrm{Tr}\sqrt{\sqrt{{\chi }_{\mathcal{E}}}{\chi }_{\pi }\sqrt{{\chi }_{\mathcal{E}}}}.$$

$$C_R(\pi ):=\min \left\{S({\chi }_{\mathcal{E}}||{\chi }_{\pi }),S({\chi }_{\pi }||{\chi }_{\mathcal{E}})\right\},$$

$$\mathrm{Minimize}2{\Vert {\mathrm{Tr}}_{2}Z\Vert }_{\infty },$$

$$\text{Subject to }Z\ge 0\text{and }Z\ge d({\chi }_{\mathcal{E}}-\Lambda (\pi )).$$

$$\partial C(\pi )=\{Z:C(\sigma )-C(\pi )\ge \text{Tr}[Z(\sigma -\pi )],\forall \sigma \in \mathcal{S}\},$$

$$\nabla C_1(\pi )$$

$$\nabla C_F(\pi )$$

$$\nabla F(\pi )$$

$$\begin{array}[c]{l}1)~{}\mathrm{Select~{}an~{}operator~{}}g_{i}{~{}\mathrm{from~{}}}\partial C(\pi_{i}),\\ 2)~{}\pi_{i+1}=\mathcal{P}_{\mathcal{S}}\left(\pi_{i}-\alpha_{i}g_{i}\right),\end{array}$$

$$\begin{array}[c]{l}1)~{}\text{{Find the smallest eigenvalue}}~{}\left|\sigma_{i}\right\rangle~{}\text{{of}}~{}\nabla C(\pi_{i}),\\ 2)~{}\pi_{i+1}=\frac{i}{i+2}\pi_{i}+\frac{2}{i+2}\left|\sigma_{i}\right\rangle\left\langle\sigma_{i}\right|.\end{array}$$

$$\nabla C_F(\pi )=-{\Lambda }^{\ast }\left[|{\chi }_U⟩⟨{\chi }_U|\right],$$

$${\mathcal{E}}_{\pi }^{\mathrm{tele}}(\rho )=\sum _i{U}_i^B⟨{\Phi }_i^{SA}|{\rho }^S\otimes {\pi }^{AB}|{\Phi }_i^{SA}⟩U_i^{B†}.$$

$${\Lambda }_{\text{tele}}(\pi )=d^{-2}\sum _i\left(U_i^{\ast }\otimes U_i\right)\pi {\left(U_i^{\ast }\otimes U_i\right)}^{†},$$

$${\mathcal{P}}_{\pi }(\rho )$$

$$=\sum _{i=1}^N\underset{{\mathbf{A}\bar{\mathbf{B}}}_{i}C}{\mathrm{Tr}}{\left[{\Pi }_i({\pi }_{\mathbf{AB}}\otimes {\rho }_C)\right]}_{B_i\to B_{\mathrm{out}}},$$

$${\Pi }_i={\tilde{\Pi }}_i+\frac{1}{N}\left(\text{openone}-\sum _k{\tilde{\Pi }}_k\right),$$

$${\tilde{\Pi }}_i$$

$${\sigma }_{\mathbf{A}C}$$

$${\pi }_{\mathbf{AB}}={\text{openone}}_A\otimes {\mathcal{E}}_B^{\otimes N}\left[\underset{k=1}{\overset{N}{⨂}}{\Phi }_{A_k{B}_k}\right]=\underset{k=1}{\overset{N}{⨂}}{\chi }_{\mathcal{E}}^{A_k{B}_k}.$$

$$C_{\diamond }^N=\Vert \mathcal{E}\circ \mathcal{I}-\mathcal{E}\circ {\mathcal{I}}_N\Vert \le \Vert \mathcal{I}-{\mathcal{I}}_N{\Vert }_{\diamond }\le 2d(d-1)N^{-1}.$$

$$\Lambda (\pi )={\sum }_{ij}{K}_{ij}\pi K_{ij}^{†},K_{ij}:=⟨e_j^i|\sqrt{{\Pi }_i}\otimes {\text{openone}}_{BD}|{\Phi }_{DC}⟩.$$

$$\mathcal{C}=\{\pi :\pi ={\sum }_k{p}_k{\rho }_{AB}^{k\otimes N},{\mathrm{Tr}}_B({\rho }_{AB}^k)=d^{-1}\text{openone}\}.$$

$${\hat{U}}_j=\exp \left(i{H}_0\otimes |0{⟩}_j{}_j⟨0|+i{H}_1\otimes |1{⟩}_j{}_j⟨1|\right).$$

$$\Lambda (\pi )={\mathrm{Tr}}_{\mathbf{R}}\left[{\hat{U}}_{A\mathbf{R}}\left({\Phi }_{BA}\otimes {\pi }_{\mathbf{R}}\right){\hat{U}}_{A\mathbf{R}}^{†}\right],$$

Full text

Introduction

Back in 1997 a seminal work by Nielsen and Chuang Nielsen and Chuang (1997) proposed a quantum version of the programmable gate array that has become a fundamental model for quantum computation Nielsen and Chuang (2000). This is a quantum processor where a fixed quantum operation is applied to an input state together with a program state. The aim of the program state is to induce the processor to apply some target quantum gate or channel Watrous (2018) to the input state. Such a desired feature of quantum programmability comes with a cost: The model cannot be universal, unless the program state is allowed to have an infinite dimension, i.e., infinite qubits Nielsen and Chuang (1997); Knill et al. (2001). Even though this limitation has been known for many years, there is still no exact characterization on how well a finite-dimensional programmable quantum processor can generate or simulate an arbitrary quantum channel. Also there is no literature on how to find the corresponding optimal program state or even to show that this state can indeed be found by some optimization procedure. Here we show the solutions to these long-standing open problems.

Here we show that the optimization of programmable quantum computers is a convex problem for which the solution can always be found by means of classical semidefinite programming (SDP) and classical gradient-based methods that are commonly employed for machine learning applications. Machine learning (ML) methods have found wide applicability across many disciplines  Bishop (2006), and we are currently witnessing the development of new hybrid areas of investigation where ML methods are interconnected with quantum information theory, such as quantum-enhanced machine learning Wittek (2014); Biamonte et al. (2017); Dunjko and Briegel (2018); Schuld et al. (2015); Ciliberto et al. (2018) (e.g., quantum neural networks, quantum annealing etc.), protocols of quantum-inspired machine learning (e.g., for recommendation systems Tang (2018a) or component analysis and supervised clustering Tang (2018b)) and classical learning methods applied to quantum computers, as explored here in this manuscript.

In our work, we quantify the error between an arbitrary target channel and its programmable simulation in terms of the diamond distance Kitaev et al. (2002); Watrous (2018) and other suitable cost functions, including the trace distance and the quantum fidelity. For all the considered cost functions, we are able to show that the minimization of the simulation error is a convex optimization problem in the space of the program states. This already solves an outstanding problem which affects various models of quantum computers (e.g., variational quantum circuits) where the optimization over classical parameters is non-convex and therefore not guaranteed to converge to a global optimum. By contrast, because our problem is proven to be convex, we can use SDP to minimize the diamond distance and always find the optimal program state for the simulation of a target channel, therefore optimizing the programmable quantum processor. Similarly, we may find suboptimal solutions by minimizing the trace distance or the quantum fidelity by means of gradient-based techniques adapted from the ML literature, such as the projected subgradient method Boyd et al. (2003) and the conjugate gradient method Jaggi (2011, 2013). We note indeed that the minimization of the $\ell_{1}$-norm, mathematically related to the quantum trace distance, is widely employed in many ML tasks Duchi et al. (2008); Liu et al. (2013), so many of those techniques can be adapted for learning program states.

With these general results in our hands, we first discuss the optimal learning of arbitrary unitaries with a generic programmable quantum processor. Then, we consider specific designs of the processor, from a shallow scheme based on the teleportation protocol, to higher-depth designs based on port-based teleportation (PBT) Ishizaka and Hiroshima (2008, 2009); Ishizaka (2015) and parametric quantum circuits (PQCs) Lloyd (1996), introducing a suitable convex reformulation of the latter. In the various cases, we benchmark the processors for the simulation of basic unitary gates (qubit rotations) and various basic channels, including the amplitude damping channel which is known to be the most difficult to simulate Pirandola et al. (2017, 2018a). For the deeper designs, we find that the optimal program states do not correspond to the Choi matrices of the target channels, which is rather counter-intuitive and unexpected.

RESULTS

We first present our main theoretical results on how to train the program state of programmable quantum processors, either via convex optimization or first-order gradient based algorithms. We then apply our general methods to study the learning of arbitrary unitaries, and the simulation of different channels via processors built either from quantum teleportation and its generalization, or from parametric quantum circuits.

Programmable quantum computing

Let us consider an arbitrary mapping from $d$-dimensional input states into $d^{\prime}$-dimensional output states, where $d^{\prime}\neq d$ in the general case. This is described by a quantum channel $\mathcal{E}$ that may represent the overall action of a quantum computation and does not need to be a unitary transformation. Any channel $\mathcal{E}$ can be simulated by means of a programmable quantum processor, which is modeled in general by a fixed completely positive trace-preserving (CPTP) map $Q$ which is applied to both the input state and a variable program state $\pi$. In this way, the processor transforms the input state by means of an approximate channel $\mathcal{E}_{\pi}$ as

[TABLE]

where $\operatorname*{Tr}_{2}$ is the partial trace over the program state. A fundamental result Nielsen and Chuang (1997) is that there is no fixed quantum “processor” $Q$ that is able to exactly simulate any quantum channel $\mathcal{E}$. In other terms, given $\mathcal{E}$, we cannot find the corresponding program $\pi$ such that $\mathcal{E}\equiv\mathcal{E}_{\pi}$. Yet simulation can be achieved in an approximate sense, where the quality of the simulation may increase for larger program dimension. In general, the open problem is to determine the optimal program state $\tilde{\pi}$ that minimizes the simulation error, that can be quantified by the cost function

[TABLE]

namely the diamond distance Kitaev et al. (2002); Watrous (2018) between the target channel $\mathcal{E}$ and its simulation $\mathcal{E}_{\pi}$. In other words,

[TABLE]

From theory Nielsen and Chuang (1997); Knill et al. (2001) we know that we cannot achieve $C_{\diamond}=0$ for arbitrary $\mathcal{E}$ unless $\pi$ and $Q$ have infinite dimensions. As a result, for any finite-dimensional realistic design of the quantum processor, finding the optimal program state $\tilde{\pi}$ is an open problem. Recall that the diamond distance is defined by $\left\|\mathcal{E}-\mathcal{E}_{\pi}\right\|_{\diamond}:=\max_{\varphi}\left\|\mathcal{I}\otimes\mathcal{E}(\varphi)-\mathcal{I}\otimes\mathcal{E}_{\pi}(\varphi)\right\|_{1}$, where $\mathcal{I}$ is the identity map and $\left\|O\right\|_{1}:=\mathrm{Tr}\sqrt{O^{\dagger}O}$ is the trace norm Nielsen and Chuang (2000).

It is important to note that this problem can be reduced to a simpler one by introducing the channel’s Choi matrix

[TABLE]

where $\Phi:=\left|\Phi\right\rangle\!\left\langle\Phi\right|$ is a $d$-dimensional maximally-entangled state. From this expression, it is clear that the Choi matrix $\chi_{\mathcal{E}_{\pi}}$ is linear in the program state $\pi$. More precisely, the Choi matrix $\chi_{\mathcal{E}_{\pi}}$ at the output of the processor $Q$ can be directly written as a CPTP linear map $\Lambda$ acting on the space of the program states $\pi$, i.e.,

[TABLE]

This map is also depicted in Fig. 1

and fully describes the action of the processor $Q$. Then, using results from Refs. Watrous (2018); Nechita et al. (2018); Fuchs and van de Graaf (1999), we may write

[TABLE]

where

[TABLE]

is the trace distance Nielsen and Chuang (2000) between target and simulated Choi matrices, and

[TABLE]

where $F(\pi)$ is Bures’ fidelity between the two Choi matrices $\chi_{\mathcal{E}}$ and $\chi_{\pi}$, i.e.,

[TABLE]

Another possible upper bound can be written using the quantum Pinsker’s inequality Pinsker (1964); Carlen and Lieb (2012). In fact, we may write $C_{1}(\pi)\leq(2\ln\sqrt{2})\sqrt{C_{R}(\pi)}$, where

[TABLE]

and $S(\rho||\sigma):=\mathrm{Tr}[\rho(\log_{2}\rho-\log_{2}\sigma)]$ is the quantum relative entropy between $\rho$ and $\sigma$. In Supplementary Note 1.3 we also introduce a cost function $C_{p}(\pi)$ based on the Schatten p-norm.

Convex optimization

One of the main problems in the optimization of reconfigurable quantum chips is that the relevant cost functions are not convex in the set of classical parameters. This problem is completely solved here thanks to the fact that the optimization of a programmable quantum processor is done with respect to a quantum state. In fact, in the methods section we prove the following

Theorem 1**.**

Consider the simulation of a target quantum channel $\mathcal{E}$ by means of a programmable quantum processor $Q$. The optimization of the cost functions $C_{\diamond}$, $C_{1}$, $C_{F}$, $C_{R}$ or $C_{p}$ is a convex problem in the space of program states $\pi$. In particular, the global minimum $\tilde{\pi}$ for $C_{\diamond}$ can always be found as a local minimum.

This convexity result is generally valid for any cost function which is convex in $\pi$. This is the case for any desired norm, not only the trace norm, but also the Frobenius norm, or any Schatten p-norm. It also applies to the relative entropy. Furthermore, the result can also be extended to any convex parametrization of the program states.

When dealing with convex optimization with respect to positive operators, the standard approach is to map the problem to a form that is solvable via semi-definite programming (SDP) Watrous (2009, 2013a). Since the optimal program is the one minimizing the cost function, it is important to write the computation of the cost function itself as a minimization. For the case of the diamond distance, this can be achieved by using the dual formulation Watrous (2009). More precisely, consider the linear map $\Omega_{\pi}:=\mathcal{E}-\mathcal{E}_{\pi}$ with Choi matrix $\chi_{\Omega_{\pi}}=\chi_{\mathcal{E}}-\chi_{\pi}=\chi_{\mathcal{E}}-\Lambda({\pi})$, and the spectral norm $\left\|O\right\|_{\infty}:=\max\{\left\|Ou\right\|:u\in\mathbb{C}^{d},\left\|u\right\|\leq 1\}$, which is the maximum eigenvalue of $\sqrt{O^{\dagger}O}$. Then, by the strong duality of the diamond norm, $C_{\diamond}(\pi)=\left\|\Omega_{\pi}\right\|_{\diamond}$ is given by the SDP Watrous (2013b)

[TABLE]

The importance of the above dual formulation is that the diamond distance is a minimization, rather than a maximization over a set of matrices. In order to find the optimal program $\tilde{\pi}$ we apply the unique minimization of Eq. (11) where $\pi$ is variable and satisfies the additional constraints $\pi\geq 0$ and $\mathrm{Tr}(\pi)=1$.

In the methods section we show that other cost functions such as $C_{1}$ and $C_{F}$ can also be written as SDPs. Correspondingly, the optimal programs $\tilde{\pi}$ can be obtained by numerical SDP solvers. Most numerical packages implement second-order algorithms such as the interior point method Vandenberghe and Boyd (1996). However, second order methods tend to be computationally heavy for large problem sizes Chao (2013); Monteiro (2003), namely when $\tilde{\pi}$ contains many qudits. In the following section we introduce first order methods, that are better suited for larger program states. It is important to remark that there also exist zeroth-order (derivative-free) methods, such as the simultaneous perturbation stochastic approximation method Spall (2000), which was utilized for a quantum problem in Zhuang and Zhang (2019). However, it is known that zeroth-order methods normally have slower convergence times Harrow and Napp (2019) compared to first-order methods.

Gradient based optimization

In machine learning applications, where a large amount of data is commonly available, there have been several works that study the minimization of suitable matrix norms for different purposes Liu et al. (2013); Duchi et al. (2008); Cai et al. (2010); Recht et al. (2010). First-order methods, are preferred for large dimensional problems, as they are less computationally intensive and require less memory. Here we show how to apply first-order (gradient-based) algorithms, which are widely employed in machine learning applications, to find the optimal quantum program.

For this purpose, we need to introduce the subgradient of the cost function $C$ at any point $\pi\in\mathcal{S}$, which is the set

[TABLE]

where $Z$ is Hermitian Nesterov (2013); Coutts et al. (2018). If $C$ is differentiable, then $\partial C(\pi)$ contains a single element: its gradient $\nabla C(\pi)$. We explicitly compute this gradient for an arbitrary programmable quantum processor (1) whose Choi matrix $\chi_{\mathcal{E}_{\pi}}\equiv\chi_{\pi}=\Lambda(\pi)$, can be written as a quantum channel $\Lambda$ that maps a generic program state to the processor’s Choi matrix. This map can be defined by its Kraus decomposition $\Lambda(\pi)=\sum_{k}A_{k}\pi A_{k}^{\dagger}$ for some operators $A_{k}$. In fact, let us call $\Lambda^{\ast}(\rho)=\sum_{k}A_{k}^{\dagger}\rho A_{k}$ the dual map, then in the methods section we prove the following

Theorem 2**.**

Consider an arbitrary quantum channel $\mathcal{E}$ with Choi matrix $\chi_{\mathcal{E}}$ which is simulated by a quantum processor $Q$ with map $\Lambda(\pi)=\chi_{\pi}$ (and dual map $\Lambda^{\ast}$). Then, we may write the following gradients for the trace distance cost $C_{1}(\pi)$ and the infidelity cost $C_{F}(\pi)$

[TABLE]

where $\lambda_{k}$ ($P_{k}$) are the eigenvalues (eigenprojectors) of the Hermitian operator $\chi_{\pi}-\chi_{\mathcal{E}}$. When $C_{1}(\pi)$ or $C_{F}(\pi)$ are not differentiable in $\pi$, then the above expressions provide an element of the subgradient $\partial C(\pi)$.

Once we have the (sub)gradient of the cost function $C$, we can solve the optimization $\min_{\pi\in\mathcal{S}}C(\pi)$ using the projected subgradient method Nesterov (2013); Boyd et al. (2003). Let $\mathcal{P}_{\mathcal{S}}$ be the projection onto the set of program states $\mathcal{S}$, namely $\mathcal{P}_{\mathcal{S}}(X)=\mathrm{argmin}_{\pi\in S}\|X-\pi\|_{2}$, that we show to be computable from the spectral decomposition of any Hermitian $X$ (see Theorem 3 in the methods section). Then, we iteratively apply the steps

[TABLE]

where $i$ is the iteration index, $\alpha_{i}$ is what is called “learning rate”, and Theorem 2 can be employed to find $g_{i}$ at each step. It is simple to show that $\pi_{i}$ converges to the optimal program state $\tilde{\pi}$ in $\mathcal{O}(\epsilon^{-2})$ steps, for any desired precision $\epsilon$ such that $|C(\pi)-C(\tilde{\pi})|\leq\epsilon$. Another approach is the conjugate gradient method Jaggi (2011); Nesterov (2013), sometimes called Frank-Wolfe algorithm. Here, we apply

[TABLE]

When the gradient of $f$ is Lipschitz continuous with constant $L$, the method converges after $\mathcal{O}(L/\epsilon)$ steps Jaggi (2013); Nesterov (2005). To justify the applicability of this method a suitable smoothening of the cost function must be employed 111The downside of the conjugate gradient method is that it necessarily requires a differentiable cost function $C$, with gradient $\nabla C$. Specifically, this may create problems for the trace distance cost $C_{1}$ which is generally non-smooth. A solution to this problem is to define the cost function in terms of the smooth trace distance $C_{\mu}(\pi)=\mathrm{Tr}\left[h_{\mu}\left(\chi_{\pi}-\chi_{\mathcal{E}}\right)\right]$ where $h_{\mu}$ is the so-called Huber penalty function $h_{\mu}(x):=x^{2}/(2\mu)\mathrm{~{}if~{}}|x|<\mu$ and $|x|-\mu/2\mathrm{~{}if~{}}|x|\geq\mu.$ This quantity satisfies $C_{\mu}(\pi)\leq C_{1}(\pi)\leq C_{\mu}(\pi)+\mu d/2$ and is a convex function over program states, with gradient $\nabla C_{\mu}(\pi)=\Lambda^{\ast}[h_{\mu}^{\prime}(\chi_{\pi}-\chi_{\mathcal{E}})]$..

Learning of arbitrary unitaries

One specific application is the simulation of quantum gates or, more generally, unitary transformations Lloyd (1996); Khaneja et al. (2005); Banchi et al. (2016); Innocenti et al. (2018); Mitarai et al. (2018). Here, the infidelity provides the most convenient cost function, as the optimal program can be found analytically. In fact, suppose we use a quantum processor with map $\Lambda$ to simulate a target unitary $U$. Because the Choi matrix of $U$ is pure $|\chi_{U}\rangle\langle\chi_{U}|$, we first note that $F(\pi)^{2}=\langle\chi_{U}|\Lambda(\pi)|\chi_{U}\rangle$ and then we see that Eq. (15) drastically simplifies to $\nabla F(\pi)=\Lambda^{\ast}\left(|\chi_{U}\rangle\langle\chi_{U}|\right)/\sqrt{4F(\pi)^{2}}$. As a result, we find

[TABLE]

where there is no dependence on $\pi$. Therefore, using the conjugate gradient method in Eq. (17), we see that the optimal program state $\tilde{\pi}$ for the infidelity cost function $C_{F}$ is a fixed point of the iteration and is equal to the maximum eigenvector of $\Lambda^{\ast}\left[|\chi_{U}\rangle\langle\chi_{U}|\right]$.

Teleportation processor

Once we have shown how to optimize a generic programmable quantum processor, we discuss some specific designs, over which we will test the optimization procedure. One possible (shallow) design for the quantum processor $Q$ is a generalized teleportation protocol Bennett et al. (1993) over an arbitrary program state $\pi$. In dimension $d$, the protocol involves a basis of $d^{2}$ maximally-entangled states $|\Phi_{i}\rangle$ and a basis $\{U_{i}\}$ of teleportation unitaries such that $\mathrm{Tr}(U_{i}^{\dagger}U_{j})=d\delta_{ij}$ Pirandola et al. (2015). An input $d$-dimensional state $\rho$ and the $A$ part of the program $\pi_{AB}$ are subject to the projector $|\Phi_{i}\rangle\left\langle\Phi_{i}\right|$. The classical outcome $i$ is communicated to the $B$ part of $\pi_{AB}$ where the correction $U_{i}^{-1}$ is applied.

The above procedure defines the teleportation channel $\mathcal{E}_{\pi}$ over $\rho$

[TABLE]

Its Choi matrix can be written as $\chi_{\pi}=\Lambda_{\text{tele}}(\pi)$, where the map of the teleportation processor is equal to

[TABLE]

which is clearly self-dual $\Lambda^{\ast}=\Lambda$. Given a target quantum channel $\mathcal{E}$ which is teleportation-covariant Pirandola et al. (2017, 2018a), namely when $[\pi,U_{i}^{\ast}\otimes U_{i}]=0$, then we know that that its simulation is perfect and the optimal program $\tilde{\pi}$ is the channel’s Choi matrix, i.e., one of the fixed points of the map $\Lambda_{\text{tele}}$. For a general channel, the optimal program $\tilde{\pi}$ can be approximated by using the cost functions in our Theorem 2 with $\Lambda$ being given in Eq. (20), or directly found by optimizing $C_{\diamond}(\pi)$.

Port-based teleportation

A deeper design is provided by a PBT processor, whose overall protocol is illustrated in Fig. 2. Here we consider a more general formulation of the original PBT protocol Ishizaka and Hiroshima (2008, 2009) where the resource entangled pairs are replaced by an arbitrary program state $\pi$. In a PBT processor, each party has $N$ systems (or ‘ports’), $\mathbf{A}=\{A_{1},\dots,A_{N}\}$ for Alice and $\mathbf{B}=\{B_{1},\dots,B_{N}\}$ for Bob. These are prepared in a program state $\pi_{\mathbf{AB}}$. To teleport an input state $\rho_{C}$, Alice performs a joint positive operator-value measurement (POVM) $\{\Pi_{i}\}$ Ishizaka and Hiroshima (2008) on system $C$ and the $\mathbf{A}$-ports. She then communicates the outcome $i$ to Bob, who discards all ports except $B_{i}$ which is the output $B_{\text{out}}$. The resulting PBT channel $\mathcal{P}_{\pi}:\mathcal{H}_{C}\mapsto\mathcal{H}_{B_{\mathrm{out}}}$ is then

[TABLE]

where $\mathbf{\bar{B}}_{i}=\mathbf{B}\backslash B_{i}=\{B_{k}:k\neq i\}$.

In the standard PBT protocol Ishizaka and Hiroshima (2008, 2009), the program state is fixed as $\pi_{\mathbf{AB}}=\bigotimes_{k=1}^{N}\Phi_{A_{k}B_{k}}$, where $|\Phi_{A_{k}B_{k}}\rangle$ are Bell states, and the following POVM is used

[TABLE]

where

[TABLE]

and $\sigma^{-1/2}$ is an operator defined only on the support of $\sigma$. The PBT protocol is formulated for $N\geq 2$ ports. However, we also include here the trivial case for $N=1$, corresponding to the process where Alice’s input is traced out and the output is the reduced state of Bob’s port, i.e., a maximally mixed state. In the limit $N\rightarrow\infty$, the standard PBT protocol approximates an identity channel $\mathcal{P}_{\pi}(\rho)\approx\rho$, with fidelity Ishizaka and Hiroshima (2008); Ishizaka (2015) $F_{\pi}=1-\mathcal{O}\left(\frac{1}{N}\right)$, so perfect simulation is possible only in the limit $N\rightarrow\infty$. Since the standard PBT-protocol provides an approximation to the identity channel, we call it $\mathcal{I}_{N}$.

From the PBT-simulation of the identity channel it is possible to approximate any general channel $\mathcal{E}$ by noting that $\mathcal{E}$ can be written as a composition $\mathcal{E}\circ\mathcal{I}$, where $\mathcal{I}$ is the identity channel. This is done by replacing the identity channel $\mathcal{I}$ with its PBT simulation $\mathcal{I}_{N}$, and then applying $\mathcal{E}$ to $B_{i}$. However, since Bob does not perform any post-processing on his systems $\mathbf{B}$, aside from discarding all ports $B_{k}$ with $k\neq i$, he can also apply first the channel $\mathcal{E}^{\otimes N}$ to all his ports and then discard all the ports $B_{k}$ with $k\neq i$. In doing so, he changes the program state to

[TABLE]

In other terms, any channel $\mathcal{E}$ can be PBT-approximated by $N$ copies of its Choi matrix $\chi_{\mathcal{E}}$ as program state. Since PBT-simulation can be decomposed as $\mathcal{E}_{\pi}=\mathcal{E}\circ\mathcal{I}_{N}$, the error $C^{N}_{\diamond}=\|\mathcal{E}-\mathcal{E}_{\pi}\|_{\diamond}$ in simulating the channel $\mathcal{E}\equiv\mathcal{E}\circ\mathcal{I}$ satisfies

[TABLE]

where we used the data processing inequality and an upper bound from  Pirandola et al. (2019). While the channel’s Choi matrix assures that $C^{N}_{\diamond}\rightarrow 0$ for large $N$, for any finite $N$ it does not represent the optimal program state. In general, for any finite $N$, finding the optimal program state $\pi_{\mathbf{AB}}$ simulating a channel $\mathcal{E}$ with PBT is an open problem, and no explicit solutions or procedures are known.

We employ our convex optimization procedures to find the optimal program state. This can be done either exactly by minimizing the diamond distance cost function $C_{\diamond}$ via SDP, or approximately, by determining the optimal program state via the minimization of the trace distance cost function $C_{1}$ via either SDP or the gradient-based techniques discussed above. For this second approach, we need to derive the map $\Lambda$ of the PBT processor, between the program state $\pi$ to output Choi matrix as in Eq. (5). To compute the Choi matrix and CP-map $\Lambda$, we consider an input maximally-entangled state $|\Phi_{DC}\rangle$ and a basis $|e_{j}^{i}\rangle$ of $\mathbf{A}\{\mathbf{B}\backslash B_{i}\}C$. Then, by using Eq. 21 and the definition $\Lambda(\pi)=\chi_{\mathcal{P}_{\pi}}=\openone_{D}\otimes\mathcal{P}_{\pi}[\Phi_{DC}]$ we find the map $\Lambda_{\mathbf{AB}\rightarrow DB_{\text{out}}}$ of a PBT processor

[TABLE]

Note that a general program state for PBT consists of $2N$ qudits, and hence the parameter space has exponential size $d^{4N}$. However, because the PBT protocol is symmetric under permutation of port labels, we show in Supplementary Note 6 that one can exploit this symmetry and reduce the number of free parameters to the binomial coeffficient $\binom{N+d^{4}-1}{d^{4}-1}$, which is polynomial in the number of ports $N$. Despite this exponential reduction, the scaling in the number of parameters still represents a practical limiting factor, even for qubits for which $\mathcal{O}\left(N^{15}\right)$. A sub-optimal strategy consists in reducing the space of program states to a convex set that we call the “Choi space” $\mathcal{C}$. Consider an arbitrary probability distribution $\{p_{k}\}$ and then define

[TABLE]

One can show (see Supplementary Note 6) that a global minimum in $\mathcal{C}$ is a global minimum in the extremal (non-convex) subspace for $p_{k}=\delta_{k,1}$ consisting of tensor-products of Choi matrices $\rho_{AB}^{\otimes N}$. Among these states, there is the $N$-copy Choi matrix of the target channel $\chi_{\mathcal{E}}^{\otimes N}=[\mathcal{I}\otimes\mathcal{E}(|\Phi\rangle\left\langle\Phi\right|)]^{\otimes N}$ which is not necessarily the optimal program, as we show below.

Parametric quantum circuits

Another deep design of quantum processor is based on PQCs Lloyd (1995, 1996). A PQC is a sequence of unitary matrices $U(t)=U_{N}(t_{N})\dots U_{2}(t_{2})U_{1}(t_{1})$, where $U_{j}(t_{j})=\exp(it_{j}H_{j})$ for some Hamiltonian $H_{j}$ and time interval $t_{j}$. The problem with PQCs is that the cost functions in the classical parameters Khaneja et al. (2005) are not convex, so that numerical algorithms are not guaranteed to converge to the global optimum. Here we fix this issue by introducing a convex formulation of PQCs where classical parameters are replaced by a quantum program. This results in a programmable PQC processor which is optimizable by our methods.

The universality of PQCs can be employed for universal channel simulation. Indeed, thanks to Stinespring’s dilation theorem, any channel can be written as a unitary evolution on a bigger space, $\mathcal{E}(\rho_{A})=\mathrm{Tr}_{R_{0}}[U(\rho_{A}\otimes\theta_{0})U^{\dagger}],$ where the system is paired to an extra register $R_{0}$ and $\theta_{0}$ belongs to $R_{0}$. In the Stinespring representation $U$ acts on system $A$ and register $R_{0}$. In Ref. Lloyd (1995) it has been shown that sequences of two unitaries, $U_{0}$ and $U_{1}$, are almost universal for simulation, i.e., any target unitary $U$ can be approximated as $U\approx\cdots U_{1}^{m_{4}}U_{0}^{m_{3}}U_{1}^{m_{2}}U_{0}^{m_{1}}$ for some integers $m_{j}$. Under suitable conditions, it takes $\mathcal{O}(d^{2}\epsilon^{-d})$ steps to approximate $U$ up to precision $\epsilon$. The choice between $U_{0}$ and $U_{1}$ is done by measuring a classical bit. We may introduce a quantum version, where the two different unitaries $U_{0}=e^{iH_{0}}$ or $U_{1}=e^{iH_{1}}$ are chosen depending on the state of qubit $R_{j}$. This results in the conditional gate

[TABLE]

Channel simulation is then obtained by replacing the unitary evolution $U$ in the Stinespring dilation via its simulation.

The result is illustrated in Fig. 3, where the program state $\pi$ is defined over $\mathbf{R}=(R_{0},\dots,R_{N})$ and each $\hat{H}_{j}$ acts on the input system $A$ and two ancillary qubits $R_{0}$ and $R_{j}$. Following the universality construction of Ref. Lloyd (1995) we show in the Supplementary Note 3.4 that the channel shown in Fig. 3 provides a universal processor. Moreover, the channel $\Lambda$ that maps any program $\pi$ to the processor’s Choi matrix is obtained as

[TABLE]

where $\hat{U}_{A\mathbf{R}}=\openone_{B}\otimes\prod_{j=1}^{N}\hat{U_{j}}_{A,R_{0},R_{j}}$, from which we can identify the optimal program $|\tilde{\pi}\rangle$ via our methods.

PQCs are not inherently monotonic. A deeper (higher $N$) design may simulate a given channel worse than a more shallow design. We can design a modified PQC that is monotonic by design, which we designate a “monotonic PQC”, by replacing the qubits in our program state with qutrits, and modifying Eq. 29 to read

[TABLE]

where $\mathbf{0}$ is a zero operator, so that gate $j$ enacts the identity channel if program qutrit $j$ is in the state $|2\rangle\langle 2|$. Then, if it were the case that a PQC with $N$ program qubits could simulate a given channel better than one with $N+m$, a monotonic PQC with $N+m$ qutrits in the program state could perform at least as well as the PQC with $N$ program qubits by setting the first $m$ qutrits to $|2\rangle\langle 2|$. This processor design is both universal and monotonic. More precisely, let $C(\mathrm{PQC}_{N})$ denote the value of a cost function $C$ for simulating a channel $\mathcal{E}$ with an $N$-gate PQC, using the optimal program state, and let $C(\mathrm{mPQC}_{N})$ denote the value of $C$ for simulating $\mathcal{E}$ with an $N$-gate monotonic PQC, again using the optimal program state. We are then guaranteed that

[TABLE]

Processor benchmarking

In order to show the performance of the various architectures, we consider the simulation of an amplitude damping channel with probability $p$. The reason is because this is the most difficult channel to simulate, with a perfect simulation only known for infinite dimension, e.g., using continuous-variable quantum operations Pirandola et al. (2017). In Figs. 4-5 we compare teleportation-based, PBT, PQC and “monotonic PQC” (mPQC) programmable processors whose program states have been optimized according to the cost functions $C_{\diamond}$ and $C_{1}$. For the PBT processor the trace distance cost $C_{1}$ is remarkably close to $C_{\diamond}$ and allows us to easily explore high depths. Note that the optimal program states differ from the naive choice of the Choi matrix of the target channel. Note too that PQC processors display non-monotonic behaviour when simulating amplitude damping channels, meaning that shallow PQC processors (e.g., for $N=4$) may perform better than deeper processors 222For the PQC processor we use the universal Hamiltonians $H_{0}=\sqrt{2}(X\otimes Y-Y\otimes X)$ and $H_{1}=(\sqrt{2}Z{+}\sqrt{3}Y{+}\sqrt{5}X)\otimes(Y{+}\sqrt{2}Z)$, where $X$, $Y$, and $Z$ are Pauli operators.. Monotonic PQC processors guarantee that deeper designs always perform at least as well as any shallower design. In Fig. 5 perfect simulation is achievable at specific values of $p$ because of our choice of the universal gates $U_{0}$ and $U_{1}$. More details are provided in the Supplementary Note 3.

Many other numerical simulations are performed in the Supplementary Note 3 where we study the convergence rate in learning a unitary operation, the exact simulation of Pauli channels, and approximate simulation of both dephasing and amplitude damping channels. In particular, we study the performance of the approximate solution when optimizing over larger, but easier to compute, cost functions such as the trace distance or the infidelity.

DISCUSSION

In this work we have considered a general finite-dimensional model of a programmable quantum processor, which is a fundamental scheme for quantum computing and also a primitive tool for other areas of quantum information. By introducing suitable cost functions, based on the diamond distance, trace distance and quantum fidelity, we have shown how to characterize the optimal performance of this processor in the simulation of an arbitrary quantum gate or channel. In fact, we have shown that the minimization of these cost functions is a convex optimization problem that can always be solved.

In particular, by minimizing the diamond distance via SDP, we can always determine the optimal program state for the simulation of an arbitrary channel. Alternatively, we may minimize the simpler but larger cost functions in terms of trace distance and quantum fidelity via gradient-based methods adapted from ML, so as to provide a very good approximation of the optimal program state. This other approach can also provide closed analytical solutions, as is the case for the simulation of arbitrary unitaries, for which the minimization of the fidelity cost function corresponds to computing an eigenvector.

We have then applied our results to various designs of programmable quantum processor, from a shallow teleportation-based scheme to deeper and asymptotically-universal designs that are based on PBT and PQCs. We have explicitly benchmarked the performances of these quantum processors by considering the simulation of unitary gates, depolarizing and amplitude damping channels, showing that the optimal program states may differ from the naive choice based on the Choi matrix of the target channel. Moreover, our results can be applied also for universal quantum measurements D’Ariano and Perinotti (2005).

A potential application of our work may be the development of “programmable” model of cloud-based quantum computation, where a client has an input state to be processed by an online quantum server which is equipped with a programmable quantum processor. The client classically informs the server about what type of computation it needs (e.g., some specified quantum algorithm) and the server generates an optimal program state which closely approximates the overall quantum channel to be applied to the input. The server then accepts the input from the client, processes it, and returns the output together with the value of a cost function quantifying how close the computation was with respect to the client’s request.

Our results may also be useful in areas beyond quantum computing, wherever channel simulation is a basic problem. For instance, this is the case when we investigate the ultimate limits of quantum communications Pirandola et al. (2018a), design optimal Hamiltonians for one-way quantum repeaters, and for all those areas of quantum sensing, hypothesis testing and metrology which are based on quantum channel simulations Pirandola et al. (2018b). Indeed the study of adaptive protocols of quantum channel discrimination (or estimation) is notoriously difficult, and their optimal performance is not completely understood. Nonetheless, these protocols can be analyzed by using simulation techniques Pirandola et al. (2019, 2018b) where the channel, encoding the unknown parameter, is replaced by an approximate simulating channel, and its parameter is mapped into the label of a program state (therefore reducing the problem from channel to state discrimination/estimation). In this regard, our theory provides the optimal solution to this basic problem, by determining the best simulating channel and the corresponding program state.

METHODS

Convexity proofs

In this section we provide a proof of Theorem 1, namely we show that the minimization of the main cost functions $C_{\diamond}$, $C_{1}$ and $C_{F}$ is a convex optimization problem in the space of the program states $\pi$. This means that we can find the optimal program state $\tilde{\pi}$ by minimizing $C_{\diamond}$ or, alternatively, sub-optimal program states can be found by minimizing either $C_{1}$ or $C_{F}$. For the sake of generality, we prove the result for all of the cost functions discussed in the previous sections. We restate Theorem 1 below for completeness:

Theorem**.**

*The minimization of the generic cost function $C=C_{\diamond}$, $C_{1}$, $C_{F}$, $C_{R}$ or $C_{p}$ for any $p>1$ is a convex optimization problem in the space of program states. *

Proof.  Let us start to show the result for the diamond distance $C_{\diamond}$. In this case, we can write the following

[TABLE]

where we use $(1)$ the linearity of $\mathcal{E}$, $(2)$  the triangle inequality and $(3)$ the property $\|xA\|_{1}=|x|\|A\|_{1}$, valid for any operator $A$ and coefficient $x$.

For any Schatten $p$-norm $C_{p}$ with $p\geq 1$, we may prove convexity following a similar reasoning. Since for any combination $\bar{\pi}:=p_{0}\pi_{0}+p_{1}\pi_{1}$, with $p_{0}+p_{1}=1$, we have $\Lambda(\bar{\pi})=p_{0}\Lambda(\pi_{0})+p_{1}\Lambda(\pi_{1})$, then by exploiting the triangle inequality, and the property $\|xA\|_{p}=|x|\|A\|_{p}$, we can show that

[TABLE]

To show the convexity of $C_{F}$, defined in Eq. (8), we note that the fidelity function $F(\rho,\sigma)$ satisfies the following concavity relation Uhlmann (1976)

[TABLE]

Due to the linearity of $\chi_{\pi}=\Lambda(\pi)$, the fidelity in Eq. (9) satisfies $F_{\bar{\pi}}^{2}\geq\sum_{k}p_{k}F_{\pi_{k}}^{2}$ for $\bar{\pi}:=\sum_{k}p_{k}\pi_{k}$. Accordingly, we get the following convexity result

[TABLE]

For the cost function $C_{R}$, the result comes from the linearity of $\Lambda(\pi)$ and the joint convexity of the relative entropy. In fact, for $\bar{\pi}:=p_{0}\pi_{0}+p_{1}\pi_{1}$, we may write

[TABLE]

with a symmetric proof for $S[\chi_{\mathcal{E}}||\Lambda(\bar{\pi})]$. This implies the convexity of $C_{R}(\pi)$ in Eq. (10). $\blacksquare$

Convex classical parametrizations

The result of the theorem 1 can certainly be extended to any convex parametrization of program states. For instance, assume that $\pi=\pi(\boldsymbol{\lambda})$, where $\boldsymbol{\lambda}=\{\lambda_{i}\}$ is a probability distribution. This means that, for $0\leq p\leq 1$ and any two parametrizations, $\boldsymbol{\lambda}$ and $\boldsymbol{\lambda}^{\prime}$, we may write

[TABLE]

Then the problem remains convex in $\boldsymbol{\lambda}$ and we may therefore find the global minimum in these parameters. It is clear that this global minimum $\boldsymbol{\tilde{\lambda}}$ identifies a program state $\pi(\boldsymbol{\tilde{\lambda}})$ which is not generally the optimal state $\tilde{\pi}$ in the entire program space $\mathcal{S}$, even though the solution may be a convenient solution for experimental applications.

Note that a possible classical parametrization consists of using classical program states, of the form

[TABLE]

where $\{\left|\varphi_{i}\right\rangle\}$ is an orthonormal basis in the program space. Convex combinations of probability distributions therefore define a convex set of classical program states

[TABLE]

Optimizing over this specific subspace corresponds to optimizing the programmable quantum processor over classical programs. It is clear that global minima in $\mathcal{S}_{\text{class}}$ and $\mathcal{S}$ are expected to be very different. For instance, $\mathcal{S}_{\text{class}}$ cannot certainly include Choi matrices which are usually very good quantum programs.

Gradient-based optimization

As discussed in the main text, the SDP formulation allows the use of powerful and accurate numerical methods, such as the interior point method. However, these algorithms are not suitable for high dimensional problems, due to their higher computational and memory requirements. Therefore, an alternative approach (useful for larger program states) consists of the optimization of the larger but easier-to-compute cost function $C=C_{1}$ (trace distance) or $C_{F}$ (infidelity), for which we can use first order methods. Indeed, according to Theorem 1, all of the proposed cost functions $C:\mathcal{S}\rightarrow\mathbb{R}$ are convex over the program space $\mathcal{S}$ and, therefore, we can solve the optimization $\min_{\pi\in\mathcal{S}}C(\pi)$ by using gradient-based algorithms.

Gradient-based convex optimization is at the heart of many popular ML techniques such as online learning in a high-dimensional feature space Duchi et al. (2008), missing value estimation problems Liu et al. (2013), text classification, image ranking, and optical character recognition Duchi et al. (2011), to name a few. In all of the above applications, “learning” corresponds to the following minimization problem: $\min_{x\in\mathcal{S}}f(x)$, where $f(x)$ is a convex function and $\mathcal{S}$ is a convex set. Quantum learning falls into this category, as the space of program states is convex due to the linearity of quantum mechanics and the fact that cost functions are typically convex in this space (see Theorem 1). Gradient-based approaches are among the most applied methods for convex optimization of non-linear, possibly non-smooth functions Nesterov (2013).

When the cost function is not differentiable we cannot formally define its gradient. Nonetheless, we can always define the subgradient $\partial C$ of $C$ as in Eq. (12), which in principle contains many points. When $C$ is not only convex but also differentiable, then $\partial C(\pi)=\{\nabla C(\pi)\}$, i.e. the subgradient contains a single element, the gradient $\nabla C$, that can be obtained via the Fréchet derivative of $C$ (for more details see Supplementary Note 4). When $C$ is not differentiable, the gradient still provides an element of the subgradient that can be used in the minimization algorithm.

In order to compute the gradient $\nabla C$, it is convenient to consider the Kraus decomposition of the processor map $\Lambda$. Let us write

[TABLE]

with Kraus operators $A_{k}$. We then define the dual map $\Lambda^{\ast}$ of the processor as the one (generally non-trace-preserving) which is given by the following decomposition

[TABLE]

With these definitions in hands, we can now prove Theorem 2, which we rewrite here for convenience.

Theorem**.**

Suppose we use a quantum processor $Q$ with map $\Lambda(\pi)=\chi_{\pi}$ in order to approximate the Choi matrix $\chi_{\mathcal{E}}$ of an arbitrary channel $\mathcal{E}$. Then, the gradients of the trace distance $C_{1}(\pi)$ and the infidelity $C_{F}(\pi)$ are given by the following analytical formulas

[TABLE]

where $\lambda_{k}$ ($P_{k}$) are the eigenvalues (eigenprojectors) of the Hermitian operator $\chi_{\pi}-\chi_{\mathcal{E}}$. When $C_{1}(\pi)$ or $C_{F}(\pi)$ are not differentiable at $\pi$, then the above expressions provide an element of the subgradient $\partial C(\pi)$.

Proof.  We prove the above theorem assuming that the functions are differentiable for program $\pi$. For non-differentiable points, the only difference is that the above analytical expressions are not unique and provide only one of the possibly infinite elements of the subgradient. Further details of this mathematical proof are given in Supplementary Note 4. Following matrix differentiation, for any function $f(A)=\operatorname*{Tr}[g(A)]$ of a matrix $A$, we may write

[TABLE]

and the gradient is $\nabla f(A)=g^{\prime}(A)$. Both the trace-distance and fidelity cost functions can be written in this form. To find the explicit gradient of the fidelity function, we first note that, by linearity, we may write

[TABLE]

and therefore the following expansion

[TABLE]

From this equation and differential calculations of the fidelity (see Supplementary Note 4.2 for details), we find

[TABLE]

where $dF=F(\pi+\delta\pi)-F(\pi)$. Then, using the cyclic property of the trace, we get

[TABLE]

Exploiting this expression in Eq. (46) we get the gradient $\nabla F(\pi)$ as in Eq. (45). The other Eq. (44) simply follows from applying the definition in Eq. (8).

For the trace distance, let us write the eigenvalue decomposition

[TABLE]

Then using the linearity of Eq. (47), the definition of a processor map of Eq. (5) and differential calculations of the trace distance (see Supplementary Note 4.3 for details), we can write

[TABLE]

From the definition of the gradient in Eq. (46), we finally get

[TABLE]

which leads to the result in Eq. (43). $\blacksquare$

The above results in Eqs. (44) and (43) can be used together with the projected subgradient method Boyd et al. (2003) or conjugate gradient algorithm Jaggi (2011, 2013) to iteratively find the optimal program state in the minimization of $\min_{\pi\in\mathcal{S}}C(\pi)$ for $C=C_{1}$ or $C_{F}$. In the following sections we present two algorithms, the projected subgradient method and the conjugate gradient method, and show how they can be adapted to our problem.

Projected subgradient methods have the advantage of simplicity and the ability to optimize non-smooth functions, but can be slower, with a convergence rate $\mathcal{O}\left(\epsilon^{-2}\right)$ for a desired accuracy $\epsilon$. Conjugate gradient methods Jaggi (2011, 2013) have a faster convergence rate $\mathcal{O}\left(\epsilon^{-1}\right)$, provided that the cost function is smooth. This convergence rate can be improved even further to $\mathcal{O}\left(\epsilon^{-1/2}\right)$ for strongly convex functions Garber and Hazan (2015) or using Nesterov’s accelerated gradient method Nesterov (2005). The technical difficulty in the adaptation of these methods for learning program states comes because the latter is a constrained optimization problem, namely at each iteration step the optimal program must be a proper quantum state, and the cost functions coming from quantum information theory are, generally, non-smooth.

Projected subgradient method

Given the space $\mathcal{S}$ of program states, let us define the projection $\mathcal{P}_{\mathcal{S}}$ onto $\mathcal{S}$ as

[TABLE]

where argmin is the argument of the minimum, namely the closest state $\pi\in\mathcal{S}$ to the operator $X$. Then, a first order algorithm to solve $\min_{\pi\in\mathcal{S}}C(\pi)$ is to apply the projected subgradient method Nesterov (2013); Boyd et al. (2003), which iteratively applies the iteration (16), which we rewrite below for convenience

[TABLE]

where $i$ is the iteration index and $\alpha_{i}$ a learning rate.

The above algorithm differs from standard gradient methods in two aspects: i) the update rule is based on the subgradient, which is defined even for non-smooth functions; ii) the operator $\pi_{i}-\alpha_{i}g_{i}$ is generally not a quantum state, so the algorithm fixes this issue by projecting that operator back to the closest quantum state, via Eq. (54). The algorithm converges to the optimal solution $\pi_{\ast}$ (approximating the optimal program $\tilde{\pi}$) as Boyd et al. (2003)

[TABLE]

where $e_{1}=\|\pi_{1}-\pi_{\ast}\|_{2}^{2}$ is the initial error (in Frobenius norm) and $G$ is such that $\|g\|_{2}^{2}\leq G$ for any $g\in\partial C$. Popular choices for the learning rate that assure convergence are $\alpha_{k}\propto 1/\sqrt{k}$ and $\alpha_{k}=a/(b+k)$ for some $a,b>0$.

In general, the projection step is the major drawback, which often limits the applicability of the projected subgradient method to practical problems. Indeed, projections like Eq. (54) require another full optimization at each iteration that might be computationally intensive. Nonetheless, we show in the following theorem that this issue does not occur in learning quantum states, because the resulting optimization can be solved analytically.

Theorem 3**.**

Let $X$ be a Hermitian operator in a $d$-dimensional Hilbert space with spectral decomposition $X=UxU^{\dagger}$, where the eigenvalues $x_{j}$ are ordered in decreasing order. Then $\mathcal{P}_{\mathcal{S}}(X)$ of Eq. (54) is given by

[TABLE]

where $\theta=\frac{1}{s}\sum_{j=1}^{s}\left(x_{j}-1\right)$ and

[TABLE]

Proof.  Any quantum (program) state can be written in the diagonal form $\pi=V\lambda V^{\dagger}$ where $V$ is a unitary matrix, and $\lambda$ is the vector of eigenvalues in decreasing order, with $\lambda_{j}\geq 0$ and $\sum_{j}\lambda_{j}=1$. To find the optimal state, it is required to find both the optimal unitary $V$ and the optimal eigenvalues $\lambda$ with the above property, i.e.,

[TABLE]

For any unitarily-invariant norm, the following inequality holds (Bhatia, 2013, Eq. IV.64)

[TABLE]

with equality when $U=V$, where $X=UxU^{\dagger}$ is a spectral decomposition of $X$ such that the $x_{j}$’s are in decreasing order. This shows that the optimal unitary in Eq. (59) is the diagonalization matrix of the operator $X$. The eigenvalues of any density operator form a probability simplex. The optimal eigenvalues $\lambda$ are then obtained thanks to Algorithm 1 from Ref. Duchi et al. (2008). $\blacksquare$

In the following section we present an alternative algorithm with faster convergence rates, but stronger requirements on the function to be optimized.

Conjugate gradient method

The conjugate gradient method Jaggi (2011); Nesterov (2013), sometimes called the Frank-Wolfe algorithm, has been developed to provide a better convergence speed and to avoid the projection step at each iteration. Although the latter can be explicitly computed for quantum states (thanks to our Theorem 3), having a faster convergence rate is important, especially with higher dimensional Hilbert spaces. The downside of this method is that it necessarily requires a differentiable cost function $C$, with gradient $\nabla C$.

In its standard form, the conjugate gradient method to approximate the solution of $\operatorname*{argmin}_{\pi\in\mathcal{S}}C(\pi)$ is defined by the following iterative rule

[TABLE]

The first step in the above iteration rule is solved by finding the smallest eigenvector $|\sigma\rangle$ of $\nabla C(\pi_{i})$. Indeed, since $\pi$ is an operator and $C(\pi)$ a scalar, the gradient $\nabla C$ is an operator with the same dimension as $\pi$. Therefore, for learning quantum programs we find the iteration (17), that we rewrite below for convenience

[TABLE]

When the gradient of $C$ is Lipschitz continuous with constant $L$, the conjugate gradient method converges after $\mathcal{O}(L/\epsilon)$ steps Jaggi (2013); Nesterov (2005). The following iteration with adaptive learning rate $\alpha_{i}$ has even faster convergence rates, provided that $C$ is strongly convex Garber and Hazan (2015):

[TABLE]

where the constant $\beta_{C}$ and norm $\|\cdot\|_{C}$ depend on $C$ Garber and Hazan (2015).

In spite of the faster convergence rate, conjugate gradient methods require smooth cost functions (so that the gradient $\nabla C$ is well defined at every point). However, cost functions based on trace distance (7) are not smooth. For instance, the trace distance in one-dimensional spaces reduces to the absolute value function $|x|$ that is non-analytic at $x=0$. When some eigenvalues are close to zero, conjugate gradient methods may display unexpected behaviors, though we have numerically observed that convergence is always obtained with a careful choice of the learning rate. In the next section we show how to formally justify the applicability of the conjugate gradient method, following Nesterov’s smoothing prescription Nesterov (2005).

Smoothing: smooth trace distance

The conjugate gradient method converges to the global optimum after $\mathcal{O}\left(\frac{L}{\epsilon}\right)$ steps, provided that the gradient of $C$ is $L$-Lipschitz continuous Nesterov (2005). However, the constant $L$ can diverge for non-smooth functions like the trace distance (7) so the convergence of the algorithm cannot be formally stated, although it may still be observed in numerical simulations. To solidify the convergence proof (see also Supplementary Note 5.2), we introduce a smooth approximation to the trace distance. This is defined by the following cost function that is differentiable at every point

[TABLE]

where $\lambda_{j}$ are the eigenvalues of $\chi_{\pi}-\chi_{\mathcal{E}}$ and $h_{\mu}$ is the so-called Huber penalty function

[TABLE]

The previous definition of the trace distance, $C_{1}$ in Eq. (7), is recovered for $\mu\rightarrow 0$ and, for any non-zero $\mu$, the $C_{\mu}$ bounds $C_{1}$ as follows

[TABLE]

where $d$ is the dimension of the program state $\pi$. In Supplementary Note 5.2 we then prove the following result

Theorem 4**.**

The smooth cost function $C_{\mu}(\pi)$ is a convex function over program states and its gradient is given by

[TABLE]

where $h_{\mu}^{\prime}$ is the derivative of $h_{\mu}$. Moreover, the gradient is $L$-Lipschitz continuous with

[TABLE]

where $d$ is the dimension of the program state.

Being Lipschitz continuous, the conjugate gradient algorithm and its variants Nesterov (2005); Garber and Hazan (2015) converge up to an accuracy $\epsilon$ after $\mathcal{O}(L/\epsilon)$ steps. In some applications, it is desirable to analyze the convergence in trace distance in the limit of large program states, namely for $d\rightarrow\infty$. The parameter $\mu$ can be chosen such that the smooth trace distance converges to the trace distance, namely $C_{\mu}\rightarrow C_{1}$ for $d\rightarrow\infty$. Indeed, given the inequality (66), a possibility is to set $\mu=\mathcal{O}(d^{-(1+\eta)})$ for some $\eta>0$ so that, from Eq. (68), the convergence to the trace norm is achieved after $\mathcal{O}(d^{2+\eta})$ steps.

**Acknowledgements. ** L.B. acknowledges support by the program “Rita Levi Montalcini” for young researchers. S.P. and J.P. acknowledge support by the EPSRC via the ‘UK Quantum Communications Hub’ (Grants EP/M013472/1 and EP/T001011/1) and S.P. acknowledges support by the European Union via the project ‘Continuous Variable Quantum Communications’ (CiViQ, no 820466).

1 More on programmable simulation

As discussed in the main text, the task we are interested in is the simulation of a channel $\mathcal{E}$ using a programmable quantum processor Nielsen and Chuang (1997) that we simply call a “quantum processor” (see Fig. 6). This is represented by a completely positive trace-preserving (CPTP) universal map $Q$ as in Eq. (1). Our goal is to find the program state $\pi$ according to (3), namely the state for which the simulation $\mathcal{E}_{\pi}$ is the closest to $\mathcal{E}$. The most appropriate definition of “closeness” between two quantum channels is via the diamond norm $C_{\diamond}(\pi):=\left\|\mathcal{E}-\mathcal{E}_{\pi}\right\|_{\diamond}\leq 2$. Recall that the diamond distance is defined by the following maximization

[TABLE]

where $\left\|O\right\|_{1}:=\mathrm{Tr}\sqrt{O^{\dagger}O}$ is the trace norm Watrous (2018). Because the trace norm is convex over mixed states, one may reduce the maximization in Eq. (S1) to bipartite pure states $\varphi=\left|\varphi\right\rangle\left\langle\varphi\right|$. In general, we therefore need to consider a min-max optimization, i.e., find $\tilde{\pi}$ and (pure) $\tilde{\varphi}$ such that

[TABLE]

1.1 Programmable quantum measurements

Programmable quantum measurements D’Ariano and Perinotti (2005) represent a particular instance of the general channel approximation problem, where the channel approximates a measurement device. Consider a POVM $\{\Pi_{j}\}$ and a programmable POVM $\{\Pi_{j}^{\pi}\}$ with program $\pi$ which is obtained by performing a fixed joint POVM $\{Q_{j}\}$ on both the input state $\rho$ and the program $\pi$. We may associate to these measurement devices two measurement channels $\mathcal{E}(\rho)=\sum_{j}p_{j}|j\rangle\langle j|$, where $p_{j}=\operatorname*{Tr}[\rho\Pi_{j}]$ is the probability of getting the outcome $j$, and similarly $\mathcal{E}_{\pi}(\rho)=\sum_{j}p_{j}^{\pi}|j\rangle\langle j|$, where

[TABLE]

is the probability of getting the outcome $j$ using the programmable measurement $\Pi_{j}^{\pi}$. The above discussion shows that programmable quantum measurements represent a particular instance of the general case that we consider for arbitrary channels, so we may optimize the program state $\pi$ with the techniques presented in our paper.

Note that we may also consider a different cost function, namely the worst-case distance between the two probability distributions given by

[TABLE]

It was shown in D’Ariano and Perinotti (2005) that there exist (fixed) universal POVMs $\{Q_{j}\}$ such that $\Pi_{j}^{\pi}$ approximates any arbitrary measurement $\Pi_{j}$ with an optimal program $\pi$. The error in the approximation, as quantified by $C_{M}$, decreases with the dimension of the program. Via the measurement channels defined above, it is easy to see that

[TABLE]

so that our theory includes the results of D’Ariano and Perinotti (2005) as a special case.

1.2 SDP minimization

We show that some of the convex cost functions that we have introduced can be explicitly evaluated via semi-definite programming (SDP). This allows us to use standard SDP algorithms for finding the optimal program.

We first fix the program state $\pi$ and show how for fixed $\pi$ it is possible to compute $C_{\diamond}(\pi)$ via semidefinite programming. Let us introduce the linear map $\Omega_{\pi}:=\mathcal{E}-\mathcal{E}_{\pi}$ with corresponding Choi matrix

[TABLE]

Thanks to the property of strong duality of the diamond norm, for any program $\pi$ we can compute the cost function $C_{\diamond}(\pi)=\|\Omega_{\pi}\|_{\diamond}$ via the following SDP Watrous (2009)

[TABLE]

where $M_{0}\geq 0$ and $M_{1}\geq 0$ in $\mathbb{C}^{d\times d^{\prime}}$, and the spectral norm $\|O\|_{\infty}$ equals the maximum singular value of $O$.

Moreover, because $\chi_{\Omega_{\pi}}$ is Hermitian, the above SDP can be simplified into

[TABLE]

Not only does this procedure compute $C_{\diamond}(\pi)$, but it also provides the upper bound $C_{\diamond}(\pi)\leq d\left\|\mathrm{Tr}_{2}\left|\chi_{\mathcal{E}}-\chi_{\pi}\right|\right\|_{\infty}$ Nechita et al. (2018). In fact, it is sufficient to choose $Z=d~{}\chi_{\Omega_{\pi}}^{+}$, where $\chi^{+}=(\chi+|\chi|)/2$ is the positive part of $\chi$. Using $\mathrm{Tr}_{2}\chi_{\Omega_{\pi}}=0$, we may write $\mathrm{Tr}_{2}Z\leq d\mathrm{Tr}_{2}\chi_{\Omega_{\pi}}^{+}=\frac{d}{2}\mathrm{Tr}_{2}|\chi_{\Omega_{\pi}}|$.

The SDP form in Eq. (S8) is particularly convenient for finding the optimal program. In fact, suppose now that $\pi$ is not fixed but we want to optimize on this state too, so as to compute the optimal program state $\tilde{\pi}$ such that $\tilde{\pi}={\rm argmin}_{\pi\in\mathcal{S}}C_{\diamond}(\pi)$. The problem is therefore mapped into the following unique minimization

[TABLE]

Unlike the min-max optimization of Eq. (S2), the above SDP is much simpler as it contains a unique minimization. Therefore, this algorithm can be used to optimize the performance of any programmable quantum processor.

Using a similar argument, and exploiting known convex programming formulations for the trace norm and the fidelity cost Watrous (2013a); Recht et al. (2010), we can also compute the optimal program states $\tilde{\pi}$ as $\tilde{\pi}_{1}={\rm argmin}_{\pi\in S}C_{1}(\pi)$ and $\tilde{\pi}_{F}={\rm argmax}_{\pi\in S}F(\pi)$. Note indeed that, clearly, the fidelity has to be maximized, rather than minimized. The optimal program $\tilde{\pi}_{1}$ and its associated cost $C_{1}(\tilde{\pi}_{1})$ can be computed with the following minimization

[TABLE]

The optimal program $\tilde{\pi}_{F}$ and its associated cost $F(\tilde{\pi}_{F})$ can be computed with the following maximization

[TABLE]

1.3 Cost functions and their relative dependence

For completeness, we report here some inequalities between the different cost functions introduced in the paper. For example, we may connect the minimization of the diamond distance $C_{\diamond}(\pi)$ to the minimization of the trace distance $C_{1}(\pi)$ via the sandwich relation Watrous (2018)

[TABLE]

While the lower bound is immediate from the definition of Eq. (S1), the upper bound can be proven using the following equivalent form of the diamond distance

[TABLE]

where the optimization is carried out over the density matrices $\rho_{0}$ and $\rho_{1}$ (Watrous, 2013b, Theorem 3.1). In fact, consider the Frobenius norm $\|A\|_{2}:=\sqrt{\operatorname*{Tr}[A^{\dagger}A]}$ and the spectral norm

[TABLE]

which satisfy the following properties Watrous (2018)

[TABLE]

Then, from Eqs. (S13), (S15) and (S16), one gets

[TABLE]

Thanks to Eq. (S12), we may avoid the maximization step in the definition of the diamond distance and simplify the original problem to approximating the Choi matrix $\chi_{\mathcal{E}}$ of the channel by varying the program state $\pi$. This is a process of learning Choi matrices as depicted in Fig. 7. Because the simpler cost function $C_{1}(\pi)$ is an upper bound, its minimization generally provides a sub-optimal solution for the program state.

Finally we may consider other cost functions in terms of any Shatten p-norm $C_{p}(\pi):=\|\chi_{\mathcal{E}}-\chi_{\pi}\|_{p}$, even though this option provides lower bounds instead of upper bounds for the trace distance. Recall that, given an operator $O$ and a real number $p\geq 1$, we may define its Schatten p-norm as Watrous (2018)

[TABLE]

where $|O|=\sqrt{O^{\dagger}O}$. For any $1\leq p\leq q\leq\infty$, one has the monotony $\left\|O\right\|_{p}\geq\left\|O\right\|_{q}$, so that $\left\|O\right\|_{\infty}\leq\ldots\leq\left\|O\right\|_{1}$. An important property is duality. For each pair of operators $A$ and $B$, and each pair of parameters $p,q\in[1,\infty]$ such that $p^{-1}+q^{-1}=1$, we may write Watrous (2018)

[TABLE]

where $\left\langle B,A\right\rangle=\mathrm{Tr}(B^{\dagger}A)$ is the Hilbert-Schmidt product, and the second inequality follows since we can arbitrarily change the sign of $B$.

2 Convergence in learning arbitrary unitaries

The simulation of quantum gates or, more generally, unitary transformations is crucial for quantum computing applications Lloyd (1996) so ML techniques have been developed for this purpose Banchi et al. (2016); Innocenti et al. (2018); Mitarai et al. (2018). In the main text we have shown that, for learning arbitrary unitaries, the fidelity cost function provides a convenient choice for which the optimal program can be found analytically. Indeed, the optimal program is always a pure state and is given by the eigenstate of $\Lambda^{\ast}\left[|\chi_{U}\rangle\langle\chi_{U}|\right]$ with the maximum eigenvalue. Here we consider the convergence of the Frank-Wolfe iteration Eq. (17) towards that state.

Let $\pi_{1}$ be the initial guess for the program state. After $k$ iterations of Eq. (17), we find the following approximation to the optimal program state

[TABLE]

where $\frac{2}{k+k^{2}}=\prod_{j=1}^{k-1}\frac{j}{j+2}$. The above equation shows that $\pi_{k}\rightarrow\tilde{\pi}_{F}$ for $k\rightarrow\infty$, with error in trace distance

[TABLE]

This example shows that the convergence rate $\mathcal{O}(\epsilon^{-1})$ of the conjugate method provides a worst case instance that can be beaten in some applications with some suitable cost functions. From Eq. (S21) we see that $\epsilon=k^{-2}$ for learning arbitrary unitaries via the minimization of $C_{F}$, meaning that convergence is obtained with the faster rate $\mathcal{O}(\epsilon^{-1/2})$. On the other hand, there are no obvious simplifications for the optimization of the trace distance, since the latter still requires the diagonalization of Eq. (51). For the trace distance, or its smooth version, only numerical approaches are feasible.

3 Applications

3.1 Learning a unitary with the teleportation processor

Here we consider the following example. Assume that the target channel is a unitary $U$, so that its Choi matrix is $\chi_{U}:=|\chi_{U}\rangle\langle\chi_{U}|$ with $|\chi_{U}\rangle=\openone\otimes U|\Phi\rangle$ and where $|\Phi\rangle$ is maximally entangled. By using Eq. (20), the fact that $\Lambda_{\rm tele}=\Lambda_{\rm tele}^{*}$ and $U^{\ast}\otimes\openone|\Phi\rangle=\openone\otimes U^{\dagger}|\Phi\rangle$, we may write the dual processor map

[TABLE]

where $V_{i}^{U}=U_{i}UU_{i}^{\dagger}$. The maximum eigenvector of $\Lambda_{\text{tele}}^{\ast}[|\chi_{U}\rangle\langle\chi_{U}|]$ represents the optimal program state $\tilde{\pi}_{F}$ for simulating the unitary $U$ via the teleportation processor (according to the fidelity cost function). In some cases, the solution is immediate. For instance, this happens when $V_{i}^{U}\propto U$ is independent of $i$. This is the case when $U$ is a teleportation unitary, because it satisfies the Weyl-Heisenberg algebra Pirandola et al. (2017). For a teleportation unitary $U$, we simply have

[TABLE]

so that the unique optimal program is $\tilde{\pi}_{F}=|\chi_{U}\rangle\langle\chi_{U}|$.

In Fig. 8 we show the convergence of the projected subgradient algorithm using the teleportation processor and target unitaries $R(\theta)=e^{i\theta X}$, for different values of $\theta$. When $\theta$ is a multiple of $\pi/2$, then the above unitary is teleportation covariant and the Frank-Wolfe algorithm converges to zero trace distance. For other values of $\theta$ perfect simulation is impossible, and we notice that the algorithm converges to a non zero value of the trace distance (7). For comparison, in Fig. 8 we also plot the value of the fidelity upper bound $\sqrt{1-F[\Lambda(\tilde{\pi}_{F}),\chi_{\mathcal{E}}]^{2}}$, where $\tilde{\pi}_{F}$ is the optimal program that maximizes the fidelity of Eq. (9), namely the eigenvector of Eq. (S22) with the maximum eigenvalue. We note that for $\theta=\pi/2^{\ell}$, the trace distance decreases for larger $\theta$. The limit case $\ell\rightarrow\infty$ is perfectly simulable as $R(0)$ is teleportation covariant.

3.2 Pauli channel simulation

Pauli channels are defined as Nielsen and Chuang (2000)

[TABLE]

where $U_{i}$ are generalized Pauli operators and $p_{i}$ are some probabilities. For $d=2$ the Pauli operators are the four Pauli matrices $I,X,Y,Z$ and in any dimension they form the Weyl-Heisenberg group Nielsen and Chuang (2000). These operators are exactly the teleportation unitaries $U_{j}$ defined in the previous section. The Choi matrix $\chi_{\mathcal{P}}$ of a Pauli channel $\mathcal{P}$ is diagonal in the Bell basis, i.e., we have

[TABLE]

where $\Phi_{i}=\openone\otimes U_{i}|\Phi\rangle$ and $|\Phi\rangle=\sum_{j=1}^{d}|jj\rangle/\sqrt{d}$.

We now consider the simulation of a Pauli channel with the teleportation quantum processor introduced in the previous section. Let

[TABLE]

be an arbitrary program state expanded in the Bell basis. For any program state, the Choi matrix of the teleportation-simulated channel is given by Eq. (20). Using standard properties of the Pauli matrices, we find

[TABLE]

namely a generic state is transformed into a Bell diagonal state. Therefore, the cost function

[TABLE]

can be minimized analytically for any Pauli channel by choosing $\pi_{ij}=p_{i}\delta_{ij}$. With this choice we find $C_{1}^{\mathrm{Pauli}}=0$, meaning that the simulation is perfect.

From theory Bowen and Bose (2001); Cope et al. (2017); Bennett et al. (1996) we know that only Pauli channels can be perfectly simulated in this way. No matter how much more general we make the states $\pi$, it is proven Bowen and Bose (2001); Cope et al. (2017) that these are the only channels we can perfectly simulate. This is true even if we apply the Pauli corrections in a probabilistic way, i.e., we assume a classical channel from the Bell outcomes to the corresponding label of the Pauli correction operator Cope et al. (2017).

3.3 PBT: Numerical examples

We first consider the simulation of an amplitude damping channel $\mathcal{E}_{\mathrm{AD}}(\rho)=\sum_{i}K_{i}^{\mathrm{AD}}\rho K_{i}^{\mathrm{AD}\dagger}$, which is defined by the Kraus operators

[TABLE]

In Fig. 9 we study the performance of the PBT simulation of the amplitude damping channel $\mathcal{E}_{\mathrm{AD}}$ for different choices of $p$. For $p=0$ the amplitude damping channel is equal to the identity channel, while for $p=1$ it is a “reset” channel sending all states to $|0\rangle$. We compare the simulation error with program states $\pi$ either made by products of the channel’s Choi matrix $\chi_{\mathcal{E}_{\mathrm{AD}}}^{\otimes N}$ as in Eq. (25) or obtained from the minimization of the trace distance cost function of Eq. (7) with the projected subgradient iteration in Eq. (16). Alternative methods, such as the conjugated gradient algorithm, perform similarly for this channel. We observe that, surprisingly, the optimal program $\tilde{\pi}_{1}$ obtained by minimizing the trace distance $C_{1}$ is always better than the natural choice $\chi_{\mathcal{E}_{\mathrm{AD}}}^{\otimes N}$.

In Fig. 10 we study the PBT simulation of the amplitude damping channel by considering the subset of program states $\pi=\chi^{\otimes N}$ which is made of tensor products of the $4\times 4$ generic Choi matrices $\chi$ (satisfying $\mathrm{Tr}_{2}\chi=\openone/2$). As discussed in Sec. 6, this is equivalent to optimizing over the Choi set $\mathcal{C}_{N}$ and it practically reduces to the convex optimization of the channel $\tilde{\Lambda}$ over the generic single-copy Choi matrix $\chi$. Moreover, $\tilde{\Lambda}$ itself can be simplified, as shown in Appendix 6, so that all of the operations depend polynomially on the number $N$ of ports. This allows us to numerically explore much larger values of $N$, even for the minimization of $C_{\diamond}$. In Fig. 10 the dotted lines correspond to the value of $C_{\diamond}$ when the program $\pi=\chi_{\mathcal{E}_{\mathrm{AD}}}^{\otimes N}$ is employed, where $\chi_{\mathcal{E}_{\mathrm{AD}}}$ is the channel’s Choi matrix. As Fig. 10 shows, the cost $C_{\diamond}$ may be significantly smaller with an optimal $\chi$, thus showing that the optimal program may be different from the channel’s Choi matrix, especially when $p$ is far from the two boundaries $p=0$ and $p=1$.

As an other example, we consider the simulation of the depolarizing channel defined by

[TABLE]

It was shown in Ishizaka and Hiroshima (2009); Pirandola et al. (2019) that PBT generates a depolarizing channel, whose depolarizing probability $p_{\rm th}$ depends on $N$. This implies that a quantum processor based on PBT can perfectly simulate a depolarizing channel when $p\geq p_{\rm th}$. In Fig. 11 we study the performance of PBT simulation of the depolarizing channel in terms of $p$. Perfect simulation is possible for $p\geq p_{\rm th}\simeq 0.71$ for $N=2$ and $p\geq p_{\rm th}=0.5$ for $N=3$, where the explicit values of $p_{\rm th}$ are taken from Ref. Pirandola et al. (2019). For $p=0$ the depolarizing channel is equal to the identity channel, while for $p=1$ it sends all states to the maximally mixed state. Again we compare the simulation error with program states either composed of copies of the channel’s Choi matrix $\chi_{\mathcal{E}_{\mathrm{dep}}}^{\otimes N}$ or obtained from the minimization of $C_{1}$ with the conjugate gradient method of Eq. (17), which performs significantly better than the projected subgradient for this channel. Also for the depolarizing channel we observe that, for any finite $N$, we obtain a lower error by optimizing over the program states instead of the naive choice $\chi_{\mathcal{E}_{\mathrm{dep}}}^{\otimes N}$.

Finally, in Fig. 12 we study the PBT simulation of a unitary gate $U_{\theta}=e^{i\theta X}$ for different values of $\theta$. Unlike the previous non-unitary channels, in Fig. 12 we observe a flat error where different unitaries have the same simulation error of the identity channel $\theta=0$. This is expected because both the trace distance and the diamond distance are invariant under unitary transformations. In general, we have the following.

Proposition 5**.**

Given a unitary $\mathcal{U}(\rho)=U\rho U^{\dagger}$ and its PBT simulation $\mathcal{U}_{\pi}$ with program $\pi$ we may write

[TABLE]

where $\mathcal{I}_{\pi}$ is the PBT simulation of the identity channel.

Proof.  In fact, we simultaneously prove

[TABLE]

where (1) comes from the fact that $||\mathcal{U}-\mathcal{U}_{\pi}||_{\diamond}=||\mathcal{U}^{-1}\mathcal{U}-\mathcal{U}^{-1}\mathcal{U}_{\pi}||_{\diamond}=||\mathcal{I}-\mathcal{U}^{-1}\mathcal{U}_{\pi}||_{\diamond}$ and $\mathcal{U}^{-1}\mathcal{U}_{\pi}$ is a possible PBT simulation of the identity $\mathcal{I}$ with program state $\mathcal{I}\otimes(\mathcal{U}^{-1})^{\otimes N}(\pi)$ once $\mathcal{U}^{-1}$ is swapped with the filtering of the ports; then (2) comes from the fact that the composition $\mathcal{U}\circ\mathcal{I}_{\pi}$ is a possible simulation of the unitary $\mathcal{U}$ with program state $\mathcal{I}\otimes\mathcal{U}^{\otimes N}(\pi)$ and we have the inequality$\ ||\mathcal{U}\circ\mathcal{I}-\mathcal{U}\circ\mathcal{I}_{\pi}||_{\diamond}\leq||\mathcal{I}-\mathcal{I}_{\pi}||_{\diamond}$. $\blacksquare$

The scaling of $||\mathcal{I}-\mathcal{I}_{\pi}||_{\diamond}$ for different values of $N$ is plotted in Fig. 13 where numerical values are obtained from SDP, while the upper bound is given by Eq. (26).

3.4 Parametric quantum circuits

We now study another design of universal quantum processor that can simulate any target quantum channel in the asymptotic limit of an arbitrarily large program state. This is based on a suitable reformulation of the PQCs, which are known to simulate any quantum computation with a limited set of quantum gates Lloyd (1995, 1996).

A PQC is composed of a sequence of unitary matrices $U_{j}(\theta_{j})$, each depending on a classical parameter $\theta$. The resulting unitary operation is then

[TABLE]

A convenient choice is $U_{j}(\theta_{j})=\exp(i\theta_{j}H_{j})$, where each elementary gate corresponds to a Schrödinger evolution with Hamiltonian $H_{j}$ for a certain time interval $\theta_{j}$. For certain choices of $H_{j}$ and suitably large $N$ the above circuit is universal Lloyd (1996), namely any unitary can be obtained with $U(\theta)$ and a suitable choice of $\theta_{j}$. The optimal parameters can be found via numerical algorithms Khaneja et al. (2005), e.g. by minimizing the cost function $C(\theta)=|\operatorname*{Tr}[U_{\rm target}^{\dagger}U(\theta)]|$. However, the above cost function is not convex, so the numerical algorithms are not guaranteed to converge to the global optimum.

As a first step, we show that the task of learning the optimal parameters in a PQC can be transformed into a convex optimization problem by using a quantum program. This allows us to use SDP and gradient-based ML methods to find the global optimum solution.

3.4.1 Convex reformulation

Consider a program state $|\pi\rangle=|\theta_{1},\dots,\theta_{N}\rangle$ composed of $N$ registers $R_{j}$, each in a separable state $|\theta_{j}\rangle$. We can transform the classical parameters in Eq. (S33) into quantum parameters via the conditional gates

[TABLE]

that act non-trivially on systems and registers $R_{j}$. If the parameters $\theta_{j}$ are continuous, then we can replace the sum with an integral.

With the above gates, we define the parametric quantum channel

[TABLE]

whose action on a generic state $|\psi\rangle$ is shown in Fig. 14. For a pure separable program $|\pi\rangle=|\theta_{1},\dots,\theta_{N}\rangle$, we obtain the standard result, i.e.,

[TABLE]

where $U(\theta)$ is defined in Eq. (S33). The parametric quantum processor $Q_{\pi}$ in Eq. (S35) is capable of simulating any parametric quantum channel, but it is more general, as it allows entangled quantum parameters and also parameters in quantum superposition.

An equivalent measurement-based protocol is obtained by performing the trace in Eq. (S36) over the basis $|\theta_{1},\dots,\theta_{N}\rangle$, so that

[TABLE]

where $U(\theta)$ is defined in Eq. (S33). In this alternative yet equivalent formulation, at a certain iteration $j$, the processor measures the qubit register $R_{j}$. Depending on the measurement outcome $\theta_{j}$, the processor then applies a different unitary $U(\theta_{j})$ on the system. However, in this formulation the program state $|\pi\rangle$ is destroyed after each channel use. From Eq. (S37) we note that $Q_{\pi}$ depends on $\pi$ only via the probability distribution $\langle\theta_{1},\dots,\theta_{N}|\pi|\theta_{1},\dots,\theta_{N}\rangle$. As such, any advantage in using quantum states can only come from the capability of quantum systems to model computationally hard probability distributions Boixo et al. (2018).

3.4.2 Universal channel simulation via PQCs

The universality of PQCs can be employed for universal channel simulation, thanks to Stinespring’s dilation theorem:

[TABLE]

where $\theta_{0}$ belongs to $R_{0}$, and $U$ acts on system $A$ and register $R_{0}$. In Ref. Lloyd (1995) it was shown that two quantum gates are universal for quantum computation. Specifically, given $U_{0}=e^{it_{0}H_{0}}$ and $U_{B}=e^{it_{1}H_{1}}$ for fixed times $t_{i}$ and Hamiltonians $H_{j}$, it is possible to write any unitary as

[TABLE]

for some integers $m_{j}$. Under suitable conditions, it was shown that with $M=\sum_{j}m_{j}=\mathcal{O}(d^{2}\epsilon^{-d})$ it is possible to approximate any unitary $U$ with a precision $\epsilon$. More precisely, the conditions are the following

1. i)

The Hamiltonians $H_{0}$ and $H_{1}$ are generators of the full Lie algebra, namely $H_{0}$, $H_{1}$ and their repeated commutators generate all the elements of su(d). 2. ii)

The eigenvalues of $U_{0}$ and $U_{1}$ have phases that are irrationally related to $\pi$.

The decomposition in Eq. (S39) is a particular case of Eq. (S33) where $\theta_{j}$ can only take binary values $\theta_{j}=0,1$. As such we can write the conditional gates of Eq. (S34) as in Eq. 29, which is rewritten below

[TABLE]

for some times $t_{j}$. Channel simulation is then obtained by replacing the unitary evolution $U$ of Eq. (S38) with the approximate form in Eq. (S39) and its simulation in Eq. (S41). The result is illustrated in Fig. 3 and described by the following channel

[TABLE]

where the program state $\pi$ is defined over $\mathbf{R}=(R_{0},\dots,R_{N})$ and each $\hat{H}_{j}$ acts on the input system $A$ and two ancillary qubits $R_{0}$ and $R_{j}$. The decomposition of Eq. (S39) assures that, with the program

[TABLE]

the product of unitaries approximates $U$ in Eq. (S38) with precision $\epsilon$. This is possible in general, provided that the program state has dimension $\mathcal{O}(d^{2}\epsilon^{-d})$. However, the channel (S41) is more general, as it allows both quantum superposition and entanglement.

The processor map $\Lambda$, written in Eq. (30) easily follows from this construction , while the (non-trace-preserving) dual channel may be written as

[TABLE]

This channel requires $2N$ quantum gates at each iteration and can be employed for the calculation of gradients, following Theorem 2. When we are interested in simulating a unitary channel $U$ via the quantum fidelity, then following the results of Section 2, the corresponding optimal program $\tilde{\pi}_{F}$ is simply the eigenvector $\Lambda^{\ast}[|\chi_{U}\rangle\langle\chi_{U}|]$ with the maximum eigenvalue, where $|\chi_{U}\rangle=\openone\otimes U|\Phi\rangle$. Note also that $\Lambda^{\ast}[|\chi_{U}\rangle\langle\chi_{U}|]=Z^{\dagger}Z$ where

[TABLE]

so the optimal program $\tilde{\pi}_{F}$ is the principal component of $Z$. Since there are quantum algorithms for principal component analysis Lloyd et al. (2014), the optimization may be efficiently performed on a quantum computer.

3.4.3 Monotonicity by design

Two unitaries are sufficient for universality, however such a design may not be monotonic as a function of $N$ when simulating a given channel. By adding the identity as a third possible unitary, as in Eq. (31), we get a processor that is monotonic by design. Note that the identity cannot be one of the initial choices of unitaries, due to the condition that the eigenvalues of $U_{0}$ and $U_{1}$ have phases that are irrationally related to $\pi$. Because we want a single program qudit to control the application of one of three unitaries, we use qutrits to control the gates.

Such a design is more powerful than the original design, because it is guaranteed to both be monotonic as a function of $N$, where $N$ is the number of controlled gates, but also to be able to simulate any channel at least as well as the original $M$-gate processor can, for any $M\leq N$. This is because a valid program state for the monotonic processor is $\pi_{M}\otimes|2\rangle\langle 2|^{\otimes(N-M)}$, where $\pi_{M}$ is any program state for the original $M$-gate processor. Such a program state would result in the same cost function that the original processor obtains using $\pi_{M}$.

Due to having a larger program state space, and more unitaries to choose from, the $N$-gate monotonic design can also potentially perform not just as well as but better than any $M$-gate processor using the original design, when simulating many channels. This comes at the cost of higher dimensionality. The number of parameters to optimise over scales with order $\mathcal{O}(3^{2N})$ for the monotonic processor, whilst it scales with order $\mathcal{O}(2^{2N})$ for the original design.

3.5 PQC: Numerical examples

As an example we study the simulation of an amplitude damping channel, with Kraus operators in Eq. (S29). A possible Stinespring dilation for this channel is obtained with $|\theta_{0}\rangle=|0\rangle$ and

[TABLE]

where the Hamiltonian is given by

[TABLE]

with $X$ and $Y$ being Pauli operators. We may construct a PQC simulation by taking

[TABLE]

for some $\alpha$ and taking $U_{1}$ to be a different unitary that makes the pair $U_{0}$, $U_{1}$ universal. Here we may choose $\alpha=\sqrt{2}$ and $U_{1}=e^{iH_{1}}$ with

[TABLE]

Results are shown in Fig. 15. Compared with the similar PBT simulation of Fig. 9, we observe that the PQC simulation (using the non-monotonic design) displays a non-monotonic behavior as a function of $N$. PBT with $N$ pairs requires a register of $2N$ qubits, while PQC requires $N+1$ qubits, namely $N$ qubits for the conditional gates and an extra auxiliary qubit coming from the Stinespring decomposition (see Fig. 3). We observe that, with a comparable yet finite register size, PQC can outperform PBT in simulating the amplitude damping channel.

In Fig. 16 we also study the PQC simulation of the depolarizing channel for different values of $p$. Although the gates $U_{0}$ and $U_{1}$ were chosen with inspiration from the Stinespring decomposition of the amplitude damping channel, those gates are universal and capable of simulating other channels. Indeed, we observe in Fig. 16 that a depolarizing channel is already well simulated with $N=4$ for all values of $p$.

4 Matrix calculus

4.1 Matrix differentiation

For a general overview of these techniques, the reader may consult Ref. Stickel (1987). Thanks to Cauchy’s theorem, a matrix function can be written as

[TABLE]

For the same reason

[TABLE]

Applying a basic rule of matrix differentiation, $d(A^{-1})=-A^{-1}(dA)A^{-1}$ we obtain

[TABLE]

Clearly, $df(A)=f^{\prime}(A)dA$ only when $[A,dA]=0$. In general $df(A)$ is a superoperator that depends on $A$ and is applied to $dA$. The explicit form is easily computed using the eigenvalue decomposition or other techniques Stickel (1987). Note that in some cases the expressions are simple. Indeed, using the cyclic invariance of the trace, we have

[TABLE]

while in general $d$Tr$[Bf(A)]\neq$Tr$[Bf^{\prime}(A)dA]$.

4.2 Differential of the quantum fidelity

The quantum fidelity can be expanded as

[TABLE]

where in the second line we have applied Eq. (S49). Taking the differential with respect to $Y$ and using the cyclic property of the trace we get

[TABLE]

where in (1) we use Eq. (S51) and the cyclic property of the trace; in (2) we use Eq. (S50) with $f(\lambda)=\sqrt{\lambda}$, so $f^{\prime}(\lambda)=\frac{1}{2}\lambda^{-1/2}$; and in (3) we use the cyclic property of the trace. See also Lemma 11 in Coutts et al. (2018).

4.3 Differential of the trace distance

The trace norm for a Hermitian operator $X$ is defined as

[TABLE]

where in the second line we applied Eq. (S49). From the spectral decomposition $X=U\lambda U^{\dagger}$, we find $t(X)=\sum_{j}|\lambda_{j}|$, so the trace distance reduces to the absolute value function for one-dimensional Hilbert spaces. The absolute value function $|\lambda|$ is differentiable at every point, except $\lambda=0$. Therefore, for any $\lambda\neq 0$, the subgradient of the absolute value function is composed of only its gradient, i.e.

[TABLE]

For $\lambda=0$ we can use the definition (12) to write

[TABLE]

which is true iff $-1\leq z\leq 1$. Therefore,

[TABLE]

The sign function in (S56) can be extended to $\lambda=0$ in multiple ways (common choices are $\mathrm{sign}(0)=-1,0,1$). From the above equation, it appears that for any extension of the sign function, provided that $\mathrm{sign}(0)\in[-1,1]$ we may write the general form

[TABLE]

which is true for any value of $\lambda$.

With the same spirit we extend the above argument to any matrix dimension, starting from the case where $X$ is an invertible operator (no zero eigenvalues). Taking the differential with respect to $X$ we find

[TABLE]

where in (1) we use Eq. (S51), the cyclic property of the trace and the identity $dX^{2}=X(dX)+(dX)X$; in (2) we use Eq. (S50) with $f(\lambda)=\sqrt{\lambda}$, so $f^{\prime}(\lambda)=\frac{1}{2}\lambda^{-1/2}$; and in (3) we use the cyclic property of the trace and the commutation of $X$ and $\sqrt{X^{2}}$. Let

[TABLE]

be the eigenvalue decomposition of $X$ with eigenvalues $\lambda_{k}$ and eigenprojectors $P_{k}$. For non-zero eigenvalues we may write

[TABLE]

and accordingly

[TABLE]

Therefore, for invertible operators we may write

[TABLE]

We now consider the general case where some eigenvalues of $X$ may be zero. We do this by generalizing Eq. (S59), namely we show that even if $\partial t(X)$ may contain multiple elements, it is always true that $\nabla t\in\partial t$, provided that $-\openone\leq\mathrm{sign}(X)\leq\openone$. Following (12) we may write, for fixed $X$ and arbitrary $Y$,

[TABLE]

where in (1) we use the property $\|X\|_{1}=\operatorname*{Tr}[{\rm sign}(X)X]$ and in (2) we use the assumption $-\openone\leq\mathrm{sign}(X)\leq\openone$. From the definition of the subgradient (12), the above equation shows that $\mathrm{sign}(X)\in\partial t(X)$, so we may always use $\nabla t(X)=\mathrm{sign}(X)$ in the projected subgradient algorithm (16).

5 Smoothing techniques

5.1 Stochastic smoothing

The conjugate gradient algorithm converges after $\mathcal{O}(c/\epsilon)$ steps Jaggi (2011, 2013), where $\epsilon$ is the desired precision and $c$ is a curvature constant that depends on the function. However, it is known that $c$ could diverge for non-smooth functions. This is the case for the trace norm, as shown in Example 0.1 in Ravi et al. (2017).

A general solution, valid for arbitrary functions, is stochastic smoothing Yousefian et al. (2012). In this approach the non-smooth function $C(\pi)$ is replaced by the average

[TABLE]

where $\sigma$ is such that $\|\sigma\|_{\infty}\leq 1$. If $|C(x)-C(y)|\leq M\|x-y\|_{\infty}$, then

[TABLE]

so that $C_{\eta}(\pi)$ provides a good approximation for $C(\pi)$. Moreover, $C_{\eta}$ is differentiable at any point, so we may apply the conjugate gradient algorithm. A modified conjugate gradient algorithm with adaptive stochastic approximation was presented in Ref. Lan (2013). At each iteration $k$ the algorithm reads

[TABLE]

where $g$ denotes any element of the subgradient $\partial C$. The above algorithm converges after $\mathcal{O}(\epsilon^{2})$ iterations. Since Eqs. (45) and (43) provide an element of the subgradient, the above algorithm can be applied to both fidelity and trace distance. However, this algorithm requires $k$ evaluations of the subgradient to perform the averages, so it may be impractical when the number of iterations get larger. In the following we study an alternative that does not require any average.

5.2 Nesterov’s smoothing

An alternative smoothing scheme is based on Nesterov’s dual formulation Nesterov (2005). Suppose that the non-smooth objective function $f$ admits a dual representation as follows

[TABLE]

for some inner product $\langle\cdot,\cdot\rangle$. Nesterov’s approximation consists of adding a strongly convex function $d$ to the dual

[TABLE]

The resulting $\mu$-approximation is smooth and satisfies

[TABLE]

The trace norm admits the dual representation Watrous (2018)

[TABLE]

where $\langle Y,X\rangle$ is the Hilbert Schmidt product. This can be regularized with any strongly convex function $d$. A convenient choice Liu et al. (2013) that enables an analytic solution is $d(X)=\frac{1}{2}\|X\|_{2}^{2}:=\frac{1}{2}\langle X,X\rangle$, so

[TABLE]

This function is smooth and its gradient is given by Liu et al. (2013)

[TABLE]

where $X=U\Sigma V^{\dagger}$ is the singular value decomposition of $X$ and $\Sigma_{\mu}$ is a diagonal matrix with diagonal entries $(\Sigma_{\mu})_{i}=\min\{\Sigma_{i}/\mu,1\}$. Plugging this into Eq. (S71) we get

[TABLE]

For a diagonalizable matrix $X$ with spectral decomposition $X=U\lambda U^{\dagger}$, the singular value decomposition is obtained with $\Sigma=|\lambda|$ and $V=U\mathrm{sign}(\lambda)$. Inserting these expressions in (S72) we find

[TABLE]

where $h_{\mu}$ is the so called Huber penalty function

[TABLE]

The gradient $\nabla t_{\mu}$ is then $h_{\mu}^{\prime}(X)\equiv Uh^{\prime}(\lambda)U^{\dagger}$, where

[TABLE]

We then find that, via the smooth trace norm $t_{\mu}$, we can define the smooth trace distance of Eq. (64) that is differentiable at every point

[TABLE]

Thanks to the inequalities in (S69), the smooth trace distance bounds the cost $C_{1}$ as

[TABLE]

where we employed the identity $\sup_{\|Y\|_{\infty}\leq 1}\|Y\|^{2}_{2}\leq d$ to get the upper bound. Moreover, we find the following

Lemma 6**.**

The smooth trace distance, defined in Eq. (64), is a convex function of $\pi$.

Proof.  From the definition and Eq. (S71) we find

[TABLE]

Now for $\bar{\pi}=p\pi_{1}+(1-p)\pi_{2}$ linearity implies $f(\bar{\pi}):=\langle Y,\Lambda(\bar{\pi})-\chi_{\mathcal{E}}\rangle=pf(\pi_{1})+(1-p)f(\pi_{2})$. Therefore

[TABLE]

showing the convexity. $\blacksquare$

Then, using the definitions from Nesterov (2005), the following theorem bounds the growth of the gradient

Theorem 7**.**

*The gradient of the smooth trace norm is Lipschitz continuous with Lipschitz constant *

[TABLE]

In particular, if the gradient is Lipschitz continuous, the smooth trace norm satisfies the following inequality for any state $\pi,\sigma$

[TABLE]

Proof.  Given the linearity of the quantum channel $\Lambda$, we can apply theorem 1 from Nesterov (2005) to find

[TABLE]

Since all eigenvalues of $y$ are less than or equal to 1, we can write $y\leq 1$ and as such

[TABLE]

$\blacksquare$

6 PBT: program state compression

The dimension of the program state grows exponentially with the number of ports $N$ as $d^{2N}$ where $d$ is the dimension of the Hilbert space. However, as also discussed in the original proposal Ishizaka and Hiroshima (2008, 2009) and more recently in Ref. Christandl et al. (2018), the resource state of PBT can be chosen with extra symmetries, so as to reduce the number of free parameters. In particular, we may consider the set of program states that are symmetric under the exchange of ports, i.e., such that rearranging any $A$ modes and the corresponding $B$ modes leaves the program state unchanged.

Let $P_{s}$ be the permutation operator swapping labels 1 to $N$ for the labels in the sequence $s$, which contains all the numbers 1 to $N$ once each in some permuted order. Namely $P_{s}$ exchanges all ports according to the rule $i\mapsto s_{i}$. Since PBT is symmetric under exchange of ports, we may write

[TABLE]

Consider then an arbitrary permutation-symmetric resource state $\pi_{\mathrm{sym}}$ as

[TABLE]

where the sum is over all possible sequences $s$ that define independent permutations and $N!$ is the total number of possible permutations. Clearly $\mathcal{P}_{\pi_{\mathrm{sym}}}=\mathcal{P}_{\pi}$, so any program state gives the same PBT channel as some symmetric program state. It therefore suffices to consider the set of symmetric program states. This is a convex set: any linear combination of symmetric states is a symmetric state.

To construct a basis of the symmetric space, we note that each element of a density matrix is the coefficient of a dyadic (of the form $\left|x\right\rangle\left\langle y\right|$). If permutation of labels maps one dyadic to another, the coefficients must be the same. This allows us to constrain our density matrix using fewer global parameters. For instance, for $d=2$ we can define the 16 parameters $n_{00,00}$, $n_{00,01}$, $n_{00,10}$, etc., corresponding to the number of ports in the dyadic of the form $\left|0_{A}0_{B}\right\rangle\left\langle 0_{A}0_{B}\right|$, $\left|0_{A}0_{B}\right\rangle\left\langle 0_{A}1_{B}\right|$, $\left|0_{A}0_{B}\right\rangle\left\langle 1_{A}0_{B}\right|$, etc. Each element of a symmetric density matrix can then be defined solely in terms of these parameters, i.e., all elements corresponding to dyadics with the same values of these parameters have the same value.

For the general qudit case, in which our program state consists of $N$ ports, each composed of two $d$-dimensional qudits, we can find the number of independent parameters from the number of independent dyadics. Each port in a dyadic can be written as $|a_{A},b_{B}\rangle\langle c_{A},d_{B}|$ where the extra indices $A$ and $B$ describe whether those states are modeling either qudit $A$ or $B$. There are $d^{4}$ different combinations of $\{a,b,c,d\}$, so we can place each qudit into one of $d^{4}$ categories based on these values. If two elements in the density matrix correspond to dyadics with the same number of ports in each category, they must take the same value. Hence, the number of independent coefficients is given by the number of ways of placing $N$ (identical) ports into $d^{4}$ (distinguishable) categories. This is exactly the binomial coefficient

[TABLE]

Consequently, exploiting permutation symmetry of the PBT protocol, we can exponentially reduce the number of parameters for the optimization over program states.

The number of parameters can be reduced even further by considering products of Choi matrices. We may focus indeed on the Choi set, first defined in Eq. (28),

[TABLE]

where each $\chi_{k}=\chi_{AB}^{k}$ is a generic Choi matrix, therefore satisfying $\mathrm{Tr}_{B}\chi_{k}=d^{-1}\openone$, and $p_{k}$ form a probability distribution. Clearly $\mathcal{C}$ is a convex set. We now show that this set can be further reduced to just considering $N=1$.

When the program state $\pi=\chi^{\otimes N}$ is directly used in the PBT protocol we find

[TABLE]

namely that the optimization can be reduced to the $\mathcal{O}(d^{4}$) dimensional space of Choi matrices $\chi$. Note that, in the above equation, we used the identity

[TABLE]

where $\mathbf{\bar{A}}_{i}=\mathbf{A}\backslash A_{i}$.

Now let $\pi$ be a linear combination of tensor products of Choi matrix states, $\chi_{k}^{\otimes N}$, each with probability $p_{k}$ as in Eq. (S86). Then we can write

[TABLE]

However, this is precisely the partial trace over the tensor product $\chi^{\prime\otimes N}$ of some other Choi matrix $\chi^{\prime}=\sum_{k}p_{k}\chi_{k}$. Hence, the program state $\pi=\sum_{k}p_{k}\chi_{k}^{\otimes N}$ simulates the same channel as the resource state $\pi^{\prime}=\left(\sum_{k}p_{k}\chi_{k}\right)^{\otimes N}$.

Therefore, the optimization over the convex set $\mathcal{C}_{N}$ can be reduced to the optimization over products of Choi matrices $\chi^{\otimes N}$. From Eq. (S89) this can be further reduced to the optimization of the quantum channel $\tilde{\Lambda}$ over the convex set of single-copy Choi matrices $\chi$

[TABLE]

which is $\mathcal{O}(d^{4})$. Using $\mathcal{C}_{1}$ drastically reduces the difficulty of numerical simulations, thus allowing the exploration of significantly larger values of $N$.

Finally, we provide an explicit expression for the reduced map $\tilde{\Lambda}$ of Eq. (S89) in the case of qubits. For $d=2$ we can rewrite PBT in a language that can be more easily formulated from representations of SU(2). For simplicity of notation, here we do not use bold letters for vectorial quantities. Let us modify the POVM in Eq. (21) as

[TABLE]

where $|\Psi^{-}\rangle=(|01\rangle-|10\rangle)/\sqrt{2}$ is a singlet state. For $\pi=\chi^{\otimes n}$ the quantum channel is simplified. In fact, since $\operatorname*{Tr}_{B}\chi=\openone/2$, we may write

[TABLE]

where $\ell$ and $\ell^{\prime}$ are multi-indices and, in defining the Kraus operators, we have separated the contributions from $\tilde{\Pi}_{i}$ and $\Delta$ (see below).

In order to express these operators, we write

[TABLE]

so that

[TABLE]

where $\vec{S}=\vec{\sigma}/2$ is a vector of spin operators, $\vec{S}_{A}=\sum_{j}\vec{S}_{A_{j}}$ and $\vec{S}_{\mathrm{tot}}=\vec{S}_{C}+\vec{S}_{A}$. The eigenvalues of $\sigma_{AC}$ are then obtained from the eigenvalues of the three commuting Casimir operators

[TABLE]

where $S_{\mathrm{tot}}=s_{A}\pm 1/2$.

Substituting the definition of $S_{\mathrm{tot}}$, we find two classes of eigenvalues

[TABLE]

with corresponding eigenvectors

[TABLE]

where $-\frac{N+1}{2}\leq M\leq\frac{N+1}{2}$, $\alpha=1,\dots,g^{[N]}(s)$ describes the degeneracy, $g^{[N]}(s)$ is the size of the degenerate subspace, and

[TABLE]

are Clebsch-Gordan coefficients.

Note that the Clebsch-Gordan coefficients define a unitary transformation between the two bases $|s_{1},m_{1};s_{2};m_{2}\rangle$ and $|S,M;s_{1},s_{2}\rangle$. From the orthogonality relations, of these coefficients we find the equalities

[TABLE]

where $\delta(S,S^{\prime},s)=1$ iff $S=S^{\prime}$ and $|s-1/2|\leq S\leq s+1/2$. The eigenvalues in Eq. (S102) are zero iff $S_{\mathrm{tot}}=S_{A}+1/2$ and $S_{A}=N/2$. These eigenvalues have degeneracy $2S_{\mathrm{tot}}+1=N+2$ and the corresponding eigenvectors are

[TABLE]

Thus, the operator $\Delta$ from Eq. (S97) may be written as

[TABLE]

To finish the calculation we need to perform the partial trace over all spins except those in port $i$. We use $s_{\bar{A}_{i}}$, $m_{\bar{A}_{i}}$ and $\alpha_{i}$ to model the state of the total spin in ports ${A_{j}}$ with $j\neq i$. These refer to the value of total spin and the projection along the $z$ axis, as well as the degeneracy. Moreover, since $S_{\bar{A}_{i}}$ commutes with both $S_{A}^{2}$ and $S_{A}^{z}$, we may select a basis for the degeneracy that explicitly contains $s_{\bar{A}_{i}}$. We may write then $\alpha=(s_{\bar{A}_{i}},\tilde{\alpha}_{i})$ where $\tilde{\alpha}_{i}$ represents some other degrees of freedom.

With the above definitions, when we insert several resolutions of the identity in Eq. (S98), we may write the Kraus operators as

[TABLE]

where each set of states $|s_{\bar{A}_{i}},m_{\bar{A}_{i}},\alpha_{i}\rangle$ represents a basis of the space corresponding to all ports $j$ with $j\neq i$. To simplify the Kraus operators we study the overlap

[TABLE]

In the last line we find that the overlap is independent of $\alpha$ and $\alpha_{i}$, though with constraints $\alpha=(s_{\bar{\imath}},\alpha_{i})$, which requires $\alpha_{i}=\alpha_{i}^{\prime}$. Therefore, different Kraus operators provide exactly the same operation and, accordingly, we can sum over these equivalent Kraus operators to reduce the number of indices. After this process, we get

[TABLE]

The Kraus operators of the reduced channel $\tilde{\Lambda}$ are obtained as $(K_{\ell}^{u}\otimes\openone_{D})(|\Psi_{CD}^{-}\rangle\otimes\openone_{AB})$. It is simple to check that the above operators define a CPTP-map.