Download Notebook

Povm and povm_typical

[1]:

import numpy as np
np.set_printoptions(linewidth=200)

Povm

Quantum measurement has two mathematical treatments. One is positive operator-valued measure (POVM), which can describe the effect of quantum measurement on the probability distribution of its measurement outcome only. The other is measurement apparatus, which can describe both of the effect on the probability distribution and states after the measurement. Current version of quara prepare a class for POVM only. A class for measurement apparatus will be added in the near future.

Povm can be represented by a list of POVM element as vecs.
Each emelemts is a linear combination of basis and represents the coefficients of this linear combination in the form of a numpy array.
Example:
Povm elements \(\Pi_0 = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ \end{bmatrix}\), \(\Pi_1 = \begin{bmatrix} 0 & 0 \\ 0 & 1 \\ \end{bmatrix}\).
Povm \(\Pi = \{ \Pi_x \}_{x=0,1}\). index \(x\) is called outcomes.
We can write \(\Pi_0= 1/\sqrt{2} \cdot I/\sqrt{2} + 0 \cdot X/\sqrt{2} + 0 \cdot Y/\sqrt{2} + 1/\sqrt{2} \cdot Z/\sqrt{2}\) and \(\Pi_1= 1/\sqrt{2} \cdot I/\sqrt{2} + 0 \cdot X/\sqrt{2} + 0 \cdot Y/\sqrt{2} - 1/\sqrt{2} \cdot Z/\sqrt{2}\).
In this case vecs of \(\Pi\) is a list of \([ 1/\sqrt{2}, 0, 0, 1/\sqrt{2} ]\) and \([ 1/\sqrt{2}, 0, 0, -1/\sqrt{2} ]\).

The methods for generating a State includes the following:

  • Generate from povm_typical module

  • Generate Povm object directly

Generate from povm_typical module by specifying CompositeSystem and povm name (ex. “z”).

[2]:

from quara.objects.composite_system_typical import generate_composite_system
from quara.objects.povm_typical import generate_povm_from_name

c_sys = generate_composite_system("qubit", 1)
povm = generate_povm_from_name("z", c_sys)
print(povm)

Type:
Povm

Dim:
2

Number of outcomes:
2

Vecs:
[[ 0.70710678  0.          0.          0.70710678]
 [ 0.70710678  0.          0.         -0.70710678]]

Generate Povm object directly using CompositeSystem and a list of numpy array.

[3]:

from quara.objects.composite_system import CompositeSystem
from quara.objects.elemental_system import ElementalSystem
from quara.objects.matrix_basis import get_normalized_pauli_basis
from quara.objects.povm import Povm

basis = get_normalized_pauli_basis(1)
e_sys = ElementalSystem(0, basis)
c_sys = CompositeSystem([e_sys])
vec1 = np.array([1, 0, 0, 1]) / np.sqrt(2)
vec2 = np.array([1, 0, 0, -1]) / np.sqrt(2)
vecs = [vec1, vec2]
povm = Povm(c_sys, vecs)
print(povm)

Type:
Povm

Dim:
2

Number of outcomes:
2

Vecs:
[[ 0.70710678  0.          0.          0.70710678]
 [ 0.70710678  0.          0.         -0.70710678]]

specific properties

The property vecs of Povm is a list of a numpy array specified by the constructor argument vecs.

[4]:

povm = Povm(c_sys, vecs)
print(f"vecs: \n{povm.vecs}")

vecs:
(array([0.70710678, 0.        , 0.        , 0.70710678]), array([ 0.70710678,  0.        ,  0.        , -0.70710678]))

The property dim of Povm is a square root of the size of element of vecs.

[5]:

print(f"dim: {povm.dim}")
print(f"square root of the size of element of vecs: {int(np.sqrt(len(vecs[0])))}")

dim: 2
square root of the size of element of vecs: 2

The property num_outcomes of Povm is the number of POVM elements.

[6]:

print(f"num_outcomes: {povm.num_outcomes}")
print(f"number of POVM elements: {len(vecs)}")

num_outcomes: 2
number of POVM elements: 2

functions to check constraints

