Computing the Classical-Quantum channel capacity: experiments on a Blahut-Arimoto type algorithm and an approximate solution for the binary inputs, two-dimensional outputs channel

TL;DR

This paper demonstrates a faster convergence of a Blahut-Arimoto type algorithm for classical-quantum channel capacity and introduces an explicit approximate solution for a specific binary-input, two-dimensional-output channel.

Contribution

It provides empirical evidence of accelerated convergence and proposes a new approximate solution for a particular class of classical-quantum channels.

Findings

Algorithm converges faster than predicted

Explicit approximate solution with minimal error

Applicable to binary-input, two-dimensional-output channels

Abstract

In our previous work, we presented a Blahut-Arimoto type algorithm for computing the discrete memoryless (DM) classical-quantum channel capacity. And the speed of convergence is analyzed. In this paper, we present numerical experiment to show that the algorithm converged much faster than the theoretical prediction. Also we present an explicit approximate solution for the binary inputs, two dimensional outputs classical-quantum channel, which has an error within $3\times10^{-4}$

Tables

Table 1. Table 1:

Input $n$	Output $m$	Accuracy	Ave iteration	Max ite ($\times {10}^2$)
2	2	${10}^{-3}$	7	0.35
		${10}^{-4}$	25	0.12
		${10}^{-5}$	64	0.25
2	5	${10}^{-3}$	11	0.27
		${10}^{-4}$	24	0.11
		${10}^{-5}$	44	0.13
2	8	${10}^{-3}$	11	0.25
		${10}^{-4}$	24	0.70
		${10}^{-5}$	39	0.11
5	2	${10}^{-3}$	125	5.0
		${10}^{-4}$	245	4.5
		${10}^{-5}$	494	23
5	5	${10}^{-3}$	95	4.5
		${10}^{-4}$	243	10
		${10}^{-5}$	456	20
5	8	${10}^{-3}$	85	0.17
		${10}^{-4}$	212	0.60
		${10}^{-5}$	414	20
8	2	${10}^{-3}$	102	4.1
		${10}^{-4}$	219	12
		${10}^{-5}$	353	23
8	5	${10}^{-3}$	114	4.4
		${10}^{-4}$	275	13
		${10}^{-5}$	456	27
8	8	${10}^{-3}$	112	3.0
		${10}^{-4}$	291	9.8
		${10}^{-5}$	582	25

Equations

$$\chi (\{p_x:{\rho }_x{\}}_{x\in \mathcal{X}})=H(\sum _x{p}_x{\rho }_x)-\sum _x{p}_{x}H({\rho }_x),$$

$$H(\rho )=-\mathrm{Tr}(\rho \log \rho ).$$

$$\underset{p}{\max }H(\sum _x{p}_x{\rho }_x)-\sum _x{p}_{x}H({\rho }_x),$$

$$\displaystyle subject\ to\ \$$

$$p\in {\Delta }_n,$$

$$\chi (t)=H(\sum _x{p}_x^t{\rho }_x)-\sum _x{p}_x^{t}H({\rho }_x^t).$$

$$\sum _x{c}_x{\rho }_x=0.$$

$$D(p^{\ast }||p^t)\le (1-\delta {)}^{t-N_0}D(p^{\ast }||p^{N_0}).$$

$$\underset{t\to \infty }{\lim }\frac{p_x^{t+1}}{p_x^t}=\exp [D({\rho }_x||{\rho }^{\ast })-C],$$

$$D({\rho }_x||{\rho }^{\ast })\le C$$

$$\underset{x}{\max }\{D({\rho }_x||{\rho }^t)\}\to C.$$

$$\underset{x}{\max }$$

$$=\sum _x{p}_x^{\ast }D({\rho }_x||{\rho }^{\ast })+D({\rho }^{\ast }||{\rho }^t)$$

$$=C+D({\rho }^{\ast }||{\rho }^t)\ge C.$$

$$\{p_1:{\rho }_1;p_2:{\rho }_2\},p_1+p_2=1,{\rho }_1,{\rho }_2\in {\mathcal{D}}^2.$$

$$\cos \frac{\alpha }{2}|0⟩+\sin \frac{\alpha }{2}{e}^{i\varphi }|1⟩,$$

$$\chi ({\lambda }_1,{\lambda }_2,\theta ,p_1)=$$

$$-[p_{1}S(\frac{1}{2}+r_1)+(1-p_1)S(\frac{1}{2}+r_2)],$$

$${\nabla }_{p_1}\chi ({\lambda }_1,{\lambda }_2,\theta ,p_1).$$

