quantum-signal-processing

name	quantum-signal-processing
description	Basic Quantum Signal Processing (QSP) demo for single-qubit signal operators, phase-factor optimization, and scalar polynomial approximation such as cos(t*x). Use the Hamiltonian-simulation QSP-HS skill instead for block-encoded Hamiltonians and e^{-iHt}.

Quantum Signal Processing (QSP)

Purpose

Quantum Signal Processing implements a polynomial transformation $P(x)$ of a scalar signal $x$ using a sequence of single-qubit signal-processing rotations and phase operators. This implementation demonstrates basic QSP applied to function approximation: approximating $\cos(t \cdot x)$ at a test point.

Use this skill when you need to:

• Demonstrate the core QSP signal/phase sequence on one qubit.
• Approximate a scalar function such as $\cos(t x)$ at a point $x \in [-1,1]$.
• Inspect phase-factor optimization and QSP circuit structure.

Do not use this skill for QSP Hamiltonian simulation. If the request mentions a Hamiltonian matrix H, time evolution $e^{-iHt}$, block encoding, Chebyshev expansions of $\cos(tH)$ and $\sin(tH)$, or the QSPHSAlgorithm class, use ../../hamiltonian-simulation/qsp/SKILL.md.

Overview

QSP constructs a single-qubit circuit interleaving:

• Phase rotations $R_z(2\phi_k)$: parameterized by a phase sequence $\Phi = (\phi_0, \phi_1, \ldots, \phi_d)$.
• Signal operator $W(x) = R_x(2\arccos(x))$: encoding the signal $x$.

The circuit implements the $(0,0)$ matrix element of: $$U_\Phi(x) = e^{i\phi_0 Z}\prod_{k=1}^d [W(x) \cdot e^{i\phi_k Z}] \approx \begin{pmatrix} P(x) & \cdot \ \cdot & \cdot \end{pmatrix}$$

The phase sequence $\Phi$ is optimized (via L-BFGS-B) to approximate the target polynomial $P(x) \approx \cos(t \cdot x)$.

Prerequisites

• Single-qubit rotations: $R_z$, $R_x$.
• Polynomial approximation concepts (Chebyshev expansion).
• Python: numpy, scipy.optimize, Circuit, Register.

Using the Provided Implementation

from unitarylab_algorithms import QSPAlgorithm

algo = QSPAlgorithm(text_mode="legacy")
result = algo.run(
    t=1.0,         # Evolution time t: targets cos(t * x)
    d=10,          # Polynomial degree d
    x=0.5,         # Test signal point x ∈ [-1, 1]
    backend='torch'
)

print(result['Absolute error'])    # Absolute error |QSP(x) - cos(t*x)|
print(result['Estimated value'])   # QSP estimated value at x
print(result['Ideal value'])       # cos(t * x) at test point
print(result['circuit_path'])      # SVG circuit diagram
print(result['plot'])              # List of saved output file dicts

Core Parameters Explained

Parameter	Type	Default	Description
t	float	required	Evolution time. The target function is $\cos(t \cdot x)$.
d	int	required	Polynomial degree $d$. Higher degree achieves smaller error for complex targets.
x	float	0.5	Test point $x \in [-1, 1]$ at which to evaluate the approximation error.
backend	str	'torch'	Simulation backend.
device	str	'cpu'	Compute device (e.g. 'cpu', 'cuda').
dtype	type	np.complex128	Numeric dtype for simulation.

Common misunderstandings:

• x must be in $[-1, 1]$ (the domain of the signal operator $W(x)$). Values outside this range are clipped internally.
• d determines the number of signal $W(x)$ applications: d signal operators and d+1 phase rotations in the circuit.
• The phase sequence _find_phases() is computed via numerical optimization (L-BFGS-B), not analytically. The result may vary across runs due to random initialization.

Return Fields

