Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself

​
Welcome to the Jungle of HHL

  • Rabbits vs.
Foxes The Lanchester model is widely used to describe the dynamics of combat between two opposing forces. Originally formulated for military applications, this model can also be adapted to other contexts, such as ecological competition between two species or market competition between two businesses. The model uses linear differential equations to capture the interactions between the populations, making it suitable for scenarios where interaction terms are purely linear. In this notebook, we explore solving this modified Lanchester model using the Harrow-Hassidim-Lloyd (HHL) algorithm, a quantum algorithm designed for efficiently solving linear systems of equations.

​
Lanchester Model

Lanchester Model is described in this link. This differential equation describes the behavior of two forces, xx and yy: dxdt=ax+bxy+cy+d\frac{d x}{d t} = a x + b x y + cy +d dydt=ey+fxy+gx+h\frac{d y}{d t} = e y + f x y + g x +h where the coefficients described in the next sections indicate the sensitivities of each side to the other and to itself.

​
Describing the Classical Model Using the Finite Difference Method

Assuming b=0b=0, f=0f=0: xi+1=dt⋅axi+dt⋅cyi+dt⋅d+xi=(dt⋅a+1)xi+dt⋅cyi+dt⋅dx_{i+1} = dt \cdot a x_i + dt\cdot c y_i + dt\cdot d + x_i = (dt\cdot a+1) x_i+ dt\cdot c y_i + dt\cdot d yi+1=dt⋅eyi+dt⋅gxi+dt⋅h+yi=(dt⋅e+1)yi+dt⋅gxi+dt⋅h  ,y_{i+1} = dt\cdot e y_i + dt\cdot g x_i + dt\cdot h + y_i = (dt\cdot e+1) y_i + dt \cdot g x_i + dt\cdot h~~, where the dot product was added explicitly to emphasize that dtdt is a variable. The first sample is the initial condition: x0=X0x_0 = X_0 y0=Y0y_0 = Y_0 Then, for any next point, we solve numerically using the previous point: −(dt⋅a+1)xi−dt⋅cyi+xi+1=dt⋅d- (dt\cdot a+1) x_i - dt \cdot c y_i + x_{i+1} = dt\cdot d −(dt⋅e+1)yi−dt⋅gxi+yi+1=dt⋅h- (dt\cdot e+1) y_i - dt\cdot g x_i + y_{i+1} = dt\cdot h The system of linear equations describing the model looks like this:

​
Building the Matrix with NumPy

import numpy as np


def diff_eq_model(a, b, c0, d, e, f, g, h, dt, N, x0, y0):
    assert b == 0, "model is currently unsupported"
    assert f == 0, "model is currently unsupported"

    A = np.identity(2 * N)

    for r in range(2 * N):
        if 1 <= r <= N - 1:
            c = r - 1
            A[r][c] = -(dt * a + 1)
            c = N + r - 1
            A[r][c] = -dt * c0
        elif N + 1 <= r:
            c = r - 1
            A[r][c] = -(dt * e + 1)
            c = r - 1 

- N
            A[r][c] = -dt * g

    b1 = np.ones((N, 1)) * dt * d
    b1[0] = x0
    b2 = np.ones((N, 1)) * dt * h
    b2[0] = y0
    b = np.concatenate([b1, b2])

    return A, b
We choose specific values for the matrix, representing a predator-prey system. Here, (xx) represents the prey population (e.g., rabbits), and (yy) represents the predator population (e.g., foxes). The coefficients have the following meanings and specific values:
  • (aa) and (ee) define the rate of natural losses or birth (due to death, disease, etc.).
For example, (a=−0.01a = -0.01) (1% natural death rate for rabbits), (e=−0.02e = -0.02) (22% natural death rate for foxes).
  • (bb) and (ff) define the rate of losses due to environmental factors (affecting both species).
For example, (b=0b = 0), (f=0f = 0) (assuming no such environmental exposure for simplicity).
  • (cc) and (gg) are losses or gains due to interactions between species (prey hunted by predators and vice versa).
For example, (c=−0.1c = -0.1) (1010% loss of rabbits due to predation), (g=0.2g = 0.2) (2020% increase in predator population due to hunting prey).
  • (dd) and (hh) are gains due to migration.
For example, (d=0.4d = 0.4) (0.4K0.4 K Rabbits/year migration rate for rabbits), (h=0.01h = 0.01) (0.01K0.01K Foxes/year migration rate for foxes). 
N = 8  # N time steps
dt = 6  # Sample every dt years

A, b = diff_eq_model(
    a=-0.01, b=0, c0=-0.1, d=0.4, e=-0.02, f=0, g=0.02, h=0.01, dt=dt, N=N, x0=30, y0=1
)


import matplotlib.pyplot as plt

