# Quantum State Discrimination

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%201%20-%20Quantum%20State%20Discrimination.ipynb)

## Introduction

The maximum probability of successfully distinguishing states $\sigma_i$, created with probability 
$p_i$, is given by

\begin{align}
    \begin{aligned}
        \text{maximize} \; &  \sum_{i = 1}^{n-1} p_i \rm Tr(E_i \rho_i) + \rm p_n Tr(E_n \rho_n) & \\ 
        \text{subject to} \; & E_j \succeq 0,& \forall 1 \leq j \leq n-1, \\ 
        \text{and} \; & I - \sum_{i = 1}^{n-1} E_i \succeq 0. &
    \end{aligned}
\end{align}

where $p_n := 1 - \sum_{i = 1}^{n-1} p_i$ and $E_n := I - \sum_{i = 1}^{n-1} E_i$.


```python
%pylab inline
```

    Populating the interactive namespace from numpy and matplotlib



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

# Install PICOS and CVXOPT in Google Colab
if IN_COLAB:
    !pip install -q picos
    !pip install -q cvxopt
```


```python
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']


### Example 1
Consider a case where $n=2$, density operators $\sigma_i$ have the form
    \begin{equation}
    \sigma_1 = 
    \begin{pmatrix}
    1 & 0 \\
    0 & 0
    \end{pmatrix}, 
    \quad
    \sigma_2 = 
    \begin{pmatrix}
    0 & 0 \\
    0 & 1
    \end{pmatrix},
    \end{equation}
  and $p_1 = q, p_2 = 1-q$ and $0 \leq q \leq 1$. In this case the SDP mentioned at the top of the notebook takes the form
  
\begin{align}
    \begin{aligned}
        \text{maximize} \; &  q \rm Tr(E_1 \sigma_1) + (1-q) \rm Tr\big( (I - E_1) \sigma_2\big) & \\ 
        \text{subject to} \; & E_1 \succeq 0, &, \\ 
        \text{and} \; & I - E_1 \succeq 0. &
    \end{aligned}
\end{align}

We solve this semidefinite program using PICOS.


```python
#Example 1

sig1 = np.array([[1.,0],[0.,0.]])
sig2 = np.array([[0.,0],[0.,1.]])

qVal = np.random.rand()
```


```python

#Constants
#----------
Sgs1 = pic.Constant("sg1", sig1)
Sgs2 = pic.Constant("sg2", sig2)

q = pic.Constant("q", qVal)

#Identity matrix
shp = np.shape(sig1)
iMat = pic.Constant('I', np.eye(shp[0],dtype='complex'))

#Variables
#----------
eVars = pic.HermitianVariable("E1", shp)

prob1 = pic.Problem()
    
#Constraint
#----------
prob1.add_constraint(eVars >> 0)
prob1.add_constraint(iMat - eVars >> 0)

#Objective
#----------
obj = q*(Sgs1 | eVars) + (1-q)*(Sgs2 | iMat - eVars)

prob1.set_objective('max',obj)

```


```python
#User readable view of the problem being composed in PICOS
print(prob1)
```

    ----------------------------------------------
    Complex Semidefinite Program
      maximize q·⟨sg1, E1⟩ + (1 - q)·⟨sg2, I - E1⟩
      over
        2×2 hermitian variable E1
      subject to
        E1 ≽ 0
        I - E1 ≽ 0
    ----------------------------------------------



```python
#Solve the problem using a solver of our choice and display progress
prob1.solve(verbosity=False,solver='cvxopt')
```




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




```python
#The solver claims to have found an optimal solution. As a bonus, the SDP solution also gives the POVM element E1
pvm1 = np.array(eVars.value)
```


```python
prb = prob1.value
print('Probability of discriminating the inputs is = ', prb)
#For an explaination of why success probability is 1 even though the input p is random, see the text below
```

    Probability of discriminating the inputs is =  0.9999999994373823


