Skip to content

Completion Solvers

qsp_proc.completion

Completion step solvers: Fejér-Wilson, BerSün, optimization, rootfinding.

list_methods()

Return supported completion method names.

Source code in src\qsp_proc\completion\__init__.py 
78
79
80
def list_methods() -> list[str]:
    """Return supported completion method names."""
    return list(_SOLVERS)

complete(poly, method='bersun', **kwargs)

Compute complementary completion Q for input Laurent polynomial P.

Dispatches to the selected completion backend while enforcing the precondition ||P||_∞ <= 1 on U(1) via dense unit-circle sampling.

Parameters:

Name Type Description Default
poly LaurentPoly

Input Laurent polynomial P.

 required
method str

Completion backend name. See :func:list_methods.

 'bersun'
**kwargs Any

Forwarded directly to the selected solver.

 {}

Returns:

Type Description
LaurentPoly

Complementary Laurent polynomial Q.

Raises:

Type Description
TypeError

If poly is not a :class:~qsp_proc.polynomials.laurent.LaurentPoly.
ValueError

If method is unknown or sampled ||P||_∞ > 1.
RuntimeError

If the method is not implemented for Laurent inputs or the solver fails.
Source code in src\qsp_proc\completion\__init__.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
def complete(poly: LaurentPoly, method: str = "bersun", **kwargs: Any) -> LaurentPoly:
    """Compute complementary completion ``Q`` for input Laurent polynomial ``P``.

    Dispatches to the selected completion backend while enforcing the precondition
    ``||P||_∞ <= 1`` on ``U(1)`` via dense unit-circle sampling.

    Args:
        poly: Input Laurent polynomial ``P``.
        method: Completion backend name. See :func:`list_methods`.
        **kwargs: Forwarded directly to the selected solver.

    Returns:
        Complementary Laurent polynomial ``Q``.

    Raises:
        TypeError: If ``poly`` is not a :class:`~qsp_proc.polynomials.laurent.LaurentPoly`.
        ValueError: If ``method`` is unknown or sampled ``||P||_∞ > 1``.
        RuntimeError: If the method is not implemented for Laurent inputs or the solver fails.
    """
    if not isinstance(poly, LaurentPoly):
        raise TypeError("complete expects a LaurentPoly input.")

    linf = _linf_on_unit_circle(poly)
    if linf > 1.0 + _NORM_TOL:
        raise ValueError(
            f"Completion requires ||P||_∞ <= 1 on U(1); estimated ||P||_∞ = {linf:.6e}."
        )

    key = _normalize_method(method)
    solver = _SOLVERS.get(key)
    if solver is None:
        supported = ", ".join(_SOLVERS)
        raise ValueError(
            f"Unknown completion method {method!r}. Supported methods: {supported}."
        )

    if key == "wilson" and linf > _WILSON_RECOMMENDED_LINF + _NORM_TOL:
        warnings.warn(
            f"Wilson method works best with ||P||_∞ <= 1/sqrt(2); got {linf:.4f}. "
            "Consider using method='bersun' or subnormalizing P.",
            stacklevel=2,
        )

    return solver(poly, **kwargs)

qsp_proc.completion.bersun

BerSün FFT completion (complementary Laurent Q from P).

See Berntson & Sunderhauf arXiv:2406.04246 and Skelton arXiv:2501.05977 §5.1.4.

laurent_at_roots_of_unity_ifft(poly, num_roots=None, *, N=1)

P(ω^k) for ω = exp(2πi/n), with n from num_roots or N * support_width.

Zero-padded IFFT on causal coeffs, times ω^(k * min_degree) for the Laurent shift.

Source code in src\qsp_proc\completion\bersun.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def laurent_at_roots_of_unity_ifft(
    poly: LaurentPoly,
    num_roots: int | None = None,
    *,
    N: int = 1,
) -> np.ndarray:
    """P(ω^k) for ω = exp(2πi/n), with n from ``num_roots`` or ``N * support_width``.

    Zero-padded IFFT on causal coeffs, times ω^(k * min_degree) for the Laurent shift.
    """
    n = _grid_length(poly, num_roots, N, "num_roots")
    coeffs = poly.coeffs.astype(np.complex128, copy=False)
    padded = np.zeros(n, dtype=np.complex128)
    padded[: coeffs.size] = coeffs

    omega = np.exp(2j * np.pi / n)
    k = np.arange(n, dtype=np.float64)
    return (omega ** (k * poly.min_degree)) * (n * np.fft.ifft(padded))

downscale_p_for_bersun(poly, epsilon)

Scale P by (1 − ε/4) as in Algorithm 2 when ‖P‖∞,𝕋 = 1 (arXiv:2406.04246).

Source code in src\qsp_proc\completion\bersun.py 
68
69
70
71
72
73
74
75
def downscale_p_for_bersun(poly: LaurentPoly, epsilon: float) -> LaurentPoly:
    """Scale P by (1 − ε/4) as in Algorithm 2 when ‖P‖∞,𝕋 = 1 (arXiv:2406.04246)."""
    if epsilon <= 0:
        raise ValueError("epsilon must be positive.")
    if epsilon >= 4:
        raise ValueError("epsilon must be < 4 so that 1 - epsilon/4 remains positive.")
    scale = 1.0 - epsilon / 4.0
    return LaurentPoly(poly.coeffs * scale, poly.min_degree)

one_minus_modulus_squared_laurent(poly)

Laurent polynomial S with S(z) = 1 - |P(z)|² on |z|=1 (via P(z)conj(P(1/z))).

Implemented as 1 minus conv(p, conj(p)[::-1]) on the causal coefficient vector p.