plt.matshow(A)
plt.title("Matrix A")
plt.show()
plt.matshow(b.transpose())
plt.title("Vector b")
plt.show()
output output

​
Classical Solution

x = np.matmul(np.linalg.inv(A), b)


plt.matshow(x.transpose())
plt.title("Solution x")
plt.show()
output 
import matplotlib.pyplot as plt

t = dt * np.array(range(N))

x_sol = x[0:N]
y_sol = x[N : 2 * N]
plt.plot(t, x_sol, label="Rabbits")
plt.plot(t, y_sol, label="Foxes")
plt.xlabel("t [years]")
plt.ylabel("Population [Thosands of individuals]")
plt.legend()
plt.show()
output Note that at some critical point, the rabbit population is so low that the fox population also starts to decrease. This is typical predator-prey model behavior.

​
Redefining the Matrix

The matrix of HHL is a canonical one, assuming the following properties:
  1. The RHS vector b⃗\vec{b} is normalized.
  2. The matrix AA is a Hermitian one.
  3. The matrix AA is of size 2n×2n2^n\times 2^n .
  4. The eigenvalues of AA are in the range (0,1)(0,1).
However, any general problem that does not follow these conditions can be resolved, as below.

​

​
1) Normalized b

As preprocessing, normalize b⃗\vec{b} and then return the normalization factor as postprocessing. 
norm_factor = np.linalg.norm(b)

b_normalized = b / norm_factor

​

​
2) Hermitian Matrix

Symmetrize the problem: (0ATA0)(x⃗0)=(0b⃗).\begin{pmatrix} 0 & A^T \\ A & 0 \end{pmatrix} \begin{pmatrix} \vec{x} \\ 0 \end{pmatrix} = \begin{pmatrix} 0 \\ \vec{b} \end{pmatrix}. This increases the number of qubits by 
def to_hermitian(A):
    N = A.shape[0]
    A_hermitian = np.concatenate(
        [
            np.concatenate([np.zeros((N, N)), A.transpose().conj()], axis=1),
            np.concatenate([A, np.zeros((N, N))], axis=1),
        ]
    )
    return A_hermitian


b_new = np.concatenate([np.zeros((2 * N, 1)), b_normalized])
plt.matshow(b_new.transpose())
plt.title("Normalized and Padded Vector b")
plt.show()

A_hermitian = to_hermitian(A)
plt.matshow(A_hermitian)
plt.title("Hermitian Matrix A")
plt.show()
output output

​

​
3) Make the Matrix AA of Size 2n×2n2^n\times 2^n

Complete the matrix dimension to the closest 2n2^n with an identity matrix. The vector b⃗\vec{b} is completed with zeros. (A00I)(b⃗0)=(x⃗0).\begin{pmatrix} A & 0 \\ 0 & I \end{pmatrix} \begin{pmatrix} \vec{b} \\ 0 \end{pmatrix} = \begin{pmatrix} \vec{x} \\ 0 \end{pmatrix}. However, our matrix is already the right size.

​

​
4) Rescaled Matrix

If the eigenvalues of AA are in the range [wmin⁡,wmax⁡][w_{\min},w_{\max}], we can employ transformations to the exponentiated matrix and then undo them to extract the results: A~=(A−wmin⁡I)(1−12m)1wmax⁡−wmin⁡.\tilde{A}=(A-w_{\min}I)\left(1-\frac{1}{2^{m}}\right)\frac{1}{w_{\max}-w_{\min}}. The eigenvalues of this matrix lie in the interval [0,1)[0,1), and are related to the eigenvalues of the original matrix via λ=(wmax⁡+wmin⁡)λ~[1/(1−12m)]+wmin⁡,\lambda = (w_{\max}+w_{\min})\tilde{\lambda}\left[1/\left(1-\frac{1}{2^{m}}\right)\right]+w_{\min}, with λ~\tilde{\lambda} being an eigenvalue of A~\tilde{A} resulting from the QPE algorithm. This relation between eigenvalues is then used for the expression inserted into the eigenvalue inversion, via the AmplitudeLoading function. 
def condition_number(A):
    w, _ = np.linalg.eig(A)
    return max(np.abs(w)) / min(np.abs(w))


QPE_RESOLUTION_SIZE = 6

assert QPE_RESOLUTION_SIZE > np.log2(
    condition_number(A)
), "condition number is too big, and QPE resolution cannot hold all eigenvalues"

w, v = np.linalg.eig(A_hermitian)
w_max = np.max(w)
w_min = np.min(w)
mat_shift = -w_min
mat_rescaling = (1 - 1 / 2**QPE_RESOLUTION_SIZE) / (
    w_max - w_min
)  # assures eigenvalues in [0,1-1/2^QPE_SIZE]
min_possible_w = (
    w_max - w_min
) / 2**QPE_RESOLUTION_SIZE  # this is the minimal eigenvalue which can be resolved by the QPE