* The probability of successful discriminating two qauntum states $\sigma_1$ and $\sigma_2$ each occuring with probability $p_1$ and $p_2$ respectively, is given by 
 $$p^* = \frac{1}{2}( 1 + || p_1 \sigma_1 - p_2\sigma_2 ||_1 ).$$
 The example above numerically demonstrates this fact in a simple case where $\sigma_1$ and $\sigma_2$ are $2 \times 2$ matrices. In this example, $p_1 = q$ and $p_2 = 1 -q$, where $q \in [0,1]$ is random. Density operators $\sigma_1$ and $\sigma_2$ are orthogonal. These orthogonal density operators have $l_1$ distance, $|| \sigma_1 - \sigma_2 ||_1 = 2$ and thus $p^* = 1$, independent of $q$. In the following example, we choose $p_1,p_2,\sigma_1$ and $\sigma_2$ randomly and confirm our numerical SDP result with the algebraic result above.

### Example 2
Consider a case where $n=2, d = 6$, density operators $\sigma_1$ and $\sigma_2$ are chosen randomly, $p_1 = q$,
$p_2 = 1-q$ and $q$ chosen randomly between 0 and 1. In this case the SDP mentioned at the top of the notebook takes the form

\begin{align}
    \begin{aligned}
        \text{maximize} \; &  q \rm Tr(E_1 \sigma_1) + (1-q) \rm Tr\big( (I - E_1) \sigma_2\big) & \\ 
        \text{subject to} \; & E_1 \succeq 0, &, \\ 
        \text{and} \; & I - E_1 \succeq 0. &
    \end{aligned}
\end{align}



```python
#Example 2
d = 6

mt1 = np.random.rand(d,d) + 1j*np.random.randn(d,d)
mt1 = np.dot(mt1,mt1.conj().T)
sig1 = mt1/np.trace(mt1)

mt2 = np.random.rand(d,d) + 1j*np.random.randn(d,d)
mt2 = np.dot(mt2,mt2.conj().T)
sig2 = mt2/np.trace(mt2)

qVal = np.random.rand()
```


```python

#Constants
#----------
Sgs1 = pic.Constant("sg1", sig1)
Sgs2 = pic.Constant("sg2", sig2)

q = pic.Constant("q", qVal)

#Identity matrix
shp = np.shape(sig1)
iMat = pic.Constant('I', np.eye(shp[0],dtype='complex'))

#Variables
#----------
eVars = pic.HermitianVariable("E1", shp)

prob2 = pic.Problem()
    
#Constraint
#----------
prob2.add_constraint(eVars >> 0)
prob2.add_constraint(iMat - eVars >> 0)

#Objective
#----------
obj = q*(Sgs1 | eVars) + (1-q)*(Sgs2 | iMat - eVars)

prob2.set_objective('max',obj.real)

```


```python
#User readable view of the problem being composed in PICOS
print(prob2)
```

    --------------------------------------------------
    Complex Semidefinite Program
      maximize Re(q·⟨sg1, E1⟩ + (1 - q)·⟨sg2, I - E1⟩)
      over
        6×6 hermitian variable E1
      subject to
        E1 ≽ 0
        I - E1 ≽ 0
    --------------------------------------------------



```python
#Solve the problem using a solver of our choice and display progress
prob2.solve(verbosity=False,solver='cvxopt')
```




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




```python
#Helstrom result computes using numpy
delOpt = qVal*sig1 - (1-qVal)*sig2
# print(np.linalg.eigh(delOpt))
(u,s,vH) = np.linalg.svd(delOpt, hermitian=False)
pStarAlg = (np.sum(s) + 1)/2
```


```python
prb2 = prob2.value
print('Probability of discriminating the input')
print('Using SDP = ', prb2)
print('Helstrom result', pStarAlg)
print('Difference between these values', abs(prb2 - pStarAlg))
```

    Probability of discriminating the input
    Using SDP =  0.9129937076702335
    Helstrom result 0.9129937078308089
    Difference between these values 1.6057544183212258e-10