Key	Type	Description
status	str	'ok' on success, 'failed' on failure.
Estimated value	complex	QSP-estimated value $\langle 0
Ideal value	float	Ideal value $\cos(t \cdot x)$ at test point x.
Absolute error	float	Absolute approximation error $
Computation time (s)	float	Wall-clock simulation execution time in seconds.
circuit	Circuit	The built QSP circuit object.
circuit_path	str	Path to saved SVG circuit diagram.
plot	list	List of saved output file dicts, each with 'format' (last 3 chars of filename) and 'filename' keys.

Implementation Architecture

QSPAlgorithm in algorithm.py organizes the algorithm into five stages with one classical optimizer helper.

run(t, d, x, backend, device, dtype) — Five Stages:

Stage	Code Action	Algorithmic Role
1 — Phase Optimization	_find_phases(t, d) — runs L-BFGS-B minimization	Computes the QSP phase sequence $\Phi$ numerically
2 — Circuit Construction	Creates Circuit(Register('q', 1)); applies qc.rz(2*phases[0], 0) as initial phase; loops d times applying qc.rx(2*theta, 0) (signal) then qc.rz(2*phases[k], 0) (phase rotation)	Builds the alternating signal/phase circuit
3 — Simulation	qc.execute(backend=backend, device=device, dtype=dtype) → final_state[0]	Runs statevector evolution of the single qubit
4 — Post-Processing	Compares qsp_val = complex(final_state[0]) to ideal_val = np.cos(t * x)	Computes absolute approximation error
5 — Export	self.save_circuit(qc) and self.save_txt()	Saves SVG circuit diagram and text result file

Helper Methods:

• _find_phases(t, d) — Fully classical numerical optimizer. Builds an internal matrix-product simulation of the QSP circuit for a grid of 2d+1 samples. Defines loss as average squared magnitude error vs. cos(t * x). Runs scipy.optimize.minimize with L-BFGS-B starting from random initial phases (np.random.randn(d+1) * 0.1). Returns (d+1) phase values.
• _build_return_dict(success, circuit_path, filepath, circuit) — Packages the result. Sets status='ok' if success=True, 'failed' otherwise. Wraps filepath into a list of {'format': filename[-3:], 'filename': filename} dicts stored under 'plot'. Merges self.output (containing 'Estimated value', 'Ideal value', 'Absolute error', 'Computation time (s)') into the returned dict.

Single-qubit circuit structure (Stage 2):

Rz(2φ₀) → Rx(2θ) → Rz(2φ₁) → Rx(2θ) → Rz(2φ₂) → ... → Rx(2θ) → Rz(2φ_d)

Total gates: 2d + 1. The theta = arccos(clip(x, -1, 1)) is a fixed parameter computed from x.

Important implementation note: The simulation in _find_phases uses explicit 2×2 matrix products (NumPy) — not the Circuit engine — so the optimization is fully classical and independent of the backend.

Data flow: (t, d) → _find_phases() → phases array → Circuit with Rz/Rx gates → execute() → final_state[0] → abs_error vs. ideal_val = np.cos(t * x) → _build_return_dict() → result dict.

Understanding the Key Quantum Components

$$W(x) = \begin{pmatrix} x & i\sqrt{1-x^2} \ i\sqrt{1-x^2} & x \end{pmatrix}$$ Implemented as $R_x(2\arccos(x))$ on the single qubit. This maps the scalar signal $x$ into a 2×2 rotation.

2. Phase Rotations $e^{i\phi_k Z}$

Implemented as $R_z(2\phi_k)$ gates. Interleaved with the signal operator, they provide the degrees of freedom to approximate arbitrary polynomials.

3. QSP Circuit Structure

Rz(2φ₀) → Rx(2θ) → Rz(2φ₁) → Rx(2θ) → ... → Rz(2φ_d)

where $\theta = \arccos(x)$. The $(0,0)$ matrix element of the product equals $P(x)$.

4. Phase Optimization

The phases $\Phi$ are found by minimizing: $$L(\Phi) = \frac{1}{\tilde{d}}\sum_{j=1}^{\tilde{d}} |\langle 0|U_\Phi(x_j)|0\rangle - e^{-i\tau x_j}|^2$$ over $2d+1$ Chebyshev sample points $x_j = \cos\left(\frac{(2j-1)\pi}{4\tilde{d}}\right)$.

