# Quantum State Fidelity

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%202%20-%20Quantum%20State%20Fidelity.ipynb)


### Introduction
Fidelity between two density operators $\rho$ and $\sigma$,
$$ F(\rho,\sigma) = || \sqrt{\rho} \sqrt{\sigma}||_1,$$
is the optimum value of the semi-definite program (SDP)

\begin{align}
    \begin{aligned}
        \text{maximize} \; & \frac{1}{2} \; \rm Tr(\Lambda + \Lambda^{\dagger}), \\
        \text{subject to} \; & 
        \begin{pmatrix}
          \rho & \Lambda \\
          \Lambda^{\dagger} & \sigma
        \end{pmatrix}
        \succeq 0, & 
    \end{aligned}
\end{align}

where $\Lambda$ is a linear operator. The SDP above has a dual formulation,

\begin{align}
    \begin{aligned}
        \text{minimize} \; & \frac{1}{2} \; \big(  \rm Tr(\rho Y) +  \rm Tr( \sigma Z) \big), \\
        \text{subject to} \; & 
        \begin{pmatrix}
        Y & -I\\
        -I & Z
        \end{pmatrix}
        \succeq 0, & 
    \end{aligned}
\end{align}

where $I$ is the Identity matrix. 


```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
Fidelity $F(\rho_2,\sigma_2)$ between $d$-dimensional random pure state density operators $\rho_2 = | \psi \rangle \langle \psi |$ and $\sigma_2 = | \phi \rangle \langle \phi |$ can be computed using both the Primal and Dual SDP formulations stated at the for of this notebook.



```python
#Example 1
d= 4

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

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

```


```python
#Primal SDP

#Constants
#----------
rhoPic1 = pic.Constant("rho", rho1)
sigmaPic1 = pic.Constant("sigma", sig1)

#Variables
#----------
shp1 = np.shape(rho1)
lmPic1 = pic.ComplexVariable("Lm", shp1)

prob1 = pic.Problem()
    
