Soluzione

Progetto 1: Canale di depolarizzazione

Soluzione

# Imports
import numpy as np
import matplotlib.pyplot as plt

from qiskit import QuantumRegister, QuantumCircuit
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime.fake_provider import FakeGuadalupe
from qiskit_experiments.library.tomography import StateTomography, StateTomographyAnalysis, MitigatedStateTomography, MitigatedTomographyAnalysis
from qiskit.quantum_info import DensityMatrix, Statevector, partial_trace, state_fidelity

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 5
      2 import numpy as np
      3 import matplotlib.pyplot as plt
----> 5 from qiskit import QuantumRegister, QuantumCircuit
      6 from qiskit_aer import AerSimulator
      7 from qiskit_ibm_runtime.fake_provider import FakeGuadalupe

ModuleNotFoundError: No module named 'qiskit'

Attività 1 (1 punto)

Crea una funzione che dia come risultato un circuito quantistico che implementi un canale di depolarizzazione con parametro $p$ su uno specifico system di qubit, usando tre qubit ausiliari ancillae = [a1, a2, a3].

def depolarizing_channel(q, p, system, ancillae):
    """Dà come risultato un QuantumCircuit che implementa il canale di depolarizzazione su q[sistema]
    
    Args:
        q (QuantumRegister): il registro da usare per il circuito
        p (float): la probabilità per il canale fra 0 e 1
        system (int): indice del qubit di sistema
        ancillae (list): lista di indici per i qubit ausilirari
        
    Returns:
        Un oggetto QuantumCircuit 
    """
    
    dc = QuantumCircuit(q)
   
    # 
    theta = 1/2 * np.arccos(1-2*p)
    
    #
    dc.ry(theta, q[ancillae[0]])
    dc.ry(theta, q[ancillae[1]])
    dc.ry(theta, q[ancillae[2]])

    dc.cx(q[ancillae[0]], q[system])
    dc.cy(q[ancillae[1]], q[system])
    dc.cz(q[ancillae[2]], q[system])

    return dc

Attività 2 (1 punto)

Scrivi uno stato iniziale (initial_state) per un circuito che prepari il system qubit in uno stato iniziale con popolazioni e coerenze non nulle (parti reali e immaginarie)

# Creiamo il circuito quantistico
q = QuantumRegister(4, name='q')

# Indice del qubit di sistema
system = 1

# Indici dei qubit ausiliari
ancillae = [0, 2, 3]

# Prepariamo il qubit in uno stato con coerenza e con popolazioni diverse 
initial_state = QuantumCircuit(q)
initial_state.u(np.pi/4, np.pi/4, 0, q[system])
initial_state.draw('mpl')

Attività 3 (4 punti)

Preparare l'esperimento di tomografia di stato con cui analizzeremo lo stato generato dal circuito e dal rumore depolarizzante. In questo caso, ricostruiremo le matrici di densità e le fidelità di stato.

Per diversi valori di $p \in [0, 1]$:
1. Combina initial_state e depolarizing_channel in un circuito.
2. Prepara ed esegui StateTomography da qiskit_experiments.library.tomography per eseguire la tomografia solo sul qubit di sistema utilizzando un simulatore.
  - Calcola anche la fedeltà rispetto allo stato del qubit di sistema in initial_state (senza il canale di depolarizzazione). Potrebbe essere necessario partial_trace da qiskit.quantum_info.
  - Calcolare le barre di errore per le fedeltà con il bootstrapping). A tal fine, utilizzare StateTomographyAnalysis da qiskit_experiments.library.tomography con alcuni argomenti e fornirlo a StateTomography.