Theory-to-Code Mapping

README / Theory Concept	Code Object or Location
Signal operator $W(x)$ = $R_x(2\theta)$	qc.rx(2 * theta, 0) where theta = arccos(x)
Phase rotations $e^{i\phi_k Z}$ = $R_z(2\phi_k)$	qc.rz(2 * phases[k], 0) in the main loop
Initial phase $\phi_0$	qc.rz(2 * phases[0], 0) before loop
Phase sequence $\Phi = (\phi_0, \ldots, \phi_d)$	phases returned by _find_phases() — (d+1) floats
Target function $\cos(t \cdot x)$	ideal_val = np.cos(t * x) in Stage 4
Loss function minimization	minimize(loss, ...) with method='L-BFGS-B' in _find_phases()
$(0,0)$ matrix element $\langle 0	U_\Phi(x)
Approximation error at test point	abs_error = float(np.abs(qsp_val - ideal_val))
Sample points $x_j = \cos((2j-1)\pi/(4\tilde d))$	x_samples = np.linspace(-1, 1, 2*d+1) (uniform, not Chebyshev nodes)

Notes on implementation vs. theory: The README describes exact Chebyshev sample points, but the code uses np.linspace(-1, 1, 2*d+1) (uniformly spaced) as the optimization grid. The loss function is identical in structure. The _find_phases uses 2×2 matrix products for efficiency, not the quantum simulator. The returned 'Absolute error' is evaluated only at x, not over the full interval — it is a spot check, not a worst-case bound.

Mathematical Deep Dive with $|P(x)| \leq 1$ for all $x \in [-1,1]$ and parity $d \bmod 2$, there exist phases $\Phi$ such that:

$$\langle 0|U_\Phi(x)|0\rangle = P(x)$$

Function approximation target: $P(x) = \cos(t \cdot x)$, approximated by a degree-$d$ truncated Chebyshev expansion. For $t$ fixed, error decreases exponentially in $d$.

Complexity: The circuit uses $d+1$ single-qubit gates plus $2d+1$ $R_z$ gates — total $O(d)$ gates.

Hands-On Example (UnitaryLab)

from unitarylab_algorithms import QSPAlgorithm
import numpy as np

algo = QSPAlgorithm(text_mode="legacy")

# Approximate cos(0.5 * x) at several test points
for x_test in [0.0, 0.3, 0.7, 1.0]:
    result = algo.run(t=0.5, d=8, x=x_test)
    ideal = np.cos(0.5 * x_test)
    print(f"x={x_test:.1f}: error={result['Absolute error']:.2e}, ideal={ideal:.4f}")

Reference Implementation (PennyLane)

The main implementation path in this project remains the UnitaryLab QSP implementation based on QSPAlgorithm, Circuit, and alternating Rz/Rx signal-processing rotations.

PennyLane is provided here only as a reference implementation for the more general Quantum Singular Value Transformation (QSVT) framework. Compared with the current UnitaryLab QSP example, PennyLane QSVT applies polynomial transformations to the singular values of a block-encoded matrix or operator.

Example A: Minimal PennyLane QSVT Run with Matrix Block Encoding

import pennylane as qml
import numpy as np

# Hermitian input matrix A
A = np.array([
    [0.1, 0.2],
    [0.2, -0.4],
])

# Polynomial P(x) = 0.2 + 0.3x^2
# Coefficients are ordered from lowest to highest power.
poly = np.array([0.2, 0.0, 0.3])

dev = qml.device("default.qubit", wires=[0, 1])

@qml.qnode(dev)
def circuit():
    qml.qsvt(
        A,
        poly,
        encoding_wires=[0, 1],
        block_encoding="embedding",
    )
    return qml.state()

result = circuit()
print(result)

Minimal Manual Implementation (UnitaryLab)

