Fidelity of quantum states #1459

Merged
merged 8 commits into from Aug 11, 2012
@@ -0,0 +1,142 @@
+{
+ "metadata": {
+ "name": "fidelity"
+ },
+ "nbformat": 2,
+ "worksheets": [
+ {
+ "cells": [
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "from sympy import *",
+ "from sympy.physics.quantum import *",
+ "from sympy.physics.quantum.density import *",
+ "from sympy.physics.quantum.spin import (",
+ " Jx, Jy, Jz, Jplus, Jminus, J2,",
+ " JxBra, JyBra, JzBra,",
+ " JxKet, JyKet, JzKet,",
+ ")",
+ "from IPython.core.display import display_pretty",
+ "from sympy.physics.quantum.operator import *",
+ "",
+ "%load_ext sympyprinting"
+ ],
+ "language": "python",
+ "outputs": [],
+ "prompt_number": 2
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "##Fidelity using some kets"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "up = JzKet(S(1)/2,S(1)/2)",
+ "down = JzKet(S(1)/2,-S(1)/2)",
+ "amp = 1/sqrt(2)",
+ "updown = (amp * up ) + (amp * down)",
+ "",
+ "# represent turns Kets into matrices",
+ "up_dm = represent(up * Dagger(up))",
+ "down_dm = represent(down * Dagger(down)) ",
+ "updown_dm = represent(updown * Dagger(updown))",
+ "updown2 = (sqrt(3)/2 )* up + (1/2)*down",
+ "",
+ "",
+ "display_pretty(fidelity(up_dm, up_dm))",
+ "display_pretty(fidelity(up_dm, down_dm)) #orthogonal states",
+ "display_pretty(fidelity(up_dm, updown_dm).evalf())",
+ "",
+ "",
+ "# alternatively, puts Kets into Density object and compute fidelity",
+ "d1 = Density( [updown, 0.25], [updown2, 0.75])",
+ "d2 = Density( [updown, 0.75], [updown2, 0.25])",
+ "display_pretty(fidelity(d1, d2))"
+ ],
+ "language": "python",
+ "outputs": [
+ {
+ "output_type": "display_data",
+ "text": [
+ "1"
+ ]
+ },
+ {
+ "output_type": "display_data",
+ "text": [
+ "0"
+ ]
+ },
+ {
+ "output_type": "display_data",
+ "text": [
+ "0.707106781186548"
+ ]
+ },
+ {
+ "output_type": "display_data",
+ "text": [
+ "0.817293551913876"
+ ]
+ }
+ ],
+ "prompt_number": 7
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "## Fidelity on states as Qubits"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "collapsed": false,
+ "input": [
+ "",
+ "from sympy.physics.quantum.qubit import Qubit",
+ "state1 = Qubit('0')",
+ "state2 = Qubit('1')",
+ "state3 = (1/sqrt(2))*state1 + (1/sqrt(2))*state2",
+ "state4 = (sqrt(S(2)/3))*state1 + (1/sqrt(3))*state2",
+ "",
+ "state1_dm = Density([state1, 1])",
+ "state2_dm = Density([state2, 1])",
+ "state3_dm = Density([state3, 1])",
+ "",
+ "# mixed qubit states in density",
+ "d1 = Density([state3, 0.70], [state4, 0.30])",
+ "d2 = Density([state3, 0.20], [state4, 0.80])",
+ "",
+ "",
+ "display_pretty(fidelity(d1, d2))",
+ "",
+ ""
+ ],
+ "language": "python",
+ "outputs": [
+ {
+ "output_type": "display_data",
+ "text": [
+ "0.996370452558227"
+ ]
+ }
+ ],
+ "prompt_number": 9
+ },
+ {
+ "cell_type": "code",
+ "collapsed": true,
+ "input": [],
+ "language": "python",
+ "outputs": []
+ }
+ ]
+ }
+ ]
+}
@@ -1,10 +1,9 @@
-from sympy import Tuple, Add, Mul, Matrix, log, expand
+from sympy import Tuple, Add, Mul, Matrix, log, expand, sqrt, Rational
from sympy.core.trace import Tr
from sympy.printing.pretty.stringpict import prettyForm
from sympy.physics.quantum.dagger import Dagger
from sympy.physics.quantum.operator import HermitianOperator, OuterProduct, Operator
from sympy.physics.quantum.represent import represent
-from sympy.physics.quantum.state import KetBase
from matrixutils import numpy_ndarray, scipy_sparse_matrix, to_numpy
from sympy.physics.quantum.tensorproduct import TensorProduct, tensor_product_simp
from sympy.core.compatibility import product
@@ -254,3 +253,60 @@ def entropy(density):
return -np.sum(eigvals*np.log(eigvals))
else:
raise ValueError("numpy.ndarray, scipy.sparse or sympy matrix expected")
+
+
+def fidelity(state1, state2):
+ """ Computes the fidelity between two quantum states
+ (http://en.wikipedia.org/wiki/Fidelity_of_quantum_states)
+
+ The arguments provided to this function should be a square matrix or a
+ Density object. If it is a square matrix, it is assumed to be diagonalizable.
+
+ Parameters:
+ ==========
+
+ state1, state2 : a density matrix or Matrix
+
+
+ Examples:
+ =========
+
+ >>> from sympy import S, sqrt
+ >>> from sympy.physics.quantum.dagger import Dagger
+ >>> from sympy.physics.quantum.spin import JzKet
+ >>> from sympy.physics.quantum.density import Density, fidelity
+ >>> from sympy.physics.quantum.represent import represent
+ >>>
+ >>> up = JzKet(S(1)/2,S(1)/2)
+ >>> down = JzKet(S(1)/2,-S(1)/2)
+ >>> amp = 1/sqrt(2)
+ >>> updown = (amp * up) + (amp * down)
+ >>>
+ >>> # represent turns Kets into matrices
+ >>> up_dm = represent(up * Dagger(up))
+ >>> down_dm = represent(down * Dagger(down))
+ >>> updown_dm = represent(updown * Dagger(updown))
+ >>>
+ >>> fidelity(up_dm, up_dm)
+ 1
+ >>> fidelity(up_dm, down_dm) #orthogonal states
+ 0
+ >>> fidelity(up_dm, updown_dm).evalf().round(3)
+ 0.707
+
+ """
+ state1 = represent(state1) if isinstance(state1, Density) else state1
+ state2 = represent(state2) if isinstance(state2, Density) else state2
+
+ if (not isinstance(state1, Matrix) or
+ not isinstance(state2, Matrix)):
+ raise ValueError("state1 and state2 must be of type Density or Matrix "
+ "received type=%s for state1 and type=%s for state2" %
+ (type(state1), type(state2)))
+
+ if ( state1.shape != state2.shape and state1.is_square):
+ raise ValueError("The dimensions of both args should be equal and the"
+ "matrix obtained should be a square matrix")
+
+ sqrt_state1 = state1**Rational(1,2)
+ return Tr((sqrt_state1 * state2 * sqrt_state1)**Rational(1, 2)).doit()
@@ -2,7 +2,7 @@
from sympy.matrices.matrices import Matrix
from sympy.core.trace import Tr
from sympy.external import import_module
-from sympy.physics.quantum.density import Density, entropy
+from sympy.physics.quantum.density import Density, entropy, fidelity
from sympy.physics.quantum.state import Ket, Bra, TimeDepKet
from sympy.physics.quantum.qubit import Qubit
from sympy.physics.quantum.qapply import qapply
@@ -201,3 +201,82 @@ def _eval_trace(self, bra, **options):
t = Tr(d)
assert t.doit() == 1
+
+
+def test_fidelity():
+ #test with kets
+ up = JzKet(S(1)/2, S(1)/2)
+ down = JzKet(S(1)/2, -S(1)/2)
+ updown = (S(1)/sqrt(2))*up + (S(1)/sqrt(2))*down
+
+ #check with matrices
+ up_dm = represent(up * Dagger(up))
+ down_dm = represent(down * Dagger(down))
+ updown_dm = represent(updown * Dagger(updown))
+
+ assert abs(fidelity(up_dm, up_dm) - 1) < 1e-3
+ assert fidelity(up_dm, down_dm) < 1e-3
+ assert abs(fidelity(up_dm, updown_dm) - (S(1)/sqrt(2))) < 1e-3
+ assert abs(fidelity(updown_dm, down_dm) - (S(1)/sqrt(2))) < 1e-3
+
+ #check with density
+ up_dm = Density([up, 1.0])
+ down_dm = Density([down, 1.0])
+ updown_dm = Density([updown, 1.0])
+
+ assert abs(fidelity(up_dm, up_dm) - 1) < 1e-3
+ assert abs(fidelity(up_dm, down_dm)) < 1e-3
+ assert abs(fidelity(up_dm, updown_dm) - (S(1)/sqrt(2))) < 1e-3
+ assert abs(fidelity(updown_dm, down_dm) - (S(1)/sqrt(2))) < 1e-3
+
+ #check mixed states with density
+ updown2 = (sqrt(3)/2)*up + (S(1)/2)*down
+ d1 = Density([updown, 0.25], [updown2, 0.75])
+ d2 = Density([updown, 0.75], [updown2, 0.25])
+ assert abs(fidelity(d1, d2) - 0.991) < 1e-3
+ assert abs(fidelity(d2, d1) - fidelity(d1, d2)) < 1e-3
+
+
+ #using qubits/density(pure states)
+ state1 = Qubit('0')
+ state2 = Qubit('1')
+ state3 = (S(1)/sqrt(2))*state1 + (S(1)/sqrt(2))*state2
+ state4 = (sqrt(S(2)/3))*state1 + (S(1)/sqrt(3))*state2
+
+ state1_dm = Density([state1, 1])
+ state2_dm = Density([state2, 1])
+ state3_dm = Density([state3, 1])
+
+ assert fidelity(state1_dm, state1_dm) == 1
+ assert fidelity(state1_dm, state2_dm) == 0
+ assert abs(fidelity(state1_dm, state3_dm) - 1/sqrt(2)) < 1e-3
+ assert abs(fidelity(state3_dm, state2_dm) - 1/sqrt(2)) < 1e-3
+
+ #using qubits/density(mixed states)
+ d1 = Density([state3, 0.70], [state4, 0.30])
+ d2 = Density([state3, 0.20], [state4, 0.80])
+ assert abs(fidelity(d1, d1) - 1) < 1e-3
+ assert abs(fidelity(d1, d2) - 0.996) < 1e-3
+ assert abs(fidelity(d1, d2) - fidelity(d2, d1)) < 1e-3
+
+ #TODO: test for invalid arguments
+ # non-square matrix
+ mat1 = [[0, 0],
+ [0, 0],
+ [0, 0]]
+
+ mat2 = [[0, 0],
+ [0, 0]]
+ raises(ValueError, lambda: fidelity(mat1, mat2))
+
+ # unequal dimensions
+ mat1 = [[0, 0],
+ [0, 0]]
+ mat2 = [[0, 0, 0],
+ [0, 0, 0],
+ [0, 0, 0]]
+ raises(ValueError, lambda: fidelity(mat1, mat2))
+
+ # unsupported data-type
+ x, y = 1, 2 #random values that is not a matrix
+ raises(ValueError, lambda: fidelity(x, y))