Standard Quantum Tomography¶
Quara supports several estimators (linear, least squares, maximum-likelihood) for standard quantum tomography of state, POVM, and gate. Here we briefly explain how to use them with examples on 1-qubit system. Quara supports much larger systems. However, performances (computation speed and constraint feasibility) of numerical optimization solver implemented are quite low, and we do not recommend to perform tomographic data-processing for systems larger than 1-qubit. Improvements on the performances are necessary tasks at the development of Quara in the near future.
[1]:
import numpy as np
# quara
from quara.objects.composite_system_typical import generate_composite_system
from quara.objects.tester_typical import (
generate_tester_states,
generate_tester_povms,
)
from quara.objects.qoperation_typical import generate_qoperation
from quara.protocol.qtomography.standard.standard_qst import StandardQst
from quara.protocol.qtomography.standard.standard_povmt import StandardPovmt
from quara.protocol.qtomography.standard.standard_qpt import StandardQpt
from quara.protocol.qtomography.standard.linear_estimator import LinearEstimator
from quara.protocol.qtomography.standard.loss_minimization_estimator import (
LossMinimizationEstimator,
)
from quara.loss_function.weighted_probability_based_squared_error import (
WeightedProbabilityBasedSquaredError,
WeightedProbabilityBasedSquaredErrorOption,
)
from quara.minimization_algorithm.projected_gradient_descent_backtracking import (
ProjectedGradientDescentBacktracking,
ProjectedGradientDescentBacktrackingOption,
)
Quantum State Tomography (1-qubit)¶
First, we consider quantum state tomography on 1-qubit system. We prepare a system.
[2]:
mode = "qubit"
num = 1
c_sys = generate_composite_system(mode=mode, num=num)
print(c_sys)
elemental_systems:
[0] 0 (system_id=4572582864)
dim: 2
basis:
(array([[0.70710678+0.j, 0. +0.j],
[0. +0.j, 0.70710678+0.j]]), array([[0. +0.j, 0.70710678+0.j],
[0.70710678+0.j, 0. +0.j]]), array([[0.+0.j , 0.-0.70710678j],
[0.+0.70710678j, 0.+0.j ]]), array([[ 0.70710678+0.j, 0. +0.j],
[ 0. +0.j, -0.70710678+0.j]]))
In order to perform data-pocessing of standard quantum tomography, we need some pre-knowledge information on the tomographic experiment. Suppose that we performed projective measurements along with x, y, and z axes. We call them a tester.
[3]:
# Testers
names = ["x", "y", "z"]
testers = generate_tester_povms(c_sys=c_sys, names = names)
for i, tester in enumerate(testers):
print("Tester ", i, ":\n", tester)
Tester 0 :
Type:
Povm
Dim:
2
Number of outcomes:
2
Vecs:
[[ 0.70710678 0.70710678 0. 0. ]
[ 0.70710678 -0.70710678 0. 0. ]]
Tester 1 :
Type:
Povm
Dim:
2
Number of outcomes:
2
Vecs:
[[ 0.70710678 0. 0.70710678 0. ]
[ 0.70710678 0. -0.70710678 0. ]]
Tester 2 :
Type:
Povm
Dim:
2
Number of outcomes:
2
Vecs:
[[ 0.70710678 0. 0. 0.70710678]
[ 0.70710678 0. 0. -0.70710678]]
We generate a StandardQst object. We set seed for generating pseudo-data. The seed is unnecessary if you do not generate pseudo-data.
[4]:
seed = 7896
qst = StandardQst(testers, on_para_eq_constraint=True, schedules="all", seed=seed)
Suppose that the state of the system to be estimated is \(|A\rangle\). We call it the true state.
[5]:
mode = "state"
name = "a"
true = generate_qoperation(mode=mode, name=name, c_sys=c_sys)
print(true)
Type:
State
Dim:
2
Vec:
[0.70710678 0.5 0.5 0. ]
The true probability distribution is given as follows.
[6]:
prob_dists = qst.calc_prob_dists(true)
print(prob_dists)
[[0.85355339 0.14644661]
[0.85355339 0.14644661]
[0.5 0.5 ]]
Suppose that we performed tomographic experiment and repeated each sub-experiment with tester element 1000 times.
[7]:
#qst.reset_seed()
num_data = 1000
empi_dists = qst.generate_empi_dists(state=true, num_sum=num_data)
for f in empi_dists:
print(f)
(1000, array([0.864, 0.136]))
(1000, array([0.844, 0.156]))
(1000, array([0.49, 0.51]))
Here empi_dists is a set of empirical distributions, which is a pair of the repetition number and relative frequencies.
When we choose a linear estimator, the estimate is calculated as follows.
[8]:
estimator = LinearEstimator()
result = estimator.calc_estimate(qtomography=qst, empi_dists=empi_dists, is_computation_time_required=True)
estimate = result.estimated_qoperation
print(estimate)
print("is estimate physical? : ", estimate.is_physical())
print("\nEigenvalues are: ", estimate.calc_eigenvalues())
Type:
State
Dim:
2
Vec:
[ 0.70710678 0.51477374 0.48648947 -0.01414214]
is estimate physical? : False
Eigenvalues are: [1.0009311329913524, -0.0009311329913525568]
An eigenvalue of the estimated density matrix is negative, which violates the requirement of positive-semidefiniteness on density matrix. This kind of violation can occur when we choose a linear estimator. In order to avoid the problem, we need to perform a constraint optimization at the data-processing.
When we choose a constraint least squares estimator, the estimate is calculated as folllows. Here we choose a projected gradient descent back-tracking as the numerical optimization algorithm.
[9]:
estimator = LossMinimizationEstimator()
loss = WeightedProbabilityBasedSquaredError()
loss_option = WeightedProbabilityBasedSquaredErrorOption("identity")
algo = ProjectedGradientDescentBacktracking()
algo_option = ProjectedGradientDescentBacktrackingOption(mode_stopping_criterion_gradient_descent="sum_absolute_difference_variable", num_history_stopping_criterion_gradient_descent=1)
result = estimator.calc_estimate(qtomography=qst, empi_dists=empi_dists, loss=loss, loss_option=loss_option, algo=algo, algo_option=algo_option, is_computation_time_required=True)
estimate = result.estimated_qoperation
print(estimate)
print("\nis estimate physical? : ", estimate.is_physical())
print("\nEigenvalues are: ", estimate.calc_eigenvalues())
Type:
State
Dim:
2
Vec:
[ 0.70710678 0.51381689 0.4855852 -0.01411585]
is estimate physical? : False
Eigenvalues are: [1.000000019575682, -1.9575682325134736e-08]
An eigenvalue of the estimated density matrix is negative as same as observed at the case of linear estimator, although the strength of the violation is reduced from O(10^{-4}) to O(10^{-8}). This is due to the low constraint feasibility of the numerical solver implemented at the current version of Quara, which must be improved in the near future.
POVM tomography (1-qubit)¶
Next, we consider POVM tomography on 1-qubit system.
[10]:
# Composite System
mode = "qubit"
num = 1
c_sys = generate_composite_system(mode=mode, num=num)
#print(c_sys)
# Testers
names = ["x0", "y0", "z0", "z1"]
testers = generate_tester_states(c_sys=c_sys, names = names)
povmt = StandardPovmt(testers, num_outcomes=2, on_para_eq_constraint=True, schedules="all")
[11]:
mode = "povm"
name = "z"
true = generate_qoperation(mode=mode, name=name, c_sys=c_sys)
print(true)
Type:
Povm
Dim:
2
Number of outcomes:
2
Vecs:
[[ 0.70710678 0. 0. 0.70710678]
[ 0.70710678 0. 0. -0.70710678]]
[12]:
prob_dists = povmt.calc_prob_dists(true)
print(prob_dists)
[[5.00000000e-01 5.00000000e-01]
[5.00000000e-01 5.00000000e-01]
[1.00000000e+00 2.22044605e-16]
[0.00000000e+00 1.00000000e+00]]
[13]:
num_data = 1000
empi_dists = povmt.generate_empi_dists(povm=true, num_sum=num_data)
for f in empi_dists:
print(f)
(1000, array([0.488, 0.512]))
(1000, array([0.516, 0.484]))
(1000, array([1., 0.]))
(1000, array([0., 1.]))
[14]:
estimator = LinearEstimator()
result = estimator.calc_estimate(qtomography=povmt, empi_dists=empi_dists, is_computation_time_required=True)
estimate = result.estimated_qoperation
print(estimate)
print("\nis estimate physical? : ", estimate.is_physical())
print("\nEigenvalues are:", estimate.calc_eigenvalues())
Type:
Povm
Dim:
2
Number of outcomes:
2
Vecs:
[[ 0.70710678 -0.01697056 0.02262742 0.70710678]
[ 0.70710678 0.01697056 -0.02262742 -0.70710678]]
is estimate physical? : False
Eigenvalues are: [[1.000399840127872, -0.0003998401278722573], [1.000399840127872, -0.0003998401278720354]]
[15]:
estimator = LossMinimizationEstimator()
loss = WeightedProbabilityBasedSquaredError(4)
loss_option = WeightedProbabilityBasedSquaredErrorOption("identity")
algo = ProjectedGradientDescentBacktracking()
algo_option = ProjectedGradientDescentBacktrackingOption()
result = estimator.calc_estimate(qtomography=povmt, empi_dists=empi_dists, loss=loss, loss_option=loss_option, algo=algo, algo_option=algo_option, is_computation_time_required=True)
estimate = result.estimated_qoperation
print(estimate)
print("\nis estimate physical? : ", estimate.is_physical())
print("\nEigenvalues are:", estimate.calc_eigenvalues())
Type:
Povm
Dim:
2
Number of outcomes:
2
Vecs:
[[ 0.7071068 -0.01694356 0.02259135 0.70654269]
[ 0.70710676 0.01694356 -0.02259135 -0.70654269]]
is estimate physical? : False
Eigenvalues are: [[1.0000000327444887, -6.646159317336142e-17], [0.9999999999999999, -3.27444887981e-08]]
Quantum Process Tomography (1-qubit)¶
Finally, we consider quantum process tomography on 1-qubit system.
[16]:
# Composite System
mode = "qubit"
num = 1
c_sys = generate_composite_system(mode=mode, num=num)
#print(c_sys)
# Testers
names_states = ["x0", "y0", "z0", "z1"]
testers_states = generate_tester_states(c_sys=c_sys, names = names_states)
names_povms = ["x", "y", "z"]
testers_povms = generate_tester_povms(c_sys=c_sys, names=names_povms)
qpt = StandardQpt(states=testers_states, povms=testers_povms, on_para_eq_constraint=True, schedules="all")
[17]:
mode = "gate"
name = "hadamard"
true = generate_qoperation(mode=mode, name=name, c_sys=c_sys)
print(true)
Type:
Gate
Dim:
2
HS:
[[ 1. 0. 0. 0.]
[ 0. 0. 0. 1.]
[ 0. 0. -1. 0.]
[ 0. 1. 0. 0.]]
[18]:
prob_dists = qpt.calc_prob_dists(true)
print(prob_dists)
[[0.5 0.5]
[0.5 0.5]
[1. 0. ]
[0.5 0.5]
[0. 1. ]
[0.5 0.5]
[1. 0. ]
[0.5 0.5]
[0.5 0.5]
[0. 1. ]
[0.5 0.5]
[0.5 0.5]]
[19]:
num_data = 1000
empi_dists = qpt.generate_empi_dists(gate=true, num_sum=num_data)
for f in empi_dists:
print(f)
(1000, array([0.507, 0.493]))
(1000, array([0.513, 0.487]))
(1000, array([1., 0.]))
(1000, array([0.47, 0.53]))
(1000, array([0., 1.]))
(1000, array([0.471, 0.529]))
(1000, array([1., 0.]))
(1000, array([0.505, 0.495]))
(1000, array([0.476, 0.524]))
(1000, array([0., 1.]))
(1000, array([0.505, 0.495]))
(1000, array([0.524, 0.476]))
[20]:
estimator = LinearEstimator()
result = estimator.calc_estimate(qtomography=qpt, empi_dists=empi_dists, is_computation_time_required=True)
estimate = result.estimated_qoperation
print(estimate)
print("\nis estimate physical? : ", estimate.is_physical())
evals, evecs = np.linalg.eigh(estimate.to_choi_matrix())
print("\nEigenvalues are:", evals)
Type:
Gate
Dim:
2
HS:
[[ 1. 0. 0. 0. ]
[ 0. 0.014 -0.06 1. ]
[ 0.01 0.016 -1.01 0. ]
[ 0. 1. -0.058 -0.048]]
is estimate physical? : False
Eigenvalues are: [-0.05842999 0.01079319 0.04147905 2.00615775]
[21]:
estimator = LossMinimizationEstimator()
loss = WeightedProbabilityBasedSquaredError()
loss_option = WeightedProbabilityBasedSquaredErrorOption("identity")
algo = ProjectedGradientDescentBacktracking()
algo_option = ProjectedGradientDescentBacktrackingOption()
result = estimator.calc_estimate(qtomography=qpt, empi_dists=empi_dists, loss=loss, loss_option=loss_option, algo=algo, algo_option=algo_option, is_computation_time_required=True)
estimate = result.estimated_qoperation
print(estimate)
print("\nis estimate physical? : ", estimate.is_physical())
evals, evecs = np.linalg.eigh(estimate.to_choi_matrix())
print("\nEigenvalues are:", evals)
Type:
Gate
Dim:
2
HS:
[[ 1. 0. 0. 0. ]
[-0.01050078 0.03961536 -0.02267407 0.98598833]
[ 0.01001083 -0.01491283 -0.98970188 -0.01249551]
[-0.00595741 0.98724549 -0.02176788 -0.03950884]]
is estimate physical? : False
Eigenvalues are: [-2.51576912e-08 -1.54360436e-08 1.73812841e-02 1.98261876e+00]
[ ]: