# ********************************************************************************
# * 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
# ********************************************************************************
"""Defines the DCQOProblem class for solving QUBOs via COLD/LCD counterdiabatic driving."""
import numpy as np
import sympy as sp
from scipy.optimize import Bounds, minimize
from qrisp import h, z
from qrisp.algorithms.cold._fast_trotterization import fast_trotterization
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.
Must be symmetric: the cost is read as $y = x^T Q x$ with off-diagonal entries
counted twice, and the Ising coefficients $J_{ij} = \frac{1}{2} Q_{ij}$,
$h_i = -\frac{1}{2} \sum_j Q_{ij}$ derived from it are only valid for
symmetric $Q$. Symmetrize an upper-triangular matrix with
``Q = 0.5 * (Q + Q.T)`` before passing it in.
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 QubitOperator, 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.sum(Q, axis=1)
J = 0.5 * Q
H_init = 1 * sum([X(i) for i in range(N)])
# QubitOperator.sum merges all terms in a single pass -- for a dense
# QUBO (O(N^2) coupling terms), Python's built-in sum() would instead
# fold pairwise via repeated __add__ calls, each copying the full
# running total, costing O(N^4) instead of O(N^2).
H_prob = QubitOperator.sum(
J[i][j] * Z(i) * Z(j) for i in range(N) for j in range(i)
) + QubitOperator.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. This minimizes the action
# S = Tr[G^2], G = d_lambda H + i[A_lam, H], for the ansatz A_lam = sum_i Y_i.
def alpha(lam):
J_sq = np.sum(J**2) - np.sum(np.diag(J)**2) # sum_{i != j} J_ij**2
nom = np.sum(h)
denom = 2 * (N * (1 - lam)**2 + np.sum((lam * h)**2) + lam**2 * J_sq)
return [-nom/denom]*N
# 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.5666, np.float64(-2.5)],
'1111': [0.0986, np.float64(-1.0)],
'1001': [0.087, np.float64(-2.0)],
'1010': [0.0846, np.float64(-2.3)],
'0111': [0.0586, np.float64(-0.6000000000000001)],
'1000': [0.031, np.float64(-1.2)],
'1110': [0.0218, np.float64(-0.8000000000000003)],
'0011': [0.0216, np.float64(-1.3)],
'0110': [0.0074, np.float64(-0.40000000000000013)],
'0001': [0.0056, np.float64(-0.8)],
'1101': [0.0052, np.float64(-0.9)],
'0101': [0.005, np.float64(-0.5)],
'0100': [0.0026, np.float64(0.3)],
'0000': [0.0022, np.float64(0.0)],
'0010': [0.0018, np.float64(-1.1)],
'1100': [0.0004, np.float64(-0.09999999999999998)]
}
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 probability around 0.57 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__( # noqa: PLR0913, PLR0917 -- public, keyword-callable API shape
self,
Q,
H_init,
H_prob,
A_lam,
agp_coeffs,
lam_func,
H_control=None,
qarg_prep=None,
):
"""Construct a DCQOProblem instance. See class docstring for parameter details."""
# 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.sum(Q, axis=1)
# Placeholder for the _precompute_timegrid function
self.lam = None
self.lamdot = None
# Random CRAB frequency offsets r_k of the current run, drawn on first use by _precompute_opt_pulses
self._crab_r = 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)
# Midpoint rule: avoids t=0/T exactly, where some lam_funcs have lamdot=0, diverging g_deriv below.
dt = T / N_steps
t_list = (np.arange(N_steps) + 0.5) * dt
# 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
# Inverse scheduling function g(lam) and its derivative (later needed for opt pulses)
# must only be calculated for COLD, not for LCD
if method == "COLD":
# g is the *normalized* time t/T (dimensionless), as in Eq. (18) of
# https://doi.org/10.1103/PRXQuantum.4.010312. Since lam[s] = lam_func(t[s], T) by
# construction, the inverse is simply t[s]/T for any monotonic lam_func -- the closed
# form of g is only needed to evaluate it off the time grid.
g = t_list / T
# lam is dimensionless but lamdot = dlam/dt is not: lam depends on t and T only through
# t/T, so lamdot scales as 1/T. The chain rule in _precompute_opt_pulses needs
# dg/dlam = (1/T) * dt/dlam, which is dimensionless.
g_deriv = 1.0 / (T * lamdot)
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. The optimized beta only fit the basis they were optimized for, so
# every later call (objective evaluations, final circuit) must see the same r_k.
if self._crab_r is None or len(self._crab_r) != N_opt:
self._crab_r = np.random.uniform(-0.5, 0.5, N_opt)
r_params = self._crab_r
else:
# Otherwise add nothing
r_params = np.zeros(N_opt)
# Pulse basis f_k = sin(omega_k * g(lam)) with principal frequencies omega_k = 2*pi*k and
# CRAB randomization k -> k*(1 + r_k), see Eq. (28) of https://doi.org/10.1103/PRXQuantum.4.010312.
for k in range(N_opt):
omega = 2 * np.pi * (k + 1) * (1 + r_params[k])
sin_matrix[:, k] = np.sin(omega * self.g)
cos_matrix[:, k] = omega * np.cos(omega * 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).
``H_init``, ``H_prob``, ``A_lam``, and (for COLD) ``H_control`` are each
Trotterized, using a fast native-gate path for Ising-type Hamiltonians and
falling back to :meth:`QubitOperator.trotterization` otherwise.
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 = fast_trotterization(self.H_init)
U2 = fast_trotterization(self.H_prob)
if isinstance(self.A_lam, QubitOperator):
# Uniform AGP coefficients
U3 = fast_trotterization(self.A_lam)
else:
# Non-uniform AGP coefficients
U3 = [fast_trotterization(A_lam) 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. Must match the midpoint grid used in _precompute_timegrid.
dt = T / N_steps
t_list = (np.arange(int(N_steps)) + 0.5) * dt
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 = fast_trotterization(self.H_init)
U2 = fast_trotterization(self.H_prob)
if isinstance(self.A_lam, QubitOperator):
# Uniform AGP coefficients
U3 = fast_trotterization(self.A_lam)
else:
# Non-uniform AGP coefficients
U3 = [fast_trotterization(A_lam) for A_lam in self.A_lam]
U4 = fast_trotterization(self.H_control)
# 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):
"""Compile the circuit created by the COLD Hamiltonian application step.
See :meth:`apply_cold_hamiltonian <qrisp.cold.DCQOProblem.apply_cold_hamiltonian>` for the step being compiled.
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( # noqa: PLR0913, PLR0917 -- public, keyword-callable API shape
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``). Required --
this parameter is positional here; :meth:`run` is what applies the ``exp_value`` default.
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
if objective == "exp_value" and exp_value_backend is None:
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)
# Per-factor eigenvalue on computational basis state `basis`, for factor
# types that are diagonal in that basis. X, Y and the ladder operators
# (C, A) are not diagonal and have no such per-state value, so a term
# containing any of them disables the fast path entirely below.
diagonal_eigenvalues = {"Z": 1 - 2 * basis, "P0": 1 - basis, "P1": basis}
costs = np.zeros(num_states)
for term, coeff in self.H_prob.terms_dict.items():
term_value = np.ones(num_states)
for idx, factor in term.factor_dict.items():
if factor not in diagonal_eigenvalues:
costs = None
break
term_value = term_value * diagonal_eigenvalues[factor][:, idx]
else:
costs += coeff * term_value
if costs is None:
break
# Expectation value of the QUBO Hamiltonian
def objective_exp(params, CRAB):
"""Expectation value of the problem Hamiltonian for the given optimization parameters."""
# 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:
if costs is None:
raise ValueError("H_prob has non-diagonal terms; fast path unavailable")
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 (sum of absolute values over all timesteps)
def objective_mag(params, CRAB):
"""Weighted magnitude of the AGP coefficients for the given optimization parameters."""
# Precompute opt pulses to be multiplied with opt params.
# Must match the midpoint grid used in _precompute_timegrid.
t_list = (np.arange(int(N_steps)) + 0.5) * (T / N_steps)
sin_matrix, cos_matrix = self._precompute_opt_pulses(N_steps, T, t_list, N_opt=len(params), CRAB=CRAB)
# The objective always scores the second-order Y/XY/ZY template of _solve_alpha_gamma_chi,
# not the operator the circuit applies (local Y or nested commutators, see apply_cold_hamiltonian).
# It is a generic proxy for the AGP size. Only uniformity is matched: a single merged A_lam
# (QubitOperator) means uniform coefficients, a list one coefficient per qubit.
uniform = isinstance(self.A_lam, QubitOperator)
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=uniform)
# Weight by lamdot: the AGP enters the circuit as dt*lamdot[s]*alpha (see
# apply_cold_hamiltonian), and alpha itself diverges as 1/lamdot through f_deriv.
# Without this factor the sum is dominated by the first and last timestep, whose
# actual contribution to the evolution is negligible, and it grows as N_steps**3.
magnitude += self.lamdot[s] * (np.sum(np.abs(alpha)) + np.sum(np.abs(gamma)) + np.sum(np.abs(chi)))
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(
f"{objective} is not a valid option as objective. "
"Valid options are 'exp_value' and 'agp_coeff_magnitude'."
)
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( # noqa: PLR0913, PLR0917 -- public, keyword-callable API shape
self,
qarg,
N_steps,
T,
method,
N_opt=None,
CRAB=False,
optimizer="COBYQA",
objective="exp_value",
bounds=(),
options={},
mes_kwargs={},
precision=0.01,
exp_value_backend=None,
):
"""Run the 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. Default is ``exp_value``, which optimizes the control
pulses against the measured energy of ``H_prob`` and is the reliable choice. The alternative
``agp_coeff_magnitude`` is experimental. It is a cheap proxy that minimizes the magnitude of the
coefficients of a generic second-order AGP template instead of running the circuit, so it does
not score the AGP operators the circuit applies. It avoids the simulation cost but can select
controls that perform worse than no control at all, so opt into it only deliberately.
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):
"""Prepare the default uniform superposition state on ``q``."""
h(q)
z(q)
return q
self.qarg_prep = qarg_prep
# Run COLD routine
if method == "COLD":
if isinstance(N_opt, (bool, np.bool_)) or not isinstance(N_opt, (int, np.integer)) or N_opt < 1:
raise ValueError(f"N_opt must be a positive integer for the COLD method (received {N_opt!r}).")
qarg1, qarg2 = qarg.duplicate(), qarg.duplicate()
# Every run optimizes over a fresh random CRAB basis
self._crab_r = None
# 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
self.apply_cold_hamiltonian(qarg, N_steps, T, opt_params, CRAB=CRAB)
# 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, prob in res_dict.items():
res_array = np.fromiter(res, dtype=int)
res_dict[res] = [prob, self.QUBO_cost(res_array)]
return res_dict