$${\nabla }_{p_1}\chi ({\lambda }_1,{\lambda }_2,\theta =0,p_1)=0,$$

$${\hat{p}}_1=\frac{\frac{1}{2}\frac{1-c}{1+c}-r_2}{r_1-r_2},wherec=2^{\frac{S({\lambda }_1)-S({\lambda }_2)}{r_1-r_2}},$$

$${\rho }_1={\lambda }_1|0⟩⟨0|+(1-{\lambda }_1)|1⟩⟨1|;$$

$${\rho }_2={\lambda }_2|{\psi }_0⟩⟨{\psi }_0|+(1-{\lambda }_2)|{\psi }_1⟩⟨{\psi }_1|,$$

$$|{\psi }_0⟩=\cos \frac{\theta }{2}|0⟩+\sin \frac{\theta }{2}|1⟩;$$

$$|{\psi }_1⟩=-\sin \frac{\theta }{2}|0⟩+\cos \frac{\theta }{2}|1⟩,$$

$$Error({\lambda }_1,{\lambda }_2)=\underset{\theta \in [0,\pi ]}{\max }|\chi ({\lambda }_1,{\lambda }_2,\theta ,{\hat{p}}_1)-{\chi }_{BA}({\lambda }_1,{\lambda }_2,\theta )|.$$

$$\underset{{\lambda }_1,{\lambda }_2\in [0.5,R]}{\max }Error({\lambda }_1,{\lambda }_2).$$

Full text

• •

1 Introduction

The classical-quantum channel [2] can be considered as consisting of an input alphabet $\mathcal{X}=\{1,2,\dots,|\mathcal{X}|\}$ and a mapping $x\rightarrow\rho_{x}$ from the input alphabet to a set of quantum states in a finite dimensional Hilbert space $\mathcal{H}$. The state of a quantum system is given by a density operator $\rho$, which is a positive semi-definite operator with trace equal to one. Let $\mathcal{D}^{m}$ denote the set of all density operators acting on a Hilbert space $\mathcal{H}$ of dimension $m$. If the resource emits a letter $x$ with probability $p_{x}$, the output would be $\rho_{x}$, so the output would form an ensemble: $\{p_{x}:\rho_{x}\}_{x\in\mathcal{X}}$.

In 1998, Holevo showed [3] that the classical capacity of the classical-quantum channel is the maximization of a quantity called the Holevo information over all input distributions. The Holevo information $\chi$ of an ensemble $\{p_{x}:\rho_{x}\}_{x\in\mathcal{X}}$ is defined as

[TABLE]

where $H(\cdot)$ is the von Neumann entropy which is defined on positive semidefinite matrices:

[TABLE]

Due to the concavity of von Neumann entropy [4], the Holevo information is always non-negative. The Holevo quantity is concave in the input distribution [4], and so the maximization of equation (1) over $p$ is a convex optimization problem. However it is not a straightforward convex optimization problem.

For discrete memoryless classical channels, the capacity can be computed efficiently by using an algorithm called Blahut-Arimoto (BA) algorithm [5][6][7]. In 1998, H. Nagaoka [8] proposed a quantum version of BA algorithm. In his work he considered the quantum-quantum channel and this problem was proved to be NP-complete [9]. And Nagaoka mentioned an algorithm concerning classical-quantum channel, however, its speed of convergence was not studied there. In 2014, Davide Sutter et al. [10] promoted an algorithm based on duality of convex programing and smoothing techniques [11] with a complexity of $O(\frac{(n\vee m)m^{3}(\log n)^{1/2}}{\varepsilon})$, where $n\vee m=\max\{n,m\}$.

The structure of this paper is organized as following: In Section 2, we summarize the results of our previous work: the problem and the algorithm. In Section 3 we show the numerical experiments of BA algorithm to see how well this algorithm performs. In Section 4, we propose an approximate solution for a special case, which is the binary inputs, two dimensional outputs channel.

