Computing Quantum Channel Capacities

TL;DR

This paper extends classical Blahut-Arimoto algorithms to quantum channels, providing efficient iterative methods to compute various quantum capacities with proven error bounds and demonstrating fast convergence through numerical experiments.

Contribution

It introduces quantum generalizations of classical capacity algorithms, enabling accurate and efficient computation of multiple quantum channel capacities.

Findings

Algorithms converge rapidly in numerical tests.

Rigorous bounds on estimation errors are established.

Methods applicable to classical-quantum and quantum channels.

Abstract

The capacity of noisy quantum channels characterizes the highest rate at which information can be reliably transmitted and it is therefore of practical as well as fundamental importance. Capacities of classical channels are computed using alternating optimization schemes, called Blahut-Arimoto algorithms. In this work, we generalize classical Blahut-Arimoto algorithms to the quantum setting. In particular, we give efficient iterative schemes to compute the capacity of channels with classical input and quantum output, the quantum capacity of less noisy channels, the thermodynamic capacity of quantum channels, as well as the entanglement-assisted capacity of quantum channels. We give rigorous a priori and a posteriori bounds on the estimation error by employing quantum entropy inequalities and demonstrate fast convergence of our algorithms in numerical experiments.

Figures

Tables

Table 1. Table 1: Summary of the Blahut-Arimoto algorithms discussed in this work, with the overall asymptotic worst case complexity for an additive ε 𝜀 \varepsilon -approximation. X 𝑋 X and Y 𝑌 Y refer to classical registers while A 𝐴 A and B 𝐵 B refer to quantum registers. The Kraus rank of ℰ ℰ \mathcal{E} , i.e., the minimal number of Kraus operators required to represent ℰ ℰ \mathcal{E} (see e.g.  tomamichel_quantum_2015 for the definition of Kraus operators), is denoted by K ≤ | A | ​ | B | 𝐾 𝐴 𝐵 K\leq|A||B| . The complexity for the classical case is the same as in the original works  Blahut72 ; Arimoto72 . For the coherent information of quantum channels, we require the channel ℰ ℰ \mathcal{E} to lie in the class of less noisy channels.

Channels	Quantity	Time Complexity
$X\to Y$	Mutual information $I\left(ℰ\right)$	$𝒪\left(\frac{\left|X\right|\left|Y\right|\log \left|X\right|}{\epsilon }\right)$
$X\to B$	Holevo quantity $\chi \left(ℰ\right)$	$𝒪\left(\frac{\left({\left|B\right|}^3+{\left|B\right|}^2\left|X\right|\right)\log \left|X\right|}{\epsilon }+\left|X\right|{\left|B\right|}^3\right)$
$A\to B$	Thermodynamic capacity $T_{\Gamma }\left(ℰ\right)$	$𝒪\left(\frac{\left({\left|A\right|}^3+{\left|A\right|}^2{\left|B\right|}^2+{\left|B\right|}^3\right)\log \left|A\right|}{\epsilon }\right)$
$A\to B$	Coherent information $I_{\mathrm{coh}}\left(ℰ\right)$	$𝒪\left(\frac{\left({\left|A\right|}^3+{\left|B\right|}^3+K^3+{\left|A\right|}^2\left({\left|B\right|}^2+K^2\right)\right)\log \left|A\right|}{\epsilon }\right)$
$A\to B$	Quantum mutual information $I\left(ℰ\right)$	$𝒪\left(\frac{\left({\left|A\right|}^3+{\left|B\right|}^3+K^3+{\left|A\right|}^2\left({\left|B\right|}^2+K^2\right)\right)\log \left|A\right|}{\epsilon }\right)$

Equations

$$I_{\mathrm{coh}}(\mathcal{E})=\underset{\rho }{\max }\underset{=I_{\mathrm{coh}}(\rho ,\mathcal{E})}{\underbrace{S(\mathcal{E}(\rho ))-S({\mathcal{E}}_c(\rho ))}},$$

$${\rho }^{(n+1)}=\frac{1}{Z^{(n+1)}}\exp \left(\log {\rho }^{(n)}+\mathcal{F}({\rho }^{(n)})\right),$$

$$n=\lceil \log \frac{N}{\epsilon }\rceil$$

$${\mathcal{E}}_p^{AD}(\rho )=A_0\rho A_0^{†}+A_1\rho A_1^{†}$$

$$\mathcal{O}\left({\epsilon }^{-1}\max \{|X|,|B|\}|B{|}^3\sqrt{\log |X|}\right)$$

