[PYTHON] Basics of Quantum Information Theory: Quantum State Tomography

\def\bra#1{\mathinner{\left\langle{#1}\right|}} \def\ket#1{\mathinner{\left|{#1}\right\rangle}} \def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}

Introduction

I have covered various topics under the title of "Basics of Quantum Information Theory", but in order to express quantum states, "density operator" has always appeared. In particular, it is a very convenient theoretical tool that must be used when considering quantum systems that are in a mixed state, so we cannot live without it anymore. It's so convenient that when you have an unknown physical system, you'll be tempted to ask, "What density operator can this quantum state be represented by?" (Please!). In fact, it can be estimated using a technique called "quantum state tomograpy". In this article, I will first explain the method, and then use the quantum calculation simulator qlazy to estimate the density operator for any quantum state. I will.

The following documents were used as references.

  1. Nielsen, Chan "Quantum Computer and Quantum Communication (3)" Ohmsha (2005)
  2. Ishizaka, Ogawa, Kawachi, Kimura, Hayashi "Introduction to Quantum Information Science" Kyoritsu Shuppan (2012)
  3. Tomita "Quantum Information Engineering" Morikita Publishing (2017)

What is quantum state tomography?

About mixed state (review)

Before estimating the density operator, let's review a little about how mixed states appear. This is to understand when this estimation method is useful. I think that there are two patterns as shown below.

[1] is a mixed state that appears when, for example, Alice has a device that can fire multiple patterns of pure states, and uses it to stochastically select pure states and fire them one after another. Imagine a photon launcher that can set the polarization angle as you like. Or, if you measure some pure state and forget the measurement result, you will get a mixed state of this pattern [^ 1].

Representing a pure state stochastic ensemble as $ \ {p_i, \ ket {\ psi_i} \} $, the density operator $ \ rho $

[^ 1]: By the way, such measurement is called "non-selective measurement". When I read a textbook on quantum information, I often come across the words "forget" and "lose measurement results." I don't think there is such a stupid experimenter, but let's think of it as Aya of words (when I was a student, I sometimes made non-selective measurements, but sweat). On the contrary, the measurement that remembers the measurement result firmly is called "selective measurement". In that case, the state of the system becomes pure.

\rho = \sum_{i} p_i \ket{\psi_i} \bra{\psi_i}  \tag{1}

Can be expressed as.

[2] is a mixed state that appears by paying attention to the partial system S of the closed system (= pure state) as a whole [^ 2]. The partial system S is generally considered to interact with the rest of the other systems (environmental system E), so it becomes an open system (= mixed state). Or, when running a quantum circuit on a quantum computer, if you look at some of the qubits, they are generally in a mixed state [^ 3].

[^ 2]: It's an interesting part of the quantum state. When the whole system of the pure state whose state is fixed is divided, the "entanglement" is unwound and it becomes a mixed state. In other words, simply splitting a system with zero entropy causes entropy in a partial system for some reason!

[^ 3]: Quantum calculation is a unitary transformation to the pure state as a whole, so when you look at the whole quantum state, it always keeps the pure state, but if you pay attention to some of the qubits, It is in a mixed state. However, if the partial qubit of interest is in an independent = untangled = product state with other qubits, it is in a pure state.

When the whole system is $ \ ket {\ Psi} $, if you write the attention to the partial system with a mathematical formula,

\rho = Tr_E (\ket{\Psi} \bra{\Psi})  \tag{2}

It will be.

Density operator estimation (1 quantum state)

Now that we have briefly reviewed the specific appearance of the mixed state, let's explain how to estimate the density operator at that time. But here is one premise. Unfortunately, if the $ \ rho $ condition appears only once and you have only one chance to measure it, it is impossible to estimate. While the number obtained by measurement is one, there are many numbers for describing the density operator, so it becomes a well-posed problem and cannot be estimated. Therefore, it is assumed that the state of $ \ rho $ can be generated repeatedly and can be measured each time. The basic idea is to estimate the density operator by collecting a large number of measurement results and integrating them.

Let's first consider the density operator for one quantum state.

Only in the case of one quantum, the density operator could be expressed as a linear sum of Pauli operators [^ 4]. This expression is the basis for explaining the estimation method, so let's take a closer look at how it was derived.

