Quantum state tomography

An essential tool in the field of quantum information, quantum state tomography allows us to gain information about an unknown state as long as we can prepare it repeatedly.

In this blog post, we are going to explore the concept of quantum state tomography. Tomography is, in the general sense, imaging by sections through the use of some kind of wave. In our case, the image we obtain is of a quantum state. And, in a sense, we are using a wave whenever we measure a qubit, which happens a lot when doing this.

Tomography of a human brain. Image available here.

I dove into this topic thanks to QOSF’s monthly challenges, and you can find my solution to this particular challenge here.

Goal

Our goal here is to reconstruct the density matrix ρ of an unknown quantum state, assuming we have an efficient way of preparing this state a large amount of times.

Single qubit

For a single qubit, doing this is quite simple. Many methods to do this can be found in the review paper Quantum State Tomography of a Single Qubit: Comparison of Methods by Roman Schmied. However, to achieve this, I took the general method presented in Maximum Likelihood, Minimum Effort by John A. Smolin, Jay M. Gambetta, and Graeme Smith, and first implemented it for a single qubit.

To start, we need to remember some facts about measurement operators. The important here is that the expectation value of some measurement operator with respect to some state ρ is equal to the trace of these two matrices multiplied (remember that here ρ is the density matrix, not the statevector).

Expectation value of M is equal to the trace of the state times the operator

Now, remember that an arbitrary density matrix for the state of a mixed or pure qubit is given as a linear combination of the Pauli matrices as follows:

General density matrix

Where each r is some real number and whenever the sum of each squared is 1 the density matrix corresponds to a pure state. And one last thing we need to remember is the following relations for the trace of Pauli matrices, where index 0 corresponds to the identity matrix and the other 3 go as usual.

Trace relations of Pauli matrices

Remember that delta here represents the Kronecker delta. Thus, we can see that all Pauli matrices except the identity have a trace of 0, and the identity has a trace of 2. Furthermore, the product of any two different Pauli matrices is also 0, and the product of any Pauli matrix with itself is 2 (since they are all unitary and thus their own inverse).

With the equations introduced above, we can deduce the following.

Relation between expectation value and r parameters of general density matrix

Note that most simplifications here were thanks to the fact that trace is a linear mapping. Also note that the same can be done for the other two Pauli matrices, the derivation goes exactly the same way. Thus, this shows that the expectation values of the Pauli matrices have a direct relation with the density matrix of the state.

We also know that the expectation values can be calculated through repeated measurement. This can be done by projecting the state into each basis, running the circuit a large amount of times, and using the resulting counts to get the expectation value. This last step is done by subtracting the times 1 was measured from the times 0 was measured and dividing by the total number of shots. This can be done in Qiskit with the following code:

input_state = np.array(input_state)
circuits = []
    
# Z basis measurement
qc = QuantumCircuit(1)
qc.initialize(input_state, 0)
qc.measure_all()
circuits.append(qc)
    
# X basis measurement
qc = QuantumCircuit(1)
qc.initialize(input_state, 0)
qc.h(0)
qc.measure_all()
circuits.append(qc)
    
# Y basis measurement
qc = QuantumCircuit(1)
qc.initialize(input_state, 0)
qc.sdg(0)
qc.h(0)
qc.measure_all()
circuits.append(qc)
    
# run circuits
backend = Aer.get_backend('qasm_simulator')
job     = execute(circuits, backend, shots=samples)
counts  = job.result().get_counts()
    
# get expectation values, ordered the same way as the circuits
expectation_vals = []
for count in counts:
    m = (count.get('0', 0) - count.get('1', 0)) / samples
    expectation_vals.append(m)

This snippet will give you the expectation value stored in the expectation_vals array, ordered as X, Y, Z. Now that we have the expectation value of each Pauli matrix and we know that these come up in the general formula of the density matrix ρ, we can calculate our desired density matrix quite easily.

You may think we are done here. The density matrix as we have it gives good fidelity results, but remember we are only working with one qubit at the moment. In this case, increasing the number of shots to simulate the circuit helps increase the fidelity. However, this may not be reasonable with multi-qubit states.

Multi-qubit states