A_rescaled = (
    A_hermitian + mat_shift * np.identity(A_hermitian.shape[0])
) * mat_rescaling

# verifying that the matrix is symmetric and has eigenvalues in [0,1)
if not np.allclose(A_rescaled, A_rescaled.T, rtol=1e-6, atol=1e-6):
    raise Exception("The matrix is not symmetric")

w_rescaled, _ = np.linalg.eig(A_rescaled)

for lam in w_rescaled:
    if lam < -1e-6 or lam >= 1:
        raise Exception("Eigenvalues are not in (0,1)")

plt.matshow(A_rescaled)
plt.title("Rescaled Matrix A")
plt.show()
output

​
Defining HHL Algorithm for the Quantum Solution

This section is based on Classiq HHL in the user guide, here and here. Note the rescaling in simple_eig_inv based on the matrix rescaling. image.png 
from classiq import *


@qfunc
def simple_eig_inv(
    gamma: float,
    delta: float,
    c_param: float,
    phase: QNum,
    indicator: Output[QBit],
):
    allocate(indicator)
    assign_amplitude_table(
        lookup_table(lambda p: np.clip(c_param / ((gamma * p) + delta), -1, 1), phase),
        phase,
        indicator,
    )


sol_classical_hermitian = np.array([s[0] for s in np.linalg.solve(A_hermitian, b_new)])
compared_sol = sol_classical_hermitian * min_possible_w
amp_compared = compared_sol / np.linalg.norm(compared_sol)
Besides HHL, we also perform swap test as in the user guide, comparing the HHL solution to a state preparation of the known solution. 
import scipy

exponentiation_A_rescaled = scipy.linalg.expm(1j * 2 * np.pi * A_rescaled).tolist()
b_list = np.concatenate(b_new).tolist()
amp_compared_list = amp_compared.tolist()


@qfunc
def main(indicator: Output[QBit], test: Output[QBit]) -> None:

    state = QArray()
    compared_state = QArray()
    rescaled_eig = QNum()
    allocate(QPE_RESOLUTION_SIZE, UNSIGNED, QPE_RESOLUTION_SIZE, rescaled_eig)
    prepare_amplitudes(b_list, 0, state)

    within_apply(
        lambda: qpe(
            unitary=lambda: unitary(exponentiation_A_rescaled, state),
            phase=rescaled_eig,
        ),
        lambda: simple_eig_inv(
            gamma=mat_rescaling ** (-1),
            delta=-mat_shift,
            c_param=min_possible_w,
            phase=rescaled_eig,
            indicator=indicator,
        ),
    )
    prepare_amplitudes(amp_compared_list, 0.0, compared_state)
    swap_test(state, compared_state, test)

    drop(rescaled_eig)
    drop(state)
    drop(compared_state)
We create the circuit with depth limitation compared to the maximal width in a simulator (25). We also set the preferences for synthesis and execution. 
MAX_WIDTH_SWAP_TEST = 25
constraints = Constraints(max_width=MAX_WIDTH_SWAP_TEST)
preferences = Preferences(
    optimization_level=0, optimization_timeout_seconds=90, transpilation_option="none"
)
NUM_SHOTS = 30000
execution_preferences = ExecutionPreferences(num_shots=NUM_SHOTS)

qmod_hhl_swap_test = create_model(main, constraints, execution_preferences, preferences)

​

​
Synthesize and Show

qprog_hhl_swap = synthesize(qmod_hhl_swap_test)
show(qprog_hhl_swap)
Output:

Quantum program link: https://platform.classiq.io/circuit/35tN4TKVCACNL14VMInDJb2alDZ

​

​
Execution and Results Analysis

The analysis is explained in the swap test user guide. We compare the measured overlap with the exact overlap using the expected probability of measuring the state ∣0⟩|0\rangle, defined as α2=12(1+∣⟨ψ1∣ψ2⟩∣2).\alpha^2 = \frac{1}{2}\left(1+|\langle \psi_1 |\psi_2 \rangle |^2\right). We extract the overlap ∣⟨ψ1∣ψ2⟩∣2=2P(qtest=∣0⟩)−1|\langle \psi_1 |\psi_2 \rangle |^2=\sqrt{2 P\left(q_{\text{test}}=|0\rangle\right)-1}. The exact overlap is computed with the dot product of the two state vectors. Note that for the sake of this demonstration we execute this circuit 100,000100,000 times to improve the precision of the probability estimate. This is usually not required. Since we are in the HHL context, we filter only the indicator==1 results. 
execution_job_id = execute(qprog_hhl_swap)