Source code in src\qsp_proc\completion\bersun.py 
78
79
80
81
82
83
84
85
86
87
88
def one_minus_modulus_squared_laurent(poly: LaurentPoly) -> LaurentPoly:
    """Laurent polynomial S with S(z) = 1 - |P(z)|² on |z|=1 (via P(z)conj(P(1/z))).

    Implemented as 1 minus conv(p, conj(p)[::-1]) on the causal coefficient vector p.
    """
    p = poly.coeffs.astype(np.complex128, copy=False)
    w = p.size
    autocorr = np.convolve(p, np.conj(p)[::-1])
    s = -autocorr
    s[w - 1] += 1.0
    return LaurentPoly(s, min_degree=-(w - 1))

fourier_coeffs_log_one_minus_modulus_squared(poly, num_samples=None, *, N=1, epsilon=0.0)

FFT coefficients of log(1-|P|²) on the n-th roots: fft(g)/n (numpy order).

Uses the same IFFT root evaluation as Algorithm 1 Step 1 (laurent_at_roots_of_unity_ifft), not direct poly(z) evaluation.

Source code in src\qsp_proc\completion\bersun.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def fourier_coeffs_log_one_minus_modulus_squared(
    poly: LaurentPoly,
    num_samples: int | None = None,
    *,
    N: int = 1,
    epsilon: float = 0.0,
) -> np.ndarray:
    """FFT coefficients of log(1-|P|²) on the n-th roots: fft(g)/n (numpy order).

    Uses the same IFFT root evaluation as Algorithm 1 Step 1 (``laurent_at_roots_of_unity_ifft``),
    not direct ``poly(z)`` evaluation.
    """
    if epsilon < 0:
        raise ValueError("epsilon must be non-negative.")

    n = _grid_length(poly, num_samples, N, "num_samples")
    p_vals = laurent_at_roots_of_unity_ifft(poly, num_roots=n)
    defect = 1.0 - np.abs(p_vals) ** 2
    if epsilon > 0:
        defect = np.maximum(defect, epsilon)
    defect_real = np.real(defect).astype(np.float64, copy=False)
    if epsilon == 0.0 and np.any(defect_real <= 0):
        raise ValueError(
            "1 - |P|² has non-positive values on the grid; "
            "pass epsilon > 0 or subnormalize P."
        )
    g = np.log(defect_real)
    return np.fft.fft(g) / n

apply_causal_projection_fourier(fourier_coeffs)

Causal / Hardy projector Π on fft(g)/n: halve DC, zero negative wrapped bins.

