"""********************************************************************************
* Copyright (c) 2026 the Qrisp authors
*
* This program and the accompanying materials are made available under the
* terms of the Eclipse Public License 2.0 which is available at
* http://www.eclipse.org/legal/epl-2.0.
*
* This Source Code may also be made available under the following Secondary
* Licenses when the conditions for such availability set forth in the Eclipse
* Public License, v. 2.0 are satisfied: GNU General Public License, version 2
* with the GNU Classpath Exception which is
* available at https://www.gnu.org/software/classpath/license.html.
*
* SPDX-License-Identifier: EPL-2.0 OR GPL-2.0 WITH Classpath-exception-2.0
********************************************************************************
"""
import numpy as np
import sympy as sp
from scipy.optimize import Bounds, minimize
from qrisp import h, z
from qrisp.algorithms.cold.AGP_params import solve_alpha_gamma_chi
from qrisp.operators import QubitOperator
[docs]
class DCQOProblem:
r"""General structure to formulate Digitized Counterdiabatic Quantum Optimization problems.
This class is used to solve |dcqo_link|
problems with the algorithms `COLD <https://doi.org/10.1103/PRXQuantum.4.010312>`_
(counterdiabatic optimized local driving) or LCD (local counterdiabatic driving).
To run the COLD algorithm on the problem, you need to specify the control Hamiltonian
``H_control`` and the inverse scheduling function ``g_func``. These are not needed for
the LCD algorithm. To learn more about counterdiabatic driving, make sure to check out the `tutorial <https://www.qrisp.eu/general/tutorial/CD.html>`_.
Parameters
----------
Q : np.array
The QUBO matrix.
H_init : :ref:`QubitOperator`
Hamiltonian, the system is at the time t=0.
H_prob : :ref:`QubitOperator`
Hamiltonian, the system evolves to for t=T.
A_lam : :ref:`QubitOperator`
Operator holding an appoximation for the adiabatic gauge potential (AGP).
agp_coeffs : callable
The parameters for the adiabatic gauge potential (AGP). If the COLD method is being used,
they must depend on the optimization pulses in ``H_control``.
lam_func : callable
A function $\lambda(t, T)$ mapping $t \in [0, T]$ to $\lambda \in [0, 1]$. This function needs to return
a `sympy <https://docs.sympy.org/>`_ expression with $t$ and $T$ as `sympy.Symbols <https://docs.sympy.org/latest/modules/core.html#sympy.core.symbol.Symbol>`_.
H_control : :ref:`QubitOperator`, optional
Hamiltonian specifying the control pulses for the COLD method. If not given, the LCD method is used automatically.
qarg_prep : callable, optional
A function receiving a :ref:`QuantumVariable` for preparing the inital state.
By default, the groundstate of the x-operator $\ket{-}^n$ is prepared.
Examples
--------
For a quick demonstration we build a DCQO problem instance for a 4x4 QUBO. We choose a first order AGP ansatz with uniform coefficients and solve it with LCD.
::
import numpy as np
import sympy as sp
from qrisp.operators.qubit import X, Y, Z
from qrisp.algorithms.cold import DCQOProblem
from qrisp import QuantumVariable
Q = np.array([
[-1.2, 0.40, 0.0, 0.0],
[ 0.40, 0.30, 0.20, 0.0],
[ 0.0, 0.20,-1.1, 0.30],
[ 0.0, 0.0, 0.30,-0.80]
])
N = Q.shape[0]
# Define QUBO problem hamiltonian
h = -0.5 * np.diag(Q) - 0.5 * np.sum(Q, axis=1)
J = 0.5 * Q
H_init = 1 * sum([X(i) for i in range(N)])
H_prob = (sum([sum([J[i][j]*Z(i)*Z(j) for j in range(i)]) for i in range(N)])
+ sum([h[i]*Z(i) for i in range(N)]))
# Create AGP
A_lam = sum([Y(i) for i in range(N)]) # uniform
# Function for uniform AGP coefficients
def alpha(lam):
A = lam * h
B = 1 - lam
nom = np.sum(A + 4*B*h)
denom = 2 * (np.sum(A**2) + N * (B**2)) + 4 * (lam**2) * np.sum(np.tril(J, -1).sum(axis=1))
alph = nom/denom
alph = [alph]*N
return alph
# Simple scheduling function 0 -> 1
def lam():
t, T = sp.symbols("t T", real=True)
lam_expr = t/T
return lam_expr
# Create problem instance
lcd_problem = DCQOProblem(Q, H_init, H_prob, A_lam, alpha, lam)
# Run problem with LCD algorithm
qarg = QuantumVariable(N)
res = lcd_problem.run(qarg, N_steps=4, T=12, method="LCD")
print(res)
::
{'1011': [0.40630593694063055, np.float64(-2.5)], '1111': [0.16247837521624783, np.float64(-0.9999999999999999)], '0111': [0.13156868431315685, np.float64(-0.6000000000000001)], '1000': [0.06881931180688193, np.float64(-1.2)], '0011': [0.05949940500594993, np.float64(-1.3)], '1010': [0.04499955000449995, np.float64(-2.3)], '1101': [0.04084959150408495, np.float64(-0.9)], '0110': [0.019769802301976978, np.float64(-0.40000000000000013)], '1100': [0.01815981840181598, np.float64(-0.09999999999999998)], '0100': [0.013679863201367985, np.float64(0.3)], '0001': [0.010399896001039988, np.float64(-0.8)], '0000': [0.007659923400765992, np.float64(0.0)], '1110': [0.006329936700632993, np.float64(-0.7999999999999999)], '0101': [0.0052899471005289946, np.float64(-0.5)], '1001': [0.0024299757002429973, np.float64(-2.0)], '0010': [0.0017599824001759982, np.float64(-1.1)]}
We get a dictionary where the key is the quantum state and the values are lists of [probability, cost].
So our most likely result is '1011' with probabilty 0.4 and the QUBO cost $x^T Q x = -2.5$.
.. |dcqo_link| raw:: html
<a href="https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.4.L042030" target="_blank">DCQO</a>
"""
def __init__(
self,
Q,
H_init,
H_prob,
A_lam,
agp_coeffs,
lam_func,
H_control=None,
qarg_prep=None,
):
# Scheduling function
self.lam_func = lam_func
# Operators
self.agp_coeffs = agp_coeffs
self.H_init = H_init
self.H_prob = H_prob
self.A_lam = A_lam
self.H_control = H_control
self.qarg_prep = qarg_prep
# Qubo characteristics
self.Q = Q
self.J = 0.5 * Q
self.h = -0.5 * np.diag(Q) - 0.5 * np.sum(Q, axis=1)
# Placeholder for the _precompute_timegrid function
self.lam = None
self.lamdot = None
def _precompute_timegrid(self, N_steps, T, method):
"""Compute lambda(t, T) and the time-derivative lambdadot(t, T)
for each timestep.
Parameters
----------
N_steps : int
Number of timesteps.
T : float
Evolution time for the simulation.
method : str
The method to solve the QUBO with (either ``LCD`` or ``COLD``).
Returns
-------
lam : array
The parametrized timefunction, specified by ``lam_func`` for t in [0, T].
lamdot : array
The time derivative of ``lam_func`` for t in [0, T].
"""
# Sympy symbols for t and T
t_sym, T_sym = sp.symbols("t T", real=True)
# Array for t values
dt = T / N_steps
t_list = np.linspace(dt, T, N_steps)
# Sympy functions for lam and lamdot
lam_func = sp.lambdify((t_sym, T_sym), self.lam_func(), "numpy")
lamdot_expr = sp.diff(self.lam_func(), t_sym)
lamdot_func = sp.lambdify((t_sym, T_sym), lamdot_expr, "numpy")
# Compute functions of values t_list
lam = lam_func(t_list, T)
lamdot = lamdot_func(t_list, T)
# Convert lamdot to a constant list to make it subscriptable by timestep
if isinstance(lamdot, float):
lamdot = np.full(N_steps, lamdot)
else:
lamdot = np.array(lamdot)
self.lam = lam
self.lamdot = lamdot
# Functions for t = g(lam) and derivative (later needed for opt pulses)
# must only be calculated for COLD, not for LCD
if method == "COLD":
g = t_list # g(lam[s]) = t[s] by definition
g_deriv = 1.0 / lamdot # dt/dlambda = 1 / (dlambda/dt)
self.g = g
self.g_deriv = g_deriv
def _precompute_opt_pulses(self, N_steps, T, t_list, N_opt, CRAB=False):
"""Precompute optimization pulses for COLD routine that will be scaled by optimized paramters.
Parameters
----------
N_steps : int
Number of timesteps.
T : float
Evolution time for the simulation.
t_list : list
Time values for each step.
N_opt : int
Number of optimization parameters in ``H_control``.
CRAB : bool
If ``True``, the CRAB optimization method is being used. The default is ``False``.
Returns
-------
sin_matrix :
Numpy or sympy (if CRAB) array holding the opt pulse for each timestep.
cos_matrix :
Numpy or sympy (if CRAB) array holding the derivative of the opt pulse for each timestep.
"""
# Precompute f (sine) and f_deriv (cosine) for each timestep as numpy arrays
sin_matrix = np.zeros((N_steps, N_opt))
cos_matrix = np.zeros((N_steps, N_opt))
if CRAB:
# Random CRAB parameters
r_params = np.random.uniform(-0.5, 0.5, N_opt)
else:
# Otherwise add nothing
r_params = np.zeros(N_opt)
for k in range(N_opt):
sin_matrix[:, k] = np.sin(np.pi * (k + 1 + r_params[k]) * t_list / T)
cos_matrix[:, k] = (
(np.pi * (k + 1 + r_params[k])) * np.cos(np.pi * (k + 1 + r_params[k]) * self.g) * self.g_deriv
)
return sin_matrix, cos_matrix
[docs]
def apply_lcd_hamiltonian(self, qarg, N_steps, T):
"""Simulate the local counterdiabatic driving (LCD) Hamiltonian on a
quantum argument via trotterization. The LCD Hamiltonian consists
of the system Hamiltonian and the adiabatic gauge potential (AGP).
Parameters
----------
qarg : :ref:`QuantumVariable`
The quantum argument to which the quantum circuit is applied.
N_steps : int
Number of steps in which the timefunction ``lambda`` is split up.
T : float
Evolution time for the simulation.
"""
self.qarg_prep(qarg)
dt = T / N_steps
# Compute time-function lamda(t, T) and the derivative lamdot(t, T)
self._precompute_timegrid(N_steps, T, "LCD")
# Trotterize Hamiltonian in different parts with each one needing different coefficients
U1 = self.H_init.trotterization()
U2 = self.H_prob.trotterization()
if isinstance(self.A_lam, QubitOperator):
# Uniform AGP coefficients
U3 = self.A_lam.trotterization()
else:
# Non-uniform AGP coefficients
U3 = [A_lam.trotterization() for A_lam in self.A_lam]
# Apply hamiltonian to qarg for each timestep
for s in range(N_steps):
# Get alpha for the timestep
coeffs = self.agp_coeffs(self.lam[s])
# print(f"Alpha = {coeffs}")
# H_0 contribution scaled by dt
U1(qarg, t=dt * (1 - self.lam[s]))
U2(qarg, t=dt * self.lam[s])
# AGP contribution scaled by dt* lambda_dot(t)
if isinstance(U3, list):
for idx, U in enumerate(U3):
U(qarg, t=dt * self.lamdot[s] * coeffs[idx])
else:
U3(qarg, t=dt * self.lamdot[s] * coeffs[0])
[docs]
def apply_cold_hamiltonian(self, qarg, N_steps, T, opt_params, CRAB=False):
"""Simulate counterdiabatic optimized local driving (COLD) Hamiltonian
on a quantumvariable via trotterization. The COLD Hamiltonian consists
of the system Hamiltonian, the adiabatic gauge potential (AGP) and
local pulses (given by ``H_control``) with optimized parameters.
Parameters
----------
qarg : :ref:`QuantumVariable`
The quantum argument to which the quantum circuit is applied.
N_steps : int
Number of steps in which the timefunction ``lambda`` is split up.
T : float
Evolution time for the simulation.
opt_params : list
Either the optimized parameters or the corresponding `sympy.Symbols <https://docs.sympy.org/latest/modules/core.html#sympy.core.symbol.Symbol>`_.
CRAB : bool, optional
If ``True``, the CRAB optimization method is being used. The default is ``False``.
"""
# Initialize qarg
self.qarg_prep(qarg)
# Compute time-function lamda(t, T) and the derivative lamdot(t, T)
self._precompute_timegrid(N_steps, T, "COLD")
# Precompute opt pulses
dt = T / N_steps
t_list = np.linspace(dt, T, int(N_steps))
sin_matrix, cos_matrix = self._precompute_opt_pulses(N_steps, T, t_list, N_opt=len(opt_params), CRAB=CRAB)
beta = opt_params
# Trotterize Hamiltonian in different parts with each one needing different coefficients
U1 = self.H_init.trotterization()
U2 = self.H_prob.trotterization()
if isinstance(self.A_lam, QubitOperator):
# Uniform AGP coefficients
U3 = self.A_lam.trotterization()
else:
# Non-uniform AGP coefficients
U3 = [A_lam.trotterization() for A_lam in self.A_lam]
U4 = self.H_control.trotterization()
# Apply hamiltonian to qarg for each timestep
for s in range(N_steps):
# Get alpha, f and f_deriv for the timestep
f = sin_matrix[s, :] @ beta
f_deriv = cos_matrix[s, :] @ beta
alpha = self.agp_coeffs(self.lam[s], f, f_deriv)
# H_init contribution scaled by dt*(1-lam)
U1(qarg, t=dt * (1 - self.lam[s]))
# H_prob contribution scaled by dt*lam
U2(qarg, t=dt * (self.lam[s]))
# AGP contribution scaled by dt*lambda_dot(t)*alpha
if isinstance(U3, list):
# Non-uniform alpha
for idx, U in enumerate(U3):
U(qarg, t=dt * (self.lamdot[s] * alpha[idx]))
else:
# Uniform alpha
U3(qarg, t=dt * (self.lamdot[s] * alpha[0]))
# Control pulse contribution scaled by opt parameters f
U4(qarg, t=dt * f)
[docs]
def compile_U_cold(self, qarg, N_opt, N_steps, T, CRAB=False):
"""Compiles the circuit that is created by the :meth:`apply_cold_hamiltonian <qrisp.cold.DCQOProblem.apply_cold_hamiltonian>` method.
Parameters
----------
qarg : :ref:`QuantumVariable`
The argument to which the COLD circuit is applied.
N_opt : int
Number of optimization parameters in ``H_control``.
N_steps : int
Number of timesteps.
T : float
Evolution time for the simulation.
CRAB : bool, optional
If ``True``, the CRAB optimization method is being used. The default is ``False``.
Returns
-------
compiled_qc : :ref:`QuantumCircuit`
The compiled quantum circuit.
"""
temp = list(qarg.qs.data)
# Initzialize parameters as symbols
params = [sp.Symbol("par_" + str(i)) for i in range(N_opt)]
self.apply_cold_hamiltonian(qarg, N_steps, T, params, CRAB=CRAB)
intended_measurements = list(qarg)
compiled_qc = qarg.qs.compile(intended_measurements=intended_measurements)
qarg.qs.data = temp
return compiled_qc
[docs]
def optimization_routine(
self,
qarg,
N_opt,
N_steps,
T,
qc,
CRAB,
optimizer,
options,
objective,
bounds,
precision=0.01,
exp_value_backend=None,
):
"""Subroutine for the optimization method used in COLD.
The initial values are set and the optimization via is conducted here.
Parameters
----------
qarg : :ref:`QuantumVariable`
The argument to which the H_prob circuit is applied.
N_opt : int
Number of optimization parameters in ``H_control``.
N_steps : int
Number of time steps for the simulation.
T : float
Evolution time for the simulation.
qc : :ref:`QuantumCircuit`
The COLD circuit that is applied before measuring the qarg.
CRAB : bool, optional
If ``True``, the CRAB optimization method is being used. The default is ``False``.
optimizer : str
Specifies the `SciPy optimization routine <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html>`_.
options : dict
A dictionary of solver options.
objective : str
The objective function to be minimized (``exp_value``, ``agp_coeff_magnitude``, ``agp_coeff_amplitude``). Default is ``exp_value``.
bounds : tuple
The parameter bounds for the optimizer. Default is (-2, 2).
precision : float, optional
Precision for expectation value calculations. Default is 0.01.
exp_value_backend : BackendLike, optional
Backend for expectation value calculations, if ``exp_value`` is used as objective function.
If provided, uses measurement-based expectation value with this backend.
If not provided, tries statevector first and falls back to measurement. Default is None.
Returns
-------
res.x: array
The optimized parameters of the problem instance.
"""
# Different objective functions: exp_value, agp coeffs magnitude, agp coeffs amplitude
# Precompute costs for statevector method
n_qubits = len(qarg)
num_states = 1 << n_qubits
bit_indices = np.arange(n_qubits - 1, -1, -1, dtype=np.uint64)
state_indices = np.arange(num_states, dtype=np.uint64)
basis = ((state_indices[:, None] >> bit_indices) & 1).astype(np.int8)
costs = np.einsum("bi,ij,bj->b", basis, self.Q, basis)
# Expectation value of the QUBO Hamiltonian
def objective_exp(params, CRAB):
# Dict to assign the optimization parameters
subs_dic = {sp.Symbol("par_" + str(i)): params[i] for i in range(len(params))}
if exp_value_backend is not None:
# Use provided backend
exp_val = self.H_prob.expectation_value(
qarg,
precision=precision,
compile=False,
subs_dic=subs_dic,
precompiled_qc=qc,
backend=exp_value_backend,
)()
else:
# Try statevector first, fallback to measurement
try:
bound_qc = qc.bind_parameters(subs_dic)
sv = bound_qc.statevector_array()
probs = np.abs(sv) ** 2
exp_val = float(np.dot(probs, costs))
except (MemoryError, ValueError, RuntimeError):
exp_val = self.H_prob.expectation_value(
qarg,
precision=precision,
compile=False,
subs_dic=subs_dic,
precompiled_qc=qc,
)()
return exp_val
# Magnitude of the AGP coefficients (coeffs are treated as uniform for simplification)
# (sum of absolute values for each timestep)
def objective_mag(params, CRAB):
# Precompute opt pulses to be multiplied with opt params
t_list = np.linspace(T / N_steps, T, int(N_steps))
sin_matrix, cos_matrix = self._precompute_opt_pulses(N_steps, T, t_list, N_opt=len(params), CRAB=CRAB)
magnitude = 0
# Iterate through lambda(t)
for s in range(N_steps):
# Get alpha, f and f_deriv for the timestep
f = sin_matrix[s, :] @ params
f_deriv = cos_matrix[s, :] @ params
alpha, gamma, chi = solve_alpha_gamma_chi(self.h, self.J, self.lam[s], f, f_deriv, uniform=True)
magnitude += np.abs(gamma[0]) + np.abs(chi[0]) + np.abs(alpha[0])
return magnitude
init_point = np.random.rand(N_opt) * np.pi / 2
# Define objective function
if objective == "exp_value":
objective = objective_exp
elif objective == "agp_coeff_magnitude":
objective = objective_mag
else:
raise ValueError("{objective} is not a valid option as objective.")
res = minimize(
objective,
init_point,
method=optimizer,
options=options,
bounds=Bounds(*bounds),
args=(CRAB),
)
return res.x, objective(res.x, CRAB)
[docs]
def QUBO_cost(self, res):
"""Returns the cost y = x^T Q x for a given binary array x.
Parameters
----------
res : np.array
The array to calculate the cost for.
Returns
-------
cost : float
The QUBO cost.
"""
cost = res @ self.Q @ res
return cost
[docs]
def run(
self,
qarg,
N_steps,
T,
method,
N_opt=None,
CRAB=False,
optimizer="COBYQA",
objective="agp_coeff_magnitude",
bounds=(),
options={},
mes_kwargs={},
precision=0.01,
exp_value_backend=None,
):
"""Run the specific DCQO problem instance with given quantum arguments, number of timesteps,
evolution time and method.
There is also the option to choose if parameter optimization via the expectation value objective function should be done via a simulator or real quantum backend.
If the user chooses a quantum backend this iterative optimization can potentially use a lot of computing time.
Parameters
----------
qarg : :ref:`QuantumVariable`
The argument to which the DCQO circuit is applied.
N_steps : int
Number of time steps for the simulation.
T : float
Evolution time for the simulation.
method : str
Method to solve the QUBO with. Either ``LCD`` or ``COLD``.
N_opt : int
Number of optimization parameters in ``H_control``.
CRAB : bool
If ``True``, the CRAB optimization method is being used. The default is ``False``.
optimizer : str, optional
Specifies the `SciPy optimization routine <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html>`_.
We set the default to ``Powell``.
options : dict
A dictionary of solver options.
objective : str
The objective function to be minimized (``exp_value``, ``agp_coeff_magnitude``). Default is ``agp_coeff_magnitude``.
bounds : tuple
The parameter bounds for the optimizer. Default is (-2, 2).
options : dict
Additional options for the Scipy solver.
mes_kwargs : dict, optional
The keyword arguments for the measurement function. Default is an empty dictionary.
precision : float, optional
Precision for expectation value calculations. Default is 0.01.
exp_value_backend : BackendLike, optional
Backend for expectation value calculations, if ``exp_value`` is used as objective function.
If provided, uses measurement-based expectation value with this backend. Default is the Qrisp statevector simulator.
Returns
-------
res_dict : dict
The optimal result after running DCQO problem for a specific problem instance. It contains the measurement results after applying the optimal DCQO circuit to the quantum argument.
"""
# If no prep for qarg is specified, use uniform superposition state
if self.qarg_prep is None:
def qarg_prep(q):
h(q)
z(q)
return q
self.qarg_prep = qarg_prep
# Run COLD routine
if method == "COLD":
qarg1, qarg2 = qarg.duplicate(), qarg.duplicate()
# If we optimize the Hamiltonian expectation value,
# compile COLD routine into a circuit for the optimization
if objective == "exp_value":
U_circuit = self.compile_U_cold(qarg1, N_opt, N_steps, T, CRAB)
else:
self._precompute_timegrid(N_steps, T, "COLD")
U_circuit = None
# Find optimal params for control pulse
opt_params, cost = None, None
# Do optimization 3 times and choose best result
for i in range(3):
opt_params_temp, cost_temp = self.optimization_routine(
qarg2,
N_opt,
N_steps,
T,
U_circuit,
CRAB,
optimizer,
options,
objective=objective,
bounds=bounds,
precision=precision,
exp_value_backend=exp_value_backend,
)
if cost is None or cost_temp < cost:
opt_params = opt_params_temp
cost = cost_temp
# Apply hamiltonian with optimal parameters
# Here we do not want the randomized parameters to be included -> CRAB=False in any case
self.apply_cold_hamiltonian(qarg, N_steps, T, opt_params, CRAB=False)
# Run LCD routine
elif method == "LCD":
self.apply_lcd_hamiltonian(qarg, N_steps, T)
else:
raise ValueError(f'"{method}" is not an option for method. Choose "LCD" or "COLD".')
# Measure qarg
if "shots" not in mes_kwargs:
mes_kwargs["shots"] = 5000
res_dict = dict(qarg.get_measurement(**mes_kwargs))
# Add qubo cost in result dict
for res in res_dict.keys():
res_array = np.fromiter(res, dtype=int)
res_dict[res] = [res_dict[res], self.QUBO_cost(res_array)]
return res_dict