Therefore, to make the density matrix more accurate in cases where we have more qubits, the paper by John A. Smolin, Jay M. Gambetta, and Graeme Smith cited above introduces an optimization algorithm to increase the fidelity of the state. The technicalities of this algorithm are out of the scope of this algorithm, but they can be found in the linked paper, which includes a good summary in Fast algorithm for Subproblem 1.

Such algorithm can be implemented in code as follows:

def mu_optimize(mu, n):
    # calculate eigenvalues of µ matrix
    eigen = np.linalg.eig(mu)
    vals  = eigen[0]
    vecs  = eigen[1].transpose()

    # order eigenvalues from largest to smallest
    eig_vals = sorted(vals, reverse=True)
    idx      = []
    for val in eig_vals:
        idx.append(np.where(vals == val)[0][0])
    eig_vecs = []
    for i in idx:
        eig_vecs.append(vecs[i])
        
    # calculate eigenvalues of the density matrix
    accumulator = 0
    lamb_vals   = [None] * len(eig_vals)
    for i in range(len(eig_vals) - 1, -1, -1):
        if eig_vals[i] + (accumulator / (i + 1)) >= 0:
            for j in range(i + 1):
                lamb_vals[j] = eig_vals[j] +  (accumulator / (i + 1))
            break
        else:
            lamb_vals[i] = 0
            accumulator  += eig_vals[i]

    # calculate density matrix
    predicted_state = np.zeros((2 ** n, 2 ** n), 'complex')
    for idx, lamb_val in enumerate(lamb_vals):
        predicted_state += lamb_vals[idx] * np.outer(eig_vecs[idx], eig_vecs[idx].conj())
    
    return predicted_state

Here, the matrix μ is the density matrix ρ we obtained in the previous section. In the paper where this optimization algorithm is found, μ is the noisy density matrix which turns into our guess for ρ after going through the optimization procedure.

Another thing to note is the number of operators we need to get the expectation value of for multi-qubit states. For the single-qubit case, we only had to do the three Pauli matrices. Howerver, for n-qubit states, we have to consider all the possible tensor products of length n formed by the Pauli matrices along with the identity matrix. Thus, performing tomography for an n-qubit system requires 3^n operators to be measured. A good discussion about this can be found in this QCSE post.

Essentialy, it boils down to the fact that the density matrix space for n dimension is spanned by all the possible 4^n tensor products of length n made from the identity and Pauli matrices. But, if 4^n tensor products are needed, why do we only need to measure 3^n operators? Well, we only need to measure 3^n operators since any experiment that involves measuring an operator that includes the identity matrix is redundant with another experiment that has any Pauli matrix instead of that identity.

For example, suppose that you measure a two-qubit state in the ZZ basis, and from these results you want to get the expectation value of ZI. This can be done as follows:

Getting expectation value of ZI through measuring ZZ

Here, p represents the probability of the qubits collapsing into the state indicated by each index. Remember that |0⟩ has an eigenvalue of +1 and |1⟩ an eigenvalue of -1. And when measuring multi-qubit states, the eigenvalue of the whole system is the eigenvalues of the subsystems multiplied. That’s why we take |11⟩ to have eigenvalue +1 when calculating ⟨ZZ⟩. But, when calculating ⟨ZI⟩, we take the eigenvalue of the qubit corresponding to the identity operator to always be +1. That’s why |11⟩ is subtracted in the equation for ⟨ZI⟩.

The procedure above can be easily extended to any number of qubits, and, given the measurement results for the 3^n Pauli operators, we are able to find the expectation value for the 4^n Pauli plus identity tensor products, enabling us to perform state tomography. For the exact code that performs this and the remaining steps for n-qubit state tomography, go to level 2 of my solution to QOSF’s tomography challenge.

Using this method for general state tomography, I was able to obtain the following results:

Results for n-qubit state tomography (1 ≤ n ≤ 5)

As you can see, the more qubits in the system, the more shots needed to achieve a desirable fidelity (although most are already above 0.96 from the beginning).

It is worth exploring other methods for quantum state tomography, as they trade computational resources for fidelity, or the other way around. Just as in many other applications of quantum computing, the tradeoff desired depends on what you are trying to achieve. For example, an interesting approach is through a variational algorithm using the SWAP test as cost function. For more on this, checkout the paper Variational Quantum Circuits for Quantum State Tomography by Yong Liu et al., or my attempt at it here (scroll down to level 3).

If you enjoyed this post, remember to follow Quantum Untangled for more similar content!