$$D(\rho \Vert \sigma )=\{\begin{array}{ll}\mathrm{Tr}\left[\rho \left(\log \rho -\log \sigma \right)\right]& \text{if σ≫ρ}\\ \infty & \text{otherwise},\end{array}$$

$$J_{\gamma }(\rho ,\sigma )=-\gamma D(\rho ||\sigma )+\mathrm{Tr}\left[\rho \mathcal{F}(\sigma )\right]\in \mathbb{R},$$

$${\mathcal{F}}_1(\rho )=\underset{\sigma \text{ with }\sigma \gg \rho }{\mathrm{argmax}}{J}_{\gamma }(\rho ,\sigma ),{\mathcal{F}}_2(\sigma )=\underset{\rho \text{ with }\sigma \gg \rho }{\mathrm{argmax}}{J}_{\gamma }(\rho ,\sigma ).$$

$$\mathrm{Tr}\left[\rho \left\{\mathcal{F}\left(\sigma \right)-\mathcal{F}\left(\rho \right)\right\}\right]\le \gamma D(\rho ||\sigma )\text{for density operators σ≫ρ,}$$

$${\mathcal{F}}_1(\rho )=\rho ,$$

$${\mathcal{F}}_2(\sigma )=\frac{1}{Z(\sigma )}\exp \left(\log \sigma +\frac{1}{\gamma }\mathcal{F}\left(\sigma \right)\right)$$

$$J_{\gamma }(\rho ,\rho )=\mathrm{Tr}\left[\rho \mathcal{F}\left(\rho \right)\right]\ge \mathrm{Tr}\left[\rho \mathcal{F}\left(\sigma \right)\right]-\gamma D(\rho ||\sigma )=J_{\gamma }(\rho ,\sigma )\text{for density operators σ≫ρ.}$$

$${\mathcal{F}}_1(\rho )=\underset{\sigma \gg \rho }{\mathrm{argmax}}{J}_{\gamma }(\rho ,\sigma )=\rho .$$

$$\underset{\rho }{\mathrm{argmax}}{J}_{\gamma }(\rho ,\sigma )=\underset{\rho }{\mathrm{argmax}}-\mathrm{Tr}\left[\rho \log \rho \right]+\mathrm{Tr}\left[\rho \left\{\log \sigma +\frac{1}{\gamma }\mathcal{F}(\sigma )\right\}\right].$$

$$0\le \mathrm{Tr}\left[\rho \left\{\mathcal{F}\left(\sigma \right)-\mathcal{F}\left(\rho \right)\right\}\right]\le \gamma D(\rho ||\sigma )\text{for density operators σ≫ρ,}$$

$$C^{\star }=\underset{\rho ,\sigma \text{ with }\sigma \gg \rho }{\max }{J}_{\gamma }(\rho ,\sigma )$$

$$|C^{\star }-C(n)|\le \frac{\gamma D\left({\rho }^{\star }\Vert {\rho }^{(1)}\right)}{n},$$

$$C(t)$$

$$=\mathrm{Tr}\left[{\rho }^{(t+1)}\left[-\gamma \log {\rho }^{(t+1)}+\gamma \log {\rho }^{(t)}+\mathcal{F}\left({\rho }^{(t)}\right)\right]\right]$$

$$=\mathrm{Tr}\left[{\rho }^{(t+1)}\left[-\gamma \left\{\log {\rho }^{(t)}+\frac{1}{\gamma }\mathcal{F}\left({\rho }^{(t)}\right)-\log Z^{(t+1)}\right\}+\gamma \log {\rho }^{(t)}+\mathcal{F}\left({\rho }^{(t)}\right)\right]\right]$$

$$=\gamma \log Z^{(t+1)}.$$

$$\mathrm{Tr}\left[{\rho }^{\star }(\log {\rho }^{(t+1)}-\log {\rho }^{(t)})\right]$$

$$=\frac{1}{\gamma }\left(C^{\star }-C(t)+\mathrm{Tr}\left[{\rho }^{\star }\left\{\mathcal{F}\left({\rho }^{(t)}\right)-\mathcal{F}\left({\rho }^{\star }\right)\right\}\right]\right)$$

$$\ge \frac{1}{\gamma }|C^{\star }-C(t)|,$$

$$\sum _{t=1}^k\mathrm{Tr}\left[{\rho }^{\star }\left(\log {\rho }^{(t+1)}-\log {\rho }^{(t)}\right)\right]$$

$$=D\left({\rho }^{\star }\Vert {\rho }^{(1)}\right)-D\left({\rho }^{\star }\Vert {\rho }^{(k+1)}\right)$$

$$\le D\left({\rho }^{\star }\Vert {\rho }^{(1)}\right),$$

$$\sum _{t=1}^{\infty }|C^{\star }-C(t)|\le \gamma D\left({\rho }^{\star }\Vert {\rho }^{(1)}\right),$$

$${\gamma }^{(t+1)}=\frac{\mathrm{Tr}\left[{\rho }^{(t)}\left\{\mathcal{F}\left({\rho }^{(t-1)}\right)-\mathcal{F}\left({\rho }^{(t)}\right)\right\}\right]}{D\left({\rho }^{(t)}\Vert {\rho }^{(t-1)}\right)}.$$

$$C(t)\le C^{\star }\le {\alpha }_{\max }^{(t)}$$

$$\underset{\rho \ne \sigma }{\inf }\frac{\mathrm{Tr}\left[\rho \left\{\mathcal{F}\left(\sigma \right)-\mathcal{F}\left(\rho \right)\right\}\right]}{D(\rho ||\sigma )}=a\text{for density operators σ≫ρ},$$

$$C^{\star }-C(t+1)$$

$$D({\rho }^{\star }||{\rho }^{(t)})-D({\rho }^{\star }||{\rho }^{(t+1)})$$

$$\ge \frac{a}{\gamma }D({\rho }^{\star }||{\rho }^{(t)})$$

$$⟹D({\rho }^{\star }||{\rho }^{(t+1)})$$

Full text

I Introduction

A fundamental problem in information theory is to compute the maximum rate at which information can be reliably transmitted over a noisy channel. For channels described by classical physics, Shannon’s seminal work shannon1948mathematical shows that every channel can be characterized by a single number — the capacity of the channel. More precisely, the channel capacity corresponds to the highest rate at which messages can be transmitted while ensuring that the error in the communication vanishes in the limit of long messages. In general, there is no closed analytical expression for the capacity of a channel but instead it is given in terms of an optimization problem. As such, efficient numerical methods to solve this optimization problem are sought after.

For channels governed by quantum physics, there exist multiple distinct settings where information must be transmitted over a noisy channel. Analogously to the classical regime, a quantum channel is given by any mapping that sends physical states to physical states. A graphical depiction is shown in Figure 1 and one may, for example, choose to transmit classical or quantum information, use pre-shared entanglement between the sender and the recipient to enhance the transmission rate, or look at minimizing the entropy gain between the output and input states of the channel. Hence, there exist many different relevant quantum channel capacities and — as with classical channels — these are typically given in terms of optimization problems which have to be solved numerically.

For concreteness, we focus in this introduction on the quantum capacity of a channel which is related to the so-called coherent information lloyd1997capacity ; shor2002quantum ; devetak2005private . For a quantum channel $\mathcal{E}$ the coherent information is defined as

[TABLE]

where the maximization is taken over input quantum states, the von Neumann entropy is denoted by $S(\rho)=-\mathrm{Tr}\left[\rho\log\rho\right]$, and $\mathcal{E}_{c}$ is the complementary channel of $\mathcal{E}$ (defined in Section III.1). As information is preserved in quantum physics, the complementary channel captures the notion that any information that is not transmitted to the recipient must have leaked to the environment. A channel is called less noisy when the information that leaks to the environment is less than the information transmitted (for a precise definition, see (65)). In this case, the considered problem becomes convex watanabe2012private and as in classical information theory, one could in principle use standard methods from convex optimization boyd2004convex to solve the resulting problem. However, the gold standard to compute channel capacities in classical information theory are so-called Blahut-Arimoto algorithms Blahut72 ; Arimoto72 . They have the following particularly strong properties:

• •

The algorithm is of an iterative form where each iteration step is of a very simple form given by a matrix exponential. Hence, no matrix derivatives have to be calculated in the update rules which is in contrast to, e.g., gradient based methods.

• •

Explicit a priori and a posteriori convergence guarantees are directly derived from entropic inequalities.

• •

The algorithm typically converges extremely fast in practice. Additionally, in certain cases, there is even exponentially fast convergence.

In our work, we generalize Blahut-Arimoto algorithms to the quantum setting and show that aforementioned desirable properties are retained. For example, we show that the coherent information of less noisy channels can be approximated by iteratively improving on an estimate of the input state $\rho$ to the channel. The update rule for each iteration is given as

[TABLE]

where $\mathcal{F}(\sigma)=\mathcal{E}_{c}^{\dagger}\log\mathcal{E}_{c}(\sigma)-\mathcal{E}^{\dagger}\log\mathcal{E}(\sigma)$, $\mathcal{E}^{\dagger}$ and $\mathcal{E}_{c}^{\dagger}$ denote the adjoint channels (see e.g. tomamichel_quantum_2015 for the definition) of $\mathcal{E}$ and $\mathcal{E}_{c}$, respectively, and $Z^{(n+1)}$ takes care of normalization. The estimate of the coherent information after the $n$-th iteration step is then given by $I_{\rm{coh}}^{(n)}=I_{\rm{coh}}(\rho^{(n+1)},\mathcal{E})$ and we prove that starting from an initial input state $\rho^{(1)}=\frac{1}{N}$ we have after

[TABLE]

Here, $N$ denotes the dimension of the input system of the channel. The number of iterations required in (3) corresponds to the worst-case and in practice, convergence is much faster for the channels we study. For example, we consider the qubit amplitude damping channel which is defined as follows

[TABLE]

Setting $p=0.3$, we compute $I_{\rm coh}(\mathcal{E}^{AD}_{0.3})$ using the generalized Blahut-Arimoto algorithm up to an additive error of $\epsilon=10^{-6}$ bits. We introduce an a posteriori error criterion that allows us to find the error in our capacity estimate at each iteration and terminate the algorithm if this error is sufficiently small. When we use this criterion, the number of iterations is indeed much smaller than that required in (3) as shown in Figure 2.

In the following, we give an overview of our results showing that quantum Blahut-Arimoto algorithms can in fact be used to estimate several fundamental variants of channel capacities in quantum information theory.

II Overview of results

II.1 Quantum settings

In addition to the coherent information of less noisy channels, other settings are of interest when one considers quantum channel capacities. For example, one may also ask the following questions:

How much classical information can be sent per use of a noisy channel with classical input and quantum output? 2. 2.

How much information can be sent over a quantum channel in the presence of physical restrictions imposed by thermodynamics? 3. 3.

How much quantum information can be sent over a quantum channel with entanglement assistance, i.e., with access to pre-shared entangled states between the two communicating parties?

These quantities are characterized by different capacities that are related to well-known quantities in quantum information theory, respectively the

Holevo quantity holevo_capacity_1998 ; schumacher_sending_1997 2. 2.

Thermodynamic capacity faist_thermodynamic_2018 3. 3.

Quantum mutual information Bennett02 .

We show that all the entropic optimization problems corresponding to these quantities as well as the coherent information discussed in the introduction can be cast into the same form and solved using an update rule similar to (2) for suitable choice of $\mathcal{F}$. This generalization of the Blahut-Arimoto algorithm to the quantum setting provides the first efficient schemes for the calculation of the thermodynamic capacity of quantum channels, the coherent information of less noisy channels, and the quantum mutual information of quantum channels.

II.2 Algorithm speed and time complexity

The number of iterations required by the generalized Blahut-Arimoto algorithm to obtain an a priori error of $\varepsilon$ is $\mathcal{O}\left(\varepsilon^{-1}\log N\right)$, where $N$ denotes the input dimension of the channel. To find the worst-case time complexity of the full algorithm that estimates the capacity up to an additive error $\varepsilon$, one can simply multiply the required number of iterations with the complexity of the update rule given in (2) for appropriate choice of $\mathcal{F}$. The full complexity of the algorithm for an $\varepsilon$-estimation of the capacity is given in Table 1.

II.3 Additional speed-ups

One can also construct an accelerated version of the generalized Blahut-Arimoto algorithm similar to the classical case Matz04 ; Naja09 . The accelerated version lowers the number of iterations by a constant factor in the worst case compared to the standard version and hence, the asymptotic results given Table 1 are not affected by this. It also gives rise to heuristics that provide a significant speed-up of the algorithm in practical instances as shown in Figures 3, 4, 5 and 6.

The a priori bound on the error corresponds to the worst case and the number of iterations required is extremely conservative for most problems. Hence, we adapt classical techniques Vontobel03 to find a posteriori error bounds on the capacity estimate at each iteration of the algorithm (see also li2019computing ). This error bound can then be used as a termination criterion for our numerics as described in Sections IV.3, V.3, VI.3 and VII.3. In fact, even if we use an alternative algorithm or heursitics, the a posteriori criterion gives us guarantees on the error of our estimate compared to the actual capacity.

We notice that for the channels and capacities we consider, the a posteriori error is achieved in significantly fewer iterations compared to the a priori requirement. One of the explanations for this speedup is that the convergence of our algorithm can be exponentially fast in certain cases. Specifically, when certain entropy inequalities used by our convergence proof are strict inequalities, one can prove (see Proposition III.6 for more details) that the error in the capacity estimate after the $t^{\rm th}$ iteration for a channel with input dimension $N$ is bounded by $(1-\delta)^{t}\log N$ for some $\delta\in(0,1]$.

II.4 Comparison with previous work

Previous works have considered other methods to approximate the solution of entropic optimization problems in quantum information theory nagaoka1998algorithms ; Nagaoka01 ; Shor03 ; Hayashi05 ; kato2009computational ; tanuma2011revisiting ; Sutter16 ; Fawzi18 ; Fawzi2019 ; Gour10 ; Girard14 . Here, we provide the first fully quantum version for Blahut-Arimoto type algorithms, give convergence guarantees and demonstrate numerically that these algorithms outperform all previously known algorithms for the specific optimization tasks. In particular, calculating the classical capacity of classical-quantum channels has been the subject of several previous works. Nagaoka nagaoka1998algorithms proposed a Blahut-Arimoto type algorithm to compute the Holevo quantity but did not provide an analysis of the convergence and the time complexity. A provable convergence bound for the estimation of the Holevo quantity was previously given in Sutter16 by using an algorithm based on convex programming duality and smoothing techniques and has a time complexity of

[TABLE]

This suggests that our generalized Blahut-Arimoto algorithm is faster for the regime $|B|\gg|X|$ but slightly slower for $|X|\gg|B|$ in the worst case. However, our numerical results suggest that Blahut-Arimoto type algorithms converge much faster in practice than the algorithm given in Sutter16 . Further, the techniques used in Sutter16 could not be generalized to fully quantum settings. For other algorithms that compute the Holevo quantity, there is often no explicit complexity analysis given. In particular, in kato2009computational ; tanuma2011revisiting , an algorithm to approximate the Holevo quantity of arbitrary quantum channels is provided. However, no convergence guarantees are given and numerics suggest that the runtime of the algorithm is quite high in practice kato_voronoi_2008 .

We compare our numerics with the best known runtimes given in Fawzi18 , whose algorithm is based on semi-definite approximations of the matrix logarithm Fawzi2019 . Our numerics suggest that for the specific optimization problems considered here, Blahut-Arimoto type algorithms converge faster (Section IV.3). We emphasize that while the advantage of Blahut-Arimoto type algorithms might be explained by their conformity to the specific structure of entropy optimization problems, the method presented in Fawzi18 is applicable for a wider range of problems.

The rest of the paper is structured as follows. We first give the general structure and convergence proof of Blahut-Arimoto alternating optimization algorithms in Section III and then discuss in more detail the classical-quantum case (incorporating the classical case) in Section IV, the thermodynamic capacity of quantum channels in Section V, the coherent information of less noisy channels in Section VI and the mutual information of quantum channels in Section VII. Numerical results are provided for the different cases in Sections IV.3, V.3, VI.3 and VII.3, respectively. Finally, we discuss some conclusions and give an outlook (Section VIII).

III Structure of quantum Blahut-Arimoto algorithms

III.1 Notation

We label Hilbert spaces (which are assumed to be finite-dimensional) with capital letters e.g. $A$, $B$, and denote their dimension by $|A|$, $|B|$, and so on. The set of density operators on a system $A$, i.e., positive semi-definite matrices $\rho_{A}$ with $\mathrm{Tr}\left[\rho_{A}\right]=1$, is denoted $\mathcal{D}(A)$. Whenever we work with a single system, we may omit the sub-index referring to the system. A quantum channel $\mathcal{E}_{A\rightarrow B}$ is a linear completely positive trace-preserving map that maps density operators on a system $A$ to density operators on a system $B$. For simplicity, we may write $\mathcal{E}$ instead of $\mathcal{E}_{A\rightarrow B}$ if the input and output systems of the channel are clear from the context. The complementary channel $\mathcal{E}^{c}_{A\to E}$ of a channel $\mathcal{E}_{A\rightarrow B}$ is defined via taking the partial trace over $B$ of the Stinespring dilation $U_{A\to BE}$ of $\mathcal{E}_{A\to B}$.

For a density operator $\rho$, the von Neumann entropy is defined as $S(\rho)=-\mathrm{Tr}\left[\rho\log\rho\right]$ and for density operators $\rho,\sigma$ the quantum relative entropy is defined as

[TABLE]

where the notation $\sigma\gg\rho$ denotes that the kernel of $\sigma$ is a subset of the kernel of $\rho$ (i.e., $\ker(\sigma)\subseteq\ker(\rho)$), and where we take the logarithm only on the support of the argument. We also work with discrete probability distributions, which we represent as vectors $\lambda=[\lambda_{1},\dots,\lambda_{m}]$ with $\sum_{i}\lambda_{i}=1$ or alternatively as diagonal matrices with entries $\lambda_{1},\dots,\lambda_{m}$ i.e., $\rho_{\lambda}=\sum_{i=1}^{m}\lambda_{i}|i\rangle\langle i|$ in some pre-defined basis. The definition of the von Neumann entropy and the quantum relative entropy then simplifies for probability distributions to the Shannon entropy and the Kullback-Leibler divergence, respectively.

III.2 Blahut-Arimoto type

In the following we study a special entropic type of alternating optimization algorithms, called Blahut-Arimoto algorithms Blahut72 ; Arimoto72 111In the broader context of numerical optimization, Blahut-Arimoto algorithms can be viewed as a proximal point method that maximizes a function iteratively with a penalty term if the current guess moves away from the previous guess Matz04 ; Naja09 .. The channel capacities we are interested in are typically given as a convex optimization problem over input states. One may write an extension function $J$ in two variables such that the maximization over both variables gives back the capacity of the channel. Performing the maximizations iteratively leads to an algorithm of the following form.

Definition III.1** (Blahut-Arimoto)**

Algorithm 1 together with the following conditions on $J$, $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ defines a quantum Blahut-Arimoto algorithm: For $\gamma>0$ and density operators $\sigma\gg\rho$,

[TABLE]

*where $\mathcal{F}$ is a Hermitian matrix valued super-operator on density operators such that $\mathrm{Tr}\left[\rho\mathcal{F}(\sigma)\right]$ is continuous in $\sigma$ for $\sigma\gg\rho$. The update rules are then given by *

[TABLE]

Now, under certain conditions we can find analytic expressions for the optimizers in (8).222We can choose any optimizer in the above optimizations if there is not a unique one.

Lemma III.2** (Update rules)**

For quantum Blahut-Arimoto algorithms with

[TABLE]

the update rules appearing in Algorithm 1 are of the form

[TABLE]

for density operators $\rho$ and $\sigma>0$.

Proof. By (9) together with the definition in (7), we find

[TABLE]

Hence, the fist update rule is given by

[TABLE]

For the second update rule we assume that $\sigma>0$ and hence the optimizer $\rho$ can be chosen from the set of density operators. The form of the optimizer follows directly from Gibbs’ variational principle (Lemma A.1). To see this, we write for $\sigma>0$,

[TABLE]

Setting $\omega=\rho$ and $H=\log\sigma+\frac{1}{\gamma}\mathcal{F}(\sigma)$ in Lemma A.1 leads to the claimed form of the second update rule.

The next theorem shows convergence of quantum Blahut-Arimoto algorithms under some technical assumptions (which will be satisfied for the applications in Sections IV-VI).

Theorem III.3** (Convergence)**

For quantum Blahut-Arimoto algorithms together with a strictly positive definite initial state $\rho^{(1)}>0$ on a system $A$ and

[TABLE]

we have that $C(n)$ of Algorithm 1 is monotonically increasing and converges for $n\rightarrow\infty$ to

[TABLE]

with the approximation error bounded as

[TABLE]

where $\rho^{\star}$ is an optimizer (possibly not unique) that achieves the capacity $C^{\star}$. If $\rho^{(1)}$ is chosen to be the maximally mixed state, the approximation error is bounded as $|C^{\star}-C(n)|\leq\frac{\gamma\log|A|}{n}$.

Proof. Let $t\in\mathbb{N}$ and let $\rho^{(t)}\in\mathcal{D}(A)$ be a density operator in the $t$-th iteration step of Algorithm 1. Note that $\rho^{(t)}>0$ for all $t$, since the exponentiation of a matrix in the update rule (11) ensures full support. First, using the update rules given in (10) and in (11), and setting $Z^{(t+1)}=Z(\rho^{(t)})$, we find

[TABLE]

Let $\rho^{\star}=\operatorname*{arg\,max}_{\rho}J_{\gamma}(\rho,\rho)$, and note that such a $\rho^{\star}$ always exists because $J_{\gamma}(\rho,\rho)$ is continuous by Definition III.1 and we optimize over the compact set of all density operators. Further, we have $C^{\star}=J_{\gamma}(\rho^{\star},\rho^{\star})=\mathrm{Tr}\left[\rho^{\star}\,\mathcal{F}(\rho^{\star})\right]$. Using $C(t)=\gamma\,\log Z^{(t+1)}$, we derive an upper bound on the additive error at the iteration step $t$

[TABLE]

where we used the assumption (15) of the lemma in the last step and that $C^{\star}\geq C(t)$, since $C^{\star}=J_{\gamma}(\rho^{\star},\rho^{\star})=\max_{\rho,\sigma\textnormal{ with }\sigma\gg\rho}J_{\gamma}(\rho,\sigma)$ is the maximum value that $J_{\gamma}$ can achieve. The sum over the additive error terms is upper bounded by a telescopic sum, which can itself be upper bounded as follows

[TABLE]

where we used the positivity of the quantum relative entropy in the last inequality. We conclude

[TABLE]

and since $D\left(\rho^{\star}\middle\|\rho^{(1)}\right)$ is finite for $\rho^{(1)}\in\mathcal{D}(A)$ with full support, $C(n)$ converges to $C^{\star}$ for $n\rightarrow\infty$. Moreover, $C(n)$ increases monotonically in $n$ by construction of the updated states as optimizers in Lemma III.2. Together with (28), this implies the error bound stated in the theorem.

The value of the parameter $\gamma$ appearing in Definition III.1 corresponds to accelerated versions of the Blahut-Arimoto algorithm and hence, $\gamma$ should be chosen as small as possible. However, to have provable convergence, $\gamma$ is lower bounded by the requirement in Theorem III.3 that (15) is satisfied. The standard algorithm uses the data processing inequality for the quantum relative entropy to find a $\gamma$ satisfying (15), whereas the accelerated algorithm uses relative entropy coefficients similar to (40) to establish tight lower bounds on $\gamma$. Note that calculating these coefficients is itself a difficult optimization problem and may be intractable for many channels. However, as a heursitic, we may still introduce an adaptive acceleration parameter $\gamma^{(t)}$ at each iteration, similar to the classical case Matz04 .

Remark III.4** (Acceleration heuristic)**

For $D(\rho^{(t)}||\rho^{(t-1)})\neq 0$, define an adaptive acceleration parameter

[TABLE]

This adaptive acceleration parameter may be used when strict lower bounds on $\gamma$ cannot be calculated. The proof of convergence from Theorem III.3 no longer holds since the capacity estimate does not necessarily increase monotonically when the adaptive acceleration method is used. Nevertheless, we find in practice that the adaptive acceleration heuristic provides a significant speed-up (see Figures 3, 4, 5 and 6).

Finally, regardless of whether we use the adaptive acceleration parameter or otherwise, the following proposition allows us to terminate the algorithm when the a posteriori error is sufficiently small (see also li2019computing ).

Proposition III.5** (Termination criteria)**

For quantum Blahut-Arimoto algorithms that satisfy (15) let us denote the largest eigenvalue of $\mathcal{F}(\rho^{(t)})$ by $\alpha^{(t)}_{\max}$ for $t\in\mathbb{N}$. Then, we have

[TABLE]

with $C^{\star}=\max_{\rho,\sigma\textnormal{ with }\sigma\gg\rho}J_{\gamma}(\rho,\sigma)$. Hence, in particular, $|C^{\star}-C(t)|\leq\alpha^{(t)}_{\max}-C(t)$ provides a bound on the a posteriori error.

Proof. For any $t\in\mathbb{N}$, we have $\mathrm{Tr}\left[\rho^{\star}(\mathcal{F}(\rho^{(t)})-\mathcal{F}(\rho^{\star}))\right]\geq 0$ from (15), where $\rho^{\star}$ is the optimizer (possibly not unique) that achieves the capacity $C^{\star}$. Therefore, there must exists at least one eigenvalue $\alpha_{i}^{(t)}$ of $\mathcal{F}(\rho^{(t)})$ that satisfies $\alpha_{i}^{(t)}\geq\mathrm{Tr}\left[\rho^{\star}\,\mathcal{F}(\rho^{\star})\right]=C^{\star}$. In particular, $\alpha_{\max}^{(t)}\geq C^{\star}$. Noting that $C(t)\leq C^{\star}$ finishes the proof.   

In particular, if $\rho^{(t)}$ converges to a full rank optimizer $\rho^{\star}$, then we may take the logarithm on both sides of (11) to show that $\lim_{t\rightarrow\infty}\alpha^{(t)}_{\max}=C^{\star}$. This shows that the a posteriori criterion is indeed achieved for some sufficiently large $t$ for full rank optimizers. Next, we show in the following proposition that under certain technical conditions, the speed of convergence is exponential.

Proposition III.6** (Exponential convergence)**

For quantum Blahut-Arimoto algorithms satisfying (15) with $a\in(0,\gamma]$ such that

[TABLE]

there exists a unique optimizer $\rho^{\star}$ with $C^{\star}=J_{\gamma}(\rho^{\star},\rho^{\star})$, and the algorithm converges exponentially fast. That is, for all $t\in\mathbb{N}$ we have

[TABLE]

where $\delta=\frac{a}{\gamma}\in(0,1]$. Moreover, when $\rho^{(1)}$ is chosen to be the maximally mixed state, we obtain $C^{\star}-C(t+1)\leq\gamma(1-\delta)^{t}\log|A|$.

Proof. We first establish the uniqueness of the optimizer. From (22) and (23) in the proof of Theorem III.3, we have

[TABLE]

Note that (35) implies that there is a unique optimizer $\rho^{\star}$ since this inequality holds for each $t$. Indeed, to see this, we start by noting that (35) implies $\lim_{t\rightarrow\infty}D(\rho^{\star}||\rho^{(t)})=0$. Then, Pinsker’s inequality (watroustheory, , Theorem 5.38) implies that also $\lim_{t\rightarrow\infty}\lVert\rho^{(t)}-\rho^{\star}\rVert_{1}=0$. Since the trace distance is a norm, we find $\lim_{t\rightarrow\infty}\rho^{(t)}=\rho^{\star}$. Since this must be true for any optimizer $\rho^{\star}$, the optimizer must be unique. From (15) and (33), we find

[TABLE]

Noting that this argument holds for all $t$ finishes the proof.

IV Holevo quantity of classical quantum channels

IV.1 Definitions

Here we consider how to compute the capacity of cq channels, as given by the Holevo quantity. A cq channel can be described by a set of input-output pairs $\{(x,\tau_{x})\}_{x\in\{1,2,\dots,N\}}$, i.e., the channel takes $x\in\{1,2,\dots,N\}$ as an input and provides the quantum states $\tau_{x}\in\mathcal{D}(B)$ as outputs. Alternatively, we may consider cq channels as quantum channels $\mathcal{E}_{X\rightarrow B}$ with $|X|=N$, defined by the completely positive trace-preserving mapping $\mathcal{E}_{X\rightarrow B}:\rho_{X}\mapsto\sum_{x}\langle x|\rho_{X}|x\rangle\,(\tau_{x})_{B}$. Restricting the quantum channel to the classical input states $|x\rangle\langle x|_{X}$ then gives us back the original cq channel. For an input distribution vector $\lambda$, the output of the cq channel corresponds to

[TABLE]

where $\lambda_{i}$ denotes the $i$-th component of the probability vector $\lambda$ and $\rho_{\lambda}=\sum_{k}\lambda_{k}|k\rangle\langle k|$. As shown by Holevo, Schumacher and Westmoreland holevo_capacity_1998 ; schumacher_sending_1997 , the classical capacity of a cq channel is given by the Holevo quantity

[TABLE]

where $I(\lambda,\mathcal{E})$ is the mutual information between the input and output registers.

IV.2 Blahut-Arimoto algorithm

The following algorithm to compute the Holevo quantity is a quantum Blahut-Arimoto algorithm and was first proposed by Nagaoka nagaoka1998algorithms . Similar to the classical case Matz04 , we will show that one can accelerate the convergence of Blahut-Arimoto type algorithms. To do so, we will use the contraction coefficient of relative entropy Ruskai99 , given for a channel $\mathcal{E}$ as

[TABLE]

Note that $0\leq\eta_{\rm con}(\mathcal{E})\leq 1$ by the nonnegativity of relative entropy and by the data processing inequality. We now define the following two variable extension of the mutual information for $\gamma\in[\eta_{\rm con}(\mathcal{E}),1]$.

[TABLE]

where $\rho_{\lambda}=\sum_{k}\lambda_{k}|k\rangle\langle k|$, $\rho_{\mu}=\sum_{k}\mu_{k}|k\rangle\langle k|$ (with $\rho_{\mu}\gg\rho_{\lambda}$) and $\tau_{\mathcal{E},i}=\mathcal{E}(|i\rangle\langle i|)$. Note that $\gamma=1$ corresponds to the standard algorithm and if it is known that $\eta_{\rm con}(\mathcal{E})<1$, then one may accelerate the algorithm. The condition $\eta_{\rm con}(\mathcal{E})<1$ implies that $\forall\rho\neq\sigma,\ D(\mathcal{E}(\rho)||\mathcal{E}(\sigma))<D(\rho||\sigma)$ i.e. the map strictly decreases the relative entropy. We then bring $J_{\gamma}$ into the form

[TABLE]

Note that $\mathrm{Tr}\left[\rho_{\lambda}\mathcal{F}(\rho_{\mu})\right]$ is continuous for all $\rho_{\mu}\gg\rho_{\lambda}$, since for $\lambda_{i}\neq 0$, we have $\mu_{i}\neq 0$ and hence $\ker{\left[\mathcal{E}(\rho_{\mu})\right]}=\ker{\left[\sum\mu_{k}\tau_{\mathcal{E},k}\right]}\subset\ker{\left[\tau_{\mathcal{E},i}\right]}$. Further, a calculation shows

[TABLE]

where we used the nonnegativity of the quantum relative entropy in the last inequality (together with $\gamma\in[\eta_{\rm con}(\mathcal{E}),1]$). From Lemma III.2 whose requirement is satisfied by (46), we find a double optimization form of the capacity $\chi(\mathcal{E})$

[TABLE]

Performing the two maximizations in $\max_{\lambda,\mu}J_{\gamma}(\lambda,\mu,\mathcal{E})$ iteratively, leads to the following algorithm (see Lemma III.2 for the form of the update rules), which provably converges to the Holevo quantity by Theorem III.3. A detailed analysis of the time complexity can be found in Appendix B.1.

Remark IV.1** (Concavity)**

The function $J_{\gamma}(\lambda,\mu,\mathcal{E})=-\gamma\,\mathrm{Tr}\left[\rho_{\lambda}\log\rho_{\lambda}\right]+\mathrm{Tr}\left[\rho_{\lambda}(\gamma\log\rho_{\mu}+\mathcal{F}(\rho_{\mu}))\right]$ is concave in $\lambda$ (for fixed $\mu$). This follows from the concavity of the von Neumann entropy and by noting that the second term is linear in $\lambda$. In contrast to the two-variable function used in the classical Blahut-Arimoto algorithm Blahut72 ; Arimoto72 , the function $J_{\gamma}$ is not concave in $\mu$. To see this, consider the classical channel given by the stochastic matrix

[TABLE]

$\lambda=(\frac{1}{2},\frac{1}{2},0),\,\mu_{1}=(\frac{1}{2},\frac{1}{2},0)$, and $\mu_{2}=(\frac{1}{2},1/4,1/4)$. Indeed, for $a=0.7$ we have $J_{1}(\lambda,a\mu_{1}+(1-a)\mu_{2},\mathcal{E})<aJ_{1}(\lambda,\mu_{1},\mathcal{E})+(1-a)J_{1}(\lambda,\mu_{2},\mathcal{E})$ showing that the function is not concave. On the other hand, for $a=0.3$ we have $J_{1}(\lambda,a\mu_{1}+(1-a)\mu_{2},\mathcal{E})>aJ_{1}(\lambda,\mu_{1},\mathcal{E})+(1-a)J_{1}(\lambda,\mu_{2},\mathcal{E})$ showing that the function is not convex either.

Arimoto Arimoto72 showed in the classical case that if the optimizer that achieves the capacity is unique and full rank, there exists an integer $N\in\mathbb{N}$ such that for all $t>N$, the convergence rate is exponential. That is, for $t>N$, there exists $0<\delta\leq 1$ such that

[TABLE]

Li and Cai li2019computing show an analogous result for cq channels. Here we provide a different criterion that yields exponential convergence from the first iteration albeit with a stronger assumption. First, similar to (40), we define the relative entropy expansion coefficient as

[TABLE]

Note that $\eta_{\rm exp}(\mathcal{E})\in[0,1]$ by the nonnegativity of relative entropy and by the data processing inequality. Using (45), we see that when the coefficient $\eta_{\rm\exp}(\mathcal{E})>0$, we get exponential convergence according to Proposition III.6 with $\delta=\frac{\eta_{\rm exp}(\mathcal{E})}{\gamma}$. A consequence of the condition $\eta_{\rm\exp}(\mathcal{E})>0$ is that $\forall\rho\neq\sigma,\ D(\mathcal{E}(\rho)||\mathcal{E}(\sigma))\neq 0$ and hence $\mathcal{E}(\rho)\neq\mathcal{E}(\sigma)$ i.e. the map is injective.

IV.3 Simulation results

We numerically compute the classical capacity of a cq channel with the Holevo quantity $\chi(\mathcal{E})$ as given in (39). We choose the ensemble of output density operators randomly using the RandomDensityMatrix package from QETLAB qetlab . Figure 3 illustrates the results for a channel with input alphabet of size $10$ and output dimension $16$. The initial guess $\rho_{\lambda}^{(1)}$ is chosen to be the maximally mixed state. For several simple choices of channels, the uniform distribution is optimal and our algorithm requires only one iteration to compute the capacity (this is the reason for choosing a random channel for the illustration here). By Proposition III.3, the capacity estimate $C(n)$ is $\varepsilon$-close to $C^{\star}$ after $\lceil\frac{\gamma\log|X|}{\varepsilon}\rceil$ iterations. Setting $\varepsilon=10^{-6}$ bits, it would require around $10^{6}$ iterations for a provable convergence. To reduce the number of iterations, we use the termination criterion given in Proposition III.5. Hence, if we observe that the maximal eigenvalue of $\mathcal{F}(\rho^{(t)})$ satisfies $\alpha^{(t)}_{\max}-C(t)\leq\varepsilon$ in the $t$-th iteration step, we also have that $C^{\star}-C(t)\leq\varepsilon$. With this termination condition, the standard algorithm converges in $248$ iterations while the adaptive accelerated algorithm converges in $17$ iterations.

We may also compare our numerics to the one given by Fawzi and Fawzi Fawzi18 . Their approach to compute the capacity of this particular cq channel requires $2.21$ seconds using CVX and the Mosek solver Lofberg2004 . Their method does not allow us to bound the error between the capacity estimate and the true capacity but the solver tolerance achieved is $1.3\times 10^{-6}$ bits. For the same channel, our standard Blahut-Arimoto algorithm takes $0.98$ seconds to achieve an a posteriori error of $10^{-6}$ bits and the adaptive accelerated algorithm takes $0.09$ seconds.333All run-times correspond to a machine with a 2.3 GHz Intel Core i5 processor with 8GB of RAM.

V Thermodynamic capacity of quantum channels

V.1 Definitions

The thermodynamic capacity quantifies the information-theoretic power of quantum channels in the presence of physical restrictions imposed by thermodynamics faist_thermodynamic_2018 (see also Navascues15 ). For a quantum channel $\mathcal{E}_{A\rightarrow B}$, relative to operators $\Gamma_{A},\Gamma_{B}>0$ it can be written as

[TABLE]

The $\Gamma$ operators are thereby typically given by Gibbs states generated by Hamiltonians (i.e., some Hermitian operators) $H_{A}$ and $H_{B}$ on the input and output systems, respectively. That is, we have the choice $\Gamma_{A}=\exp\left(-\beta H_{A}\right)$ for a fixed inverse temperature parameter $\beta\geq 0$ and similarly for the output system $B$. However, for the sake of our algorithm we leave $\Gamma_{A},\Gamma_{B}>0$ in general form and refer to faist_thermodynamic_2018 for a discussion of the thermodynamic capacity and its properties. An interesting special case is $\Gamma_{A}=1_{A}$ and $\Gamma_{B}=1_{B}$, for which we get $T_{1}(\mathcal{E})=-G(\mathcal{E})$, where the minimal entropy gain of the quantum channel is given by Alicki04 ; Holevo11 ; Holevo2011

[TABLE]

We refer to Buscemi16 for a discussion of the minimal entropy gain for finite-dimensional quantum channels. The thermodynamic capacity is also related to the completely bounded minimal conditional entropy Devetak06

[TABLE]

with $\rho_{AR}$ a purification of $\rho_{A}$. Namely, we have $S_{\text{CB},\min}(\mathcal{E})=-T_{1}(\mathcal{E}^{c})$ with $\mathcal{E}^{c}_{A\to E}$ the complementary channel of $\mathcal{E}_{A\to B}$. The completely bounded minimal conditional entropy plays a role in entanglement theory as shown in (Devetak06, , Section 5).

V.2 Blahut-Arimoto algorithm

The thermodynamic capacity can be approximated using a Blahut-Arimoto algorithm of the structure given in Algorithm 1. We define the following two variable extension of $K(\rho_{A},\mathcal{E})$ for $\sigma\gg\rho$ and $\gamma\in[1-\eta_{\rm exp}(\mathcal{E}),1]$, where $K(\rho_{A},\mathcal{E})$ is defined in (51) and $\eta_{\rm exp}(\mathcal{E})$ is defined in (50).

[TABLE]

where we omitted the system indices for simplicity. Note that $\gamma=1$ corresponds to the standard algorithm and if it is known that $\eta_{\rm exp}(\mathcal{E})>0$, then one may accelerate the algorithm. With a short calculation, one can bring $J_{\gamma}$ into the form

[TABLE]

Note that $\mathrm{Tr}\left[\rho\mathcal{F}(\sigma)\right]$ is continuous for $\sigma\gg\rho$. Moreover, another short calculation leads to

[TABLE]

where we used the data processing inequality for the quantum relative entropy (together with $\gamma\geq 1-\eta_{\rm exp}(\mathcal{E})$). From Lemma III.2 whose requirement is satisfied by (59), we find a double optimization form of the thermodynamic capacity $T_{\Gamma}(\mathcal{E})$

[TABLE]

Performing the two maximizations in $\max_{\rho,\sigma}J_{\gamma}(\rho,\sigma,\mathcal{E})$ iteratively, leads to the following algorithm (see Lemma III.2 for the form of the update rules), which provably converges to the thermodynamic capacity by Theorem III.3. A detailed analysis of the time complexity can be found in Appendix B.2.

When the coefficient $\eta_{\rm con}(\mathcal{E})<1$ (from (40)), we get exponential convergence. To see this, note

[TABLE]

Using (58), we see that we have exponential convergence according to Proposition III.6 with $\delta=\frac{1-\eta_{\rm con}(\mathcal{E})}{\gamma}$ when the coefficient $\eta_{\rm con}(\mathcal{E})<1$.

V.3 Simulation results

For unital channels, the thermodynamic capacity is zero and the maximizer is the maximally mixed state, which we have chosen as our initial guess $\rho_{\lambda}^{(1)}$. Here, we consider the non-unital qubit amplitude damping channel which we recall from (4)

[TABLE]

As in Section IV.3, we choose an additive error $\varepsilon=10^{-6}$ bits. We use the termination criterion given in Proposition III.5, which significantly reduces the required number of iterations of the algorithm. The results are shown in Figure 4. The standard and the accelerated algorithm converge in $25$ and $5$ iterations, respectively.

VI Coherent information of less noisy channels

VI.1 Definitions

For a quantum channel $\mathcal{E}$, the coherent information of the channel, $I_{\rm{coh}}(\mathcal{E})$, is defined as the maximum of the coherent information $I_{\rm{coh}}(\rho,\mathcal{E})=S(\mathcal{E}(\rho))-S(\mathcal{E}_{c}(\rho))$, where the maximization is taken over input states $\rho$. That is, we have

[TABLE]

A channel $\mathcal{E}$ is called less noisy if the private capacity of its complementary channel $\mathcal{E}_{c}$ is zero watanabe2012private . In (watanabe2012private, , Proposition 4), it is shown that a channel is less noisy if and only if it is less divergence contracting, which implies that for all density operators $\rho$ and $\sigma$, we have

[TABLE]

The quantum capacity $Q(\mathcal{E})$ of the channel $\mathcal{E}$ is the regularized coherent information of a channel devetak2005private and corresponds to the rate at which Alice can transmit quantum information to Bob. In watanabe2012private it is shown that the coherent information of a channel is additive if the channel is less noisy. Hence, for less noisy channels, we have $Q(\mathcal{E})=I_{\rm{coh}}(\mathcal{E})$. We shall only consider less noisy channels henceforth and we therefore drop the distinction between the regularized and single letter formulas.

VI.2 Blahut-Arimoto algorithm

The coherent information of a quantum channel can be approximated using a Blahut-Arimoto algorithm of the structure given in Algorithm 1. We may also accelerate the algorithm and to do so, we introduce the following coefficient.

[TABLE]

We define the following two variable extension of $I_{\rm{coh}}(\rho,\mathcal{E})$ for $\sigma\gg\rho$ and $\gamma\in[\zeta_{\rm con},1]$.

[TABLE]

With a short calculation, one can bring $J_{\gamma}$ into the following form

[TABLE]

Note that $\gamma=1$ corresponds to the standard algorithm and if $\zeta_{\rm con}(\mathcal{E})<1$, then one may accelerate the algorithm. We also note that $\mathrm{Tr}\left[\rho\mathcal{F}(\sigma)\right]$ is continuous for $\sigma\gg\rho$. Moreover, another short calculation leads to

[TABLE]

where we have used (65) along with $\gamma\geq\zeta_{\rm con}(\mathcal{E})$. From Lemma III.2, whose requirement is satisfied by (71), we find a double optimization form of the coherent information $I_{\rm{coh}}(\mathcal{E})$

[TABLE]

Performing the two maximizations in $\max_{\rho,\sigma}J_{\gamma}(\rho,\sigma,\mathcal{E})$ iteratively, leads to the following algorithm (see Lemma III.2 for the form of the update rules), which provably converges to the coherent information by Theorem III.3. A detailed analysis of the time complexity is given in Appendix B.3.

We also have exponential convergence if the channel is strictly less noisy i.e. $\inf D(\mathcal{E}(\rho)||\mathcal{E}(\sigma))-D(\mathcal{E}_{c}(\rho)||\mathcal{E}_{c}(\sigma))>0$ for all $\sigma\gg\rho$ and $\sigma\neq\rho$. To see this, we define

[TABLE]

In general, $\zeta_{\rm exp}(\mathcal{E})\in[-1,1]$ but for less noisy channels, we have that $\zeta_{\rm exp}(\mathcal{E})\in[0,1]$ due to (65). Comparing it to (66), we also note that $\zeta_{\rm exp}(\mathcal{E})=-\zeta_{\rm con}(\mathcal{E}_{c})$. Using (70), we see that when the coefficient $\zeta_{\rm exp}(\mathcal{E})>0$, we get exponential convergence according to Proposition III.6 with $\delta=\frac{\zeta_{\rm exp}(\mathcal{E})}{\gamma}$. Indeed, if the channel is strictly less noisy, the condition $\zeta_{\rm exp}(\mathcal{E})>0$ follows.

VI.3 Simulation results

Once again we consider the qubit amplitude damping channel whose form is given in (63). We choose an additive error $\varepsilon=10^{-6}$ bits and use the termination criterion given in Proposition III.5 which significantly reduces the number of iterations of the algorithm required to achieve an estimate of the capacity with additive error smaller than $\varepsilon$. The results are shown in Figure 5. The standard algorithm and the adaptive accelerated algorithm converge in $24$ and $5$ iterations, respectively.

VII Mutual information of quantum channels

VII.1 Definitions

The mutual information $I(\mathcal{E})$ of a quantum channel $\mathcal{E}$ is defined as the maximum over input states $\rho$ of the mutual information $I(\rho,\mathcal{E})=S(\rho)+S(\mathcal{E}(\rho))-S(\mathcal{E}_{c}(\rho))$, i.e.,

[TABLE]

The entanglement assisted classical capacity $C(\mathcal{E})$ of the channel $\mathcal{E}$ is the regularized mutual information of a channel Bennett02 . It is also known that the mutual information of a channel is additive wildebook17 and hence we have $C(\mathcal{E})=I(\mathcal{E})$.

VII.2 Blahut-Arimoto algorithm

The mutual information of a quantum channel can be approximated using a Blahut-Arimoto algorithm of the structure given in Algorithm 1. First, recall the definition of $\zeta_{\rm con}$ given in (66). For any $\gamma\in[1+\zeta_{\rm con}(\mathcal{E}),2]$, we define the following two variable extension of $I(\rho,\mathcal{E})$ for $\sigma\gg\rho$

[TABLE]

Note that $\gamma=2$ corresponds to the standard algorithm. If it is known that $\zeta_{\rm con}(\mathcal{E})<1$, then one may accelerate the algorithm. With a short calculation, one can bring $J_{\gamma}$ into the form

[TABLE]

Note that $\mathrm{Tr}\left[\rho\mathcal{F}(\sigma)\right]$ is continuous for $\sigma\gg\rho$. Moreover, a short calculation leads to

[TABLE]

where we have used the data processing inequality and the nonnegativity of relative entropy along with $\gamma\geq 1+\zeta_{\rm con}(\mathcal{E})$. From Lemma III.2, whose requirement is satisfied by (79), we find a double optimization form of the coherent information $I(\mathcal{E})$

[TABLE]

Performing the two maximizations in $\max_{\rho,\sigma}J_{\gamma}(\rho,\sigma,\mathcal{E})$ iteratively, leads to the following algorithm (see Lemma III.2 for the form of the update rules), which provably converges to the mutual information of the channel $\mathcal{E}$ by Theorem III.3. A detailed analysis of the time complexity is given in Appendix B.4.

We also have exponential convergence under certain conditions. Recall the definition $\zeta_{\rm exp}(\mathcal{E})$ from (73). Using (78), we see that when the coefficient $\zeta_{\rm exp}(\mathcal{E})>-1$ (or equivalently, $\zeta_{\rm con}(\mathcal{E}_{c})<1$), we get exponential convergence according to Proposition III.6 with $\delta=\frac{1+\zeta_{\rm exp}(\mathcal{E})}{\gamma}$.

VII.3 Simulation results

Once again we consider the qubit amplitude damping channel which has the form as in (63). We choose an additive error $\varepsilon=10^{-6}$ bits and use the termination criterion given in Proposition III.5, which significantly reduces the number of iterations of the algorithm required to achieve an estimate of the capacity with additive error smaller than $\varepsilon$. The results are shown in Figure 6. The standard and the adaptive accelerated algorithm converge to the capacity in $12$ and $4$ iterations, respectively.

VIII Conclusion

We give an analytical and numerical analysis of alternating optimization algorithms of Blahut-Arimoto type for computing various channel capacities in quantum information theory. We note that our algorithms are of zeroth-order and do not need to take into account matrix valued derivatives. Thus, they are rather straightforward and computationally inexpensive to implement — as demonstrated in our numerical examples. It remains open if alternating optimization algorithms, in particular of Blahut-Arimoto type, can also be given for other convex optimization problems in terms of quantum entropy. Starting from the analogous classical settings this could, e.g., be quantum rate distortion functions Datta13 , quantum information bottleneck functions hirche18 , or certain quantum network capacities Winter01 . For example, in the case of the classical information bottleneck, a Blahut-Arimoto type algorithm is used for the estimation of the information bottleneck functions Tishby99 and helped to gain insight into the black box of deep neural networks tishby_deep_2015 ; shwartz-ziv_opening_2017 . The quantum information bottleneck was investigated in Watanabe16 ; salek_quantum_2019 ; hirche18 .

In contrast to classical Shannon theory, in quantum Shannon theory exact capacity formulas are often not known wildebook17 . Moreover, known upper and lower bounds are not always in the form of convex optimization problems. For example, the complexity of determining the classical capacity of general entanglement breaking channels is NP-complete Beigi07 (see also Harrow13 for hardness of approximation results). For such cases, one might rather aim for numerical tools that do well for the average case in practically relevant examples. Alternating optimization algorithms offer an interesting option in this direction. In fact, Nagaoka explored a version of his quantum Blahut-Arimoto algorithm to study the classical capacity of general quantum channels Nagaoka01 .

Acknowledgements.

We thank Joseph Renes and Marco Tomamichel for discussions. RI acknowledges support from the Swiss National Science Foundation through SNSF project No. 200020-165843 and through the National Centre of Competence in Research Quantum Science and Technology (QSIT).

Note added.

During finalization of the first version of our work we became aware of the related works Cai19 ; li2019computing by Li and Cai.

Appendix A Gibbs’ variational principle

The well-known Gibbs’ variational principle is given in the following lemma (which is used to show the form of the optimizer in Lemma III.2).

Lemma A.1** (Gibbs’ principle)**

Let $\omega$ a density operator and $H$ a Hermitian matrix on the same space. Then, we have

[TABLE]

with equality if and only if $\omega=\frac{\exp\left(H\right)}{\mathrm{Tr}\left[\exp\left(H\right)\right]}$.

Appendix B Time complexity

The time complexity of Blahut-Arimoto type algorithms is essentially given by the required number of iteration steps times the complexity of one iteration step. In the following, we give a detailed analysis for the computation of the Holevo quantity, the thermodynamic capacity, the coherent information and the mutual information.

B.1 Holevo quantity

The time complexity of Algorithm 2 is determined by the required number of iterations and the time complexity for applying the following update rule

[TABLE]

Let us sketch how to apply the update rule efficiently to derive an upper bound on the time complexity. We use the following:

• •

the channel $\mathcal{E}_{X\rightarrow B}$ is given as a lookup table, where one can access the elements $\tau_{\mathcal{E},i}$ in constant time. The application of $\mathcal{E}$ to a diagonal density operator $\rho=\sum_{i=1}^{|X|}\rho_{ii}|i\rangle\langle i|$ can then be calculated as $\mathcal{E}(\rho)=\sum_{i=1}^{|X|}\rho_{ii}\tau_{\mathcal{E},i}$ with time complexity $\mathcal{O}\left(|X||B|^{2}\right)$,

• •

the complexity of calculating the matrix logarithm of $\rho_{A}$ is $\mathcal{O}(|A|^{3})$,

• •

the complexity of calculating $\mathrm{Tr}\left[\rho_{A}\sigma_{A}\right]$ is $\mathcal{O}(|A|^{2})$.

The algorithm proceeds then with the following calculations:

The terms $c_{i}=\mathrm{Tr}\left[\tau_{\mathcal{E},i}\log\tau_{\mathcal{E},i}\right]$ can be computed once for all $i\in\{1,2,\dots,|X|\}$ at the start of the algorithm with complexity $\mathcal{O}(|X||B|^{3})$. 2. 2.

In each iteration step $t$ we have to compute:

1. (a)

the matrix $H^{(t)}=\log\mathcal{E}(\rho_{\lambda^{(t)}})$ requiring time $\mathcal{O}(|X||B|^{2}+|B|^{3})$, 2. (b)

real numbers $d^{(t)}_{i}=\mathrm{Tr}\left[\tau_{\mathcal{E},i}\,H^{(t)}\right]$ for all $i\in\{1,2,\dots,|X|\}$ requiring time $\mathcal{O}(|X||B|^{2})$, 3. (c)

the unnormalized coefficients $\tilde{\lambda}^{(t+1)}_{i}={\lambda}^{(t)}_{i}\exp\left(\frac{1}{\gamma}(c_{i}-d_{i}^{(t)})\right)$ requiring time $\mathcal{O}(|X|)$, 4. (d)

the normalization coefficient $Z^{(t+1)}=\sum_{i=1}^{|X|}\tilde{\lambda}^{(t+1)}_{i}$ requiring time $\mathcal{O}(|X|)$.

We conclude that the complexity for one iteration step is $\mathcal{O}(|B|^{3}+|X||B|^{2})$.

The number of required iterations to get an $\varepsilon$-approximation to the capacity is of order $\mathcal{O}(\log|X|/\varepsilon)$ and hence the required time for all the iteration steps is $\mathcal{O}\left(\left(|B|^{3}+|B|^{2}|X|\right)\log|X|/\varepsilon\right)$. We conclude that the time complexity of Algorithm 2 is given by

[TABLE]

B.2 Thermodynamic capacity

The time complexity of Algorithm 3 is determined by the required number of iterations and the complexity of applying the following update rule

[TABLE]

where $\mathcal{F}(\sigma)=-\log\sigma+\mathcal{E}^{\dagger}(\log\mathcal{E}(\sigma))+\log\Gamma_{A}-\mathcal{E}^{\dagger}(\log\Gamma_{B})$ and $Z^{(t+1)}=\mathrm{Tr}\left[\exp\left(\log\rho^{(t)}+\frac{1}{\gamma}\mathcal{F}(\rho^{(t)})\right)\right]$ normalizes the state. To calculate the complexity of applying the update rule, we use that

• •

the quantum channel $\mathcal{E}_{A\rightarrow B}$ and its adjoint $\mathcal{E}_{B\rightarrow A}^{\dagger}$ are both given as lookup tables whose $(i,j)$-th entry is given by $\mathcal{E}(|i\rangle\langle j|_{A})$ or $\mathcal{E}^{\dagger}(|i\rangle\langle j|_{B})$, respectively. We assume constant time access to the entries of the table. Hence, the application of $\mathcal{E}$ (or $\mathcal{E}^{\dagger}$) to a density operator $\rho_{A}$ has time complexity $\mathcal{O}\left(|A|^{2}|B|^{2}\right)$. Indeed, the channel application can be calculated as $\mathcal{E}(\rho)=\sum_{i,j=1}^{|A|}\rho_{ij}\tau_{\mathcal{E},i,j}$ with $\tau_{\mathcal{E},i,j}=\mathcal{E}(|i\rangle\langle j|_{A})$,

• •

the complexity of calculating the matrix logarithm and exponential of $\rho_{A}$ is $\mathcal{O}(|A|^{3})$,

An iteration step of the form given in (85) is then found to have complexity $\mathcal{O}(|A|^{3}+|A|^{2}|B|^{2}+|B|^{3})$. The number of required iterations to get an $\varepsilon$ approximation to the capacity is of order $\mathcal{O}(\log|A|/\varepsilon)$ and hence we conclude that the time complexity of Algorithm 3 is given by

[TABLE]

B.3 Coherent information

The time complexity of Algorithm 4 is determined by the required number of iterations and the complexity of applying the following update rule

[TABLE]

where $\mathcal{F}(\sigma)=\mathcal{E}_{c}^{\dagger}\log\mathcal{E}_{c}(\sigma)-\mathcal{E}^{\dagger}\log\mathcal{E}(\sigma)$ and $Z^{(t+1)}=\mathrm{Tr}\left[\exp\left(\log\rho^{(t)}+\frac{1}{\gamma}\mathcal{F}(\rho^{(t)}))\right)\right]$ normalizes the state. To calculate the complexity of applying the update rule, we use that

