# ********************************************************************************
# * 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
# ********************************************************************************
"""Implements the Dalzell kernel-reflection algorithm for solving quantum linear systems via QSVT."""
from collections.abc import Callable
import numpy as np
import numpy.typing as npt
from numpy.polynomial import Chebyshev
from qrisp.algorithms.gqsp.qsvt import QSVT
from qrisp.block_encodings import BlockEncoding
[docs]
def dalzell_inversion(A: BlockEncoding, prep_b: Callable, t: float, eps: float, kappa: float) -> BlockEncoding:
r"""Performs the `Dalzell quantum algorithm <https://arxiv.org/pdf/2406.12086>`_
to solve the Quantum Linear System Problem (QSLP) $A\vec{x}=\vec{b}$, using kernel reflection.
When applied to a state $\ket{0}$, the algorithm prepares a state $\tilde{x}\propto A^{-1}\ket{b}$
within target precision $\epsilon$ of the ideal solution $\ket{x}$.
See also `Dalzell's talk <https://www.youtube.com/watch?v=OwqhdCioj4Y>`_ for more details.
.. warning::
The returned BlockEncoding must be applied to operands in state $\ket{0}$ (and not in state $\ket{b}$).
Instead of solving the system $Ax=b$ directly, we consider the augmented system $A_t x_t = b'$, i.e.,
.. math::
A_t = \begin{pmatrix} A & 0 \\ 0 & t^{-1} \end{pmatrix} \begin{pmatrix} x \\ t \end{pmatrix} = \begin{pmatrix} b \\ 1 \end{pmatrix} = b'
Equivalently, we can write the augmented system as $A\ket{x_t}=\ket{b'}$ where
.. math::
A_t = \ket{0}_a\bra{0}_a \otimes A + t^{-1} \ket{1}_a\bra{1}_a \otimes \ket{0}_s\bra{0}_s,\\
\ket{b'} \propto \ket{0}_a\ket{b}_s + \ket{1}_a\ket{0}_s, \quad \ket{x_t} \propto \|x\|_2\ket{0}_a\ket{x}_s + t\ket{1}_a\ket{0}_s
Here, $\ket{\cdot}_a$ denotes the state of the (1-qubit) ancilla variable, and $\ket{\cdot}_s$ denotes the state of the system variable.
If $A_t x_t = b'$, then
.. math::
G_tx_t=(\mathbb{I} - \ket{b'}\bra{b'})A_t x_t = (\mathbb{I} - \ket{b'}\bra{b'})\ket{b'} = 0
Hence, the solution state $\ket{x_t}$ is a kernel state (sigular value 0) of the operator $G_t=(\mathbb{I} - \ket{b'}\bra{b'})A_t$.
The **kernel reflection** operator $K(G_t)$ reflects about the solution state $\ket{x_t}$ to the augmented system,
and can be implemented via QSVT using a polynomial approximation of the kernel reflection function
$K(\sigma)\colon [0,1] \to [-1,1]$ which maps the sigular value 0 to 1 and all other singular values to -1.
The algorithm consists of the following steps:
- Introduce extra basis state $\ket{e_n}=\ket{1}_a\ket{0}_s$ and prepare it as the initial state.
- Choose $t\simeq \|x\|_2$ and form the augmented linear system $A_tx_t=b'$.
- Apply the kernel reflection operator $K(G_t)$, which **reflects about the solution state to the augmented system**.
- Project the resulting state onto the $\ket{0}_a$ subspace to obtain $\ket{x}$.
The projection step can be implemented by post-selecting on the ancilla qubit being in state $\ket{0}$ after applying the kernel reflection.
The success probability of the algorithm depends on the choice of $t$.
If $t$ is chosen to be on the order of $\|x\|_2$, the algorithm succeeds with constant probability and only requires a constant number of repetitions.
.. image:: /_static/chebyshev_kernel_reflection.png
:align: center
Parameters
----------
A : BlockEncoding
The block-encoded matrix to be inverted. It is assumed that
the eigenvalues of $A/\alpha$ lie within $D_{\kappa} = [-1, -1/\kappa] \cup [1/\kappa, 1]$.
prep_b : Callable
A function ``prep_b(*operands)`` preparing the right hand side $\ket{b}$.
t : float
An estimate $t$ for the norm $\|x\|_2=\|(A/\alpha)^{-1}b\|_2$ for normalized $\|b\|_2=1$.
The success probability depends on the the ratio $t/\|x\|_2$.
The optimal choice is $t=\|x\|_2$. The estimate must lie in the interval $[1, \kappa]$. See below for more details.
eps : float
The target precision $\epsilon$.
kappa : float
An upper bound for the condition number $\kappa$ of $A$.
This value defines the "gap" around zero such that all sigular values (eigenvalues) of $A/\alpha$
lie in the domain $D_{\kappa}$.
Returns
-------
BlockEncoding
A new BlockEncoding instance representing the kernel reflection operator.
Notes
-----
- **Complexity**: The query complexitiy of the algorithm is determined by the polynomial degree and scales as $\mathcal O(\kappa\log(1/\epsilon))$.
The success probability of the algorithm depends on the choice of $t$ relative to $\|x\|_2$:
- If the norm is known up to a constant factor, the algorithm succeeds with constant probability, hence only requires a constant number of repetitions.
- If the norm is unknown, `methods for norm estimation <https://arxiv.org/pdf/2406.12086>`_ can be used to find a suitable $t$ with only a logarithmic overhead in the overall complexity.
For further insights, see `Constant Factor Analysis of Optimal Quantum Linear Solvers in Practice <https://arxiv.org/abs/2604.22185>`_,
and make sure to use :meth:`.resources() <qrisp.block_encodings.BlockEncoding.resources>` to compare the resources of different inversion methods in practice.
- The polynomial is applied to operator $G_t=(\mathbb I - \ket{b'}\bra{b'})A_t$.
Each application of $G_t$ requires 1 call to the block-encoding oracle for $A$
and 2 calls to the state preparation oracle for $b$.
Hence, in contrast to other :meth:`inversion` methods that only prepare the initial state once,
the overall circuit complexity scales multiplicatively with the complexity of the state preparation oracle.
Examples
--------
::
import numpy as np
A = np.array([[ 0.78, -0.01, -0.16, -0.1 ],
[-0.01, 0.57, -0.03, 0.08],
[-0.16, -0.03, 0.69, -0.15],
[-0.1 , 0.08, -0.15, 0.88]])
b = np.array([1, 1, 1, 1])
kappa = np.linalg.cond(A)
print("Condition number of A: ", kappa)
# Condition number of A: 2.020268873491503
::
from qrisp import *
from qrisp.block_encodings import BlockEncoding
from qrisp.gqsp import dalzell_inversion
BA = BlockEncoding.from_array(A)
def prep_b(operand):
prepare(operand, b)
BA_inv = dalzell_inversion(BA, prep_b, 2.0, 0.01, 3.0)
@terminal_sampling
def main():
# Prepare operand variable in state |0>
operand = QuantumFloat(2)
ancillas = BA_inv.apply(operand)
return operand, *ancillas
res_dict = main()
# Post-selection on ancillas being in |0> state
filtered_dict = {k[0]: p for k, p in res_dict.items() \
if all(x == 0 for x in k[1:])}
success_prob = sum(filtered_dict.values())
print("Success probability:", success_prob)
filtered_dict = {k: p / success_prob for k, p in filtered_dict.items()}
amps = np.sqrt([filtered_dict.get(i, 0) for i in range(len(b))])
::
c = (np.linalg.inv(A) @ b) / np.linalg.norm(np.linalg.inv(A) @ b)
print("QUANTUM SIMULATION\n", amps, "\nCLASSICAL SOLUTION\n", c)
# QUANTUM SIMULATION
# [0.51718522 0.43564501 0.60746181 0.41680095]
# CLASSICAL SOLUTION
# [0.51643917 0.4380885 0.60597747 0.41732523]
Alternatively, use ``apply_rus`` to directly obtain the solution state without post-selection:
::
@terminal_sampling
def main():
operand = BA_inv.apply_rus(lambda: QuantumFloat(2))()
return operand
res_dict = main()
amps = np.sqrt([res_dict.get(i, 0) for i in range(len(b))])
print(amps)
# [0.51816163 0.43295659 0.60721322 0.4187472 ]
"""
from qrisp import QuantumBool, control, h, x
def prep_b_ext(*args):
ext = args[0]
ops_A = args[1:]
h(ext)
with control(ext[0], ctrl_state=0):
prep_b(*ops_A)
# Define BlockEncoding of A_t
# A_t = |0><0| ⊗ A + (1 / t) |1><1| ⊗ |0><0|
P0 = BlockEncoding.from_projector(0, 0)
P10 = BlockEncoding.from_projector((1, 0), (1, 0))
A_t = P0.kron(A) + (A.alpha / t) * P10
# Define BlockEncoding of the kernel projector I - |b_t><b_t|
P = BlockEncoding.from_projector(prep_b_ext, kernel=True, num_ops=2)
# Define BlockEncoding of G_t = (I - |b_t><b_t|) A_t
G_t = P @ A_t
# Kernel Reflection
# Rescale kappa to account for larger normalization factor of A_t
kappa_t = kappa * (1.0 + 1.0 / t)
p = _kernel_reflection_cheb(1.0 / kappa_t, eps)
KR = QSVT(G_t, p, kind="Chebyshev", parity="even", rescale=False)
def new_unitary(*args):
anc_ext = args[0]
ancs_ = args[1 : 1 + KR.num_ancs]
args_ = args[1 + KR.num_ancs :]
x(anc_ext)
KR.unitary(*ancs_, anc_ext, *args_)
new_anc_templates = [QuantumBool().template()] + KR._anc_templates
return BlockEncoding(KR.alpha, new_anc_templates, new_unitary, num_ops=KR.num_ops - 1)
def _kernel_reflection_cheb(delta: float, eps: float = 1e-3) -> npt.NDArray[np.float64]:
r"""Constructs the Chebyshev polynomial for the Kernel Reflection
Polynomial $K_{\delta, \ell}(x)$ from `Dalzell (2024) <https://arxiv.org/pdf/2406.12086>`_ Eq 62.
Parameters
----------
delta : float
The spectral gap threshold ($0 < \delta < 1$).
eps : float
The maximum allowed error on the interval $[\delta, 1]$. Defaults to 1e-3.
Returns
-------
ndarray
1-D array containing the coefficients of the Chebyshev series representing the smooth, bounded
approximation of the kernel reflection, ordered from lowest order term to highest.
"""
# 1. Calculate the required degree parameter ell for the
# Kernel Reflection Polynomial to achieve a target error eps.
# Formula 51 (Dalzell 2024, arXiv:2406.12086):
ell = int(np.ceil((1 / (2 * delta)) * np.log(2 / eps)))
# 2. Calculate the Kernel Reflection Polynomial.
# Formula 62 (Dalzell 2024, arXiv:2406.12086):
# 2.1 Define the inner quadratic mapping z(x) in the Chebyshev basis:
# z(x) = (1 + delta^2 - 2x^2) / (1 - delta^2)
# Using the identity x^2 = 0.5*T_0(x) + 0.5*T_2(x), write z(x) exactly as:
c_0 = (delta**2) / (1 - delta**2)
c_2 = -1.0 / (1 - delta**2)
# Instantiate z(x) as a Chebyshev series object:
z_cheb = Chebyshev([c_0, 0.0, c_2])
# 2.2 Define the outer Chebyshev polynomial T_ell(z).
coeffs_T_ell = [0.0] * ell + [1.0]
T_ell = Chebyshev(coeffs_T_ell)
# 2.3 Compose the polynomials: T_ell_z(x) = T_ell(z(x)).
T_ell_z = T_ell(z_cheb)
# 2.4 Calculate the normalization constant.
# Evaluate T_ell at the scalar value z_0 = (1 + delta^2) / (1 - delta^2).
z_0 = (1 + delta**2) / (1 - delta**2)
T_ell_z0 = T_ell(z_0)
# 2.5 Construct the final normalized Kernel Reflection Polynomial.
K_cheb = -1.0 + 2.0 * (T_ell_z + 1.0) / (T_ell_z0 + 1.0)
return K_cheb.coef
