Hello! If applicable, put your complete code example down below. Make sure that your code:
- is 100% self-contained — someone can copy-paste exactly what is here and run it to
reproduce the behaviour you are observing - includes comments
# Put code here
import torch
torch.cuda.is_available()
!pip install pennylane
# Pennylane
import pennylane as qml
from pennylane import numpy as np
# Plotting
import matplotlib.pyplot as plt
"""**Setting of the main hyper-parameters of the model**"""
# Set random seed for reproducibility
np.random.seed(0)
"""## Hyperparameters"""
n_qubits = 8 # Number of system qubits
n_shots = 10 ** 6 # Number of quantum measurements
tot_qubits = n_qubits + 1 # Addition of an ancillary qubit
ancilla_idx = n_qubits # Index of the ancillary qubit (last position)
steps = 100 # Number of optimization steps
eta = 9e-3 # Adjusted learning rate (originally 0.8)
q_delta = 1e-4 # Initial spread of random quantum weights
n_layers = 8 # Number of layers in the variational circuit
"""**Circuits of the quantum linear problem**
We now define the unitary operations associated to the simple example presented in the introduction. Since we want to implement a Hadamard test, we need the unitary operations $A_{j}$ to be controlled by the state of an ancillary qubit.
"""
n = 256
A = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1) - np.eye(n, k=2) - np.eye(n, k=-2)
LCU = qml.pauli_decompose(A)
LCU_coeffs, LCU_ops = LCU.terms()
# print(f"LCU decomposition:\n {LCU}")
# print(f"Coefficients:\n {LCU_coeffs}")
# print(f"Unitaries:\n {LCU_ops}")
c = np.array(LCU_coeffs)
# Enhanced U_b preparation
def U_b():
"""Proper state preparation for |b> = [1,1,1,1,1,1,1,1]/sqrt(8)"""
for idx in range(n_qubits):
qml.Hadamard(wires=idx)
# Enhanced CA implementation
# def CA(idx):
# """Controlled versions of the unitary components"""
# if idx == 0:
# None # Identity operation
# elif idx == 1:
# qml.CNOT(wires=[ancilla_idx, 2])
# elif idx == 2:
# qml.CNOT(wires=[ancilla_idx, 1])
# elif idx == 3:
# qml.CNOT(wires=[ancilla_idx, 1])
# qml.CNOT(wires=[ancilla_idx, 2])
# elif idx == 4:
# qml.CY(wires=[ancilla_idx, 1])
# qml.CY(wires=[ancilla_idx, 2])
# elif idx == 5:
# qml.CNOT(wires=[ancilla_idx, 0])
# qml.CNOT(wires=[ancilla_idx, 1])
# elif idx == 6:
# qml.CNOT(wires=[ancilla_idx, 0])
# qml.CNOT(wires=[ancilla_idx, 1])
# qml.CNOT(wires=[ancilla_idx, 2])
# elif idx == 7:
# qml.CNOT(wires=[ancilla_idx, 0])
# qml.CY(wires=[ancilla_idx, 1])
# qml.CY(wires=[ancilla_idx, 2])
# elif idx == 8:
# qml.CY(wires=[ancilla_idx, 0])
# qml.CNOT(wires=[ancilla_idx, 1])
# qml.CY(wires=[ancilla_idx, 2])
# elif idx == 9:
# qml.CY(wires=[ancilla_idx, 0])
# qml.CY(wires=[ancilla_idx, 1])
# elif idx == 10:
# qml.CY(wires=[ancilla_idx, 0])
# qml.CY(wires=[ancilla_idx, 1])
# qml.CNOT(wires=[ancilla_idx, 2])
def CA(idx):
"""Automatically parse and build controlled unitaries from Pauli strings"""
if idx >= len(LCU_ops):
raise ValueError(f"Index {idx} out of range for LCU operations")
pauli_term = LCU_ops[idx]
pauli_str = str(pauli_term)
# Parse the Pauli string
if '⊗' in pauli_str:
terms = [term.strip() for term in pauli_str.split('⊗')]
elif '@' in pauli_str:
terms = [term.strip() for term in pauli_str.split('@')]
else:
terms = [pauli_str.strip()]
for term in terms:
# Extract Pauli operator and qubit index
if '(' in term and ')' in term:
pauli_name = term.split('(')[0].strip()
qubit_idx = int(term.split('(')[1].split(')')[0])
else:
# Handle formats like "X0" or "I"
pauli_name = term[0] if term[0] in ['I', 'X', 'Y', 'Z'] else term
qubit_idx = int(term[1:]) if len(term) > 1 else 0
# Apply controlled Pauli gate
if pauli_name == 'X':
qml.CNOT(wires=[ancilla_idx, qubit_idx])
elif pauli_name == 'Y':
qml.CY(wires=[ancilla_idx, qubit_idx])
elif pauli_name == 'Z':
qml.CZ(wires=[ancilla_idx, qubit_idx])
# Identity (I) requires no operation
"""**Variational quantum circuit**"""
# # Enhanced variational circuit
# def variational_block(weights):
# """More powerful ansatz"""
# # Layer 1: Hadamards + Entanglement
# for idx in range(n_qubits):
# qml.Hadamard(wires=idx)
# for i in range(n_qubits-1):
# qml.CNOT(wires=[i, i+1])
# # Layer 2: Parameterized rotations
# for idx, w in enumerate(weights[:n_qubits]):
# qml.RY(w, wires=idx)
# for i in range(n_qubits-1):
# qml.CNOT(wires=[i, i+1])
# # Layer 3: Final rotations
# for idx, w in enumerate(weights[n_qubits:2*n_qubits]):
# qml.RZ(w, wires=idx)
def variational_block(weights):
"""
Enhanced variational circuit with StronglyEntanglingLayers.
Args:
weights: Array of shape (n_layers, n_qubits, 3) for rotation parameters
"""
# Initial Hadamard layer for superposition
for idx in range(n_qubits):
qml.Hadamard(wires=idx)
# Strongly entangling layers for variational optimization
qml.StronglyEntanglingLayers(weights, wires=range(n_qubits))
# def variational_block(weights):
# # Apply RY rotations with the first set of parameters
# for i in range(n_qubits):
# qml.RY(weights[i], wires=i)
# # Apply CNOT gates with adjacent qubits (cyclically connected) to create entanglement
# for i in range(n_qubits):
# qml.CNOT(wires=[(i - 1) % n_qubits, (i) % n_qubits])
"""**Hadamard test**"""
dev_mu = qml.device("lightning.qubit", wires=tot_qubits)
@qml.qnode(dev_mu, interface="autograd")
def local_hadamard_test(weights, l=None, lp=None, j=None, part=None):
# First Hadamard gate applied to the ancillary qubit.
qml.Hadamard(wires=ancilla_idx)
# For estimating the imaginary part of the coefficient "mu", we must add a "-i"
# phase gate.
if part == "Im" or part == "im":
qml.PhaseShift(-np.pi / 2, wires=ancilla_idx)
# Variational circuit generating a guess for the solution vector |x>
variational_block(weights)
# Controlled application of the unitary component A_l of the problem matrix A.
CA(l)
# Adjoint of the unitary U_b associated to the problem vector |b>.
# In this specific example Adjoint(U_b) = U_b.
U_b()
# Controlled Z operator at position j. If j = -1, apply the identity.
if j != -1:
qml.CZ(wires=[ancilla_idx, j])
# Unitary U_b associated to the problem vector |b>.
U_b()
# Controlled application of Adjoint(A_lp).
# In this specific example Adjoint(A_lp) = A_lp.
CA(lp)
# Second Hadamard gate applied to the ancillary qubit.
qml.Hadamard(wires=ancilla_idx)
# Expectation value of Z for the ancillary qubit.
return qml.expval(qml.PauliZ(wires=ancilla_idx))
def mu(weights, l=None, lp=None, j=None):
"""Generates the coefficients to compute the "local" cost function C_L."""
mu_real = local_hadamard_test(weights, l=l, lp=lp, j=j, part="Re")
mu_imag = local_hadamard_test(weights, l=l, lp=lp, j=j, part="Im")
return mu_real + 1.0j * mu_imag
"""**Local cost function**"""
def psi_norm(weights):
"""Returns the normalization constant <psi|psi>, where |psi> = A |x>."""
norm = 0.0
for l in range(0, len(c)):
for lp in range(0, len(c)):
norm = norm + c[l] * np.conj(c[lp]) * mu(weights, l, lp, -1)
return abs(norm)
def cost_loc(weights):
"""Local version of the cost function. Tends to zero when A|x> is proportional to |b>."""
mu_sum = 0.0
for l in range(0, len(c)):
for lp in range(0, len(c)):
for j in range(0, n_qubits):
mu_sum = mu_sum + c[l] * np.conj(c[lp]) * mu(weights, l, lp, j)
mu_sum = abs(mu_sum)
# Cost function C_L
return 0.5 - 0.5 * mu_sum / (n_qubits * psi_norm(weights))
# def cost_loc(weights):
# mu_sum = 0.0
# for l in range(len(c)):
# for lp in range(len(c)):
# mu_sum += c[l] * np.conj(c[lp]) * mu(weights, l, lp, -1)
# return 1 - abs(mu_sum)/(psi_norm(weights))
"""**Variational optimization**"""
# np.random.seed(rng_seed)
w = q_delta * np.random.randn(n_layers, n_qubits, 3, requires_grad=True)
# w = q_delta * np.random.randn(n_qubits, requires_grad=True)
# Optimization
# opt = qml.GradientDescentOptimizer(eta)
opt = qml.AdamOptimizer(eta)
cost_history = []
for it in range(steps):
w, cost = opt.step_and_cost(cost_loc, w)
cost_history.append(cost)
if it % 10 == 0:
print(f"Step {it:3d} Cost: {cost:.8f}")
import seaborn as sns
import matplotlib.pyplot as plt
plt.style.use("seaborn-v0_8")
plt.plot(cost_history, "g")
plt.ylabel("Cost function")
plt.xlabel("Optimization steps")
plt.show()
"""**Comparison of quantum and classical results**"""
# Id = np.identity(2)
# Z = np.array([[1, 0], [0, -1]])
# X = np.array([[0, 1], [1, 0]])
# A_0 = np.identity(8)
# A_1 = np.kron(np.kron(X, Z), Id)
# A_2 = np.kron(np.kron(X, Id), Id)
# A_num = c[0] * A_0 + c[1] * A_1 + c[2] * A_2
# b = np.ones(8) / np.sqrt(8)
# Id = np.identity(2)
# Z = np.array([[1, 0], [0, -1]])
# X = np.array([[0, 1], [1, 0]])
# Y = np.array([[0, -1j], [1j, 0]])
# A_0 = np.identity(8)
# A_1 = np.kron(np.kron(Id, Id), X)
# A_2 = np.kron(np.kron(Id, X), Id)
# A_3 = np.kron(Id, np.kron(X, X))
# A_4 = np.kron(Id, np.kron(Y, Y))
# A_5 = np.kron(np.kron(X, X), Id)
# A_6 = np.kron(np.kron(X, X), X)
# A_7 = np.kron(np.kron(X, Y), Y)
# A_8 = np.kron(np.kron(Y, X), Y)
# A_9 = np.kron(np.kron(Y, Y), Id)
# A_10 = np.kron(np.kron(Y, Y), X)
# A_num = c[0] * A_0 + c[1] * A_1 + c[2] * A_2 + c[3] * A_3 + c[4] * A_4 + c[5] * A_5 + c[6] * A_6 + c[7] * A_7 + c[8] * A_8 + c[9] * A_9 + c[10] * A_10
A_num = qml.matrix(LCU)
b = np.ones(n) / np.sqrt(n)
# b = np.array(b_normalized)
print("A = \n", A_num)
print("b = \n", b)
A_inv = np.linalg.inv(A_num)
x = np.dot(A_inv, b)
c_probs = (x / np.linalg.norm(x)) ** 2
"""**Preparation of the quantum solution**"""
dev_x = qml.device("lightning.qubit", wires=n_qubits, shots=n_shots)
@qml.qnode(dev_x, interface="autograd")
def prepare_and_sample(weights):
# Variational circuit generating a guess for the solution vector |x>
variational_block(weights)
# We assume that the system is measured in the computational basis.
# then sampling the device will give us a value of 0 or 1 for each qubit (n_qubits)
# this will be repeated for the total number of shots provided (n_shots)
return qml.sample()
raw_samples = prepare_and_sample(w)
# convert the raw samples (bit strings) into integers and count them
samples = []
for sam in raw_samples:
samples.append(int("".join(str(bs) for bs in sam), base=2))
q_probs = np.bincount(samples) / n_shots
print("x_n^2 =\n", c_probs)
print("|<x|n>|^2=\n", q_probs)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 4))
ax1.bar(np.arange(0, 2 ** n_qubits), c_probs, color="blue")
ax1.set_xlim(-0.5, 2 ** n_qubits - 0.5)
ax1.set_xlabel("Vector space basis")
ax1.set_title("Classical probabilities")
ax2.bar(np.arange(0, 2 ** n_qubits), q_probs, color="green")
ax2.set_xlim(-0.5, 2 ** n_qubits - 0.5)
ax2.set_xlabel("Hilbert space basis")
ax2.set_title("Quantum probabilities")
plt.show()
vqls_512_512_5diagonal.py (11.2 KB)
If you want help with diagnosing an error, please put the full error message below:
# Put full error message here
A matrix with dimensions larger than 32 x 32 will run out of my processor RAM.
And, finally, make sure to include the versions of your packages. Specifically, show us the output of qml.about().