[^ 4]: In the Previous article that covered "quantum channels", I explained the Bloch display of the density operator.

The Pauli operator was defined as a matrix like this:

\sigma_0 = I = 
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix}, \space
\sigma_1 =  
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}, \space
\sigma_2 =
\begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix}, \space
\sigma_3 =
\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}  \tag{3}

As is clear from this definition, the Pauli operator is the Hermitian operator. And,

\frac{1}{2} Tr(\sigma_i \sigma_j) = \delta_{ij}  \tag{4}

It has the property of. By the way, $ Tr (AB) $ that appears here is a kind of inner product when the matrices $ A and B $ are regarded as vectors, and is sometimes called "Hilbert-Schmidt inner product" (later mentioned). Please remember that it will come).

Noting this property of the Pauli operator, any Hermitian operator $ A $ defined on a two-dimensional Hilbert space is $ \ {u_i \} \ space (i = 0,1,2,3). You can see that $ can be expressed as an arbitrary real number by the following linear sum.

A = u_0 \sigma_0 + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 = \sum_{i=0}^{3} u_i \sigma_i \tag{5}

It's easy, so let's prove it a little.

[Proof]

First, whether the right-hand side of equation (3) is Hermitian.

\begin{align}
A &=
u_0
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix} +
u_1
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix} +
u_2
\begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix} +
u_3
\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}  \\
&= 
\begin{pmatrix}
u_0 & 0 \\
0 & u_0
\end{pmatrix} +
\begin{pmatrix}
0 & u_1 \\
u_1 & 0
\end{pmatrix} +
\begin{pmatrix}
0 & -i u_2 \\
i u_2 & 0
\end{pmatrix} +
\begin{pmatrix}
u_3 & 0 \\
0 & -u_3
\end{pmatrix}  \\
&=
\begin{pmatrix}
u_0 + u_3 & u_1 - iu_2 \\
u_1 + iu_2 & u_0 - u_3
\end{pmatrix}  \tag{6}
\end{align}

And obviously $ A ^ {\ dagger} = A $, so $ A $ is Hermitian.

Next is the opposite proof. That is, whether any Hermitian operator can be expressed as in equation (5).

Assuming that $ A $ is Hermitian, using any real number $ a, b, c, d $,

A =
\begin{pmatrix}
a & b-ic \\
b+ic & d
\end{pmatrix}  \tag{7}

It can be expressed as. No matter what $ a, b, c, d $ you bring here

\begin{align}
a &= u_0 + u_3 \\
b &= u_1 \\
c &= u_2 \\
d &= u_0 - u_3 \tag{8}
\end{align}

Since there are always real numbers $ u_0, u_1, u_2, u_3 $ that satisfy, equation (7)

\begin{align}
A &=
\begin{pmatrix}
u_0 + u_3 & u_1 - iu_2 \\
u_1 + iu_2 & u_0 - u_3
\end{pmatrix} \\
&=
u_0
\begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix} +
u_1
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix} +
u_2
\begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix} +
u_3
\begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}  \\
&= u_0 \sigma_0 + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 \\
&= \sum_{i=0}^{3} u_i \sigma_i \tag{9}
\end{align}

Then, equation (5) was proved.

(End of proof)

Now, let's say something a little differently about what equation (5) means. That is, the Elmeet operator defined in the two-dimensional Hilbert space can be regarded as a vector in the linear space stretched by the four Pauli operators when they are orthogonal bases in the sense of the Hilbert-Schmidt product. You can say that.

Hmm? What does "vector" mean when you say "operator (= matrix)"? ?? You may hear the question, but please soften your head and follow along. Roughly speaking, if you can define the inner product, have an orthogonal basis, and have something that can be represented by its linear sum, you can say anything as a "vector".

So, the one-quantum density operator $ \ rho $ was exactly the Hermitian operator in the two-dimensional Hilbert space, so it can be expanded like the equation (5) just proved.

\rho = u_0 \sigma_0 + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 = \sum_{i=0}^{3} u_i \sigma_i \tag{10}