$\bf{Notations}$: The logarithm $\log(\cdot)$ is with basis $2$ unless specially specified. The space of all Hermitian operators of dimension $m$ is denoted by $H^{m}$. The set of all density matrices of dimension $m$ is denoted by $\mathcal{D}^{m}=\{\rho\in H^{m}:\rho\geq 0,\mathop{\mathrm{Tr}}\rho=1\}$. Each letter $x\in\mathcal{X}$ is mapped to a density matrix $\rho_{x}$ so the classical-quantum channel can be represented as a set of density matrices $\{\rho_{x}\}_{x\in\mathcal{X}}$. The set of all probability distributions of length $n$ is denoted by $\Delta_{n}=\{p:p_{x}\geq 0,\sum_{x=1}^{n}p_{x}=1\}$. The von Neumann entropy of a density matrix $\rho$ is denoted by $H(\rho)=-\mathop{\mathrm{Tr}}[\rho\log\rho]$. The relative entropy between $p,q\in\Delta_{n}$, if ${supp}(p)\subset{supp}(q)$, is denoted by $D(p||q)=\sum_{x}p_{x}(\log p_{x}-\log q_{x})$ and $+\infty$ otherwise. The relative entropy between $\rho,\sigma\in\mathcal{D}^{m}$, if ${supp}(\rho)\subset{supp}(\sigma)$, is denoted by $D(\rho||\sigma)=\mathop{\mathrm{Tr}}[\rho(\log\rho-\log\sigma)]$ and $+\infty$ otherwise.

2 The BA algorithm and its convergence

We want to compute the capacity of a given classical-quantum channel $\{x\rightarrow\rho_{x}\}_{x\in\mathcal{X}}$, which can be written as

[TABLE]

where $\rho_{x}\in\mathcal{D}^{m}$, $|s\rangle\geq 0,S>0$. We denote the maximal value of equation (3) as $C(S)$. In this optimization problem, we are to maximize the Holevo quantity with respect to the input distribution $\{p_{x}\}_{x\in\mathcal{X}}$. Practically, the preparation of different signal state $x$ has different cost, which is represented by $|s\rangle$. And we would like to bound the expected cost of the resource within some quantity, which is represented by the inequality constraint in equation (4)

In this paper we only present the algorithm which does not concern the inequality constraint (4). (The algorithms, with or without inequality constraint (4), are essentially the same.) The algorithm we presented at [1] can be summarized as following: (we take the natural logarithm for conveniece)

And we have the following theorems: (For detailed derivation and proof, please refer to [1].)

Definition 1**.**

Define the value of Holevo quantity (1) of a certain channel after $t$ iterations as

[TABLE]

Theorem 1**.**

$\chi(t)$* converges to the capacity of the classical-quantum channel from below as $t\rightarrow\infty$; and the input distribution $p^{t}$ converges to an optimal distribution $p^{*}$ as $t\rightarrow\infty$.*

Theorem 2**.**

To reach $\epsilon$ accuracy to the capacity, the BA algorithm needs an iteration complexity less than $\frac{\log n}{\epsilon}$.

If the classical-quantum channel has some extra feature, we have a better convergence performance:

Assumption 1**.**

The channel matrices $\{\rho_{x}\}_{x\in\mathcal{X}}$ are linearly independent in the complex matrix space, i.e. there doesn’t exist a vector $|c\rangle$ such that

[TABLE]

Theorem 3**.**

Under Assumption 1, the optimal solution $p^{*}$ is unique. And $p^{t}$ converges to $p^{*}$ at a geometric speed, i.e. there exist $N_{0}$ and $\delta>0$, where $N$ and $\delta$ are independent, such that for any $t>N_{0}$, we have

[TABLE]

3 Numerical experiments on BA algorithm

We will study the relations between iteration complexity and $n,m$ (i.e. the input size and output dimension) when the algorithm reaches certain accuracy. Due to we don’t know the true capacity of a certain channel, we will use the following theorem to bound the error of the algorithm.

Theorem 4**.**

With the iteration procedure in the BA algorithm 1, $\max_{x}\{D(\rho_{x}||\rho^{t})\}$ converges to $F(\lambda)$ from above.

Proof.

Followed from Algorithm 1, and Theorem 1, with some calculation, we have

[TABLE]

where $\rho^{*}=\sum_{x}p^{*}_{x}\rho_{x}$, $p^{*}$ is the optimal distribution that $p^{t}$ converges to, and $C$ is the capacity of the channel. The limit above is $1$ if $p^{*}_{x}>0$ and does not exceed $1$ if $p^{*}_{x}=0$. So

[TABLE]

for every $x\in\mathcal{X}$, with equality if $p^{*}_{x}>0$. This proves

[TABLE]

And for any $p^{t}$ and any optimal distribution $p^{*}$, we have

[TABLE]

The first equality requires some calculation and the second equality follows since $p^{*}$ is an optimal distribution. And the last inequality follows from the non-negativity of mutual information. This means $\max_{x}\{D(\rho_{x}||\rho^{t})\}$ converges to $C$ from above. ∎

