# Quantum Channel Capacity

Date : December 15, 2021

This notebook contains material supporting a paper, currently titled *Five Starter Pieces: Quantum Information Science via Semi-definite Programs*, by Vikesh Siddhu (vsiddhu@protonmail.com) and Sridhar Tayur (stayur@cmu.edu). The paper is available on this **[arXiv](http://arxiv.org/abs/2112.08276)** link. The arXiv paper is released there is under the **[arXiv.org perpetual, non-exclusive license](https://arxiv.org/licenses/nonexclusive-distrib/1.0/license.html)**, and this code is released under the **[MIT license](https://opensource.org/licenses/MIT)**.

This notebook depends upon various packages including [numpy](https://numpy.org/) >= 1.19.5, [picos](https://picos-api.gitlab.io/picos/index.html) >= 2.2.55, and [cvxopt](http://cvxopt.org/) >= 1.2.5.
    
[](https://colab.research.google.com/github/vsiddhu/SDP-Quantum-OR/blob/master/Notebook%205%20-%20Quantum%20Channel%20Capacity.ipynb)

### Introduction
Any linear map $\mathcal{M}:\hat{ \mathcal{H}_a} \mapsto \hat{ \mathcal{H}_b}$ has a Choi-Jamiolkowski representation
$$
\mathcal{J}(\mathcal{M})_{ba} = \mathcal{M} \otimes \mathcal{I}_a (\gamma)
$$
where $\gamma = \sum_{ij} | ii \rangle \langle  jj|$.

A quantum channel $\mathcal{B}:\hat{ \mathcal{H}_a} \mapsto \hat{ \mathcal{H}_b}$ with complementary channel
$\mathcal{B}^c:\hat{ \mathcal{H}_a} \mapsto \hat{ \mathcal{H}_c}$ is $\epsilon$ degradable when there is another quantum channel $\mathcal{C}:\hat{ \mathcal{H}_b} \mapsto \hat{ \mathcal{H}_c}$ such that,

$$
|| \mathcal{C} \circ \mathcal{B} - \mathcal{B}^c ||_{\diamond} \leq \epsilon, 
$$

where $|| \mathcal{M} ||_{\diamond}$ is the diamond norm of a map. Then the smallest $\epsilon$ for which $\mathcal{B}$ is $\epsilon$ degradable is given by twice the optimum value of the semi-definite program

 \begin{align}
 \begin{aligned}
   \text{mininimize} \; & \mu &\\
   \text{subject to} \; & \rm Tr_{c}(Z_{ca}) \preceq \mu I_a, &\\
  & \rm Tr_{c}(J_{cb}(\mathcal{C})) = I_b, &\\
  & Z_{ca} \succeq \mathcal{J}_{ca}(\mathcal{B}^c) - \mathcal{J}_{ca}(\mathcal{C} \circ \mathcal{B}), &\\
  & Z_{ca} \succeq 0, & \\
  \text{and} \; & \mathcal{J}_{cb}(\mathcal{C}) \succeq 0 &,
  \end{aligned}
  \end{align}
  
  where $I_a$ and $I_b$ are the identity on $\mathcal{H}_a$ and $\mathcal{H}_b$ respectively.
In what follows, consider various examples of this SDP. 

If $\mathcal{B}$ is $\epsilon$-degradable and has channel coherent information $\mathcal{Q}^{(1)}(\mathcal{B})$,
then 

\begin{equation}
|\mathcal{Q}(\mathcal{B}) - \mathcal{Q}^{(1)}(\mathcal{B})| 
    \leq \epsilon \log(d_c - 1)/2 + \epsilon \log(d_c) + h(\epsilon/2) + (1 + \epsilon/2)h(2/(2+ \epsilon))
\end{equation}
where $h(x) = -\big( x \log x +(1-x) \log (1+x) \big)$ is the binary entropy function.


```python
%pylab inline
```

    Populating the interactive namespace from numpy and matplotlib



```python
import HelperFunctionsforNotebook5 as crp
import picos as pic
import cvxopt as cvx
```


```python
print('Solvers supported on this installation of picos:', pic.solvers.all_solvers().keys())

```

    Solvers supported on this installation of picos: dict_keys(['cplex', 'cvxopt', 'ecos', 'glpk', 'gurobi', 'mosek', 'mskfsn', 'scip', 'smcp'])



```python
print('Solvers available to picos on this machine :', pic.solvers.available_solvers())
```

    Solvers available to picos on this machine : ['cvxopt', 'mosek', 'mskfsn']



```python
# For Google Colab use, commands installing packages
try:
  import google.colab
  IN_COLAB = True
except:
  IN_COLAB = False

# Let's start with Pyomo
if IN_COLAB:
    !pip install -q picos
    !pip install -q cvxopt
```

### Example 1
Compute smallest $\epsilon$ for which the qubit dephasing channel 
\begin{equation}
    \mathcal{F}_{q}(\rho) = (1-q) \rho + q Z \rho Z
\end{equation}
$0 \leq q \leq 1$ is $\epsilon$ degradable. 


```python
q = np.random.rand()

#Kraus Operators for qubit dephasing channel
k1 = np.sqrt(1-q)*np.array([[1.,0.],[0.,1.]])
k2 = np.sqrt(q)*np.array([[1.,0.],[0.,-1.]])
krsLstF = [k1,k2]

(dc,db,da) = np.shape(krsLstF)

#Choi-Jamiolkowski representations of channel and its complement 
(jbaF, jcaF) = crp.krausToChoiJ(krsLstF)

```


```python
#Constants
#----------
JbaFPic = pic.Constant("J(B)_ba", jbaF)
JcaFPic = pic.Constant("J(B)_ca", jcaF)

iMatA = pic.Constant('Ia', np.eye(da))
iMatB = pic.Constant('Ib', np.eye(db))
shpCA = (dc*da,dc*da)
shpCB = (dc*db,dc*db)


#Variables
#----------
ZPic = pic.HermitianVariable("Zca", shpCA)
JPic = pic.HermitianVariable("Jcb", shpCB)
mu = pic.RealVariable("mu")

```


```python
prob1P = pic.Problem()

#Constraint
#----------
prob1P.add_constraint(ZPic >> 0)
prob1P.add_constraint(JPic >> 0)

ZaPic = pic.partial_trace(ZPic,subsystems=(0), dimensions=(dc,da))
JbPic = pic.partial_trace(JPic,subsystems=(0), dimensions=(dc,db))

prob1P.add_constraint(ZaPic << mu*iMatA)
prob1P.add_constraint(JbPic == iMatB)

JcaPic = crp.choiJOfChanInSeriesPic(JbaFPic,JPic,da,db,dc)
prob1P.add_constraint(ZPic >> JcaFPic-JcaPic)
```




    <4×4 Complex LMI Constraint: Zca ≽ J(B)_ca - shuffled(shuffled(Jcb,ikjl,(2,2,2,2),C)·shuffled(J(B)_ba,ikjl,(2,2,2,2),C),ikjl,(2,2,2,2),C)>




```python
#Objective
#----------
prob1P.set_objective('min',mu)

#User readable view of the problem being composed in PICOS'
print(prob1P)
```

    --------------------------------------------------------------------------------------------------
    Complex Semidefinite Program
      minimize mu
      over
        4×4 hermitian variable Jcb, Zca
        1×1 real variable mu
      subject to
        Zca ≽ 0
        Jcb ≽ 0
        Zca.{tr([2×2])⊗[2×2]} ≼ mu·Ia
        Jcb.{tr([2×2])⊗[2×2]} = Ib
        Zca ≽ J(B)_ca -
          shuffled(shuffled(Jcb,ikjl,(2,2,2,2),C)·shuffled(J(B)_ba,ikjl,(2,2,2,2),C),ikjl,(2,2,2,2),C)
    --------------------------------------------------------------------------------------------------



```python
#Solve the problem using cvxopt as a solver
prob1P.solve(verbosity=False,solver='cvxopt')
```




    <primal feasible solution pair (claimed optimal) from cvxopt>




```python
#Solver claims to have found optimal saolution
mu1P =  prob1P.value
eps1 = 2*mu1P
print(eps1)
```

    6.20984358392016e-10



```python
print('Least epsilon for which the dephasing channel is epsilon degradable')
print('Using SDP = ', eps1)
print('Difference between the SDP value and the algberaic values', abs(eps1 - 0))
```

    Least epsilon for which the dephasing channel is epsilon degradable
    Using SDP =  6.20984358392016e-10
    Difference between the SDP value and the algberaic values 6.20984358392016e-10


The qubit dephasing channel, 

$$\mathcal{F}_{q}(\rho) = (1-q) \rho + q Z \rho Z, $$

is degradable for any $0 \leq q \leq 1$. This degradability implies that 
$\mathcal{F}_{q}$ is $\epsilon$ degradable for $\epsilon = 0$, as confirmed by
the numerics above. Given $\epsilon = 0$, the bound

\begin{equation}
|\mathcal{Q}(\mathcal{F}_q) - \mathcal{Q}^{(1)}(\mathcal{F}_q)| 
    \leq \epsilon \log(d_c - 1)/2 + \epsilon \log(d_c) + h(\epsilon/2) + (1 + \epsilon/2)h(2/(2+ \epsilon))
\end{equation}

trivially evaluates to 
\begin{equation}
   |\mathcal{Q}(\mathcal{F}_q) - \mathcal{Q}^{(1)}(\mathcal{F}_q)| =0
\end{equation}

Using the above equation, along with the expression for the channel coherent information, we get 
\begin{equation}
    \mathcal{Q}(\mathcal{F}_q) = \mathcal{Q}^{(1)}(\mathcal{F}_{q}) = 1 - h(q).
\end{equation}

### Example 2
Compute smallest $\epsilon$ for which the qubit erasure channel 
\begin{equation}
    \mathcal{E}_{p}(\rho) = (1-p) \rho + p \rm Tr(\rho) |e\rangle\langle e|
\end{equation}
with erasure probability $0 \leq p \leq 1$ is $\epsilon$ degradable. 


```python
#Consider a random erasure probability less than half
p = np.random.rand()/2

#Kraus Operators for Erasure
k1 = np.sqrt(p)*np.array([[0.,0.],[0.,0.],[1.,0.]])
k2 = np.sqrt(p)*np.array([[0.,0.],[0.,0.],[0.,1.]])
k3 = np.sqrt(1-p)*np.array([[1.,0.],[0.,1.],[0.,0.]])
krsLstE = [k1,k2,k3]

(dc,db,da) = np.shape(krsLstE)

#Choi-Jamiolkowski representations of channel and its complement 
(jbaE, jcaE) = crp.krausToChoiJ(krsLstE)
```


```python
#Constants
#----------
JbaEPic = pic.Constant("J(B)_ba", jbaE)
JcaEPic = pic.Constant("J(B)_ca", jcaE)
iMatA = pic.Constant('Ia', np.eye(da))
iMatB = pic.Constant('Ib', np.eye(db))
shpCA = (dc*da,dc*da)
shpCB = (dc*db,dc*db)

#Variables
#----------
ZPic = pic.HermitianVariable("Zca", shpCA)
JPic = pic.HermitianVariable("Jcb", shpCB)
mu = pic.RealVariable("mu")

```


```python
prob2P = pic.Problem()
```


```python
#Constraint
#----------
prob2P.add_constraint(ZPic >> 0)
prob2P.add_constraint(JPic >> 0)

ZaPic = pic.partial_trace(ZPic,subsystems=(0), dimensions=(dc,da))
JbPic = pic.partial_trace(JPic,subsystems=(0), dimensions=(dc,db))

prob2P.add_constraint(ZaPic << mu*iMatA)
prob2P.add_constraint(JbPic == iMatB)

JcaPic = crp.choiJOfChanInSeriesPic(JbaEPic,JPic,da,db,dc)
prob2P.add_constraint(ZPic >> JcaEPic-JcaPic)

```




    <6×6 Complex LMI Constraint: Zca ≽ J(B)_ca - reshaped(shuffled(shuffled(Jcb,ikjl,(3,3,3,3),C)·reshaped(shuffled(J(B)_ba,ikjl,(3,2,3,2),C)ᵀ, 4×9)ᵀ,ikjl,(3,3,2,2),C)ᵀ, 6×6)ᵀ>




```python
#Objective
#----------
prob2P.set_objective('min',mu)

#User readable view of the problem being composed in PICOS'
print(prob2P)
```

    ----------------------------------------------------------------------------------------------------
    Complex Semidefinite Program
      minimize mu
      over
        6×6 hermitian variable Zca
        9×9 hermitian variable Jcb
        1×1 real variable mu
      subject to
        Zca ≽ 0
        Jcb ≽ 0
        Zca.{tr([3×3])⊗[2×2]} ≼ mu·Ia
        Jcb.{tr([3×3])⊗[3×3]} = Ib
        Zca ≽ J(B)_ca -
          reshaped(shuffled(shuffled(Jcb,ikjl,(3,3,3,3),C)·reshaped(shuffled(J(B)_ba,ikjl,(3,2,3,2),C)ᵀ,
          4×9)ᵀ,ikjl,(3,3,2,2),C)ᵀ, 6×6)ᵀ
    ----------------------------------------------------------------------------------------------------



```python
#Solve the problem using cvxopt as a solver
prob2P.solve(verbosity=False,solver='cvxopt')
```




    <primal feasible solution pair (claimed optimal) from cvxopt>




```python
#Solver claims to have found optimal saolution
mu2P =  prob2P.value
eps2 = 2*mu2P
```


```python
print('Using SDP = ', eps2)
print('Bound on difference between capacity and coherent information =', crp.boundFunc(eps2,dc))
```

    Using SDP =  1.2914823168265619e-09
    Bound on difference between capacity and coherent information = 4.3982699098789774e-08


The qubit erasure channel $\mathcal{E}_{p}$ has coherent information 

\begin{equation}
    \mathcal{Q}^{(1)}(\mathcal{E}_{p}) = \max(1 - 2p,0).
\end{equation}

The channel is degradable for erasure probability $0 \leq p \leq 1/2$. This degradability implies that $\mathcal{E}_{p}$ is $\epsilon$ degradable for $\epsilon = 0$, as confirmed by the numerics above when $0 \leq p \leq 1/2$. For values of $0 \leq p \leq 1/2$, the bound

\begin{equation}
|\mathcal{Q}(\mathcal{E}_p) - \mathcal{Q}^{(1)}(\mathcal{E}_{p})| 
    \leq \epsilon \log(d_c - 1)/2 + \epsilon \log(d_c) + h(\epsilon/2) + (1 + \epsilon/2)h(2/(2+ \epsilon))
\end{equation}

trivially evaluates to 
\begin{equation}
   |\mathcal{Q}(\mathcal{E}_p) - \mathcal{Q}^{(1)}(\mathcal{E}_p)| =0
\end{equation}
For $1/2 \leq p \leq 1$, $\epsilon$ is non-zero,the bound is also non-zero. It turns out that even for $1/2 \leq p \leq 1$, $\mathcal{Q}(\mathcal{E}_p) = \mathcal{Q}^{(1)}(\mathcal{E}_p)$, due to anti-degradability of the
erasure channel for erasure probability $1/2 \leq p \leq 1$.

### Example 3
Compute smallest $\epsilon$ for which the qubit depolarizing channel 
\begin{equation}
    \Delta(\rho) = (1-p) \rho + \frac{p}{3} (X \rho X + Y \rho Y  + Z \rho Z),
\end{equation}
$0 \leq p \leq 1$ is $\epsilon$ degradable. The depolarizing channel above can also be written using the parameter $\lambda = 1 - 4p/3$ as
\begin{equation}
    \Delta(\rho) = \lambda \rho + \frac{(1-\lambda)}{2} \rm Tr(\rho) I.
\end{equation}



```python
#Example 3
p = np.random.rand()/10

#Kraus Operators for Depolarizing Channel
k1 = np.sqrt(1-p)*np.array([[1.,0.],[0.,1.]])
k2 = np.sqrt(p/3)*np.array([[0.,1.],[1.,0.]])
k3 = np.sqrt(p/3)*np.array([[0.,-1j],[1j,0.]])
k4 = np.sqrt(p/3)*np.array([[1.,0.],[0.,-1.]])

krsLstD = [k1,k2,k3,k4]
(dc,db,da) = np.shape(krsLstD)

(jbaD, jcaD) = crp.krausToChoiJ(krsLstD)
```


```python
#Constants
#----------
JbaDPic = pic.Constant("J(B)_ba", jbaD)
JcaDPic = pic.Constant("J(B)_ca", jcaD)
iMatA = pic.Constant('Ia', np.eye(da))
iMatB = pic.Constant('Ib', np.eye(db))

#Variables
#----------
shpCA = (dc*da,dc*da)
shpCB = (dc*db,dc*db)

ZPic = pic.HermitianVariable("Zca", shpCA)
JPic = pic.HermitianVariable("Jcb", shpCB)
mu = pic.RealVariable("mu")
```


```python
prob3P = pic.Problem()
```


```python
#Constraint
#----------
prob3P.add_constraint(ZPic >> 0)
prob3P.add_constraint(JPic >> 0)

ZaPic = pic.partial_trace(ZPic,subsystems=(0), dimensions=(dc,da))
JbPic = pic.partial_trace(JPic,subsystems=(0), dimensions=(dc,db))

prob3P.add_constraint(ZaPic << mu*iMatA)
prob3P.add_constraint(JbPic == iMatB)

JcaPic = crp.choiJOfChanInSeriesPic(JbaDPic,JPic,da,db,dc)
prob3P.add_constraint(ZPic >> JcaDPic-JcaPic)

```




    <8×8 Complex LMI Constraint: Zca ≽ J(B)_ca - reshaped(shuffled(reshaped(shuffled(Jcb,ikjl,(4,2,4,2),C)ᵀ, 4×16)ᵀ·shuffled(J(B)_ba,ikjl,(2,2,2,2),C),ikjl,(4,4,2,2),C)ᵀ, 8×8)ᵀ>




```python
#Objective
#----------
prob3P.set_objective('min',mu)

#User readable view of the problem being composed in PICOS'
print(prob3P)
```

    -------------------------------------------------------------------
    Complex Semidefinite Program
      minimize mu
      over
        8×8 hermitian variable Jcb, Zca
        1×1 real variable mu
      subject to
        Zca ≽ 0
        Jcb ≽ 0
        Zca.{tr([4×4])⊗[2×2]} ≼ mu·Ia
        Jcb.{tr([4×4])⊗[2×2]} = Ib
        Zca ≽ J(B)_ca -
          reshaped(shuffled(reshaped(shuffled(Jcb,ikjl,(4,2,4,2),C)ᵀ,
          4×16)ᵀ·shuffled(J(B)_ba,ikjl,(2,2,2,2),C),ikjl,(4,4,2,2),C)ᵀ,
          8×8)ᵀ
    -------------------------------------------------------------------



```python
#Solve the problem using cvxopt as a solver
prob3P.solve(verbosity=False,solver='cvxopt')
```




    <primal feasible solution pair (claimed optimal) from cvxopt>




```python
#Solver claims to have found optimal solution
mu3P =  prob3P.value
eps3 = 2*mu3P
```


```python
print('Least epsilon for which the depolarizing channel is epsilon degradable')
print('Using the SDP = ', eps3)
print('Bound on difference between capacity and coherent information =', crp.boundFunc(eps3,dc))
```

    Least epsilon for which the depolarizing channel is epsilon degradable
    Using the SDP =  0.0030743301525582175
    Bound on difference between capacity and coherent information = 0.041751518313526206


The qubit depolarizing channel is not degradable and the least value of $\epsilon$ for which the channel is $\epsilon$ degradable is found numerically, as demonstrated above.

The channel coherent information 
\begin{equation}
    \mathcal{Q}^{(1)}(\Delta) = \max(1 - h(p) - p \log 3, 0),
\end{equation}
along with $\epsilon$ computed using the SDP can be used to obtain a bound,

\begin{equation}
|\mathcal{Q}(\Delta) - \mathcal{Q}^{(1)}(\Delta)| 
    \leq \epsilon \log(d_c - 1)/2 + \epsilon \log(d_c) + h(\epsilon/2) + (1 + \epsilon/2)h(2/(2+ \epsilon)),
\end{equation}
which is usually not tight, except when $p = 0$ or $p = 1$.

566d181c3b596beb0c29498e3542882fd6ee705a
ipynb
Jupyter Notebook
Notebook 5 - Quantum Channel Capacity.ipynb
vsiddhu/SDP-Quantum-OR
32864fca86ad33e40258ccbde99bd315f7d5ff67
MIT
2021-12-16T08:21:09.000Z
2022-03-29T18:07:37.000Z
Notebook 5 - Quantum Channel Capacity.ipynb
vsiddhu/SDP-Quantum-OR
32864fca86ad33e40258ccbde99bd315f7d5ff67
MIT
Notebook 5 - Quantum Channel Capacity.ipynb
vsiddhu/SDP-Quantum-OR
32864fca86ad33e40258ccbde99bd315f7d5ff67
MIT
2021-12-19T06:42:39.000Z
2022-01-22T08:54:14.000Z
Qwen/Qwen-72B
1. YES
2. YES