However, the density operator has a constraint of $ Tr (\ rho) = 1 $, and $ Tr (\ sigma_1) = Tr (\ sigma_2) = Tr (\ sigma_3) = 0 $, so $ u_0 = Tr (\ rho) = 1 $. That is, one variable is reduced,

\rho = I + u_1 \sigma_1 + u_2 \sigma_2 + u_3 \sigma_3 = I + \sum_{i=1}^{3} u_i \sigma_i \tag{11}

Can be written. Here, paying attention to the nature of Pauli operator in Eq. (4),

u_i = \frac{1}{2} Tr(\sigma_i \rho)  \tag{12}

Therefore, equation (11) is, after all,

\begin{align}
\rho &= \frac{1}{2} (Tr(\rho) I + Tr(\sigma_1 \rho) \sigma_1 + Tr(\sigma_2 \rho) \sigma_2+ Tr(\sigma_3 \rho) \sigma_3) \\
&= \frac{1}{2} \sum_{i=0}^{3} Tr(\sigma_i \rho) \sigma_i  \tag{13}
\end{align}

It will be.

Now it's time to talk about how to estimate the density operator.

Take a closer look at equation (13). $ Tr (\ sigma_i \ rho) $ for $ i = 1,2,3 $ was the expected value when measuring the observable $ \ sigma_i $. $ \ sigma_1, \ sigma_2, \ sigma_3 $ are like spins in the X, Y, and Z axes, respectively. Therefore, for the flying quantum, the measurement in the X-axis direction is performed many times, the measured value (+1, -1) is recorded to obtain the expected value, and then the measurement in the Y-axis direction is performed. Is executed many times to obtain the expected value in the same manner, and then the measurement in the Z-axis direction is performed many times to obtain the expected value, which is substituted into Eq. (13). $ Tr (\ sigma_0 \ rho) $ $ \ sigma_0 $ is an identity matrix, and the trace of $ \ rho $ is 1, so it doesn't matter whether it exists or not, this term should be set to 1. So, equation (13) can be calculated perfectly.

I have just imagined the quantum state physically emitted by Alice, but as another example, you may want to estimate the density operator for one qubit in a quantum computer. It is basically the same, but you need to be careful only where you switch the measurement direction of the qubit each time. If you are using a simulator or cloud service that only supports Z-direction measurements, you need to bite a quantum gate that changes the measurement direction before measuring in the Z-direction.

Measurement in the X direction

--H--M--

The measurement in the Y direction is

--S+--H--M--

Z direction measurement

--M--

I can go. If you imagine a Bloch sphere, you can see that it is simply rotating the axis.

Density operator estimation (2 quantum states)

The discussion so far can be extended to the density operator of two quantum states.

I explained that the density operator in the case of one quantum can be thought of as a vector in a linear space based on four Pauli operators. In the case of a two-quantum state, we can bring another linear space to create a Cartesian product space and define a vector on it. The basis of this Cartesian product is defined by the tensor product of the four Pauli operators. In other words

\{\sigma_i \otimes \sigma_j \} \space (i,j = 0,1,2,3)  \tag{14}

16 4X4 matrices are the basis of this Cartesian product space. As you can see, this basis is

\frac{1}{2^2} Tr[(\sigma_i \otimes \sigma_j)(\sigma_k \otimes \sigma_l)] = \delta_{ik} \delta_{jl} \tag{15}

Since the orthogonal condition is satisfied, the density operator of this two-quantum state is the same as before.

\rho = \sum_{i}^{3} \sum_{j}^{3} u_{ij} \sigma_i \otimes \sigma_j \tag{16}

It can be expanded on the basis using Pauli operator like. Note that equation (15), this $ u_ {ij} $ is

u_{ij} = \frac{1}{2^2} Tr[(\sigma_i \otimes \sigma_j) \rho]  \tag{17}

Since it can be calculated, equation (16) is

\rho = \frac{1}{2^2} \sum_{i=0}^{3} \sum_{j=0}^{3} Tr[(\sigma_i \otimes \sigma_j) \rho] \cdot \sigma_i \otimes \sigma_j  \tag{18}

It will be.

