> ## Documentation Index
> Fetch the complete documentation index at: https://docs.classiq.io/llms.txt
> Use this file to discover all available pages before exploring further.
# Matrix Inversion with Quantum Singular Value Transform (QSVT)
Open this notebook in GitHub to run it yourself
> **Quantum Matrix Inversion (or Quantum Linear Solver) via QSVT** provides a framework for solving linear systems exponentially faster than known classical methods under certain conditions [\[1\]](#ref-grand). The algorithm employs a polynomial transformation of singular values to block-encode the inverse of a matrix, which can then be applied on a prepared quantum state to solve a linear equation. Given an efficient routine to embed the classical matrix as a quantum function (block-encoding), this algorithm gives a clean and optimal way to implement matrix inversion compared to other quantum methods. Quantum linear solvers based on block-encoding have many applications, e.g., in [Plasma Physics](https://github.com/Classiq/classiq-library/blob/main/applications/plasma/vlasov_ampere/vlasov_ampere.ipynb) and [Fluid Dynamics](https://github.com/Classiq/classiq-library/blob/main/applications/CFD/QLS_for_hybrid_solvers/qls_qsvt.ipynb).
>
> * **Input:** An $N \times N$ matrix $A$ with condition number $\kappa$ (the ratio between the maximal and minimal singular values), and a vector $|b\rangle$ prepared as a quantum state.
> * **Promise:** The matrix $A$ is well-conditioned (invertible, with bounded spectrum) and can be efficiently block-encoded in a unitary $U_A$ (see technical note at the end of this demo).
> * **Output:** A quantum state proportional to $|x\rangle = A^{-1} |b\rangle$, encoding the solution of the linear system $A x = b$.
>
> **Complexity:** Using QSVT, matrix inversion requires $O(\kappa \log(\kappa / \epsilon))$ queries of the block-encoding unitary, where $\epsilon$ is the precision parameter. Typically, an efficient block-encoding scales as $\mathrm{poly}(\log N)$, resulting in an overall $\tilde{O}(\kappa \,\mathrm{polylog}(N, 1 / \epsilon))$ complexity for the QSVT linear solver. This improves exponentially in $\epsilon$ upon the original Harrow-Hassidim-Lloyd ([HHL](https://github.com/Classiq/classiq-library/blob/main/algorithms/quantum_linear_solvers/hhl/hhl.ipynb)) algorithm, and exponentially in $N$ compared to classical linear solvers, while matching known lower bounds for quantum linear solvers \[2].
>
> (As with all quantum linear solvers, extracting the full solution state $|x\rangle$ removes the speedup; the advantage arises when only global properties such as expectation values, overlaps, or samples are required.)
>
> ***
>
> **Keywords:** Linear systems, Quantum Singular Value Transformation (QSVT), Block encoding, Matrix inversion, Exponential speedup, Oracle/Query complexity.
In this demo we will restrict ourselves to matrices of size $N\times N$, with $N=2^n$. A generalization to arbitrary $N$ can be done by completing the matrix into a larger one, or modifying the QSVT routine to an arbitrary $N$. We emphasize that compared to the basic HHL algorithm, the QSVT approach *is not restricted to Hermitian matrices*.
The QSVT algorithm is based on the idea of block-encoding a matrix into a larger unitary:
$$
U_{(A,s)} = \begin{pmatrix}
A/s & * \\
* & *
\end{pmatrix}, \tag{1}
$$
where $s$ is some scaling factor, which is sometimes necessary for this unitarity completion.
Given an access to $U_{(A,s)}$, matrix inversion with QSVT implements a block-encoding for the inverse of $A$
$$
U_{(A^{-1},s')} = \begin{pmatrix}
A^{-1}/s' & * \\
* & *
\end{pmatrix}.
$$
Subsequently, we can use this unitary to solve a linear system, $A|x\rangle = |b\rangle$, by applying it on an initially prepared state $|b\rangle_n|0\rangle_{\rm block}$
(the state $|0\rangle_{\rm block}$ identifies the block in which the matrix is encoded). Up to post-selection we have:
$$
U_{(A^{-1},s')} |b\rangle_n|0\rangle_{\rm block} = \frac{A^{-1}}{s'}|b\rangle_n|0\rangle_{\rm block} + \text{ garbage}.
$$
In more detail, Singular Value Decomposition (SVD) of a matrix $A$ is defined by
$$
A = W\Sigma V^\dagger,
$$
where $W$ and $V$ are unitary matrices, and $\Sigma$ is a diagonal matrix that contains the singular values of $A$ (The square root of the eigenvalues of $A A^{\dagger}$). A singular value polynomial transform refers to
$$
\mathrm{Poly}^{\rm (SV)}(A) = \left\{\begin{array}{l l}
W\mathrm{Poly}(\Sigma) V^\dagger, & \text{ for odd polynomial} \\
V\mathrm{Poly}(\Sigma) V^\dagger, & \text{ for even polynomial}
\end{array}
\right . \,\,,
$$
where $\mathrm{Poly(\sigma)}$ is a polynomial with a well-defined parity.
The QSVT routine allows to block-encode such singular value polynomial transforms, based on the Quantum Signal Processing (QSP) approach.
For the case of matrix inversion, we have the identities,
$$
\begin{aligned}
A^\dagger &=V^\dagger \Sigma W,\\
A^{-1} &= V^\dagger \Sigma^{-1}W,
\end{aligned}
$$
from which we can deduce that:
**Applying a QSVT of an odd polynomial that approximates $\mathrm{Poly(\sigma)}\sim 1/\sigma$, on the block-encoding of $A^{\dagger}$, gives the matrix inversion of $A$.**
Below we demonstrate how to implement a quantum linear solver for a specific problem.
The `BlockEncoding` class provides a high-level interface.
Given a block-encoding of the matrix $A$ and its effective condition number $(\sigma_{\min}/s)^{-1}$, we can easily obtain the inverse block-encoding using the `inverse` method, which handles the QSVT polynomial approximation internally.
## Example: Inverting a Random Matrix
```python theme={null}
!pip install -qq "classiq[qsp]"
```
#
## Setting an Example
We start by defining a specific problem: a matrix and its condition number. We take a simple use case: a random $4 \times 4$ matrix $A$ that we want to invert.
```python theme={null}
import numpy as np
# Generate a random matrix A
np.random.seed(4)
A_dim = 4 # 4x4 matrix (2 qubits)
A = np.random.rand(A_dim, A_dim)
# Normalize to ensure singular values are bounded
A = A / np.linalg.norm(A, ord=2)
print("Matrix A:")
print(A)
```
**Output:**
```
Matrix A:
[[0.39182404 0.2217292 0.39411516 0.28963128]
[0.28270785 0.08755579 0.3955698 0.00252439]
[0.10250415 0.17617013 0.3157927 0.08009863]
[0.34967018 0.39845723 0.06638609 0.24202955]]
```
Let's verify the matrix properties:
```python theme={null}
svd = np.linalg.svd(A)[1]
print(f"Singular values: {svd}")
print(f"Max singular value: {max(svd):.4f}")
print(f"Min singular value: {min(svd):.4f}")
print(f"Condition number: {max(svd)/min(svd):.4f}")
```
**Output:**
```
Singular values: [
1. 0.37943797 0.13596402 0.12337668]
Max singular value: 1.0000
Min singular value: 0.1234
Condition number: 8.1053
```
#
## Block-Encoding
The `BlockEncoding` class provides a convenient `from_matrix` method that automatically decomposes any matrix into a block encoding via Pauli decomposition and LCU (Linear Combination of Unitaries).
```python theme={null}
from classiq import *
from classiq.applications.block_encoding import BlockEncoding
# Create block encoding directly from matrix A
block_encoding = BlockEncoding.from_matrix(A)
print(f"Block encoding properties:")
print(f" Data size: {block_encoding.data_size} qubits")
print(f" Block size: {block_encoding.block_size} qubits")
print(f" Scaling factor (alpha): {block_encoding.alpha:.4f}")
```
**Output:**
```
Block encoding properties:
Data size: 2 qubits
Block size: 4 qubits
Scaling factor (alpha): 1.5894
```
The `from_matrix` method automatically:
* Decomposes the matrix into Pauli terms
* Constructs an LCU block encoding
* Computes the scaling factor $\alpha$ (Pauli 1-norm)
We evaluate the condition number of $A$, $\kappa = \sigma_{\max}/\sigma_{\min} = 1/\sigma_{\min}$ (using $\sigma_{\max}=1$ after normalization), and the *effective* condition number of the block-encoded operator $A/\alpha$, $\kappa_{\rm eff} = \alpha\,\kappa = \alpha/\sigma_{\min}$, which controls the QSVT polynomial degree. In real applications, this value needs to be approximated through prior knowledge of the problem or classical methods.
For our small example, we compute it directly.
```python theme={null}
# Bare condition number of A (sigma_max = 1 after normalization).
kappa = 1.0 / min(svd)
# Effective condition number of the block-encoded operator A / alpha.
kappa_eff = block_encoding.alpha * kappa
print(f"Condition number: {kappa:.4f}")
print(f"Effective condition number (of A / alpha): {kappa_eff:.4f}")
```
**Output:**
```
Condition number: 8.1053
Effective condition number (of A / alpha): 12.8823
```
#
## Building the Inverse by Hand with QSVT
To make the algorithm explicit, we first build the inverse block-encoding ourselves with the QSVT machinery, rather than calling the high-level `inverse` method (which we revisit at the end). We need a polynomial approximation of $1/x$.
Since $1/x$ diverges at $x=0$ it cannot be bounded by $1$ on $[-1,1]$, so we approximate it only on the range that contains the singular values of the block-encoded operator $A/\alpha$, namely $[\sigma_{\min}/\alpha, \sigma_{\max}/\alpha]$.
The target function is
$$
f(x) = \frac{\mathrm{SCALE}}{\kappa_{\rm eff}\, x}, \qquad x\in[\sigma_{\min}/\alpha, \sigma_{\max}/\alpha],
$$
with a factor $\mathrm{SCALE}<1$ ensuring $|f(x)|\le 1$ (which eases the QSVT angle finding).
According to Ref. \[2], a degree of $\sim \kappa_{\rm eff} \log(\kappa_{\rm eff}/\epsilon)$ suffices to approximate $f$ with error $\epsilon$. We use `qsp_approximate` to fit $f$ as a Chebyshev sum bounded by $1$ on $[-1,1]$, and `qsvt_phases` to obtain the corresponding QSVT phases.
```python theme={null}
from classiq.applications.qsp import qsp_approximate, qsvt_phases
EPS = 1e-4
SCALE = 0.25
degree = int(kappa_eff * np.log(kappa_eff / EPS))
if degree % 2 == 0: # inversion requires an odd polynomial
degree += 1
# Singular values of the block-encoded operator A / alpha.
be_svd = svd / block_encoding.alpha
def target_function(x):
return SCALE / (kappa_eff * x)
pcoefs, approx_error = qsp_approximate(
target_function,
degree=degree,
parity=1,
interval=[min(be_svd), max(be_svd)],
plot=True,
)
inversion_phases = qsvt_phases(pcoefs)
print(f"Polynomial degree: {degree}")
print(f"Max approximation error on the interval: {approx_error}")
```
**Output:**
```
/Users/roiedann/projects/classiq-library/.venv/lib/python3.11/site-packages/classiq/applications/qsp/qsp.py:203: UserWarning: Solution may be inaccurate.
Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
prob.solve()
```
**Output:**
```
Polynomial degree: 151
Max approximation error on the interval: 5.243390751308308e-07
```
#
## Solving a Linear Equation
Next, we solve an instance of a linear equation $A\vec{x} = \vec{b}$.
For this, we define some arbitrary $\vec{b}$ vector.
```python theme={null}
b = np.arange(A_dim)
b_norm = np.linalg.norm(b)
b_normalized = (b / b_norm).tolist()
print(b)
```
**Output:**
```
[0 1 2 3]
```
Now we build the quantum model. We prepare the initial state $|b\rangle$ and then call `qsvt_inversion`, driving it with the block encoding: the projector marks its block register being zero, its unitary is applied through `block_encoding.unitary`, and the extra QSVT auxiliary is the last block qubit.
```python theme={null}
@qfunc
def projector(block: QNum, res: QBit) -> None:
res ^= block == 0
@qfunc
def main(
data: Output[QNum[block_encoding.data_size]],
block: Output[QArray[QBit, block_encoding.block_size + 1]],
) -> None:
allocate(data)
allocate(block)
inplace_prepare_amplitudes(b_normalized, 0, data)
# QSVT inversion driven by the block encoding: reflect about its block being
# zero, apply its unitary, and use the last block qubit as the QSVT auxiliary.
qsvt_inversion(
phase_seq=inversion_phases,
block_encoding_cnot=lambda aux: projector(
block[0 : block_encoding.block_size], aux
),
u=lambda: block_encoding.unitary(data, block[0 : block_encoding.block_size]),
aux=block[block_encoding.block_size],
)
```
Let us synthesize, visualize, and execute.
For the execution, we choose a statevector simulator, as we are considering a small problem for demonstrating the algorithm. In addition, later on we would like to verify our result against the expected classical one.
```python theme={null}
qprog = synthesize(main)
show(qprog)
# Calculate state vector (filter in post-processing)
result = calculate_state_vector(qprog)
```
**Output:**
```
Submitting state-vector job to simulator
```
**Output:**
```
Quantum program link: https://platform.classiq.io/circuit/3K2wIOU7x80TMdTsjnpIy09uRfh
```
If our matrix $A$ is block-encoded with a block variable of size $m$, then the QSVT routine block encodes the inverse of $A$ on $m+1$ qubits, where the additional qubit is the QSVT auxiliary. Therefore, we need to post-select on `block` being zero, if we want to obtain a state $\sim A^{-1}|b\rangle$.
```python theme={null}
# Filter to block == all zeros (post-select on success)
zero_block = [0] * (block_encoding.block_size + 1)
df = result[result["block"].apply(lambda x: x == zero_block)].sort_values("data")
df
```
| | data | block | amplitude | magnitude | phase | probability | bitstring |
| - | - | - | - | - | - | - | - |
| 56 | 0 | \[0, 0, 0, 0, 0] | -0.004095+0.020587j | 0.02 | 0.56π | 0.000441 | 0000000 |
| 10 | 1 | \[0, 0, 0, 0, 0] | 0.024985-0.125607j | 0.13 | -0.44π | 0.016401 | 0000001 |
| 63 | 2 | \[0, 0, 0, 0, 0] | 0.001560-0.007844j | 0.01 | -0.44π | 0.000064 | 0000010 |
| 24 | 3 | \[0, 0, 0, 0, 0] | -0.015714+0.079001j | 0.08 | 0.56π | 0.006488 | 0000011 |
#
## Comparing to the Expected Result
Finally, we verify our quantum solution, by comparing it to the expected classical one.
```python theme={null}
expected_x = np.linalg.inv(A) @ b
```
After post-selecting the block on zero, the QSVT routine returns
$$
|x\rangle \propto p\!\left(A/\alpha\right)\frac{\vec b}{|\vec b|}, \qquad
p\!\left(A/\alpha\right) \approx \frac{\mathrm{SCALE}}{\kappa_{\rm eff}}\left(\frac{A}{\alpha}\right)^{-1} = \frac{\mathrm{SCALE}\,\alpha}{\kappa_{\rm eff}}\,A^{-1}.
$$
Collecting the prefactors — the $\alpha$ cancels since $\kappa_{\rm eff}=\alpha\kappa$ — the solution is recovered as
$$
\vec x = \frac{|\vec b|\,\kappa}{\mathrm{SCALE}}\,\mathrm{Re}\!\left(\text{amplitudes}\right).
$$
```python theme={null}
global_phase = np.angle(df.amplitude.iloc[0])
prefactor = b_norm * kappa / SCALE
computed_x = prefactor * np.real(np.array(df.amplitude) / np.exp(1j * global_phase))
```
Let's compare it with the expected solution:
*Comment: even after removing the global phase to ensure a real solution, the quantum solution can be obtained up to a sign.*
```python theme={null}
import matplotlib.pyplot as plt
plt.plot(expected_x, "o", label="classical")
plt.plot(computed_x, ".", label="quantum")
plt.xlabel(r"$i$", fontsize=16)
plt.ylabel(r"$x_i$", fontsize=16)
plt.legend();
```
```python theme={null}
assert (
min(
np.linalg.norm(computed_x - expected_x),
np.linalg.norm(-computed_x - expected_x),
)
< 0.05
)
```
## One-Line Inversion with `BlockEncoding.inverse`
The implementation above of the inverse (finding the QSVT phases and wiring up `qsvt_inversion`) is packaged in the `BlockEncoding.inverse` method. It takes the (bare) condition number and a target precision, forms the effective condition number $\kappa_{\rm eff}=\alpha\kappa$ internally, and returns a ready-to-use `BlockEncoding` of $A^{-1}$.
The same linear system is then solved with a single call, reproducing the by-hand result.
```python theme={null}
be_inverse = block_encoding.inverse(kappa=kappa, eps=1e-2)
# The entry point must be named `main`; this shadows the by-hand `main` above,
# but `qprog` was already synthesized from it and is unaffected.
@qfunc
def main(
data: Output[QNum[block_encoding.data_size]],
block: Output[QArray[QBit, be_inverse.block_size]],
) -> None:
allocate(data)
allocate(block)
inplace_prepare_amplitudes(b_normalized, 0, data)
be_inverse.unitary(data, block)
qprog_direct = synthesize(main)
result_direct = calculate_state_vector(qprog_direct)
zero_block_direct = [0] * be_inverse.block_size
df_direct = result_direct[
result_direct["block"].apply(lambda x: x == zero_block_direct)
].sort_values("data")
gp_direct = np.angle(df_direct.amplitude.iloc[0])
computed_x_direct = (
b_norm
* be_inverse.alpha
* np.real(np.array(df_direct.amplitude) / np.exp(1j * gp_direct))
)
print("by-hand :", np.round(computed_x, 3))
print("direct :", np.round(computed_x_direct, 3))
print("classical :", np.round(expected_x, 3))
```
**Output:**
```
Submitting state-vector job to simulator
```
**Output:**
```
by-hand : [ 2.546 -15.536 -0.97 9.771]
direct : [ 2.546 -15.536 -0.971 9.771]
classical : [-2.545 15.536 0.971 -9.771]
```
## References
\[1]: [Martyn JM, Rossi ZM, Tan AK, Chuang IL. Grand unification of quantum algorithms. PRX Quantum 2, 040203. (2021).](https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.040203)
\[2]: [A. M. Childs, R. Kothari, and R. D. Somma. Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. SIAM Journal on Computing 46, 1920-1950 (2017).](https://arxiv.org/abs/1511.02306)
## A Note on Efficient Block-Encoding, Condition Number, and Complexity
An efficient block-encoding of a matrix $A$ means that the resources to construct the unitary $U_{A,s}$, as well as the scaling factor $s$, scales poly-logarithmically with the matrix dimension.
This technical note concerns the latter point.
First, we note that an optimal value for the scaling factor $s$ is the maximal singular value $\sigma_{\max}$, namely $s\geq \sigma_{\max}$.
This comes directly from the requirement that $||A/s||\leq 1$.
For that reason, the effective condition number used throughout the demo, $\kappa = s/\sigma_{\min}$, is bounded from below by the actual condition number $\sigma_{\max}/\sigma_{\min}$.
Next, note that if $s$ scales linearly with the problem size, then the exponential speedup is lost.
This is because the polynomial degree then scales as $O(N)$. In addition, the amplitude of the encoded solution scales the same (although it also depends on the input vector $|b\rangle$), which requires an additional amplitude amplification routine that adds a factor of $O(\sqrt{N})$ to the overall resources.
This documentation is built and hosted on [Mintlify](https://mintlify.com), a developer documentation platform.