BlockEncoding.pseudo_inv#

BlockEncoding.pseudo_inv(eps: float, theta: float, delta: float | None = None) → BlockEncoding#

Returns a BlockEncoding approximating the threshold matrix pseudoinverse of the operator.

Approximates the pseudoinverse of a matrix by ignoring singular values below a specified threshold. This regularizes ill-conditioned linear systems where tiny singular values would otherwise amplify noise and cause large numerical errors. By discarding these negligible values, this method trades a small amount of exact accuracy for a significantly more stable and reliable solution.

Let \(A\) be a matrix with Singular Value Decomposition

\[A = U\Sigma V^{\dagger} = \sum_i \sigma_i u_iv_i^{\dagger}\]

where \(\sigma_i\) are the singular values, \(u_i\) are the left singular vectors, and \(v_i\) are the right singular vectors.

Given a threshold \(\theta>0\), define a scalar thresholding function \(f_{\theta}(\sigma)\) that acts on the singular values:

\[\begin{split}f_{\theta}(\sigma)=\begin{cases} \frac{1}{\sigma} & \sigma\geq\theta \\ 0 & \sigma<\theta \end{cases}\end{split}\]

The threshold pseudo inverse \(A_{\theta}^{+}\) is constructed by applying this function to the singular values and recombining the matrix:

\[A_{\theta}^{+} = V f_{\theta}(\Sigma)U^{\dagger} = \sum_{\sigma_i\geq\theta}\frac{1}{\sigma_i}v_iu_i^{\dagger}\]

For a block-encoded matrix \(A\) with normalization factor \(\alpha\), this function returns a BlockEncoding of an operator \(\tilde{A}_{\alpha\cdot\theta}^{+}\) such that \(\|\tilde{A}_{\alpha\cdot\theta}^{+} - A_{\alpha\cdot\theta}^{+}\| \leq \epsilon\).

The pseudo inverse is implemented via Quantum Singular Value Transform (QSVT) using a polynomial approximation of \(1/x\) over the domain \(D_{\theta} = [-1, -\theta] \cup [\theta, 1]\), and a smoothed rectangle filter over the domain \(D_{\theta}' = [-\theta, \theta]\).

../../../_images/chebyshev_pseudo_inversion.png

Parameters:

epsfloat: The target precision \(\epsilon\).
thetafloat: This threshold value defines the boundaries of the “gap” around zero \([-\theta, \theta]\subset [-1,1]\) where the function \(1/x\) is not approximated.
deltafloat, optional: The width of the transition region. The function will smoothly decay from 1 to 0 over the intervals \([-\theta, -\theta + 2\delta]\) and \([\theta - 2\delta, \theta]\). Defaults to \(\delta = \theta/4\).

Returns:

BlockEncoding: A new BlockEncoding instance representing an approximation of the pseudo inverse \(A_{\alpha\cdot\theta}^{+}\).

Notes

Complexity: The polynomial degree scales as \(\mathcal{O}(\log(1/(\epsilon\theta))/\theta)\).
It is assumed that the relevant singular values of \(A/\alpha\) lie within \(D_{\theta}\).

References

Childs et al. (2017) Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.
Gilyén et al. (2019) Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.

Examples

First, define a matrix \(A\) and a right-hand side vector \(\vec{b}\).

import numpy as np

N = 4
A = 1.4 * np.eye(N, k=1) + 1.1 * np.eye(N, k=-1)
A[0, N-1] = 1.1
A[N-1, 0] = 1.4

b = np.array([0, 1, 1, 1])

print(A)
#[[0.  1.4 0.  1.1]
# [1.1 0.  1.4 0. ]
# [0.  1.1 0.  1.4]
# [1.4 0.  1.1 0. ]]

_, S, _ = np.linalg.svd(A)
print("Singular values of A: ", S)
# Singular values of A:  [2.5 2.5 0.3 0.3]

Generate a block-encoding of \(A\) and use pseudo_inv() to find a block-encoding approximating \(A_{1}^{+}\).

from qrisp import QuantumFloat, prepare, terminal_sampling
from qrisp.block_encodings import BlockEncoding

# Define efficient block-encoding.
def f0(x): x+=1
def f1(x): x-=1
BE = BlockEncoding.from_lcu(np.array([1.4, 1.1]), [f0, f1])
# alpha = 1.4 + 1.1 = 2.5

# Choose threshold theta > 0.3 / 2.5 = 0.12 
# to cut off smallest singular values.
BE_inv = BE.pseudo_inv(eps=0.01, theta=0.4, delta=0.1)

# Prepares operand variable in state |b>
def prep_b():
    operand = QuantumFloat(2)
    prepare(operand, b)
    return operand

@terminal_sampling
def main():
    operand = prep_b()
    ancillas = BE_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_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))])
print(amps)

Finally, compare the quantum simulation result with the classical solution:

def threshold_pseudoinverse(A, threshold):
    U, S, Vh = np.linalg.svd(A, full_matrices=False)
    S_inv = np.zeros_like(S)
    valid_mask = S >= threshold
    S_inv[valid_mask] = 1.0 / S[valid_mask]
    return (Vh.conj().T * S_inv) @ U.conj().T

c = (threshold_pseudoinverse(A, 0.4 * 2.5) @ b) 
c = c / np.linalg.norm(c)

print("QUANTUM SIMULATION\n", amps, "\nCLASSICAL SOLUTION\n", c)
# QUANTUM SIMULATION
# [0.63244232 0.31179432 0.63244232 0.32065201] 
# CLASSICAL SOLUTION
# [0.63245553 0.31622777 0.63245553 0.31622777]