result = execution_job_id.result_value()
fidelity_basic = np.sqrt(
    result.counts_of_multiple_outputs(["indicator", "test"])[("1", "0")]
    * 2
    / (result.counts_of_output("indicator")["1"])
    - 1
)

print("Fidelity between basic HHL and classical solutions:", fidelity_basic)
Output:

Fidelity between basic HHL and classical solutions: 0.9526849928765404
In the IDE we can also check that among the indicator=1 results, the bar of test=0 is much higher than test=1. 
execution_job_id.open_in_ide()

​
State Vector Simulation

We can also run the state vector of the HHL result and examine it without a swap test. Extracting the full information is exponentially hard, and used here just for educational purposes. Extension of this work can extract some information from the state vector, such as the last element, which is the most important. We can also pad the last solution with an extension, and therefore measure the last point with high probability. 
class TimeIndexAndGroup(QStruct):
    rabbits: QBit
    time_index: QNum


@qfunc
def main(
    indicator: Output[QBit],
    time_index_and_group: Output[TimeIndexAndGroup],
    rescaled_eig: Output[QNum],
) -> None:
    allocate(QPE_RESOLUTION_SIZE, False, QPE_RESOLUTION_SIZE, rescaled_eig)
    prepare_amplitudes(b_list, 0, time_index_and_group)
    within_apply(
        lambda: qpe(
            unitary=lambda: unitary(exponentiation_A_rescaled, time_index_and_group),
            phase=rescaled_eig,
        ),
        lambda: simple_eig_inv(
            gamma=mat_rescaling ** (-1),
            delta=-mat_shift,
            c_param=min_possible_w,
            phase=rescaled_eig,
            indicator=indicator,
        ),
    )


MAX_WIDTH_BASIC = 18
constraints = Constraints(max_width=MAX_WIDTH_BASIC)
preferences = Preferences(
    optimization_level=0, optimization_timeout_seconds=90, transpilation_option="none"
)
backend_preferences = ClassiqBackendPreferences(
    backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR
)
execution_preferences = ExecutionPreferences(
    num_shots=1, backend_preferences=backend_preferences
)

qmod_hhl_basic = create_model(
    main,
    constraints=constraints,
    preferences=preferences,
    execution_preferences=execution_preferences,
)
qprog_hhl_basic = synthesize(qmod_hhl_basic)
show(qprog_hhl_basic)
Output:

Quantum program link: https://platform.classiq.io/circuit/35tNJgp5JBnQjy397u9i5ZqNFvs
  


execution_job_id = execute(qprog_hhl_basic)
result = execution_job_id.result_value()
# statevector = result.state_vector
execution_job_id.open_in_ide()
We filter the results, then extract the amplitudes of the state vector and compare them to the classical solution. We filter only the indicator=1 and phase=0 results. 
filtered_hhl_statevector = dict()
for sample in result.parsed_state_vector:
    if sample.state["indicator"] == 1 and sample.state["rescaled_eig"] == 0:
        filtered_hhl_statevector[sample.state["time_index_and_group"]] = (
            sample.amplitude
        )

states = sorted(filtered_hhl_statevector)
raw_qsol = np.array([filtered_hhl_statevector[s] for s in states])
Let’s compare the raw solution to the Hermitian matrix with the classical one. 
qsol_hermitian = raw_qsol / (min_possible_w)

plt.plot(np.real(states), np.real(qsol_hermitian))
plt.plot(np.real(states), np.real(sol_classical_hermitian))
plt.xlabel("states")
plt.ylabel("amplitude")
plt.legend(["hhl result", "original result"])
plt.show()
output Now we can compare the two solutions in the time domain, after putting back the normalization factor. 
plt.plot(t, norm_factor * np.real(qsol_hermitian[0:N]), label="Rabbits HHL")
plt.plot(
    t, norm_factor * np.real(sol_classical_hermitian[0:N]), label="Rabbits Classical"
)

plt.plot(t, norm_factor * np.real(qsol_hermitian[N : 2 * N]), label="Foxes HHL")
plt.plot(
    t,
    norm_factor * np.real(sol_classical_hermitian[N : 2 * N]),
    label="Foxes Classical",
)
plt.xlabel("time [years]")
plt.ylabel("Population [Thousands of individuals]")
plt.legend()
plt.show()
output The fidelity of the state vector solution can be calculated as the overlap of the two solutions. 
fidelity = (
    np.abs(
        np.dot(
            sol_classical_hermitian / np.linalg.norm(sol_classical_hermitian),
            qsol_hermitian / np.linalg.norm(qsol_hermitian),
        )
    )
    ** 2
)
print("Statevector Solution Fidelity:", fidelity)
Output:

Statevector Solution Fidelity: 0.9996462667876236