3. Raccogli le matrici di densità, le fedeltà e gli errori delle fedeltà.
Descrivi brevemente con parole tue cosa fa StateTomography. Cosa bisogna misurare? Come ci assicuriamo che la matrice di densità sia fisica? (Tomografia dello stato quantistico nella documentazione di qiskit-experiments e Sistemi quantistici aperti con Qiskit potrebbero essere d'aiuto. )

# Per esempio, consideriamo 10 valori di p spaziati ugualmente
p_values = np.linspace(0, 1, 10)

# Creiamo una lista di circuiti per ciascun valore di p
tomography_circuits = []

for p in p_values:
    circ = initial_state.compose(depolarizing_channel(q, p, system, ancillae))
    tomography_circuits.append(circ)

# Calcoliamo la fedeltà rispetto a questo stato
psi = partial_trace(Statevector(initial_state), [q for q in range(initial_state.num_qubits) if q != system])

# Facciamo boostrapping per gli errori di fedeltà
analysis = StateTomographyAnalysis()
analysis.set_options(target_bootstrap_samples=20)

# Facciamo tomografia di stato nella simulazione e salviamo i risultati 
backend = AerSimulator()
# backend = FakeGuadalupe()

# backend = AerSimulator(noise_model=noise_model)
rhos = []
fidelities = []
errors = []
for circ in tomography_circuits:
    st = StateTomography(circ, backend=backend, measurement_indices=[system], analysis=analysis, target=psi)
    data = st.run()
    rhos.append(data.analysis_results('state').value)
    fidelities.append(data.analysis_results('state_fidelity').value.n)
    errors.append(data.analysis_results('state_fidelity').value.s)

La tomografia dello stato quantistico misura lo stato preparato utilizzando una base di misura tomograficamente completa. Grazie a questi risultati è possibile ricostruire la matrice di densità completa dello stato. Tuttavia, utilizzando direttamente questi risultati si potrebbero ottenere stati non fisici; pertanto, Qiskit ricorre al metodo della massima verosimiglianza per individuare una matrice di densità fisica che massimizzi le probabilità in linea con i risultati.

Attività 4 (4 punti)

Calcola numericamente la matrice di densità esatta del qubit del sistema dopo il canale di depolarizzazione in funzione di $p$.
Traccia i valori di $\rho_{00}$, $\rho_{11}$, $\mathrm{Re}(\rho_{01})$, $\mathrm{Im}(\rho_{01})$ in funzione di $p$ e confrontali con la previsione analitica.
Trova numericamente le fedeltà esatte del qubit del sistema dopo il canale di depolarizzazione in funzione di $p$.
Traccia sia le fedeltà esatte che quelle simulate del qubit del sistema rispetto allo stato iniziale del sistema in funzione di $p$. Aggiungi ai grafici le barre di errore calcolate nell'attività 3.

Fatti salvi gli errori statistici dovuti al numero finito di shot, i punti simulati dovrebbero essere vicini alla previsione analitica. Per la fedeltà, le barre di errore coprono 1 deviazione standard (~68%).

tomo_rhos = np.zeros((2,2,len(p_values)), dtype=complex)

for (i, res) in enumerate(rhos):
    tomo_rhos[:,:,i] = res

# Risultati simulati
plt.plot(p_values, np.real(tomo_rhos[0,1,:]),"C0*", label='Re $\\rho_{01}$')
plt.plot(p_values, np.imag(tomo_rhos[0,1,:]),"C1*", label='Im $\\rho_{01}$')
plt.plot(p_values, np.real(tomo_rhos[0,0,:]),"C2x", label='$\\rho_{00}$')
plt.plot(p_values, np.real(tomo_rhos[1,1,:]),"C3x", label='$\\rho_{11}$')

# Predizione teorica

# Otteniamo l'operatore densità dello stato iniziale
rho0 = partial_trace(Statevector(initial_state), ancillae).data

plt.plot(p_values, np.real(rho0[0,1])*(1-p_values), "C0", linewidth=.5)
plt.plot(p_values, np.imag(rho0[0,1])*(1-p_values), "C1", linewidth=.5)
plt.plot(p_values, 0.5*p_values + np.real(rho0[0,0])*(1-p_values), "C2", linewidth=.5)
plt.plot(p_values, 0.5*p_values + np.real(rho0[1,1])*(1-p_values), "C3", linewidth=.5)

plt.xlabel('p')
plt.ylabel('$\\rho_{xx}$')
plt.legend()

plt.title("SIMULATION Depol. channel. Full tomo. $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$")

Text(0.5, 1.0, 'SIMULATION Depol. channel. Full tomo. $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$')

Le espressioni esatte degli elementi della matrice di densità relative all'azione del canale di depolarizzazione derivano dall'azione che viene applicata alla matrice di densità. Con probabilità p, lo stato rimane intatto; con probabilità (1-p), lo stato subisce un errore. Partiamo dall'espressione della mappa fornita nel progetto.

$$ \begin{align*} \mathcal{E} (\rho) &= \left[1-\frac 3 4 p \right] \rho + \frac{p}{4} \sum_i \sigma_i \rho \sigma_i \\ &= \left[1-\frac 3 4 p \right] \rho + \frac{p}{4} (\sigma_x \rho \sigma_x + \sigma_y \rho \sigma_y + \sigma_z \rho \sigma_z) \\ &= \left[1- p \right] \rho + \frac{p}{4} (\sigma_0 \rho \sigma_0 + \sigma_x \rho \sigma_x + \sigma_y \rho \sigma_y + \sigma_z \rho \sigma_z) \\ &= \left[1- p \right] \rho + \frac{p}{4} (2 I) \\ \end{align*} $$

L'ultimo passaggio viene da: $$ \begin{align*} \sigma_0 \rho \sigma_0 + \sigma_x \rho \sigma_x + \sigma_y \rho \sigma_y + \sigma_z \rho \sigma_z&= \begin{bmatrix} \rho_{00} & \rho_{01} \\ \rho_{10} & \rho_{11} \end{bmatrix}+ \sigma_x \begin{bmatrix} \rho_{00} & \rho_{01} \\ \rho_{10} & \rho_{11} \end{bmatrix} \sigma_x + \sigma_y \begin{bmatrix} \rho_{00} & \rho_{01} \\ \rho_{10} & \rho_{11} \end{bmatrix} \sigma_y + \sigma_z \begin{bmatrix} \rho_{00} & \rho_{01} \\ \rho_{10} & \rho_{11} \end{bmatrix}\sigma_z \\ &= \begin{bmatrix} \rho_{00} & \rho_{01} \\ \rho_{10} & \rho_{11} \end{bmatrix}+\begin{bmatrix} \rho_{11} & \rho_{10} \\ \rho_{01} & \rho_{00} \end{bmatrix}+\begin{bmatrix} \rho_{11} & -\rho_{10} \\ \rho_{01} & \rho_{00} \end{bmatrix} + \begin{bmatrix} \rho_{00} & -\rho_{01} \\ -\rho_{10} & -\rho_{11} \end{bmatrix} \\ &= \begin{bmatrix} 2(\rho_{00} + \rho_{11}) & 0 \\ 0 & 2(\rho_{00} + \rho_{11}) \end{bmatrix} \\ &= 2I \end{align*} $$ visto che per una matrice di densità $\rho$, $\rho_{00} + \rho_{11} = \mathrm{Tr}[\rho] = 1$.

Quindi, puoi scrivere il canale di depolarizzazione come: $$\mathcal{E} (\rho) = (1-p) \rho + p \frac{I}{2}$$ Molto probabilmente vedrai più spesso questa forma, ad esempio su Wikipedia. Le due forme sono equivalenti per i sistemi a singolo qubit.

Se guardiamo agli elementi della matrice densità dopo l'azione del canale $\rho'$, otteniamo $$\begin{align*} \rho'_{00} &=(1-p)\rho_{00} + 0.5p \\ \rho'_{11} &=(1-p)\rho_{11} + 0.5p \\ \rho'_{01} &=(1-p)\rho_{01}. \end{align*} $$

# Risultati simulati
exact_fidelities = []
for circ in tomography_circuits:
    exact_fidelities.append(state_fidelity(psi, partial_trace(DensityMatrix(circ), ancillae)))

plt.plot(p_values, fidelities, "C0", label="Estimated fidelity")
plt.fill_between(p_values, np.array(fidelities) - np.array(errors), np.array(fidelities) + np.array(errors), alpha=.3)
plt.plot(p_values, exact_fidelities, "C1", label="Exact fidelity")

plt.xlabel('p')
plt.ylabel('Fidelity')
plt.legend()

plt.title("SIMULATION Depol. channel. Fidelity $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$")

Text(0.5, 1.0, 'SIMULATION Depol. channel. Fidelity $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$')

Attività facoltativa

Esegui tutte le attività su un dispositivo reale con mitigazione del rumore e confronta i risultati con quelli della simulazione. A tal fine puoi utilizzare MitigatedStateTomography invece di StateTomography. Si noti che MitigatedTomographyAnalysis funziona in modo leggermente diverso da StateTomographyAnalysis. Inoltre, è possibile ottenere il backend del dispositivo reale utilizzando il codice seguente (le primitive come Sampler sono relativamente nuove, quindi non sono ancora supportate in qiskit-experiments). Infine, assicurati di aver creato e salvato il tuo account IBM Quantum!

# Per esempio, consideriamo 10 valori di p spaziati ugualmente
p_values = np.linspace(0, 1, 10)

# Creiamo una lista di circuiti per ciascun valore di p
tomography_circuits = []

for p in p_values:
    circ = initial_state.compose(depolarizing_channel(q, p, system, ancillae))
    tomography_circuits.append(circ)

# Calcoliamo la fedeltà rispetto a questo stato
psi = partial_trace(Statevector(initial_state), ancillae)

# Facciamo boostrapping per gli errori di fedeltà
state_tomo = StateTomographyAnalysis()
state_tomo.set_options(target_bootstrap_samples=20)

analysis = MitigatedTomographyAnalysis(tomography_analysis=state_tomo)
analysis.set_options(target=psi, unmitigated_fit=True)

# Facciamo tomografia di stato nella simulazione e salviamo i risultati 
backend = FakeGuadalupe()

rhos = []
fidelities = []
errors = []

rhos_mit = []
fidelities_mit = []
errors_mit = []

for circ in tomography_circuits:
    st = MitigatedStateTomography(circ, measurement_indices=[system])
    data = st.run(backend=backend, analysis=analysis)
    rhos_mit.append(data.analysis_results('state')[0].value)
    fidelities_mit.append(data.analysis_results('state_fidelity')[0].value.n)
    errors_mit.append(data.analysis_results('state_fidelity')[0].value.s)
    rhos.append(data.analysis_results('state')[1].value)
    fidelities.append(data.analysis_results('state_fidelity')[1].value.n)
    errors.append(data.analysis_results('state_fidelity')[1].value.s)

Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.
Adding a job from a backend (AerSimulator('aer_simulator'
             noise_model=<NoiseModel on ['reset', 'x', 'measure', 'id', 'sx', 'cx']>)) that is different than the current backend (fake_guadalupe). The new backend will be used, but service is not changed if one already exists.

tomo_rhos = np.zeros((2,2,len(p_values)), dtype=complex)

for (i, res) in enumerate(rhos):
    tomo_rhos[:,:,i] = res

tomo_rhos_mit = np.zeros((2,2,len(p_values)), dtype=complex)

for (i, res) in enumerate(rhos_mit):
    tomo_rhos_mit[:,:,i] = res

# Risultati simulati
plt.plot(p_values, np.real(tomo_rhos[0,1,:]), label='Re $\\rho_{01}$')
plt.plot(p_values, np.imag(tomo_rhos[0,1,:]), label='Im $\\rho_{01}$')
plt.plot(p_values, np.real(tomo_rhos[0,0,:]), label='$\\rho_{00}$')
plt.plot(p_values, np.real(tomo_rhos[1,1,:]), label='$\\rho_{11}$')

plt.plot(p_values, np.real(tomo_rhos_mit[0,1,:]),label='Re $\\rho_{01} mitigated$')
plt.plot(p_values, np.imag(tomo_rhos_mit[0,1,:]), label='Im $\\rho_{01} mitigated$')
plt.plot(p_values, np.real(tomo_rhos_mit[0,0,:]), label='$\\rho_{00} mitigated$')
plt.plot(p_values, np.real(tomo_rhos_mit[1,1,:]), label='$\\rho_{11} mitigated$')

# Predizione teorica

# Otteniamo l'operatore densità dello stato iniziale
rho0 = partial_trace(Statevector(initial_state), ancillae).data

plt.plot(p_values, np.real(rho0[0,1])*(1-p_values), "C0", linewidth=.5)
plt.plot(p_values, np.imag(rho0[0,1])*(1-p_values), "C1", linewidth=.5)
plt.plot(p_values, 0.5*p_values + np.real(rho0[0,0])*(1-p_values), "C2", linewidth=.5)
plt.plot(p_values, 0.5*p_values + np.real(rho0[1,1])*(1-p_values), "C3", linewidth=.5)

plt.xlabel('p')
plt.ylabel('$\\rho_{xx}$')
plt.legend()

plt.title("SIMULATION Depol. channel. Full tomo. Noisy. $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$")

Text(0.5, 1.0, 'SIMULATION Depol. channel. Full tomo. Noisy. $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$')

# Risultati simulati
exact_fidelities = []
for circ in tomography_circuits:
    exact_fidelities.append(state_fidelity(psi, partial_trace(DensityMatrix(circ), ancillae)))

plt.plot(p_values, fidelities, "C0", label="Estimated fidelity")
plt.fill_between(p_values, np.array(fidelities) - np.array(errors), np.array(fidelities) + np.array(errors), alpha=.3)
plt.plot(p_values, fidelities_mit, color="tab:orange", label="Estimated fidelity mitigated")
plt.fill_between(p_values, np.array(fidelities_mit) - np.array(errors_mit), np.array(fidelities_mit) + np.array(errors_mit), color="tab:orange", alpha=.3)
plt.plot(p_values, exact_fidelities, "C1", label="Exact fidelity")

plt.xlabel('p')
plt.ylabel('Fidelity')
plt.legend()

plt.title("SIMULATION Depol. channel. Fidelity. Mitigated. $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$")

Text(0.5, 1.0, 'SIMULATION Depol. channel. Fidelity. Mitigated. $|\\psi_0\\rangle = U_3(\\pi/4,\\pi/4,0)|0\\rangle$')

Sebbene la riduzione dell'errore di lettura sia di qualche aiuto, notiamo che i risultati continuano a non corrispondere alle precisioni previste, motivo per cui devono esserci altre fonti di errore!