Quantum Simulation of Linear Kinetic Plasma Models (Qiskit)

View on GitHub

Open this notebook in GitHub to run it yourself

**This notebook implements the same model that appears under vlasov_ampere.ipynb example, using qiskit. This is the code used to collect the benchmarking data in the paper**. In particular, it contains: Definition of a qiskit code for the block-encoding, verification of the block-encoding functionality, and definition of QSVT step with qiskit, operating on the block-encoding unitary. The results were obtained using Qiskit version 2.1.

!pip install -qq  qiskit==2.1.1

import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, transpile
from qiskit.circuit.library import StatePreparation
from qiskit.circuit.library.standard_gates import XGate
from qiskit.synthesis import synth_qft_full
from sympy import fwht

Adding inplace adder and assign_amplitudes with Qiskit

These are missing in Qiskit library.

# Creating an inplace adder


class DraperQFTAdderConstant(QuantumCircuit):
    def __init__(
        self, num_state_qubits: int, constant: int, name: str = "DraperQFTAdderConst"
    ) -> None:
        # Create the quantum register
        qr_a = QuantumRegister(num_state_qubits, name="a")
        super().__init__(qr_a, name=name)

        # Apply the QFT
        self.append(synth_qft_full(num_state_qubits, do_swaps=False).to_gate(), qr_a)

        # Add the constant by applying controlled rotations
        for qubit in range(num_state_qubits):
            angle = (constant % (2 ** (qubit + 1))) * np.pi / (2**qubit)
            self.p(angle, qr_a[qubit])

        # Apply the inverse QFT
        self.append(
            synth_qft_full(num_state_qubits, do_swaps=False).inverse().to_gate(), qr_a
        )

# Creating assign amplitudes


def get_graycode(size: int, i: int) -> int:
    if i == 2**size:
        return get_graycode(size, 0)
    return i ^ (i >> 1)


def get_controller(size: int, i: int) -> int:
    return (get_graycode(size, i) ^ get_graycode(size, i + 1)).bit_length() - 1


def get_graycode_angles_wh(size, angles):
    transformed_angles = fwht(np.array(angles) / 2**size)
    return [transformed_angles[get_graycode(size, j)] for j in range(2**size)]


def assign_amplitudes_qs(amps, x, ind):
    qc = QuantumCircuit(x, ind, name="assign amplitudes")
    size = len(amps).bit_length() - 1
    angles = 2 * np.arcsin(amps)
    gray_code_angles = np.array(get_graycode_angles_wh(size, angles)).astype(float)
    gray_code_controllers = [get_controller(size, i) for i in range(2**size)]

    for i in range(2**size):
        qc.ry(gray_code_angles[i], ind)
        qc.cx(x[gray_code_controllers[i]], ind)

    return qc.to_gate()

Setting problem parameters

# plasma parameters
Temperature = 1
N = 1  # simulate homogenous plasma - n(x) = 1

# source parameters
DS = 3
X_0 = 50
W_0 = 0.8

# grid  parameters
X_MIN = 0
X_MAX = 100

V_MIN = -4
V_MAX = 4

N_X = 3
N_V = 3

DX = (X_MAX 

- X_MIN) / (2**N_X - 1)
DV = (V_MAX 

- V_MIN) / (2**N_V - 1)

Functions for getting quantum circuits (“gates”) to construct the block encoding

Circuits for the Advective part

BC_VALUES = 0.5 * np.array([-3, 3, -1, 0])
BC_NORM = np.linalg.norm(BC_VALUES)
BC_VALUES = BC_VALUES / BC_NORM
lcu_angle = 2 * np.arccos(np.sqrt(BC_NORM / (1 + BC_NORM)))

def get_derivative_dirichlet_be(x, flag, lcu):
    qc = QuantumCircuit(x, flag, lcu, name="derivative_dirichlet_be")
    qc.h(lcu[0])
    qc.z(lcu[0])
    qc.append(
        DraperQFTAdderConstant(num_state_qubits=len(x) + 1, constant=-1)
        .to_gate()
        .control(1, ctrl_state=0),
        [lcu] + x[:] + [flag],
    )
    qc.append(
        DraperQFTAdderConstant(num_state_qubits=len(x) + 1, constant=1)
        .to_gate()
        .control(1, ctrl_state=1),
        [lcu] + x[:] + [flag],
    )
    qc.h(lcu[0])

    return qc.to_gate()