Source code in src\qsp_proc\completion\bersun.py 
121
122
123
124
125
126
127
128
129
130
131
def apply_causal_projection_fourier(fourier_coeffs: np.ndarray) -> np.ndarray:
    """Causal / Hardy projector Π on fft(g)/n: halve DC, zero negative wrapped bins."""
    a = np.asarray(fourier_coeffs)
    if a.ndim != 1:
        raise ValueError("fourier_coeffs must be a 1D array.")
    out = a.astype(np.complex128, copy=True)
    out[0] *= 0.5
    n = out.size
    if n > 1:
        out[n // 2 + 1 :] = 0.0
    return out

complementary_q_on_roots_from_projected_log(projected_log_fourier)

Q(ω^k) ≈ exp(n * ifft(Π ĝ)) with ĝ = fft(g)/n from the log-defect step.

Source code in src\qsp_proc\completion\bersun.py 
134
135
136
137
138
139
140
def complementary_q_on_roots_from_projected_log(projected_log_fourier: np.ndarray) -> np.ndarray:
    """Q(ω^k) ≈ exp(n * ifft(Π ĝ)) with ĝ = fft(g)/n from the log-defect step."""
    a = np.asarray(projected_log_fourier)
    if a.ndim != 1:
        raise ValueError("projected_log_fourier must be a 1D array.")
    n = a.size
    return np.exp(n * np.fft.ifft(a.astype(np.complex128, copy=False)))

laurent_q_from_unit_root_samples(q_on_roots, *, trim_tol=0.0, target_degree=None)

Causal Laurent Q from samples Q(ω^k): coefficients are fft(samples)/n.

Optional target_degree caps the recovered degree (Algorithm 1, Step 6).

Source code in src\qsp_proc\completion\bersun.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
def laurent_q_from_unit_root_samples(
    q_on_roots: np.ndarray,
    *,
    trim_tol: float = 0.0,
    target_degree: int | None = None,
) -> LaurentPoly:
    """Causal Laurent Q from samples Q(ω^k): coefficients are fft(samples)/n.

    Optional ``target_degree`` caps the recovered degree (Algorithm 1, Step 6).
    """
    a = np.asarray(q_on_roots)
    if a.ndim != 1:
        raise ValueError("q_on_roots must be a 1D array.")
    if a.size < 1:
        raise ValueError("q_on_roots must be non-empty.")
    if target_degree is not None and target_degree < 0:
        raise ValueError("target_degree must be non-negative.")
    n = a.size
    coeffs = np.fft.fft(a.astype(np.complex128, copy=False)) / n
    if target_degree is not None:
        coeffs = coeffs[: min(target_degree + 1, coeffs.size)]
    return LaurentPoly(coeffs, min_degree=0, trim_tol=trim_tol)

bersun_complementary_laurent(poly, num_samples=None, *, N=1, epsilon=0.0, trim_tol=0.0, target_degree=None)

Run BerSün Algorithm 1: complementary causal Q from Laurent P.

By default, truncates recovered coefficients to degree max_degree - min_degree (coefficient span), matching Step 6 for the usual P support.

Source code in src\qsp_proc\completion\bersun.py 
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def bersun_complementary_laurent(
    poly: LaurentPoly,
    num_samples: int | None = None,
    *,
    N: int = 1,
    epsilon: float = 0.0,
    trim_tol: float = 0.0,
    target_degree: int | None = None,
) -> LaurentPoly:
    """Run BerSün Algorithm 1: complementary causal Q from Laurent P.

    By default, truncates recovered coefficients to degree ``max_degree - min_degree``
    (coefficient span), matching Step 6 for the usual P support.
    """
    cap = (
        (poly.max_degree - poly.min_degree)
        if target_degree is None
        else target_degree
    )
    hat = fourier_coeffs_log_one_minus_modulus_squared(
        poly, num_samples, N=N, epsilon=epsilon
    )
    proj = apply_causal_projection_fourier(hat)
    on_roots = complementary_q_on_roots_from_projected_log(proj)
    return laurent_q_from_unit_root_samples(
        on_roots, trim_tol=trim_tol, target_degree=cap
    )

qsp_proc.completion.wilson

Fejér-Wilson completion method (Section 4).

fejer_problem_from_poly(p)

Construct the Fejér-factorization input Laurent polynomial.

Builds :math:F(z) = 1 - P(z)P(1/z)^*, which satisfies :math:F(e^{i\theta}) = 1 - |P(e^{i\theta})|^2 on :math:U(1).

Parameters:

Name Type Description Default
p LaurentPoly

Input Laurent polynomial :math:P.

 required

Returns:

Type Description
LaurentPoly

Self-adjoint Laurent polynomial :math:F.
Source code in src\qsp_proc\completion\wilson.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def fejer_problem_from_poly(p: LaurentPoly) -> LaurentPoly:
    """Construct the Fejér-factorization input Laurent polynomial.

    Builds
    :math:`F(z) = 1 - P(z)P(1/z)^*`, which satisfies
    :math:`F(e^{i\\theta}) = 1 - |P(e^{i\\theta})|^2` on :math:`U(1)`.

    Args:
        p: Input Laurent polynomial :math:`P`.

    Returns:
        Self-adjoint Laurent polynomial :math:`F`.
    """
    f_coeffs = _f_coeffs_one_minus_modulus_sq(p.coeffs)
    f_min = p.min_degree - p.max_degree
    f_poly = LaurentPoly(f_coeffs, min_degree=f_min)

    num_samples = max(128, 8 * (p.max_degree - p.min_degree + 1))
    theta = np.linspace(0.0, 2.0 * np.pi, num_samples, endpoint=False, dtype=np.float64)
    f_vals = f_poly(np.exp(1j * theta))
    min_real = float(np.min(np.real(f_vals)))
    if min_real < -1e-10:
        warnings.warn(
            (
                "fejer_problem_from_poly: sampled F(e^{iθ}) has negative values; "
                f"minimum real part = {min_real:.3e}. "
                "Wilson completion may fail or require a scaled input polynomial."
            ),
            RuntimeWarning,
            stacklevel=2,
        )

    return f_poly

wilson_iterate(f_coeffs, gamma_current)

One Wilson (Newton on :math:U(1)) step for Fejér spectral factorization.

At roots :math:z_k, uses :math:\gamma \leftarrow \tfrac12(\gamma + F/\overline{\gamma}), then recovers causal coefficients via fft(samples)/L truncated to length n+1. This is the FFT-based equivalent of the Toeplitz linear-system update in Skelton (2025), Eq. 42, and runs in :math:O(n \log n) per iteration.

Parameters:

Name Type Description Default
f_coeffs ndarray

Laurent coefficients of :math:F, degrees -n … n (length 2n+1).

 required
gamma_current ndarray

Causal coefficients of the current :math:\gamma (length n+1).

 required

Returns:

Type Description
ndarray

Updated causal coefficients (length n+1).
Source code in src\qsp_proc\completion\wilson.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def wilson_iterate(f_coeffs: np.ndarray, gamma_current: np.ndarray) -> np.ndarray:
    """One Wilson (Newton on :math:`U(1)`) step for Fejér spectral factorization.

    At roots :math:`z_k`, uses
    :math:`\\gamma \\leftarrow \\tfrac12(\\gamma + F/\\overline{\\gamma})`, then recovers
    causal coefficients via ``fft(samples)/L`` truncated to length ``n+1``.
    This is the FFT-based equivalent of the Toeplitz linear-system update in
    Skelton (2025), Eq. 42, and runs in :math:`O(n \\log n)` per iteration.

    Args:
        f_coeffs: Laurent coefficients of :math:`F`, degrees ``-n`` … ``n`` (length ``2n+1``).
        gamma_current: Causal coefficients of the current :math:`\\gamma` (length ``n+1``).

    Returns:
        Updated causal coefficients (length ``n+1``).
    """
    f = np.asarray(f_coeffs, dtype=np.complex128)
    gamma = np.asarray(gamma_current, dtype=np.complex128)

    if f.ndim != 1 or gamma.ndim != 1:
        raise ValueError("f_coeffs and gamma_current must be 1D arrays.")
    if gamma.size < 1:
        raise ValueError("gamma_current must be non-empty.")

    n = gamma.size - 1
    if f.size != 2 * n + 1:
        raise ValueError(
            f"Expected len(f_coeffs) == 2*len(gamma_current)-1 ({2 * n + 1}), got {f.size}."
        )

    fft_len = _fft_len_for_degree(n)
    f_vals = _laurent_on_unit_roots(f, min_degree=-n, fft_len=fft_len)
    g_vals = _causal_on_unit_roots(gamma, fft_len)
    denom = np.conj(g_vals)

    if np.any(np.abs(denom) < 1e-14):
        raise ValueError(
            "gamma_current has near-zeros on the FFT grid; F/γ* is unstable."
        )

    next_on_roots = 0.5 * (g_vals + f_vals / denom)
    coeffs = np.fft.fft(next_on_roots) / fft_len
    return np.asarray(coeffs[: n + 1], dtype=np.complex128)

wilson_completion(poly, max_iter=50, tol=1e-13, initial_guess=None)

Run Wilson iteration to compute a causal complementary factor :math:\gamma.

Finds causal :math:\gamma with :math:F(z)=1-|P(z)|^2 \approx |\gamma(z)|^2 on :math:U(1) (equivalently the Fejér Laurent factorization).

Parameters:

Name Type Description Default
poly LaurentPoly

Input Laurent polynomial :math:P.

 required
max_iter int

Maximum Wilson iterations.

 50
tol float

Stop when max coefficient change falls below this.

 1e-13
initial_guess ndarray | None

Optional length-n+1 starting vector for :math:\gamma.

 None

Returns:

Type Description
LaurentPoly

Completed causal Laurent polynomial :math:\gamma.

Raises:

Type Description
ValueError

Invalid arguments or unusable default initial guess.
RuntimeError

No convergence or post-validation failure.
Source code in src\qsp_proc\completion\wilson.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
def wilson_completion(
    poly: LaurentPoly,
    max_iter: int = 50,
    tol: float = 1e-13,
    initial_guess: np.ndarray | None = None,
) -> LaurentPoly:
    """Run Wilson iteration to compute a causal complementary factor :math:`\\gamma`.

    Finds causal :math:`\\gamma` with :math:`F(z)=1-|P(z)|^2 \\approx |\\gamma(z)|^2` on
    :math:`U(1)` (equivalently the Fejér Laurent factorization).

    Args:
        poly: Input Laurent polynomial :math:`P`.
        max_iter: Maximum Wilson iterations.
        tol: Stop when max coefficient change falls below this.
        initial_guess: Optional length-``n+1`` starting vector for :math:`\\gamma`.

    Returns:
        Completed causal Laurent polynomial :math:`\\gamma`.

    Raises:
        ValueError: Invalid arguments or unusable default initial guess.
        RuntimeError: No convergence or post-validation failure.
    """
    if max_iter < 1:
        raise ValueError("max_iter must be >= 1.")
    if tol <= 0:
        raise ValueError("tol must be positive.")

    f_coeffs = _f_coeffs_one_minus_modulus_sq(poly.coeffs)
    n = poly.coeffs.size - 1

    if initial_guess is None:
        f0 = np.real_if_close(f_coeffs[n]).item()
        if np.real(f0) < 0 and np.real(f0) > -tol:
            f0 = 0.0
        if np.real(f0) < 0:
            raise ValueError(
                "Cannot form default Wilson initial guess: DC coefficient F_0 is negative."
            )
        gamma = np.zeros(n + 1, dtype=np.complex128)
        gamma[0] = np.sqrt(np.real(f0))
    else:
        gamma = np.asarray(initial_guess, dtype=np.complex128)
        if gamma.ndim != 1:
            raise ValueError("initial_guess must be a 1D array.")
        if gamma.size != n + 1:
            raise ValueError(
                f"initial_guess must have length n+1={n + 1}, got {gamma.size}."
            )

    last_err = np.inf
    converged = False
    for it in range(1, max_iter + 1):
        gamma_next = wilson_iterate(f_coeffs, gamma)
        last_err = float(np.max(np.abs(gamma_next - gamma)))
        logger.info("Wilson iteration %d: l_inf update error = %.3e", it, last_err)
        gamma = gamma_next
        if last_err < tol:
            converged = True
            break

    if not converged:
        raise RuntimeError(
            f"Wilson iteration did not converge in {max_iter} steps; "
            f"last l_inf error = {last_err:.3e}."
        )

    fft_len = _fft_len_for_degree(n)
    f_vals = _laurent_on_unit_roots(f_coeffs, min_degree=-n, fft_len=fft_len)
    g_vals = _causal_on_unit_roots(gamma, fft_len)
    validation_err = float(np.max(np.abs(f_vals - np.abs(g_vals) ** 2)))
    validation_threshold = max(100.0 * tol, 1e-10)
    logger.info("Wilson validation error on sampled unit circle: %.3e", validation_err)
    if validation_err >= validation_threshold:
        raise RuntimeError(
            "Wilson completion failed validation: "
            f"||F - |γ|²||_inf = {validation_err:.3e} >= {validation_threshold:.3e}."
        )

    return LaurentPoly(gamma, min_degree=0)

qsp_proc.completion.optimization

PyTorch-based optimization for completion (Section 5.2).

build_qsp_unitary(phi, x)

Build batched QSP unitaries from phase sequence and sample points.

For each sample :math:x_j, constructs :math:U_\\Phi(x_j) = R(\\phi_0)\\prod_{k=1}^{d}[W(x_j)R(\\phi_k)], where

:math:W(x)=\\begin{bmatrix}x & i\\sqrt{1-x^2}\\\\ i\\sqrt{1-x^2} & x\\end{bmatrix}, :math:R(\\phi)=\\mathrm{diag}(e^{i\\phi}, e^{-i\\phi}).

Parameters:

Name Type Description Default
phi Tensor

Phase vector of length d+1 as a 1D tensor.

 required
x Tensor

Real sample points of length m as a 1D tensor.

 required

Returns:

Type Description
Tensor

Batched 2x2 unitaries with shape (m, 2, 2) and dtype torch.complex128.

Raises:

Type Description
ValueError

If phi or x is not 1D, or if phi is empty.
Source code in src\qsp_proc\completion\optimization.py 
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
def build_qsp_unitary(phi: torch.Tensor, x: torch.Tensor) -> torch.Tensor:
    r"""Build batched QSP unitaries from phase sequence and sample points.

    For each sample :math:`x_j`, constructs
    :math:`U_\\Phi(x_j) = R(\\phi_0)\\prod_{k=1}^{d}[W(x_j)R(\\phi_k)]`, where

    :math:`W(x)=\\begin{bmatrix}x & i\\sqrt{1-x^2}\\\\ i\\sqrt{1-x^2} & x\\end{bmatrix}`,
    :math:`R(\\phi)=\\mathrm{diag}(e^{i\\phi}, e^{-i\\phi})`.

    Args:
        phi: Phase vector of length ``d+1`` as a 1D tensor.
        x: Real sample points of length ``m`` as a 1D tensor.

    Returns:
        Batched 2x2 unitaries with shape ``(m, 2, 2)`` and dtype ``torch.complex128``.

    Raises:
        ValueError: If ``phi`` or ``x`` is not 1D, or if ``phi`` is empty.
    """
    if phi.ndim != 1:
        raise ValueError("phi must be a 1D tensor.")
    if x.ndim != 1:
        raise ValueError("x must be a 1D tensor.")
    if phi.numel() < 1:
        raise ValueError("phi must have length at least 1.")

    device = phi.device
    phi_f = phi.to(device=device, dtype=torch.float64)
    x_f = x.to(device=device, dtype=torch.float64)
    m = x_f.numel()

    phi_c = phi_f.to(torch.complex128)
    phase_pos = torch.exp(1j * phi_c)
    phase_neg = torch.exp(-1j * phi_c)

    x_c = x_f.to(torch.complex128)
    root = torch.sqrt(torch.clamp(1.0 - x_f * x_f, min=0.0)).to(torch.complex128)
    iroot = 1j * root
    # W(x) for all samples, shape (m, 2, 2).
    w = torch.stack(
        (torch.stack((x_c, iroot), dim=-1), torch.stack((iroot, x_c), dim=-1)),
        dim=-2,
    )

    u = _qsp_rotation(phase_pos[0], phase_neg[0]).unsqueeze(0).expand(m, -1, -1)
    for j in range(1, phi_f.numel()):
        rj = _qsp_rotation(phase_pos[j], phase_neg[j])
        u = torch.matmul(u, torch.matmul(w, rj))

    return u

symmetric_qsp_loss(phi, target_fn, chebyshev_nodes, parity)

Compute the symmetric-QSP least-squares objective (Skelton 2025, Eq. 51-52).

Reconstructs the full symmetric phase sequence :math:\\Phi from reduced phases :math:\\hat\\Phi, evaluates :math:U_\\Phi(x_j) on collocation nodes, and returns mean-squared error between :math:\\Re\\,U_{00}(x_j) and target values.

Parameters:

Name Type Description Default
phi Tensor

Reduced phase vector :math:\\hat\\Phi (1D).

 required
target_fn Callable[..., object]

Target function f to evaluate at collocation points.

 required
chebyshev_nodes Tensor

Collocation nodes (1D).

 required
parity int

Degree parity indicator (even/odd) used for symmetry reconstruction.

 required

Returns:

Type Description
Tensor

Scalar MSE loss tensor.

Raises:

Type Description
ValueError

If inputs have invalid shapes or unsupported parity.
Source code in src\qsp_proc\completion\optimization.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
def symmetric_qsp_loss(
    phi: torch.Tensor,
    target_fn: Callable[..., object],
    chebyshev_nodes: torch.Tensor,
    parity: int,
) -> torch.Tensor:
    r"""Compute the symmetric-QSP least-squares objective (Skelton 2025, Eq. 51-52).

    Reconstructs the full symmetric phase sequence :math:`\\Phi` from reduced phases
    :math:`\\hat\\Phi`, evaluates :math:`U_\\Phi(x_j)` on collocation nodes, and returns
    mean-squared error between :math:`\\Re\\,U_{00}(x_j)` and target values.

    Args:
        phi: Reduced phase vector :math:`\\hat\\Phi` (1D).
        target_fn: Target function ``f`` to evaluate at collocation points.
        chebyshev_nodes: Collocation nodes (1D).
        parity: Degree parity indicator (even/odd) used for symmetry reconstruction.

    Returns:
        Scalar MSE loss tensor.

    Raises:
        ValueError: If inputs have invalid shapes or unsupported parity.
    """
    if phi.ndim != 1:
        raise ValueError("phi must be a 1D tensor.")
    if chebyshev_nodes.ndim != 1:
        raise ValueError("chebyshev_nodes must be a 1D tensor.")
    if phi.numel() < 1:
        raise ValueError("phi must be non-empty.")
    if chebyshev_nodes.numel() < 1:
        raise ValueError("chebyshev_nodes must be non-empty.")
    if not isinstance(parity, (int, np.integer)):
        raise ValueError("parity must be an integer parity indicator.")

    device = phi.device
    phi_reduced = phi.to(device=device, dtype=torch.float64)
    degree_odd = _degree_is_odd(parity)
    phi_full = _full_phases_from_reduced_torch(phi_reduced, degree_is_odd=degree_odd)

    x = chebyshev_nodes.to(device=device, dtype=torch.float64)
    pred = torch.real(build_qsp_unitary(phi_full, x)[:, 0, 0])

    target = _evaluate_target(target_fn, x, device)
    if target.numel() != x.numel():
        raise ValueError(
            f"target_fn must return one value per node: expected {x.numel()}, got {target.numel()}."
        )

    target = torch.real(target.to(dtype=torch.complex128)).to(dtype=torch.float64)
    return torch.mean((pred - target) ** 2)

optimize_symmetric_qsp(target_coeffs, parity, max_iter=1000, tol=1e-12, lr=0.1, phase_padding=False)

Optimize symmetric QSP phases with PyTorch LBFGS.

Uses a reduced symmetric phase vector :math:\\hat\\phi, initializes :math:\\hat\\phi=(\\pi/4,0,\\dots,0), and optimizes with LBFGS (strong Wolfe line search). Returns the full symmetric phase vector :math:\\Phi as a NumPy array.

Parameters:

Name Type Description Default
target_coeffs ndarray

Target polynomial coefficients in Chebyshev basis.

 required
parity int

Target polynomial parity indicator.

 required
max_iter int

Maximum number of optimizer iterations per attempt.

 1000
tol float

Convergence tolerance on loss.

 1e-12
lr float

LBFGS learning rate.

 0.1
phase_padding bool

If True, do one restart with padded phases (prepend/append :math:\\pi/4) to approximately double degree if first attempt does not converge.

 False

Returns:

Type Description
ndarray

Optimized full symmetric phase array :math:\\Phi.
Source code in src\qsp_proc\completion\optimization.py 
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
def optimize_symmetric_qsp(
    target_coeffs: np.ndarray,
    parity: int,
    max_iter: int = 1000,
    tol: float = 1e-12,
    lr: float = 0.1,
    phase_padding: bool = False,
) -> np.ndarray:
    r"""Optimize symmetric QSP phases with PyTorch LBFGS.

    Uses a reduced symmetric phase vector :math:`\\hat\\phi`, initializes
    :math:`\\hat\\phi=(\\pi/4,0,\\dots,0)`, and optimizes with LBFGS
    (strong Wolfe line search). Returns the full symmetric phase vector
    :math:`\\Phi` as a NumPy array.

    Args:
        target_coeffs: Target polynomial coefficients in Chebyshev basis.
        parity: Target polynomial parity indicator.
        max_iter: Maximum number of optimizer iterations per attempt.
        tol: Convergence tolerance on loss.
        lr: LBFGS learning rate.
        phase_padding: If True, do one restart with padded phases (prepend/append
            :math:`\\pi/4`) to approximately double degree if first attempt
            does not converge.

    Returns:
        Optimized full symmetric phase array :math:`\\Phi`.
    """
    coeffs = np.asarray(target_coeffs, dtype=np.float64)
    if coeffs.ndim != 1 or coeffs.size < 1:
        raise ValueError("target_coeffs must be a non-empty 1D array.")
    if max_iter < 1:
        raise ValueError("max_iter must be >= 1.")
    if tol <= 0:
        raise ValueError("tol must be positive.")
    if lr <= 0:
        raise ValueError("lr must be positive.")
    if not isinstance(parity, (int, np.integer)):
        raise ValueError("parity must be an integer parity indicator.")

    degree_odd = _degree_is_odd(parity)
    reduced_len = (coeffs.size + 1) // 2

    init = np.zeros(reduced_len, dtype=np.float64)
    init[0] = np.pi / 4.0

    best_full = _full_phases_from_reduced_numpy(init, degree_is_odd=degree_odd)
    converged = False
    current_init = init
    attempts = 2 if phase_padding else 1

    for attempt in range(attempts):
        full_len = 2 * current_init.size if degree_odd else 2 * current_init.size - 1
        num_nodes = max(64, 2 * full_len)
        nodes_np = _chebyshev_first_kind_nodes(num_nodes)
        target_np = np.polynomial.chebyshev.chebval(nodes_np, coeffs)

        nodes = torch.tensor(nodes_np, dtype=torch.float64)
        targets = torch.tensor(target_np, dtype=torch.float64)

        def target_fn(_: torch.Tensor) -> torch.Tensor:
            return targets

        phases = torch.nn.Parameter(torch.tensor(current_init, dtype=torch.float64))
        optimizer = torch.optim.LBFGS(
            [phases],
            lr=lr,
            max_iter=1,
            line_search_fn="strong_wolfe",
        )

        last_loss = float("inf")
        for it in range(1, max_iter + 1):

            def closure() -> torch.Tensor:
                optimizer.zero_grad()
                loss = symmetric_qsp_loss(phases, target_fn, nodes, parity)
                loss.backward()
                return loss

            loss_tensor = optimizer.step(closure)
            last_loss = float(loss_tensor.detach().cpu().item())

            if it % 50 == 0:
                logger.info(
                    "Symmetric-QSP optimize (attempt %d) iter %d: loss = %.3e",
                    attempt + 1,
                    it,
                    last_loss,
                )

            if last_loss < tol:
                converged = True
                break

        best_full = _full_phases_from_reduced_torch(
            phases.detach(),
            degree_is_odd=degree_odd,
        ).numpy()
        if converged:
            break

        if phase_padding and attempt == 0:
            padded = best_full.copy()
            target_len = 2 * best_full.size
            while padded.size < target_len:
                padded = np.concatenate(([np.pi / 4.0], padded, [np.pi / 4.0]))
            current_init = _reduced_from_full_numpy(padded, degree_is_odd=degree_odd)

    if not converged:
        logger.info(
            "Symmetric-QSP optimization reached max_iter without hitting tol; "
            "returning last iterate."
        )

    return best_full

gqsp_completion_loss(q_coeffs, p_coeffs)

G-QSP completion loss from Skelton 2025 Eq. (54).

Squared Frobenius / :math:\ell_2 norm of :math:p \star \mathrm{reverse}(p^*) + q \star \mathrm{reverse}(q^*) - \delta, with :math:\delta the Kronecker delta at lag zero.

Source code in src\qsp_proc\completion\optimization.py 
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
def gqsp_completion_loss(q_coeffs: torch.Tensor, p_coeffs: torch.Tensor) -> torch.Tensor:
    r"""G-QSP completion loss from Skelton 2025 Eq. (54).

    Squared Frobenius / :math:`\ell_2` norm of
    :math:`p \star \mathrm{reverse}(p^*) + q \star \mathrm{reverse}(q^*) - \delta`,
    with :math:`\delta` the Kronecker delta at lag zero.
    """
    _gqsp_validate_pair(q_coeffs, p_coeffs)
    p = p_coeffs.to(dtype=torch.complex128)
    q = q_coeffs.to(dtype=torch.complex128, device=p.device)

    autocorr_sum = _laurent_autocorrelation(p) + _laurent_autocorrelation(q)
    delta = torch.zeros_like(autocorr_sum)
    delta[autocorr_sum.numel() // 2] = 1.0
    residual = autocorr_sum - delta
    return torch.linalg.vector_norm(residual, ord=2) ** 2

optimize_gqsp_completion(p_coeffs, max_iter=2000, tol=1e-10, *, num_restarts=8, q_init=None)

Optimize complementary Q coefficients for G-QSP completion.

Tries several deterministic random starts (seeds 0 … num_restarts-1): a single draw can trap LBFGS in a poor basin, while another seed often reaches very small Skelton loss. If q_init is given, the first run starts from those coefficients instead of a random draw (useful for high-degree P where random restarts stall).

Note

In line with Skelton 2025, this optimization method commonly plateaus near ~1e-8 loss precision, typically below FFT/root-based methods.

Source code in src\qsp_proc\completion\optimization.py 
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
def optimize_gqsp_completion(
    p_coeffs: np.ndarray,
    max_iter: int = 2000,
    tol: float = 1e-10,
    *,
    num_restarts: int = 8,
    q_init: np.ndarray | None = None,
) -> np.ndarray:
    """Optimize complementary ``Q`` coefficients for G-QSP completion.

    Tries several deterministic random starts (seeds ``0 … num_restarts-1``): a single
    draw can trap LBFGS in a poor basin, while another seed often reaches very small
    Skelton loss. If ``q_init`` is given, the first run starts from those coefficients
    instead of a random draw (useful for high-degree ``P`` where random restarts stall).

    Note:
        In line with Skelton 2025, this optimization method commonly plateaus near
        ~1e-8 loss precision, typically below FFT/root-based methods.
    """
    if max_iter < 1:
        raise ValueError("max_iter must be >= 1.")
    if tol <= 0.0:
        raise ValueError("tol must be positive.")
    if num_restarts < 1:
        raise ValueError("num_restarts must be >= 1.")

    p_np = np.asarray(p_coeffs, dtype=np.complex128)
    if p_np.ndim != 1 or p_np.size == 0:
        raise ValueError("p_coeffs must be a non-empty 1D array.")

    p_fixed = torch.tensor(p_np, dtype=torch.complex128)
    n = p_np.size

    if q_init is not None:
        q0_provided = np.asarray(q_init, dtype=np.complex128)
        if q0_provided.shape != (n,):
            raise ValueError(f"q_init must have shape ({n},), got {q0_provided.shape}.")

    best_loss = float("inf")
    best_q: np.ndarray | None = None

    for restart in range(num_restarts):
        if restart == 0 and q_init is not None:
            q0 = np.asarray(q_init, dtype=np.complex128)
        else:
            rng = np.random.default_rng(restart)
            q0 = _gqsp_random_q_init(rng, n)
        q_try, loss_try, ok = _lbfgs_gqsp_completion_run(
            p_fixed,
            q0,
            max_iter=max_iter,
            tol=tol,
            restart_label=restart + 1,
        )
        if loss_try < best_loss:
            best_loss = loss_try
            best_q = q_try
        if ok and best_loss < tol:
            break

    if best_loss >= tol:
        logger.info(
            "G-QSP completion finished without tol; best loss = %.3e over %d restart(s) "
            "(often plateaus near ~1e-8).",
            best_loss,
            num_restarts,
        )

    if best_q is None:
        raise RuntimeError("G-QSP completion produced no iterate; check p_coeffs and num_restarts.")
    return best_q

qsp_proc.completion.rootfinding

Haah-style root-finding completion (Skelton 2025 Sec. 5.1.1; Haah, Quantum 3:190, 2019).

Factors the palindromic polynomial h(z)=z^n(1-|P|^2) on the circle (Fejér data). The u=x^2 Chebyshev / Appendix-A sqrt(1-x^2) * W(u) description in the guide matches this outer factorization after x = cos(θ) when the defect is even in θ; the code uses the standard complex z-domain root split (Wilson/Fejér compatible).

partition_reciprocal_roots(roots, tol=1e-10)

Group roots into reciprocal pairs (r, s) with r * s ≈ 1.

Each root appears in exactly one pair. A root with r^2 ≈ 1 (typically ±1) is paired with itself. Matching is greedy: for each unused root, take the best unused partner among later indices; otherwise accept a self-pair if r^2 ≈ 1.

Parameters:

Name Type Description Default
roots Sequence[complex] | ndarray

Root values (e.g. from numpy.roots on a palindromic polynomial).

 required
tol float

Tolerance for |r s - 1|, scaled by 1 + |r||s|.

 1e-10

Returns:

Type Description
list[tuple[complex, complex]]

List of pairs (r, s). Order is not canonical.

Raises:

Type Description
ValueError

If some root has no reciprocal partner under the tolerance.
Source code in src\qsp_proc\completion\rootfinding.py 
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def partition_reciprocal_roots(
    roots: Sequence[complex] | np.ndarray,
    tol: float = 1e-10,
) -> list[tuple[complex, complex]]:
    """Group roots into reciprocal pairs ``(r, s)`` with ``r * s ≈ 1``.

    Each root appears in exactly one pair. A root with ``r^2 ≈ 1`` (typically ``±1``) is paired
    with itself. Matching is greedy: for each unused root, take the best unused partner among
    later indices; otherwise accept a self-pair if ``r^2 ≈ 1``.

    Args:
        roots: Root values (e.g. from ``numpy.roots`` on a palindromic polynomial).
        tol: Tolerance for ``|r s - 1|``, scaled by ``1 + |r||s|``.

    Returns:
        List of pairs ``(r, s)``. Order is not canonical.

    Raises:
        ValueError: If some root has no reciprocal partner under the tolerance.
    """
    rts = np.asarray(roots, dtype=np.complex128).ravel()
    if rts.size == 0:
        return []

    n = rts.size
    used = np.zeros(n, dtype=bool)
    pairs: list[tuple[complex, complex]] = []

    for i in range(n):
        if used[i]:
            continue
        r = complex(rts[i])

        best_j: int | None = None
        best_res = np.inf
        for j in range(i + 1, n):
            if used[j]:
                continue
            s = complex(rts[j])
            res = abs(r * s - 1.0)
            if res <= _scaled_reciprocal_tol(r, s, tol) and res < best_res:
                best_res = res
                best_j = j

        if best_j is not None:
            s = complex(rts[best_j])
            pairs.append((r, s))
            used[i] = used[best_j] = True
            continue

        if _is_reciprocal_pair(r, r, tol):
            pairs.append((r, r))
            used[i] = True
            continue

        raise ValueError(f"partition_reciprocal_roots: no reciprocal partner for root {r!r}")

    return pairs

haah_completion(poly, *, root_tol=1e-10, num_validate=4096, validate_tol=1e-05)

Complete P with causal Q via Haah root finding (Skelton Sec. 5.1.1).

Builds h(z)=z^n F(z) with F(z)=1-P(z)\overline{P(1/z)} (degree 2n, n=len(P.coeffs)-1). Roots of h pair as (r,1/r); the outer factor \gamma(z)=K\prod_{s\in S}(z-s) takes |s|>1 when possible (unit circle: break by \operatorname{Im}s). Scale K so |P|^2+|Q|^2\approx 1 on U(1).

Parameters:

Name Type Description Default
poly LaurentPoly

Laurent P with |P|\le 1 on U(1) (so F\ge 0 there).

 required
root_tol float

Reciprocal pairing tolerance for :func:partition_reciprocal_roots.

 1e-10
num_validate int

Number of circle samples for scaling and validation.

 4096
validate_tol float

Max l_\infty error allowed for |P|^2+|Q|^2-1.

 1e-05

Returns:

Type Description
LaurentPoly

Causal complementary Q (\min\deg=0).

Raises:

Type Description
TypeError

If poly is not a :class:~qsp_proc.polynomials.laurent.LaurentPoly.
ValueError

If the Fejér DC term is negative or pairing does not cover all roots.
RuntimeError

If scaling or validation fails numerically.
Source code in src\qsp_proc\completion\rootfinding.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
def haah_completion(
    poly: LaurentPoly,
    *,
    root_tol: float = 1e-10,
    num_validate: int = 4096,
    validate_tol: float = 1e-5,
) -> LaurentPoly:
    r"""Complete ``P`` with causal ``Q`` via Haah root finding (Skelton Sec. 5.1.1).

    Builds ``h(z)=z^n F(z)`` with ``F(z)=1-P(z)\overline{P(1/z)}`` (degree ``2n``,
    ``n=len(P.coeffs)-1``). Roots of ``h`` pair as ``(r,1/r)``; the outer factor
    ``\gamma(z)=K\prod_{s\in S}(z-s)`` takes ``|s|>1`` when possible (unit circle: break by
    ``\operatorname{Im}s``). Scale ``K`` so ``|P|^2+|Q|^2\approx 1`` on ``U(1)``.

    Args:
        poly: Laurent ``P`` with ``|P|\le 1`` on ``U(1)`` (so ``F\ge 0`` there).
        root_tol: Reciprocal pairing tolerance for :func:`partition_reciprocal_roots`.
        num_validate: Number of circle samples for scaling and validation.
        validate_tol: Max ``l_\infty`` error allowed for ``|P|^2+|Q|^2-1``.

    Returns:
        Causal complementary ``Q`` (``\min\deg=0``).

    Raises:
        TypeError: If ``poly`` is not a :class:`~qsp_proc.polynomials.laurent.LaurentPoly`.
        ValueError: If the Fejér DC term is negative or pairing does not cover all roots.
        RuntimeError: If scaling or validation fails numerically.
    """
    if not isinstance(poly, LaurentPoly):
        raise TypeError("haah_completion expects a LaurentPoly.")

    p_coeffs = poly.coeffs.astype(np.complex128, copy=False)
    span = p_coeffs.size - 1
    h_asc = np.asarray(_f_coeffs_one_minus_modulus_sq(p_coeffs), dtype=np.complex128)
    deg_h = h_asc.size - 1

    if deg_h > 500:
        logger.warning(
            "Haah root-finding on degree-%d h(z)=z^n(1-|P|^2) is fragile; "
            "prefer Wilson/BerSün or optimization.",
            deg_h,
        )

    defect = one_minus_modulus_squared_laurent(poly)
    dc = float(np.real_if_close(h_asc[span]).item())
    if -_DC_NEG_CLIP < dc < 0.0:
        dc = 0.0
    if dc < 0.0:
        raise ValueError(
            "haah_completion: Fejér DC F(1) is negative; input likely has |P|>1 on U(1)."
        )

    if deg_h < 1:
        return LaurentPoly(np.array([np.sqrt(max(dc, 0.0))], dtype=np.complex128), min_degree=0)

    # np.roots expects highest-degree first.
    roots = np.roots(h_asc[::-1])
    pairs = partition_reciprocal_roots(roots, tol=root_tol)
    if len(pairs) != roots.size // 2:
        raise ValueError(
            "haah_completion: reciprocal pairing incomplete "
            f"(got {len(pairs)} pairs for {roots.size} roots)."
        )

    outer_roots = np.array(
        [_outer_root_from_pair(pr, root_tol) for pr in pairs],
        dtype=np.complex128,
    )
    gamma_asc = npp.polyfromroots(outer_roots)

    z = _unit_circle_grid(num_validate)
    f_on_circle = np.clip(np.real(defect(z)), 0.0, np.inf)
    g_on_circle = npp.polyval(z, gamma_asc)
    # Heuristic: median ratio is more robust than analytic K = sqrt(F₀/γ₀²)
    # when root-finding noise corrupts individual γ evaluations.
    scale = _median_positive_scale(f_on_circle, np.abs(g_on_circle) ** 2)
    q = LaurentPoly(gamma_asc * scale, min_degree=0)

    err = float(np.max(np.abs(1.0 - np.abs(poly(z)) ** 2 - np.abs(q(z)) ** 2)))
    logger.info("Haah completion validation: ||1-|P|²-|Q|²||_inf ≈ %.3e", err)
    if err > validate_tol:
        raise RuntimeError(
            f"haah_completion: complementarity check failed "
            f"(l_inf {err:.3e} > {validate_tol:.3e})."
        )

    return q