To estimate the density operator, it is necessary to perform a two-quantum measurement while switching each direction. There are 4X4 = 16 combinations of Pauli operators, so the expected value is calculated for each case. Specifically, first, for one combination of Pauli operators, obtain the measured values (+1, -1) of each of the two flying quanta. Record the product as a measurement of the entire two quanta. Repeat this many times to get the expected value. Now that we have the expected value for one combination of Pauli operators, we will do this for all combinations of Pauli operators (16 patterns). Since 16 expected values are obtained, the estimation of the density operator is completed by substituting them into equation (17). However, there is one case where all Pauli operators are unit operators, so it is not necessary to calculate the expected value. Leave it at 1 without saying whether it is present or not. that's all.

Density operator estimation (n quantum state)

Furthermore, the same discussion is OK for the density operator of the n quantum state. The density operator is

\rho = \frac{1}{2^n} \sum_{\mu_1=0}^{3} \sum_{\mu_2=0}^{3} \cdots \sum_{\mu_{n}=0}^{3} Tr[(\sigma_{\mu_1} \otimes \sigma_{\mu_2} \otimes \cdots \otimes \sigma_{\mu_{n}})\rho] \cdot \sigma_{\mu_1} \otimes \sigma_{\mu_2} \otimes \cdots \otimes \sigma_{\mu_{n}} \tag{19}

It can be expressed as It looks a little complicated, but if you look closely, you'll understand it.

The estimation method of the density operator is

Tr[(\sigma_{\mu_1} \otimes \sigma_{\mu_2} \otimes \cdots \otimes \sigma_{\mu_{n}})\rho]   \tag{20}

You just have to find the expected value. Since it is n quantum, the combination of Pauli operators is 4 to the nth power. For each combination, the direction corresponding to the Pauli operator assigned to each qubit is measured, the measured value of n quanta (the product of the measured values of each quantum) is obtained, and the expected value is obtained. Is good.

The method of estimating as described above is "quantum state tomography".

Estimated simulation

Implementation

Now, I would like to create an arbitrary quantum state ensemble and estimate the density operator by quantum state tomography [^ 5]. The quantum calculation simulator qlazy is equipped with a function that allows you to directly obtain the density operator from the quantum state ensemble by matrix calculation, even if you do not do such a tricky thing. It has been. However, since this time we are studying quantum state tomography, let's forget about it and estimate the density operator only by measuring the quantum state ensemble. Use the density operator calculated with qlazy as the reference value (true value) to measure the accuracy of the estimation. Accuracy is evaluated by fidelity (fidelity is also included as a function of qlazy).

[^ 5]: This simulation is an estimation of the mixed state that appears from the quantum state ensemble. However, I think that the estimation of the mixed state that appears by focusing on the pure state subsystem can also be done by slightly changing the program (?).

Now let's look at the whole Python code.

import random
import math
import numpy as np
from scipy.stats import unitary_group
from qlazypy import QState,DensOp

def random_qstate_ensemble(qubit_num, mimxed_num):

    qstate = []
    dim = 2**qubit_num
    for _ in range(mixed_num):
        vec_ini = np.zeros(dim)
        vec_ini[0] = 1.0
        mat = unitary_group.rvs(dim)
        vec = np.dot(mat, vec_ini)
        qstate.append(QState(vector=vec))
        
    prob = [random.random() for _ in range(mixed_num)]
    total = sum(prob)
    prob = [p/total for p in prob]

    return (qstate,prob)

def get_pauli_index(index, total):
    # ex) index = 9 = 1001 -> [10,01] -> [2,1] --reverse-> [1,2] = pauli_index
    #     '1' means X for 0th-qubit, '2' means Y for 1st-qubit

    pauli_index = [0]*int(math.log2(total)/2)
    count = 0
    while index > 0:
        pauli_index[count] = index%4
        index = index // 4
        count += 1
    pauli_index.reverse()

    return pauli_index

def make_pauli_product(index, total, pauli_mat):

    pauli_index = get_pauli_index(index, total)
    pauli_prod_dim = 2**len(pauli_index) 
    pauli_prod = np.array([1.0])
    for pid in pauli_index:
        pauli_prod = np.kron(pauli_prod, pauli_mat[pid])

    return pauli_prod
    