• •

the application of the quantum channel $\mathcal{E}_{A\rightarrow B}$ and of its adjoint have time complexity $\mathcal{O}\left(|A|^{2}|B|^{2}\right)$, where the application of the channel $(\mathcal{E}_{c})_{A\rightarrow E}$ and of its adjoint have time complexity $\mathcal{O}\left(|A|^{2}|E|^{2}\right)$ (see Appendix B.2 for the argumentation),

• •

the complexity of calculating the matrix logarithm and exponential of a density operator $\rho_{A}$ is $\mathcal{O}(|A|^{3})$.

An iteration step of the form given in (87) is then found to have complexity $\mathcal{O}(|A|^{3}+|B|^{3}+|E|^{3}+|A|^{2}(|B|^{2}+|E|^{2}))$. By Stinespring’s factorization theorem, the dimension of the environment can be bounded by $|E|\leq K$, where $K\leq|A||B|$ denotes the Kraus rank of $\mathcal{E}$, i.e., the minimal number of Kraus operators required to represent $\mathcal{E}$. Hence, the complexity of one iteration step is $\mathcal{O}(|A|^{3}+|B|^{3}+K^{3}+|A|^{2}(|B|^{2}+K^{2}))$. The number of required iterations to get an $\varepsilon$ approximation to the capacity is of order $\mathcal{O}(\log|A|/\varepsilon)$ and hence we conclude that the time complexity of Algorithm 4 is given by

[TABLE]

B.4 Mutual information

The time complexity of Algorithm 5 is determined by the required number of iterations and the complexity of applying the following update rule

[TABLE]

where $\mathcal{F}(\sigma)=-\log(\sigma)+\mathcal{E}_{c}^{\dagger}\log\mathcal{E}_{c}(\sigma)-\mathcal{E}^{\dagger}\log\mathcal{E}(\sigma)$ and $Z^{(t+1)}=\mathrm{Tr}\left[\exp\left(\log\rho^{(t)}+\frac{1}{\gamma}\mathcal{F}(\rho^{(t)}))\right)\right]$ normalizes the state. Note that the complexity of the update rule is the same as in Appendix B.3 since the additional term $\log(\sigma)$ in $\mathcal{F}{(\sigma)}$ does not change the asymptotic complexity. We conclude that the time complexity of Algorithm 5 is given by

[TABLE]