So our accuracy criterion is: for a given classical-quantum channel, we run the BA algorithm (with no input constraint), until $[\max_{x}\{D(\rho_{x}||\rho^{t})\}-\chi(t)]$ is less than $10^{-k}$, and record the number of iteration. At this time, the accuracy is of order $10^{-(k+1)}$ since $\max_{x}\{D(\rho_{x}||\rho^{t})\}$ and $\chi(t)$ converges to the true capacity from above and below respectively.

We do the following numerical experiments: for given values of input size $n$ , output dimension $m$ and accuracy, we generate $200$ classical-quantum channels randomly, and record the numbers of iterations then calculate the average number of iterations and find the maximum number of iterations in these $200$ experiments. The results are shown in Table 1. Note that the accuracy $10^{-k}$ in Table 1 means we run the BA algorithm until $[\max_{x}\{D(\rho_{x}||\rho^{t})\}-\chi(t)]$ is less than $10^{-k}$, and the error between the true capacity and the computed value is of order $10^{-(k+1)}$.

We can see from the table that the iteration complexity scales good as accuracy increases. Notice that the datas of $n=5$ and $n=8$ are vary similar which means the iteration complexity also scales very good as the input size $n$ increases. We can also see for given input size $n$ and accuracy, the output dimension has vary little influence on iteration complexity, which means the iteration complexity also scales good as the output dimension $m$ increases. Compared with our theoretical analysis of iteration complexity in Theorem 2: to reach $\epsilon$ accuracy, we need $\frac{\log n}{\epsilon}$ number of iterations. numerical experiments show that both the average and maximum number of iterations are far smaller than $\frac{logn}{\epsilon}$ to reach $\epsilon$ accuracy, no matter whether the output quantum states are linearly independent (cases in $(n,m)=(5,2),(8,2)$) or not . The reason of this phenomenon is that the inequalities we used in the proof of Theorem 2 are quite loose. So Theorem 2 only provide a very loose upper bound on iteration complexity. In [10], only the binary two dimensional case is put into numerical experiments and the iteration complexity is of order $1/\varepsilon$ to reach $\varepsilon$ accuracy. So compared with the results in [10], our algorithm is much better both theoretically and practically (Please refer to [1] for the theoretical analysis).

4 An approximate solution of $p$ in binary two dimensional case

In this section we provide an approximate optimal input distribution for the case of binary inputs, two dimensional outputs channel:

[TABLE]

4.1 Use Bloch sphere to get an approximate solution

Any two-dimensional density matrix can be represented as a point in the Bloch sphere [4], as shown in the following:

Any density matrix can be represented as a vector in the Bloch sphere starting from the origin. Suppose $\rho_{1},\rho_{2}$ can be represented as $\bm{r_{1}},\bm{r_{2}}$ respectively, then the two eigenvalues would be $0.5\pm r_{1}/2$ and $0.5\pm r_{2}/2$ respectively. Extend $\bm{r_{1}}$ we get two intersections on the surface of the Bloch sphere, then these two intersections represents the two eigenvectors of $\rho_{1}$ (the points on the surface of the sphere represent pure state and the interior points represent mixed states.) A probabilistic combination of $\rho_{1},\rho_{2}$ can be represented as $p_{1}\rho_{1}+p_{2}\rho_{2}=p_{1}\bm{r_{1}}+p_{2}\bm{r_{2}}$ ([4] Exercise 4.4.13). And any point on the surface of Bloch sphere can be represented as

[TABLE]

where $\alpha$ is the angle to the $Z$ axis and $\phi$ is the angle of the $X$ axis to the projection of the point on the $X-Y$ plane.

By symmetry, it is obvious that the Holevo quantity is only related to $r_{1},r_{2},\theta,p_{1}$, where $\theta$ is the angle between $\bm{r_{1}},\bm{r_{2}}$. One interesting result is that the angle $\theta$ has very little influence on $p^{*}$, where $p^{*}$ is the optimal distribution that maximizes Holevo quantity. If we know $\lambda_{1},\lambda_{2}$ (the bigger eigenvalues of $\rho_{1},\rho_{2}$ respectively), $\theta$ and $p_{1}$, then the Holevo quantity can be written as

[TABLE]

where $S(\cdot)$ is the binary entropy ($S(x)=-(x\log x+(1-x)\log(1-x))$ and $r_{i}=\lambda_{i}-\frac{1}{2},i=1,2$.

Using Cosine theorem to calculate $||p_{1}\bm{r_{1}}+(1-p_{1})\bm{r_{2}}||_{2})$, the gradient of $\chi$ w.r.t. $p_{1}$ can be calculated directly, denoted as

[TABLE]