def make_densop(expect, qubit_num, pauli_mat):

    dim = 2**qubit_num
    measure_num = len(expect)
    matrix = np.zeros((dim,dim))
    for i in range(measure_num):
        pauli_prod = make_pauli_product(i,measure_num,pauli_mat)
        matrix = matrix + expect[i] * pauli_prod
    matrix = matrix / (2.0**qubit_num)
    
    return DensOp(matrix=matrix)

def calc_expect(prb):

    N = len(prb)
    val = np.zeros(N)
    for index in range(N):
        bin_list = [int(s) for s in list(format(index, 'b'))]
        val[index] = (-1)**sum(bin_list)  # odd -> -1, even -> +1

    return np.dot(prb, val)

def estimate_densop(prob,qstate,shots):

    pauli_mat = [np.eye(2),                   # = I
                 np.array([[0,1],[1,0]]),     # = X
                 np.array([[0,-1j],[1j,0]]),  # = Y
                 np.array([[1,0],[0,-1]])]    # = Z
    
    mixed_num = len(prob)
    qubit_num = qstate[0].qubit_num
    measure_num = 4**qubit_num

    for i in range(mixed_num):

        expect = {}
        for index in range(measure_num):

            pauli_index = get_pauli_index(index,measure_num)

            if index == 0:
                expect[index] = 1.0
                continue
                    
            qs = qstate[i].clone()

            for qid in range(len(pauli_index)):

                if pauli_index[qid] == 0:
                    continue
                elif pauli_index[qid] == 1:
                    qs.h(qid)
                elif pauli_index[qid] == 2:
                    qs.s_dg(qid).h(qid)
                else:
                    pass

            frq = qs.m(shots=shots).frq
            prb = np.array([f/shots for f in frq])
            expect[index] = calc_expect(prb)

            qs.free()

        de_tmp = make_densop(expect, qubit_num, pauli_mat)
        
        if i == 0:
            de_tmp.mul(prob[i])
            densop_est = de_tmp.clone()
        else:
            de_tmp.mul(prob[i])
            densop_est.add(de_tmp)

    de_tmp.free()

    return densop_est
    
if __name__ == '__main__':

    # settings
    qubit_num = 2
    mixed_num = 4
    shots = 100

    # quantum state ensemble (original)
    qstate, prob = random_qstate_ensemble(qubit_num, mixed_num)

    # density operator (original)
    densop_ori = DensOp(qstate=qstate, prob=prob)

    # density operator estimation only from quantum state ensemble
    # (quantum state tomography)
    densop_est = estimate_densop(prob,qstate,shots)

    print("** density operator (original)")
    densop_ori.show()
    print("** density operator (estimated)")
    densop_est.show()
    print("** fidelity =",densop_ori.fidelity(densop_est))

    for q in qstate:
        q.free()
    densop_ori.free()
    densop_est.free()

I will explain what you are doing. Look at the main processing section.

# settings
qubit_num = 2
mixed_num = 4
shots = 1000

This is a setting for simulation. qubit_num is the number of assumed quanta, mixed_num is the number of pure states to be mixed, and shots is the number of measurements.

# quantum state ensemble (original)
qstate, prob = random_qstate_ensemble(qubit_num, mixed_num)

Here, we are randomly creating a quantum state ensemble. qstate is a list of four quantum states. prob is a list of probabilities that determine the proportion of the four pure states to mix. Both the pure state and the probability are randomly determined. See the function definition for details.

# density operator (original)
densop_ori = DensOp(qstate=qstate, prob=prob)

This is where we are using qlazy to calculate the true value of the density operator we are thinking about. Finally, it is used as a reference value for accuracy evaluation.

# density operator estimation only from quantum state ensemble
# (quantum state tomography)
densop_est = estimate_densop(prob,qstate,shots)

This is the part of the density operator estimation that is the key to this time, so I will explain it in a little more detail. Look inside the function estimate_densop.

def estimate_densop(prob,qstate,shots):

    pauli_mat = [np.eye(2),                   # = I
                 np.array([[0,1],[1,0]]),     # = X
                 np.array([[0,-1j],[1j,0]]),  # = Y
                 np.array([[1,0],[0,-1]])]    # = Z
    
    mixed_num = len(prob)
    qubit_num = qstate[0].qubit_num
    measure_num = 4**qubit_num
    ...