#Constraint
#----------
prob1.add_constraint(((rhoPic1 & lmPic1) // (lmPic1.H & sigmaPic1)) >> 0)

#Objective
#----------
obj1 = pic.trace(lmPic1 + lmPic1.H)*0.5

prob1.set_objective('max',obj1)

```


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

    -----------------------------
    Complex Semidefinite Program
      maximize tr(Lm + Lmᴴ)·0.5
      over
        4×4 complex variable Lm
      subject to
        [rho, Lm; Lmᴴ, sigma] ≽ 0
    -----------------------------



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




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




```python
#Solver claims to have found optimal solution
fid1 =  prob1.value
```


```python
#Dual SDP

#Constants
#----------
rhoD1 = pic.Constant("rho", rho1)
sigmaD1 = pic.Constant("sigma", sig1)

shp1 = np.shape(rho1)
iMat = pic.Constant('I', np.eye(shp1[0]))

#Variables
#----------
AD = pic.HermitianVariable("A", shp1)
BD = pic.HermitianVariable("B", shp1)

prob1D = pic.Problem()
    
#Constraint
#----------
prob1D.add_constraint(((AD & -iMat) // (-iMat & BD)) >> 0)

#Objective
#----------
obj1D = .5*(rhoD1 | AD).real + .5*(sigmaD1 | BD).real

prob1D.set_objective('min',obj1D)
```


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

    ------------------------------------------------
    Complex Semidefinite Program
      minimize 0.5·Re(⟨rho, A⟩) + 0.5·Re(⟨sigma, B⟩)
      over
        4×4 hermitian variable A, B
      subject to
        [A, -I; -I, B] ≽ 0
    ------------------------------------------------



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




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




```python
#Solver claims to have found optimal solution
fid1D =  prob1D.value

```

$F(\rho_2,\sigma_2)$ between pure state density operators $\rho_2 = | \psi \rangle \langle \psi |$ and $\sigma_2 = | \phi \rangle \langle \phi |$ is simply given by

\begin{equation}
F(\rho_2,\sigma_2) = \sqrt{ \rm Tr (\rho_2 \sigma_2) } = |\langle \psi | \phi \rangle|.
\end{equation}


```python
#Compute fidelity using the formula above
fid1Alg = np.sqrt(np.trace(np.dot(rho1,sig1)).real)
```


```python
print('Fidelity between two random pure states')
print('Using Primal SDP = ', fid1)
print('Using DualSDP = ', fid1D)
print('Using Numpy', fid1Alg)
diffMx1 = max(abs(fid1 - fid1D),abs(fid1 - fid1Alg),abs(fid1D - fid1Alg))
print('Maximum difference between any pair of values above =', diffMx1)    
```

    Fidelity between two random pure states
    Using Primal SDP =  0.18015900161366682
    Using DualSDP =  0.18015900229617843
    Using Numpy 0.1801590015472587
    Maximum difference between any pair of values above = 7.489197317855911e-10


### Example 2
Fidelity $F(\rho_2,\sigma_2)$ between a $d$-dimensional random pure state density operators $\rho_2 = | \psi \rangle \langle \psi |$ and a random mixed state $\sigma_2$. 



```python
#Example 2

d = 3

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

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


```python
#Primal SDP

#Constants
#----------
rho2Pic = pic.Constant("rho2", rho2)
sigma2Pic = pic.Constant("sigma2", sig2)

#Variables
#----------
shp2 = np.shape(rho2)
lm2Pic = pic.ComplexVariable("Lm2", shp2)

prob2 = pic.Problem()
    
#Constraint
#----------
prob2.add_constraint(((rho2Pic & lm2Pic) // (lm2Pic.H & sigma2Pic)) >> 0)

#Objective
#----------
obj2 = pic.trace(lm2Pic + lm2Pic.H)*0.5

prob2.set_objective('max',obj2)

```


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

    ---------------------------------
    Complex Semidefinite Program
      maximize tr(Lm2 + Lm2ᴴ)·0.5
      over
        3×3 complex variable Lm2
      subject to
        [rho2, Lm2; Lm2ᴴ, sigma2] ≽ 0
    ---------------------------------



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




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




```python
#Solver claims to have found optimal solution
fid2 =  prob2.value
```


```python
#Constants
#----------
rhoD2 = pic.Constant("rho2", rho2)
sigmaD2 = pic.Constant("sigma2", sig2)

shp2 = np.shape(rho2)
iMat = pic.Constant('I', np.eye(shp2[0]))

#Variables
#----------
A2D = pic.HermitianVariable("A2", shp2)
B2D = pic.HermitianVariable("B2", shp2)

prob2D = pic.Problem()
    
#Constraint
#----------
prob2D.add_constraint(((A2D & -iMat) // (-iMat & B2D)) >> 0)

#Objective
#----------
obj2D = .5*(rhoD2 | A2D).real + .5*(sigmaD2 | B2D).real

prob2D.set_objective('min',obj2D)
```


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

    ----------------------------------------------------
    Complex Semidefinite Program
      minimize 0.5·Re(⟨rho2, A2⟩) + 0.5·Re(⟨sigma2, B2⟩)
      over
        3×3 hermitian variable A2, B2
      subject to
        [A2, -I; -I, B2] ≽ 0
    ----------------------------------------------------



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




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




```python
#Solver claims to have found optimal solution
fid2D =  prob2D.value
```

$F(\rho_3,\sigma_3)$ between the pure state density operator $\rho_3 = | \psi \rangle \langle \psi |$ and mixed state $\sigma_3$ is simply given by

\begin{equation}
F(\rho_3,\sigma_3) = \sqrt{ \langle \psi|\sigma_3|\psi \rangle} = \sqrt{ \rm Tr (\rho_3 \sigma_3) }
\end{equation}


```python
#Compute fidelity using the formula above
fid2Alg = np.sqrt(np.trace(np.dot(rho2,sig2)).real)
```


```python
print('Fidelity between two random states, one pure and another mixed')
print('Using Primal SDP = ', fid2)
print('Using DualSDP = ', fid2D)
print('Using Numpy', fid2Alg)
diffMx1 = max(abs(fid2 - fid2D),abs(fid2 - fid2Alg),abs(fid2D - fid2Alg))
print('Maximum difference between any pair of values above =', diffMx1)
```

    Fidelity between two random states, one pure and another mixed
    Using Primal SDP =  0.40275233796048393
    Using DualSDP =  0.4028013483991454
    Using Numpy 0.40272188430512085
    Maximum difference between any pair of values above = 7.946409402453947e-05


### Example 3
Fidelity $F(\rho_3,\sigma_3)$ between two $d$-dimensional random density operators $\rho_3$ and $\sigma_3$ is calculated using both primal and dual SDP formulations stated at the top of this notebook.


```python
#Example 3 using the primal formulation
d = 10

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

mt2E3 = np.random.rand(d,d) + 1j*np.random.randn(d,d)
mt2E3 = np.dot(mt2E3,mt2E3.conj().T)
sig3 = mt2E3/np.trace(mt2E3)
```


```python
#Constants
#----------
rho3Pic = pic.Constant("rho3", rho3)
sigma3Pic = pic.Constant("sigma3", sig3)

#Variables
#----------
shp3 = np.shape(rho3)
lm3Pic = pic.ComplexVariable("Lm3", shp3)

prob3 = pic.Problem()
    
#Constraint
#----------
prob3.add_constraint(((rho3Pic & lm3Pic) // (lm3Pic.H & sigma3Pic)) >> 0)

#Objective
#----------
obj3 = pic.trace(lm3Pic + lm3Pic.H)*0.5

prob3.set_objective('max',obj3)
```


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

    ---------------------------------
    Complex Semidefinite Program
      maximize tr(Lm3 + Lm3ᴴ)·0.5
      over
        10×10 complex variable Lm3
      subject to
        [rho3, Lm3; Lm3ᴴ, sigma3] ≽ 0
    ---------------------------------



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




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




```python
#Solver claims to have found optimal solution
fid3 =  prob3.value
```


```python
#Solve the same problem using the Dual Formulation
```


```python
#Constants
#----------
iMat = pic.Constant('I', np.eye(shp3[0]))

#Variables
#----------
A3D = pic.HermitianVariable("A3", shp3)
B3D = pic.HermitianVariable("B3", shp3)

prob3D = pic.Problem()
    
#Constraint
#----------
prob3D.add_constraint(((A3D & -iMat) // (-iMat & B3D)) >> 0)

#Objective
#----------
obj3D = .5*(rho3Pic | A3D).real + .5*(sigma3Pic | B3D).real

prob3D.set_objective('min',obj3D)
```


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

    ----------------------------------------------------
    Complex Semidefinite Program
      minimize 0.5·Re(⟨rho3, A3⟩) + 0.5·Re(⟨sigma3, B3⟩)
      over
        10×10 hermitian variable A3, B3
      subject to
        [A3, -I; -I, B3] ≽ 0
    ----------------------------------------------------



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




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




```python
#Solver claims to have found optimal solution
fid3D =  prob3D.value
```

The Fidelity between two mixed states $\rho_3$ and $\sigma_3$ is given by
$$ F(\rho_3,\sigma_4) = || \sqrt{\rho_3} \sqrt{\sigma_3}||_1,$$



```python
#Calculate the Fidelity using the usual route via NumPy Libraries
diag1,U = np.linalg.eigh(rho3)
sqRho3 = np.dot(np.dot(U, np.diag(np.sqrt(diag1))), U.conj().T)  # Square root of rho3.
diag2,V = np.linalg.eigh(sig3)
sqSig3 = np.dot(np.dot(V, np.diag(np.sqrt(diag2))), V.conj().T)  # Square root of sig3.
fid3Alg = sum(np.linalg.svd( np.dot(sqRho3, sqSig3) )[1] )  # Trace-norm of sqrt(P)·sqrt(Q).

```


```python
print('Fidelity between two random mixed states')
print('Using Primal SDP = ', fid3)
print('Using DualSDP = ', fid3D)
print('Using Numpy', fid3Alg)
diffMx1 = max(abs(fid3 - fid3D),abs(fid3 - fid3Alg),abs(fid3D - fid3Alg))
print('Maximum difference between any pair of values above =', diffMx1)    
```

    Fidelity between two random mixed states
    Using Primal SDP =  0.7599657404167818
    Using DualSDP =  0.7599659619497257
    Using Numpy 0.7599656910322277
    Maximum difference between any pair of values above = 2.709174979909079e-07



```python

```

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