def get_derivative_boundary_min_be(x, flag):
    qc = QuantumCircuit(x, flag, name="derivative_boundary_min_be")
    qc.append(StatePreparation(BC_VALUES).inverse(), x[0:2])
    qc.append(XGate().control(len(x), ctrl_state=0), x[:] + [flag])
    qc.x(flag)

    return qc.to_gate()


def get_derivative_boundary_max_be(x, flag):
    qc = QuantumCircuit(x, flag, name="derivative_boundary_max_be")
    qc.x(x)
    qc.ry(2 * np.pi, flag)
    qc.append(get_derivative_boundary_min_be(x, flag), x[:] + [flag])
    qc.x(x)

    return qc.to_gate()


def get_derivative_boundaries_be(x, flag):
    qc_start = get_derivative_boundary_min_be(x[0:-1], flag)
    qc_end = get_derivative_boundary_max_be(x[0:-1], flag)
    qc = QuantumCircuit(x, flag, name="derivative_boundaries_be")
    qc.append(qc_start.control(1, ctrl_state=0), [x[-1]] + x[0:-1] + [flag])
    qc.append(qc_end.control(1, ctrl_state=1), [x[-1]] + x[0:-1] + [flag])

    return qc.to_gate()


def get_derivative_be(x, flag, lcu1, lcu2):
    qc = QuantumCircuit(x, flag, lcu1, lcu2, name="derivative_be")
    qc.ry(lcu_angle, lcu1)
    qc.append(
        get_derivative_dirichlet_be(x, flag, lcu2).control(1, ctrl_state=1),
        [lcu1] + x[:] + [flag] + [lcu2],
    )
    qc.append(
        get_derivative_boundaries_be(x, flag).control(1, ctrl_state=0),
        [lcu1] + x[:] + [flag],
    )
    qc.ry(-lcu_angle, lcu1)

    return qc.to_gate()


def get_be_zeta(x, v, flag):
    temp1 = QuantumRegister(1, "temp1")
    temp2 = QuantumRegister(1, "temp2")

    qc = QuantumCircuit(x, v, flag, temp1, temp2, name="zeta_bc matrix")
    qc.append(XGate().control(len(x), ctrl_state=2 ** len(x) - 1), x[:] + [temp1])
    qc.cx(v[-1], temp2)
    qc.ccx(temp1, temp2, flag)
    qc.cx(v[-1], temp2)
    qc.append(XGate().control(len(x), ctrl_state=2 ** len(x) - 1), x[:] + [temp1])
    qc.append(XGate().control(len(x), ctrl_state=0), x[:] + [temp1])
    qc.append(XGate().control(1, ctrl_state=0), [v[-1]] + [temp2])
    qc.ccx(temp1, temp2, flag)
    qc.append(XGate().control(1, ctrl_state=0), [v[-1]] + [temp2])
    qc.append(XGate().control(len(x), ctrl_state=0), x[:] + [temp1])

    return qc.to_gate()