The first pauli_mat is setting Pauli matrices. mixed_num and qubit_num are calculated again from the arguments. The number of Pauli operator combinations can be obtained from qubit_num with 4 ** qubit_num (this is stored in the variable measure_num).

The next for loop is

for i in range(mixed_num):

    expect = {}
    for index in range(measure_num):

        pauli_index = get_pauli_index(index,measure_num)
	...

It is doubled like. The outer for loop is for estimating the density operator for each pure state that is mixed. Later, we will weight the probabilities and add them together to make the final density operator [^ 5].

[^ 5]: If you want to get closer to a more realistic experiment, it would have been better to implement this part so that one of the pure states stochastically flies and measures once, and so on. However, I thought that the code would be long and difficult to understand, so I skipped the implementation a little.

The inner for loop is for calculating the expected value for the number of Pauli operator combinations. First, initialize the dictionary expect for storing the expected value corresponding to each combination as an empty dictionary. I repeat for the number of measure_num, but the combination number is in the variable index. Since this is just an integer value, it is awkward to control the measurement direction of each qubit. So, the function get_pauli_index converts it to a list of Pauli operator numbers. Convert the integer value to quaternary and return each digit as a list. See the function definition for details.

Next is the contents of that for loop.

if index == 0:
    expect[index] = 1.0
    continue
        
qs = qstate[i].clone()

for qid in range(len(pauli_index)):

    if pauli_index[qid] == 0:
        continue
    elif pauli_index[qid] == 1:
        qs.h(qid)
    elif pauli_index[qid] == 2:
        qs.s_dg(qid).h(qid)
    else:
        pass

frq = qs.m(shots=shots).frq
prb = np.array([f/shots for f in frq])
expect[index] = calc_expect(prb)

qs.free()

If index = 0, the expected value is set to 1 regardless of the presence or absence. If not,

qs = qstate[i].clone()

So, copy the i-th pure state and store it in qs as temporary [^ 6].

[^ 6]: Since it is a simulator, I copy it easily, but when doing it in an actual experiment or a real quantum computer, this part generates the same thing from the quantum generator many times, or the same quantum circuit It is necessary to perform the operation of getting the state through many times.

Next, the for loop came out again.

for qid in range(len(pauli_index)):
...

This is for the quantum gate to be applied before the measurement because each qubit is measured in order in the measurement direction corresponding to the Pauli operator number. The contents are as you can see above. When the loop is over

frq = qs.m(shots=shots).frq
prb = np.array([f/shots for f in frq])
expect[index] = calc_expect(prb)

qs.free()

It is said. Now you are doing the actual measurement. Since the shot was 1000, I will measure 1000 times [^ 7]. As a result, we get a frequency list, so we normalize it so that the sum is 1 for probability. Pass it to the function calc_expect to calculate the expected value. Since we are assuming 2 quantum states now, the probability is 4. Each of them represents a particular set of Pauli operators. For a particular set of Pauli operators, the two-quantum measurements can be calculated as the product of each (which can be +1 or -1). Once the four probabilities and four measurements are determined, the expected value can be calculated. The function calc_expect does that. See the function definition for details. After using qs, free the memory with qs.free () [^ 8].

[^ 7]: Again, the implementation is skipped by the processing unique to the simulator. In a realistic experiment, qs really has to be generated and measured many times.

[^ 8]: This process is important when using the same quantum state in a loop many times. qlazy uses the C language library inside, and memory allocation is also done in the C world. Therefore, the specification is that the quantum state and density operator should be explicitly released in memory. I left it to Python's \ _ \ _ del \ _ \ _ and tried a specification that does not explicitly release it, but I do not know when the memory will be released on the Python side, so I am suspicious. Since there were cases where it was done, we have made this specification.

de_tmp = make_densop(expect, qubit_num, pauli_mat)

if i == 0:
    de_tmp.mul(prob[i])
    densop_est = de_tmp.clone()
else:
    de_tmp.mul(prob[i])
    densop_est.add(de_tmp)

Once all four expected values are found, the function make_densop once calculates the density operator for this pure state (see the function definition for details on make_densop). Add the weights of the probabilities by the number of mixed numbers (4). The if statement is that process. After exiting the outer loop, we have the final densop_est, so we will return it.

Return to the main processing section