The is_eq_constraint_satisfied() function returns True, if and only if \(\sum_x \Pi_x = I\), where \(\Pi_x\) are Povm elements.

[7]:

from operator import add
from functools import reduce

print(f"is_eq_constraint_satisfied(): {povm.is_eq_constraint_satisfied()}")
print(f"sum of matrices: \n{reduce(add, povm.matrices())}")

is_eq_constraint_satisfied(): True
sum of matrices:
[[1.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

The is_ineq_constraint_satisfied() function returns True, if and only if all \(\Pi_x\) are positive semidifinite matrices.

[8]:

print(f"is_eq_constraint_satisfied(): {povm.is_eq_constraint_satisfied()}")
print(f"matrix(0): \n{povm.matrix(0)}")
print(f"matrix(1): \n{povm.matrix(1)}")

is_eq_constraint_satisfied(): True
matrix(0):
[[1.+0.j 0.+0.j]
 [0.+0.j 0.+0.j]]
matrix(1):
[[0.+0.j 0.+0.j]
 [0.+0.j 1.+0.j]]

projection functions

calc_proj_eq_constraint() function calculates the projection of Povm on equal constraint.
Let num_outcomes = \(m\), dim = \(d\), and \(\tilde{I} = I/\sqrt{d}\). Then \(\sum_{x=0}^{m-1} \Pi_x = I = \sqrt{d} \cdot \tilde{I}\).
Therefore \(\sum_{x=0}^{m-1} |\Pi_x\rangle\rangle = |\tilde{I}\rangle\rangle = [ \sqrt{d}, 0, \dots, 0 ]^T\). Thus, the last element of Povm \(\Pi_{m-1}\) can be determined, depending on other elements, as follows:

  • \(|\Pi_x\rangle\rangle = [a_{x,0}, \dots, a_{x,d^2-1}]^T\)

  • \(\bar{a_\alpha} := \frac{1}{m} \sum_{x=0}^{m-1} (|\Pi_x\rangle\rangle)_\alpha\), where \(\alpha = 0, \dots, d^2-1\).

  • \(c_\alpha := \begin{cases} \sqrt{d}/m & (\alpha = 0) \\ 0 & (\alpha = 1, \dots, d^2-1) \end{cases}\)

  • \(a_{x,\alpha}^\prime := a_{x,\alpha} - \bar{a_\alpha} + c_\alpha\)

  • The projection of Povm \(|\tilde{\Pi}_x\rangle\rangle\) is \([ a_{x,0}^\prime, \dots, a_{x,d^2-1}^\prime ]^T\).

This function replaces the last element of vecs with \(1/\sqrt{d}\), where \(d\) is dim.

[9]:

vec1 = np.array([1.0, 2.0, 3.0, 4.0])
vec2 = np.array([5.0, 6.0, 7.0, 8.0])
vecs = [vec1, vec2]

[10]:

povm = Povm(c_sys, vecs, is_physicality_required=False)
proj_povm = povm.calc_proj_eq_constraint()
print(f"vecs: \n{proj_povm.vecs}")

vecs:
(array([-1.29289322, -2.        , -2.        , -2.        ]), array([2.70710678, 2.        , 2.        , 2.        ]))

calc_proj_ineq_constraint() function calculates the projection of Povm \(\{|\Pi_x\rangle\rangle\}_x\) on inequal constraint as follows:

  • For each \(x\), calculate the following:

  • Executes singular value decomposition on the elements of Povm \(\Pi_x\), \(\Pi_x = U \Lambda U^{\dagger}\), where \(\Lambda = \text{diag}[\lambda_0, \dots , \lambda_{d-1}]\), and \(\lambda_{i} \in \mathbb{R}\).

  • \(\lambda^{\prime}_{i} := \begin{cases} \lambda_{i} & (\lambda_{i} \geq 0) \\ 0 & (\lambda_{i} < 0) \end{cases}\)

  • \(\Lambda^{\prime} = \text{diag}[\lambda^{\prime}_0, \dots , \lambda^{\prime}_{d-1}]\)

  • \(\Pi_x^{\prime} = U \Lambda^{\prime} U^{\dagger}\)

  • The projection of Povm is \(\{|\Pi_x^{\prime} \rangle\rangle\}_x\).

[11]:

povm = Povm(c_sys, vecs, is_physicality_required=False)
proj_povm = povm.calc_proj_ineq_constraint()
print(f"vecs: \n{proj_povm.vecs}")

vecs:
(array([3.1925824 , 1.18569534, 1.77854301, 2.37139068]), array([8.60327781, 4.22884788, 4.93365586, 5.63846384]))

functions to transform parameters

to_stacked_vector() function returns a one-dimensional numpy array of all variables. This is equal to a concatenated vector of elements of vecs.

[12]:

print(f"to_stacked_vector(): {povm.to_stacked_vector()}")
print(f"vecs: \n{povm.vecs}")

to_stacked_vector(): [1. 2. 3. 4. 5. 6. 7. 8.]
vecs:
(array([1., 2., 3., 4.]), array([5., 6., 7., 8.]))
If Povm \(\Pi = \{ \Pi_0, \dots , \Pi_{m-1} \}\) and on_para_eq_constraint is True, then the last element of Povm \(\Pi_{m-1}\) is equal to \(I - \sum_{x=0}^{m-2} \Pi_x\). Thus, Povm is characterized by vecs[0],…, vecs[m-2].
Therefore, to_var() function returns a vector combining elements of of vecs[0],…, vecs[m-2], where on_para_eq_constraint is True.
[13]:

vec1 = np.array([1, 0, 0, 1]) / np.sqrt(2)
vec2 = np.array([1, 0, 0, -1]) / np.sqrt(2)
vecs = [vec1, vec2]

[14]:

# on_para_eq_constraint=True
povm = Povm(c_sys, vecs, on_para_eq_constraint=True)
print(f"to_var(): {povm.to_var()}")

to_var(): [0.70710678 0.         0.         0.70710678]

[15]:

# on_para_eq_constraint=False
povm = Povm(c_sys, vecs, on_para_eq_constraint=False)
print(f"to_var(): {povm.to_var()}")

to_var(): [ 0.70710678  0.          0.          0.70710678  0.70710678  0.          0.         -0.70710678]

functions to generate special objects

[16]:

zero_povm = povm.generate_zero_obj()
print(f"zero: \n{zero_povm.vecs}")
origin_povm = povm.generate_origin_obj()
print(f"origin: \n{origin_povm.vecs}")

zero:
(array([0., 0., 0., 0.]), array([0., 0., 0., 0.]))
origin:
(array([0.70710678, 0.        , 0.        , 0.        ]), array([0.70710678, 0.        , 0.        , 0.        ]))

supports arithmetic operations

[17]:

vec11 = np.array([1.0, 2.0, 3.0, 4.0])
vec12 = np.array([5.0, 6.0, 7.0, 8.0])
povm1 = Povm(c_sys, [vec11, vec12], is_physicality_required=False)
vec21 = np.array([9.0, 10.0, 11.0, 12.0])
vec22 = np.array([13.0, 14.0, 15.0, 16.0])
povm2 = Povm(c_sys, [vec21, vec22], is_physicality_required=False)

print(f"sum: \n{(povm1 + povm2).vecs}")
print(f"subtraction: \n{(povm1 - povm2).vecs}")
print(f"right multiplication: \n{(2 * povm1).vecs}")
print(f"left multiplication: \n{(povm1 * 2).vecs}")
print(f"division: \n{(povm1 / 2).vecs}")

sum:
(array([10., 12., 14., 16.]), array([18., 20., 22., 24.]))
subtraction:
(array([-8., -8., -8., -8.]), array([-8., -8., -8., -8.]))
right multiplication:
(array([2., 4., 6., 8.]), array([10., 12., 14., 16.]))
left multiplication:
(array([2., 4., 6., 8.]), array([10., 12., 14., 16.]))
division:
(array([0.5, 1. , 1.5, 2. ]), array([2.5, 3. , 3.5, 4. ]))

calc_gradient functions

Calculates gradient of Povm with variable index.

[18]:

grad_povm = povm.calc_gradient(0)
print(f"vecs: {grad_povm.vecs}")

vecs: (array([1., 0., 0., 0.]), array([0., 0., 0., 0.]))

convert_basis function

Returns vecs converted to the specified basis.

[19]:

from quara.objects.matrix_basis import get_comp_basis

povm = generate_povm_from_name("z", c_sys)
converted_vecs = povm.convert_basis(get_comp_basis())
print(f"vecs: {converted_vecs}")

vecs: [array([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]), array([0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j])]

generate_mprocess function

Generates MProcess from this Povm.

[20]:

print(povm.generate_mprocess())

Type:
MProcess

Dim:
2

HSs:
[array([[0.5, 0. , 0. , 0.5],
       [0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. ],
       [0.5, 0. , 0. , 0.5]]), array([[ 0.5,  0. ,  0. , -0.5],
       [ 0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ],
       [-0.5,  0. ,  0. ,  0.5]])]

ModeSampling:
False

some utility functions

[21]:

print(f"vec(0): \n{povm.vec(0)}")
print(f"matrices(): \n{povm.matrices()}")
print(f"matrix(0): \n{povm.matrix(0)}")
print(f"is_hermitian(): \n{povm.is_hermitian()}")
print(f"is_positive_semidefinite(): \n{povm.is_positive_semidefinite()}")
print(f"is_identity_sum(): \n{povm.is_identity_sum()}")
print(f"calc_eigenvalues(): \n{povm.calc_eigenvalues()}")

vec(0):
[0.70710678 0.         0.         0.70710678]
matrices():
[matrix([[1.+0.j, 0.+0.j],
        [0.+0.j, 0.+0.j]]), matrix([[0.+0.j, 0.+0.j],
        [0.+0.j, 1.+0.j]])]
matrix(0):
[[1.+0.j 0.+0.j]
 [0.+0.j 0.+0.j]]
is_hermitian():
True
is_positive_semidefinite():
True
is_identity_sum():
True
calc_eigenvalues():
[[0.9999999999999998, 0.0], [0.9999999999999998, 0.0]]

povm_typical

generate_povm_object_from_povm_name_object_name() function in povm_typical module can easily generate objects related to Povm.
The generate_povm_object_from_povm_name_object_name() function has the following arguments:

  • The string that can be specified for povm_name can be checked by executing the get_povm_names() function. The tensor product of povm_name “a”, “b” is written “a_b”.

  • object_name can be the following string:

  • “pure_state_vectors” - list of vector of pure states.

  • “matrices” - matrices of vector of pure states.

  • “vectors” - list of vectorized matrices.

  • “povm” - Povm object.

  • c_sys - CompositeSystem of objects related to Povm. Specify when object_name is “povm”.

  • basis - MatrixBasis of objects related to Povm. Specify when object_name is “vectors”.

  • is_physicality_required - Whether the generated object is physicality required, by default True.

[22]:

from quara.objects.povm_typical import (
    get_povm_names,
    generate_povm_object_from_povm_name_object_name,
)

#get_povm_names()

object_name = “pure_state_vectors”

[23]:

vecs = generate_povm_object_from_povm_name_object_name("z", "pure_state_vectors")
print(vecs)

[array([1.+0.j, 0.+0.j]), array([0.+0.j, 1.+0.j])]

object_name = “matrices”

[24]:

matrices = generate_povm_object_from_povm_name_object_name("z", "matrices")
print(matrices)

[array([[1.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j]]), array([[0.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j]])]

object_name = “vectors”

[25]:

basis = get_normalized_pauli_basis(1)
vectors = generate_povm_object_from_povm_name_object_name("z", "vectors", basis=basis)
print(vectors)

[array([0.70710678, 0.        , 0.        , 0.70710678]), array([ 0.70710678,  0.        ,  0.        , -0.70710678])]

object_name = “povm”

[26]:

c_sys = generate_composite_system("qubit", 1)
povm = generate_povm_object_from_povm_name_object_name("z", "povm", c_sys=c_sys)
print(povm)

Type:
Povm

Dim:
2

Number of outcomes:
2

Vecs:
[[ 0.70710678  0.          0.          0.70710678]
 [ 0.70710678  0.          0.         -0.70710678]]

Download Notebook