### Example 3
Consider a case where $n=3$, the density operators $\sigma_i$ have the form
    \begin{equation}
    \sigma_1 = 
    \begin{pmatrix}
    1 & 0 \\
    0 & 0
    \end{pmatrix}
    \quad
    \sigma_2 = 
    \frac{1}{2}\begin{pmatrix}
    1 & 1 \\
    1 & 1
    \end{pmatrix}
    \quad
    \sigma_3 = 
    \frac{1}{2}\begin{pmatrix}
    1 & 0 \\
    0 & 1
    \end{pmatrix}
    \end{equation}
 $p_1 = \frac{q}{2}, p_2 = \frac{q}{2}, p_3 = 1-q$, and $0 \leq q \leq 1$. In this case the SDP mentioned at the top of the notebook takes the form
\begin{align}
    \begin{aligned}
        \text{maximize} \; &  \frac{q}{2} \rm Tr(E_1 \sigma_1)+ \frac{q}{2} Tr(E_2 \sigma_2) + (1-q) \rm Tr\big( (I -E_1 - E_2) \sigma_3\big) & \\ 
        \text{subject to} \; & E_1 \succeq 0, &, \\ 
         & E_1 \succeq 0. & \\
        \text{and} \; & I - E_1 - E_2 \succeq 0. &
    \end{aligned}
\end{align}




```python
#Example 3

sig1 = np.array([[1.,0.],[0.,0.]])
sig2 = np.array([[0.,0.],[0.,1.]])
sig3 = np.array([[1.,0.],[0.,1.]])/2

qVal = np.random.rand()
```


```python
#Constants
#----------
Sgs1 = pic.Constant("sg1", sig1)
Sgs2 = pic.Constant("sg2", sig2)
Sgs3 = pic.Constant("sg2", sig3)

q = pic.Constant("q", qVal)

#Identity matrix
shp = np.shape(sig1)
iMat = pic.Constant('I', np.eye(shp[0],dtype='complex'))

#Variables
#----------
eVars1 = pic.HermitianVariable("E1", shp)
eVars2 = pic.HermitianVariable("E2", shp)

prob3 = pic.Problem()
    
#Constraint
#----------
prob3.add_constraint(eVars1 >> 0)
prob3.add_constraint(eVars2 >> 0)
prob3.add_constraint(iMat - eVars1 -eVars2 >> 0)

#Objective
#----------
obj3 = (q/2)*(Sgs1 | eVars1) + (q/2)*(Sgs2 | eVars2) + (1-q)*(Sgs3 | iMat - eVars1 - eVars2)

prob3.set_objective('max',obj3)

```


```python
#User readable view of the problem being composed in PICOS
print(prob3)
```

    ---------------------------------------------------------------------
    Complex Semidefinite Program
      maximize q/2·⟨sg1, E1⟩ + q/2·⟨sg2, E2⟩ + (1 - q)·⟨sg2, I - E1 - E2⟩
      over
        2×2 hermitian variable E1, E2
      subject to
        E1 ≽ 0
        E2 ≽ 0
        I - E1 - E2 ≽ 0
    ---------------------------------------------------------------------



```python
#Let us solve the problem using a solver of our choice
prob3.solve(verbosity=False,solver='cvxopt')
```




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




```python
#The solver claims to have found an optimal solution. As a bonus, the SDP solution also gives the POVM element E
pov1 = np.array(eVars1.value)
pov2 = np.array(eVars2.value)
```

Using the structure of the SDP being solved, one can show its maximum value to be $\max(q, 1-q)$. The numerics above match this algebraic solution


```python
prb3 = prob3.value
print('Probability of discriminating the input given by SDP = ', prb3)
print('Algebraic solution = ', max(qVal,1-qVal))

```

    Probability of discriminating the input given by SDP =  0.7133864729887238
    Algebraic solution =  0.7133864734236646


83d755e43be07c11842dae5b5964f95dfcd67fed
ipynb
Jupyter Notebook
Notebook 1 - Quantum State Discrimination.ipynb
vsiddhu/SDP-Quantum-OR
32864fca86ad33e40258ccbde99bd315f7d5ff67
MIT
2021-12-16T08:21:09.000Z
2022-03-29T18:07:37.000Z
Notebook 1 - Quantum State Discrimination.ipynb
vsiddhu/SDP-Quantum-OR
32864fca86ad33e40258ccbde99bd315f7d5ff67
MIT
Notebook 1 - Quantum State Discrimination.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