print("** density operator (original)")
densop_ori.show()
print("** density operator (estimated)")
densop_est.show()
print("** fidelity =",densop_ori.fidelity(densop_est))

Now, the elements of the real density operator (densop_ori) and the elements of the estimated density operator (densop_est) are displayed in order, and finally the fidelity of both is calculated and displayed.

result

The result of the simulation is shown.

The above program was for estimating the two-quantum state, but the first is the estimation of the one-quantum state. The parameters are

qubit_num = 1
mixed_num = 2
shots = 100

It was made. When I ran it

** density operator (original)
elm[0][0] = +0.3270-0.0000*i : 0.1069 |++
elm[0][1] = +0.1508-0.2138*i : 0.0684 |++
elm[1][0] = +0.1508+0.2138*i : 0.0684 |++
elm[1][1] = +0.6730+0.0000*i : 0.4530 |++++++
** density operator (estimated)
elm[0][0] = +0.3224+0.0000*i : 0.1040 |++
elm[0][1] = +0.1584-0.1887*i : 0.0607 |++
elm[1][0] = +0.1584+0.1887*i : 0.0607 |++
elm[1][1] = +0.6776+0.0000*i : 0.4591 |++++++
** fidelity = 0.9996155234243419

have become. Since the fidelity is 0.9996 ..., I think we were able to estimate it with considerable accuracy.

Next is the estimation of the two-quantum state. The parameters are

qubit_num = 2
mixed_num = 4
shots = 1000

I tried it. Since the element of the density operator becomes long, if only the fidelity is shown,

** fidelity = 0.9814099144563044

And the accuracy dropped a little. However, for the time being, I think that the purpose of studying quantum state tomography has been achieved, so I will leave it as good around here [^ 9].

[^ 9]: Quantum Information Engineering has the following description. "In the case of one qubit, it is relatively easy like this, but in order to determine the state of multiple qubits, it is necessary to determine the component of the density operator that increases in the order of the square of the number of qubits. Therefore, it is necessary. In addition, if the state happens to be a pure state, the rank of the density matrix will be 1 (leaving one diagonal element and all other elements will be 0), which is unnecessary. It seems that it is quite difficult to estimate the parameters. " When you actually simulate it, you can realize this kind of difficulty.

in conclusion

This time, we have taken up "quantum state tomography" for estimating the state, but there is also a method called "quantum process tomography" that estimates the "process" instead of the "state". It is a method of estimating the Klaus operator that represents changes in the physical system, and basically seems to be a method of repeatedly using this "state" tomography. However, I haven't studied yet, so please take this opportunity again.

Next time, I'm thinking about topics related to quantum communication. I would like to take up "data compression", which is a method for communicating by putting the quantum state itself on the quantum communication channel instead of classical information (although the schedule is undecided).

that's all

Recommended Posts

Basics of Quantum Information Theory: Quantum State Tomography
Basics of Quantum Information Theory: Entropy (2)
Basics of Quantum Information Theory: Data Compression (1)
Basics of Quantum Information Theory: Horebaud Limits
Basics of Quantum Information Theory: Trace Distance
Basics of Quantum Information Theory: Data Compression (2)
Basics of Quantum Information Theory: Topological Toric Code
Basics of Quantum Information Theory: Fault Tolerant Quantum Computation
Basics of Quantum Information Theory: Quantum Error Correction (Shor's Code)
Basics of Quantum Information Theory: Quantum Error Correction (CSS Code)
Basics of Quantum Information Theory: Quantum Error Correction (Stabilizer Code: 4)
Basics of Quantum Information Theory: Quantum Error Correction (Classical Linear Code)
Basics of Quantum Information Theory: Universal Quantum Calculation by Toric Code (1)
Basics of Quantum Information Theory: Logical Operation by Toric Code (Brading)
Read "Basics of Quantum Annealing" Day 5
Quantum computer implementation of 3-state quantum walk
Read "Basics of Quantum Annealing" Day 6
Basics of Tableau Basics (Visualization Using Geographic Information)
Basics of Python ①
Basics of python ①
Basics of Python scraping basics
# 4 [python] Basics of functions
Basics of network programs?
Basics of Perceptron Foundation
Basics of regression analysis
Basics of python: Output