Density-matrix simulation¶

This notebook shows how to get the density matrix corresponding to a noisy circuit.

We first define a quantum circuit:

In [1]:
from qat.lang.AQASM import Program, H
from qat.lang.AQASM.qftarith import QFT
prog = Program()
reg = prog.qalloc(4)
prog.apply(H, reg[1])
prog.apply(QFT(4), reg)
qft_circuit = prog.to_circ(inline=True) #convert program to circuit
qft_circuit.display()
No description has been provided for this image

We then create a hardware model with perfect gates and T1/T2 noise:

In [2]:
from qat.hardware import DefaultGatesSpecification, HardwareModel
from qat.quops import ParametricPureDephasing, ParametricAmplitudeDamping

ibm_gates_spec = DefaultGatesSpecification(gate_times={"H":60, "X":120, "Y":120, "S":1,"T":1,
                                                       'D-T':1, "Z":1, "RZ":lambda angle : 1,
                                                       "PH": lambda angle : 1, "CNOT":386,
                                                       "SWAP":150, "C-PH":lambda angle :150, "C-T":150, "C-S":150})
T1, T2 = 44000, 38900 #nanosecs
hardware_model = HardwareModel(ibm_gates_spec, 
                               idle_noise=[ParametricAmplitudeDamping(T_1=T1),
                                           ParametricPureDephasing(T_phi=1/(1/T2 - 1/(2*T1)))])

qft_circuit.display(hardware_model=hardware_model)

/tmp/ipykernel_41913/1413445602.py:13: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  qft_circuit.display(hardware_model=hardware_model)
No description has been provided for this image

Finally, we compute the density matrix and check that the final state is mixed, i.e $\mathrm{Tr}(\rho^2) < \mathrm{Tr}(\rho) = 1$

Note: if no hardware model is specified, we assume ideal gates.

In [3]:
from qat.qpus import NoisyQProc
from qat.noisy import compute_density_matrix
qpu = NoisyQProc(hardware_model=hardware_model)
rho = compute_density_matrix(qft_circuit, qpu)

import numpy as np
print("Tr(rho) = ", np.trace(rho).real)
print("Tr(rho^2) = ", np.trace(np.dot(rho, rho)).real)

Tr(rho) =  0.9999999999999989
Tr(rho^2) =  0.9794366320481338
In [4]:
from qat.noisy.plot import plot_density_matrix
%matplotlib inline
plot_density_matrix(rho)
No description has been provided for this image

Note: because of the high memory overhead of storing the density matrix, this method is limited to a small number of qbits.

Initializing the simulation with a density matrix (or a pure state)¶

One may want to start the computation with a state other than $|00\dots\rangle$. For this, one uses a StatePreparation gate in the circuit, and one stores the information about which state (be it pure or not) it corresponds to in the gate specification:

In [5]:
import numpy as np
from qat.core.circuit_builder.matrix_util import get_predef_generator
from qat.linalg import StatePreparation

# here we want to start from a pure state
psi0 = np.array([0.2, 0., 0.3, 0.5])
psi0 /= np.linalg.norm(psi0)


gate_durations = {"H": 60, "X": 120, "Y": 120, "S": 1, "T": 1,
                  "D-T": 1, "Z": 1, "RZ": lambda angle: 1,
                  "PH": lambda angle: 1, "C-PH": lambda angle: 1,
                  "CNOT": 386, "SWAP": 100, "STATE_PREPARATION": 10,
                  "RX": lambda angle: 1}

pred_generator = get_predef_generator()
pred_generator["STATE_PREPARATION"] = psi0
ibm_gates_spec = DefaultGatesSpecification(gate_durations, predef_generator=pred_generator)
hardware_model = HardwareModel(ibm_gates_spec)
qpu = NoisyQProc(hardware_model=hardware_model)    
    
nqbits = 2
prog = Program()
reg = prog.qalloc(nqbits)
# if one uses NoisyQProc, the value inside StatePreparation is not used
# it is superseded by that in the HW model
# one can thus put any normalized vector
dummy = np.random.rand(4)
dummy /= np.linalg.norm(dummy)
prog.apply(StatePreparation(dummy), reg)
# nonetheless, if one wants to initialize in a pure state and then compare
# a noisy evolution and a perfect evolution, one may want to put the "true"
# psi0 in the StatePreparation
# prog.apply(StatePreparation(psi0), reg)
circ = prog.to_circ()

rho_final = compute_density_matrix(circ, qpu)

print("rho = ", rho_final)
print("rho0 = ", np.outer(psi0, psi0.conj()))

rho =  [[0.10526316 0.         0.15789474 0.26315789]
 [0.         0.         0.         0.        ]
 [0.15789474 0.         0.23684211 0.39473684]
 [0.26315789 0.         0.39473684 0.65789474]]
rho0 =  [[0.10526316 0.         0.15789474 0.26315789]
 [0.         0.         0.         0.        ]
 [0.15789474 0.         0.23684211 0.39473684]
 [0.26315789 0.         0.39473684 0.65789474]]

Same, but with a mixed state as a starting point:

In [6]:
# here we want to start from a mixed state
rho_0 = 0.25 * np.identity(4)


gate_durations = {"H": 60, "X": 120, "Y": 120, "S": 1, "T": 1,
                  "D-T": 1, "Z": 1, "RZ": lambda angle: 1,
                  "PH": lambda angle: 1, "C-PH": lambda angle: 1,
                  "CNOT": 386, "SWAP": 100, "STATE_PREPARATION": 10,
                  "RX": lambda angle: 1}

pred_generator = get_predef_generator()
pred_generator["STATE_PREPARATION"] = rho_0
ibm_gates_spec = DefaultGatesSpecification(gate_durations, predef_generator=pred_generator)
hardware_model = HardwareModel(ibm_gates_spec)
qpu = NoisyQProc(hardware_model=hardware_model)    
    
nqbits = 2
prog = Program()
reg = prog.qalloc(nqbits)
dummy = np.random.rand(4)
dummy /= np.linalg.norm(dummy)
prog.apply(StatePreparation(dummy), reg)
circ = prog.to_circ()

rho_final = compute_density_matrix(circ, qpu)

print("rho = ", rho_final)

rho =  [[0.25 0.   0.   0.  ]
 [0.   0.25 0.   0.  ]
 [0.   0.   0.25 0.  ]
 [0.   0.   0.   0.25]]