def get_v_be(v, ind):
    qc = QuantumCircuit(v, ind, name="v dot matrix")
    N_V = len(v)
    amplitudes = np.linspace(-1, 1 - 2 ** (-N_V + 1), 2**N_V)
    amplitudes = np.roll(amplitudes, len(amplitudes) // 2)  # adjust to signed number
    qc.append(assign_amplitudes_qs(amplitudes, v, ind), v[:] + [ind])
    qc.x(ind)

    return qc.to_gate()


# The upper left matrix: 6 block terms, 2 aux
def get_advective_be(x, v, e, flags):
    temps = QuantumRegister(2, "temps")
    qc = QuantumCircuit(x, v, e, flags, temps, name="Full v grad_x")
    qc.append(
        get_derivative_be(x, [flags[0]], [flags[1]], [flags[2]]), x[:] + flags[0:3]
    )
    qc.append(get_v_be(v, [flags[3]]), v[:] + [flags[3]])
    qc.append(get_be_zeta(x, v, [flags[4]]), x[:] + v[:] + [flags[4]] + temps[:])
    qc.cx(e, flags[5])

    return qc.to_gate()

Circuits for off-diagonal terms

v_amplitudes = np.linspace(-1, 1 - 2 ** (-N_V + 1), 2**N_V) * V_MAX
v_amplitudes = np.roll(v_amplitudes, len(v_amplitudes) // 2)  # adjust to signed number
v_H_amplitudes = (
    v_amplitudes
    * np.exp(-(v_amplitudes**2) / (2 * Temperature))
    / (np.sqrt(2 * np.pi * Temperature))
)


def get_force_term_be(v, flag, v_H_amplitudes):
    qc = QuantumCircuit(v, flag, name="force_term")
    qc.append(XGate().control(len(v), ctrl_state=0), v[:] + [flag])
    qc.x(flag)
    qc.append(StatePreparation(-v_H_amplitudes / np.linalg.norm(v_H_amplitudes)), v[:])

    return qc.to_gate()


def get_current_term_be(v, flag, v_amplitudes):
    qc = QuantumCircuit(v, flag, name="current_term")
    qc.append(
        StatePreparation(v_amplitudes / np.linalg.norm(v_amplitudes)).inverse(), v[:]
    )
    qc.append(XGate().control(len(v), ctrl_state=0), v[:] + [flag])
    qc.x(flag)

    return qc.to_gate()


def get_equalize_amplitude(e, ind, ratio):
    qc = QuantumCircuit(e, ind, name="equalize_amplitude")
    amplitudes = np.array([ratio, 1])
    amplitudes /= max(amplitudes)
    qc.append(assign_amplitudes_qs(amplitudes, e, ind), e[:] + [ind])
    qc.x(ind)

    return qc.to_gate()


def get_off_diag_be(e, v, flag, ind, v_H_amplitudes, v_amplitudes, ratio):
    qc = QuantumCircuit(e, v, flag, ind, name="off diagonal")
    qc.x(e)
    qc.append(
        get_force_term_be(v, flag, v_H_amplitudes).control(1, ctrl_state=0),
        [e] + v[:] + [flag],
    )
    qc.append(
        get_current_term_be(v, flag, v_amplitudes).control(1, ctrl_state=1),
        [e] + v[:] + [flag],
    )
    qc.append(get_equalize_amplitude(e, ind, ratio), [e] + [ind])

    return qc.to_gate()

Block-encoding of the full matrix

DERIVATIVE_TERM_FACTOR = 2 * (1 + BC_NORM)
ADVECTIVE_TERM_FACTOR = DERIVATIVE_TERM_FACTOR * (V_MAX / (2 * DX))

FORCE_TERM_FACTOR = np.linalg.norm(v_H_amplitudes)
CURRENT_TERM_FACTOR = np.linalg.norm(v_amplitudes) * DV

OFF_DIAG_FACTOR = max(FORCE_TERM_FACTOR, CURRENT_TERM_FACTOR)

# Defining a function for 1j factor


def get_phase():
    q = QuantumRegister(1, "q")
    qc = QuantumCircuit(q, name="1j phase")
    qc.p(np.pi, q)
    qc.rz(-np.pi, q)

    return qc.to_gate()

data_size = N_V + N_X + 1
data_qs = QuantumRegister(data_size, "data_qs")
block_qs = QuantumRegister(2 + 6 + 2, "block_qs")

lcu_amps = np.sqrt(
    np.array([ADVECTIVE_TERM_FACTOR, OFF_DIAG_FACTOR, W_0, 0]).astype("float")
)
lcu_amps /= np.linalg.norm(lcu_amps)


qc = QuantumCircuit(data_qs, block_qs)

qc.append(StatePreparation(lcu_amps), block_qs[0:2])
qc.append(
    get_advective_be(
        data_qs[N_V : data_size - 1],
        data_qs[0:N_V],
        [data_qs[data_size - 1]],
        block_qs[0 + 2 : 6 + 2],
    ).control(2, ctrl_state=0),
    block_qs[0:2]
    + data_qs[N_V : data_size - 1]
    + data_qs[0:N_V]
    + [data_qs[data_size - 1]]
    + block_qs[2:],
)
qc.append(
    get_off_diag_be(
        [data_qs[data_size - 1]],
        data_qs[0:N_V],
        [block_qs[0 + 2]],
        [block_qs[3 + 2]],
        v_H_amplitudes,
        v_amplitudes,
        FORCE_TERM_FACTOR / CURRENT_TERM_FACTOR,
    ).control(2, ctrl_state=1),
    block_qs[0:2]
    + [data_qs[data_size - 1]]
    + data_qs[0:N_V]
    + [block_qs[0 + 2]]
    + [block_qs[3 + 2]],
)
qc.append(get_phase().control(2, ctrl_state=2), block_qs[0:2] + [block_qs[2]])
inv_amps = lcu_amps
inv_amps[0] *= -1  # minus sign for the advective term
qc.append(StatePreparation(inv_amps).inverse(), block_qs[0:2])

be_qc = qc.to_gate()

# Scaling factor
BE_NORM_FACTOR = ADVECTIVE_TERM_FACTOR + OFF_DIAG_FACTOR + W_0

Verification

For the given problem size, N_X, N_V= 3, 3, we construct the classical matrix

A+i\omega

that we block-encoded via a quantum circuit

U_{A+i\omega}

. Then for some random initial state

|\psi\rangle

we verify that indeed:

U_{A+i\omega} |\psi\rangle |0\rangle_{\rm block} = (A+i\omega)|\psi\rangle|0\rangle_{\rm block}+{\rm garbage}

Classical matrix

import numpy as np


def get_advective_mat(nx, nv, dx, v_max, cyclic=False):
    dx_mat = np.diag(np.ones(2**nx - 1), k=1) - np.diag(np.ones(2**nx - 1), k=-1)
    boundary = np.pad([-3, 4, -1], (0, dx_mat.shape[1] - 3))
    dx_mat[0, :] = boundary
    dx_mat[-1, :] = -np.flip(boundary)

    v_amplitudes = np.linspace(-1, 1 - 2 ** (-nv + 1), 2**nv) * v_max
    v_amplitudes = np.roll(v_amplitudes, len(v_amplitudes) // 2)

    x_values = np.arange(2**nx)
    max_x = 2**nx - 1
    xi = (
        1
        - np.kron(x_values == 0, v_amplitudes > 0)
        - np.kron(x_values == max_x, v_amplitudes <= 0)
    )

    advective = np.kron(dx_mat / (2 * dx), np.diag(v_amplitudes))
    advective[xi == 0] = 0

    return advective


def get_off_diag_mat(nx, nv, dv, v_max, temp):
    v_amplitudes = np.linspace(-1, 1 - 2 ** (-nv + 1), 2**nv) * v_max
    v_amplitudes = np.roll(v_amplitudes, len(v_amplitudes) // 2)
    v_H_amplitudes = (
        v_amplitudes
        * np.exp(-(v_amplitudes**2) / (2 * temp))
        / (np.sqrt(2 * np.pi * temp))
    )

    vec_col = v_H_amplitudes
    vec_row = v_amplitudes * dv

    n = len(vec_col)
    dim_x = 2**nx
    dim_v = 2**nv
    size = dim_x * dim_v

    A = np.zeros((size, size))
    B = np.zeros((size, size))
    for i in range(dim_x):
        A[i * dim_v : (i + 1) * dim_v, i * dim_v] = vec_col
        B[i * dim_v, i * dim_v : (i + 1) * dim_v] = vec_row

    # Compose off-diagonal block matrix
    top = np.hstack((np.zeros_like(A), -A))
    bottom = np.hstack((B, np.zeros_like(B)))
    off_diag = np.vstack((top, bottom))
    return off_diag


def get_block_encoding(nx, nv, v_max, x_max, w0, temp):
    dx = x_max / (2**nx - 1)
    dv = (2 * v_max) / (2**nv - 1)

    mat_size = 2 ** (nx + nv + 1)
    advective = get_advective_mat(nx, nv, dx, v_max)
    off_diag = get_off_diag_mat(nx, nv, dv, v_max, temp)

    return 1j * w0 * np.eye(mat_size) - np.kron(np.diag([1, 0]), advective) + off_diag


mat_cl = get_block_encoding(N_X, N_V, V_MAX, X_MAX, W_0, Temperature)

# Random initial state
init_data = np.random.rand(2**7)
init_data /= np.linalg.norm(init_data)

block_qs = QuantumRegister(10, "block_qs")
data_qs = QuantumRegister(data_size, "data_qs")

qc_test = QuantumCircuit(data_qs, block_qs)
qc_test.append(StatePreparation(init_data), data_qs[:])
qc_test.append(be_qc, data_qs[:] + block_qs[:])

tqc_test = transpile(qc_test, basis_gates=["u", "cx"], optimization_level=3)
print(f"width: {tqc_test.width()}")
print(f"depth: {tqc_test.depth()}")
print(f"cx: {tqc_test.count_ops()['cx']}")

Output:

width: 17
  depth: 237580
  cx: 124273
  

from qiskit.quantum_info import Statevector

final_sv = Statevector.from_instruction(qc_test)
psi_tensor = final_sv.data.reshape([2] * qc_test.width())

projected = psi_tensor[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, :]
reduced_sv = Statevector(projected.flatten())


cl_vector = mat_cl @ init_data

We compare the projected state on

|0\rangle_{\rm block}

with the expected classical result (in absolute value)

import matplotlib.pyplot as plt

plt.plot(np.abs(cl_vector), "o", label="expected state")
plt.plot(np.abs(reduced_sv) * BE_NORM_FACTOR, ".", label="projected quantum state")
plt.legend();

We also compute the norm for the difference between the results

l2_norm = np.linalg.norm(cl_vector - np.array(reduced_sv) * BE_NORM_FACTOR)
print(f"The L2-norm: {l2_norm}")
assert l2_norm < 1e-10

Output:


The L2-norm: 1.5013973106103977e-11

QSVT step

def get_reflect_around_zero(size):
    qc = QuantumCircuit(size)
    qc.x(0)
    qc.h(0)
    qc.mcx(
        control_qubits=[k for k in range(1, size)],
        ctrl_state="0" * (size - 1),
        target_qubit=[0],
    )
    qc.h(0)
    qc.x(0)

    return qc


def apply_projector_controlled_phase(qc, phase, block_reg, aux_reg):
    qc.append(XGate().control(len(block_reg), ctrl_state=0), block_reg[:] + aux_reg[:])
    qc.rz(phase, aux_reg)

    qc.append(XGate().control(len(block_reg), ctrl_state=0), block_reg[:] + aux_reg[:])


def apply_qsvt_step(qc, phase1, phase2, u, data, block, aux, qsvt_aux):

    qc.append(u, data[:] + block[:] + aux[:])
    apply_projector_controlled_phase(qc, phase1, block, qsvt_aux)
    qc.append(u.inverse(), data[:] + block[:] + aux[:])
    apply_projector_controlled_phase(qc, phase2, block, qsvt_aux)

def get_qsvt_circuit():

    data = QuantumRegister(N_X + N_V + 1, "data")
    block = QuantumRegister(8, "block")
    aux = QuantumRegister(2, "aux")
    qsvt_aux = QuantumRegister(1, "qsvt_aux")

    qsvt_cir = QuantumCircuit(data, block, qsvt_aux, aux)

    apply_qsvt_step(
        qsvt_cir,
        0.1,  # dummy angles
        0.2,  # dummy angles
        be_qc,
        data,
        block,
        aux,
        qsvt_aux,
    )

    return qsvt_cir


qc_qsvt = get_qsvt_circuit()

tqc_qsvt = transpile(qc_qsvt, basis_gates=["u", "cx"], optimization_level=3)
print(f"width: {tqc_qsvt.width()}")
print(f"depth: {tqc_qsvt.depth()}")
print(f"cx: {tqc_qsvt.count_ops()['cx']}")

Output:

width: 18
  depth: 475033
  cx: 248500
  

Algorithms

Applications

Tutorials

Quantum Simulation of Linear Kinetic Plasma Models (Qiskit)

View on GitHub

Adding inplace adder and assign_amplitudes with Qiskit

Setting problem parameters

Functions for getting quantum circuits (“gates”) to construct the block encoding

Circuits for the Advective part

Circuits for off-diagonal terms

Block-encoding of the full matrix

Verification

Classical matrix

QSVT step

View on GitHub

​Adding inplace adder and assign_amplitudes with Qiskit

​Setting problem parameters

​Functions for getting quantum circuits (“gates”) to construct the block encoding

​

​Circuits for the Advective part

​

​Circuits for off-diagonal terms

​Block-encoding of the full matrix

​Verification

​

​Classical matrix

​QSVT step

Adding inplace adder and assign_amplitudes with Qiskit

Setting problem parameters

Functions for getting quantum circuits (“gates”) to construct the block encoding

Circuits for the Advective part

Circuits for off-diagonal terms

Block-encoding of the full matrix

Verification

Classical matrix

QSVT step