The following Python skeleton reconstructs the key components: phase optimization and the alternating QSP circuit.

Step 1: Simulate QSP circuit classically for optimization

# Simplified reconstruction — mirrors QSPAlgorithm._find_phases() internals
import numpy as np
from scipy.optimize import minimize

def qsp_matrix(x: float, phases: np.ndarray) -> np.ndarray:
    """Compute the QSP circuit product matrix for signal x and phase vector phases."""
    theta = np.arccos(np.clip(x, -1.0, 1.0))
    # Signal operator W(x) = Rx(2*arccos(x))
    W = np.array([[x, 1j*np.sqrt(1-x**2)],
                  [1j*np.sqrt(1-x**2), x]])
    # Initial phase rotation Rz(2*phi_0)
    mat = np.array([[np.exp(1j*phases[0]), 0], [0, np.exp(-1j*phases[0])]])
    for k in range(1, len(phases)):
        # Pauli Rz gate: e^{i phi_k Z}
        Rz_k = np.array([[np.exp(1j*phases[k]), 0], [0, np.exp(-1j*phases[k])]])
        mat = Rz_k @ W @ mat
    return mat

def find_qsp_phases(tau: float, d: int, seed: int = 42) -> np.ndarray:
    """Find QSP phases that approximate exp(-i*tau*x) on [-1,1]."""
    np.random.seed(seed)
    x_samples = np.linspace(-1, 1, 2*d + 1)
    target_vals = np.exp(-1j * tau * x_samples)

    def loss(phi):
        total = 0.0
        for xi, ti in zip(x_samples, target_vals):
            mat = qsp_matrix(xi, phi)
            total += abs(mat[0, 0] - ti) ** 2
        return total / len(x_samples)

    res = minimize(loss, x0=np.random.randn(d+1)*0.1, method='L-BFGS-B',
                   options={'maxiter': 400})
    return res.x

Step 2: Build and execute the QSP circuit

# Simplified reconstruction — mirrors QSPAlgorithm.run() Stage 2 and Stage 3
from unitarylab.core import Circuit, Register
import numpy as np

def build_qsp_circuit(phases: np.ndarray, x_value: float, backend: str = 'torch') -> Circuit:
    """Construct the QSP circuit: Rz(2φ₀) → [Rx(2θ) → Rz(2φₖ)] × d."""
    d = len(phases) - 1
    theta = float(np.arccos(np.clip(x_value, -1.0, 1.0)))
    qc = Circuit(Register('q', 1))
    # Initial phase rotation
    qc.rz(2 * phases[0], 0)
    for k in range(1, d + 1):
        qc.rx(2 * theta, 0)      # signal operator W(x) = Rx(2*arccos(x))
        qc.rz(2 * phases[k], 0)  # phase rotation
    return qc

# Full QSP pipeline
def qsp_approximate(t: float, d: int, x: float, backend: str = 'torch'):
    phases = find_qsp_phases(t, d)
    circ = build_qsp_circuit(phases, x, backend)
    state = circ.execute(backend=backend).state
    qsp_val = complex(state[0]) if hasattr(state, '__getitem__') else 0.0
    ideal = np.cos(t * x)
    return qsp_val, abs(qsp_val - ideal)

Debugging Tips

1. Large Absolute error for small d: $\cos(t \cdot x)$ requires degree $\sim O(t + \log(1/\epsilon))$ for $\epsilon$-approximation. Increase d.
2. Optimizer not converging: The L-BFGS-B optimizer uses default iteration limits. Modify _find_phases() to pass options={'maxiter': N} to minimize() if needed.
3. x outside $[-1, 1]$: Internally clipped by np.clip. Ensure your signal values are in the valid domain.
4. Error non-zero at x=0: $\cos(t \cdot 0) = 1$; if error is large here with small t, the optimizer may have gotten stuck. Re-run to get a different random initialization.
5. Random phase initialization: _find_phases() uses np.random.randn(d+1)*0.1 as starting phases. Results may vary across runs; simply re-run algo.run() to try a different initialization.