If we can find a $\hat{p}_{1}$ such that $\nabla_{p_{1}}\chi(\lambda_{1},\lambda_{2},\theta,p_{1})|_{p_{1}=\hat{p}_{1}}=0$, then this $\hat{p}_{1}$ is the optimal solution (because $\chi(\lambda_{1},\lambda_{2},\theta,p_{1})$ is concave in $p_{1}$). However, we cannot solve the equation $\nabla_{p_{1}}\chi(\lambda_{1},\lambda_{2},\theta,p_{1})=0$ w.r.t. $p_{1}$ when $\theta\neq 0$. Now that $\theta$ has little influence on the optimal distribution $p^{*}$, let $\theta=0$ (this is actually the classical case), and let

[TABLE]

the above equation is easy to solve and we get a solution $\hat{p}_{1}$:

[TABLE]

where we assume $r_{1}\neq r_{2}$. (It can be easily seen from the Bloch sphere that if $r_{1}=r_{2}$, the optimal distribution would be $\{\frac{1}{2},\frac{1}{2}\}$.)

This $\hat{p}_{1}$ can be used as an approximate optimal solution. Next we need numerical experiments to see how accurate $\hat{p}_{1}$ is.

4.2 Numerical experiments on the approximated solution $\hat{p}_{1}$

It is obvious that the maximum of Holevo quantity only depends on $r_{1},r_{2}$ and $\theta$, so without loss of generality, we let $\rho_{1}$ be on the $Z$ axis and $\rho_{2}$ be on the $X-Z$ plane:

[TABLE]

where

[TABLE]

which means the angle between $\rho_{1}$ and $\rho_{2}$ (i.e. $\bm{r_{1}}$ and $\bm{r_{2}}$) is $\theta$.

In the numerical experiments, we let $\lambda_{1},\lambda_{2}$ range from $0.5$ to $1$, and $\theta$ ranges from [math] to $\pi$. For each value of $(\lambda_{1},\lambda_{2},\theta)$, we substitute $(\lambda_{1},\lambda_{2})$ into (21) to compute $\hat{p}_{1}$. Then substitute $(\lambda_{1},\lambda_{2},\theta,\hat{p}_{1})$ into (17) to get the approximate maximum of Holevo quantity over $p_{1}$: $\chi(\lambda_{1},\lambda_{2},\theta,\hat{p}_{1})$. To see how accurate this approximate maximum is, we need BA algorithm to provide an accurate maximum. The termination criterion for the iteration process of BA algorithm is, stopping when $[\max_{x}\{D(\rho_{x}||\rho^{t})\}-f(p^{t},p^{t})]$ is less than $10^{-6}$, then the BA algorithm outputs a value of Holevo quantity $\chi_{BA}(\lambda_{1},\lambda_{2},\theta)$. And we can compute the error of $\chi(\lambda_{1},\lambda_{2},\theta,\hat{p}_{1})$ then take the maximum over $\theta\in[0,\pi]$

[TABLE]

FIG. 1 is the numerical result, which is a plot of $(\lambda_{1},\lambda_{2},{Error}(\lambda_{1},\lambda_{2}))$.

From Figure 1 we can see that if $\lambda_{1},\lambda_{2}$ are not “too big”, the error can be upper bounded by $10^{-3}$. To see this more directly, we take the maximum of ${Error}(\lambda_{1},\lambda_{2})$ for different ranges of $\lambda_{1},\lambda_{2}$:

[TABLE]

Figure. 2 is a plot of $(R,\max_{\lambda_{1},\lambda_{2}\in[0.5,R]}{Error}(\lambda_{1},\lambda_{2}))$.

From Figure. 2 we can see that if $\lambda_{1},\lambda_{2}<0.95$, the error of approximate maximum of Holevo quantity can be upper bounded by $3\times 10^{-4}$. So we can conclude that when the bigger eigenvalues of $\rho_{1},\rho_{2}$ are not too big (no bigger than $0.95$), (21) can make the error of the maximum of Holevo quantity smaller than $3\times 10^{-4}$.

The approximate solution is an interesting phenomenon. The reason why the angle $\theta$ has such little influence on the maximum of Holevo quantity is unclear.

5 Conclusion

In this paper, we used numerical experiments to show that the BA algorithm works much faster than the theoretical analysis. We also presented an explicit approximate solution for the binary inputs, two dimensional outputs channel, which has an error within $3\times 10^{-4}$ when the eigenvalues of the two channel matrices are not bigger than $0.95$. However, the theoretical reason behind this approximation is not clear.

$\bf{References}$