Skip to content

Decomposition

qsp_proc.decomposition

Public API for decomposition routines and validation helpers.

DecompositionResult dataclass

Output of recursive carving: constant residue, projectors, and carve metadata.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
11
12
13
14
15
16
17
18
@dataclass
class DecompositionResult:
    """Output of recursive carving: constant residue, projectors, and carve metadata."""

    e0: np.ndarray
    projectors: list[np.ndarray]
    convention: str
    carve_sides: list[bool] | None = None

MatrixLaurentPoly

Two-by-two matrix-valued Laurent polynomial F(z) = sum_k C_k z^k.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 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
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
class MatrixLaurentPoly:
    """Two-by-two matrix-valued Laurent polynomial ``F(z) = sum_k C_k z^k``."""

    def __init__(self, coeffs: np.ndarray, min_degree: int) -> None:
        """Initialize from dense ``(num_coeffs, 2, 2)`` blocks and minimum Laurent exponent."""
        if not isinstance(min_degree, Integral):
            raise TypeError("min_degree must be an integer.")

        coeffs_arr = np.asarray(coeffs, dtype=np.complex128)
        if coeffs_arr.ndim != 3 or coeffs_arr.shape[1:] != (2, 2):
            raise ValueError("coeffs must have shape (num_coeffs, 2, 2).")
        if coeffs_arr.shape[0] <= 0:
            raise ValueError("coeffs must contain at least one 2x2 coefficient matrix.")

        self._coeffs = coeffs_arr
        self._min_degree = int(min_degree)

    @property
    def coeffs(self) -> np.ndarray:
        """Dense coefficient array (copy).

        Same encapsulation pattern as :class:`~qsp_proc.polynomials.laurent.LaurentPoly`.
        """
        return self._coeffs.copy()

    @property
    def min_degree(self) -> int:
        """Minimum Laurent exponent with a stored block."""
        return self._min_degree

    @property
    def max_degree(self) -> int:
        """Maximum Laurent exponent with a stored block."""
        return self._min_degree + self._coeffs.shape[0] - 1

    @property
    def exponents(self) -> np.ndarray:
        """Stored Laurent exponents as a dense integer range."""
        return np.arange(self.min_degree, self.max_degree + 1, dtype=np.int64)

    @property
    def degree(self) -> int:
        """Half-width degree ``max(|min_degree|, |max_degree|)``."""
        return max(abs(self.min_degree), abs(self.max_degree))

    def coefficient(self, k: int) -> np.ndarray:
        """Return coefficient matrix ``C_k``, or a ``2x2`` zero matrix if outside support."""
        idx = int(k) - self._min_degree
        if 0 <= idx < self._coeffs.shape[0]:
            return self._coeffs[idx]
        return np.zeros((2, 2), dtype=np.complex128)

    def evaluate(self, z: np.ndarray | complex) -> np.ndarray:
        """Evaluate at scalar or batched ``z`` with trailing output shape ``(2, 2)``."""
        z_arr = np.asarray(z, dtype=np.complex128)
        powers = z_arr[..., None] ** self.exponents
        values = np.asarray(
            np.einsum("...k,kij->...ij", powers, self._coeffs, optimize=True),
            dtype=np.complex128,
        )
        if z_arr.ndim == 0:
            return values.reshape(2, 2)
        return values

    __call__ = evaluate

    def __matmul__(self, other: object) -> MatrixLaurentPoly:
        """Multiply two matrix-valued Laurent polynomials via coefficient convolution."""
        if not isinstance(other, MatrixLaurentPoly):
            return NotImplemented

        left = self._coeffs
        right = other._coeffs
        out_len = left.shape[0] + right.shape[0] - 1
        out = np.zeros((out_len, 2, 2), dtype=np.complex128)

        for i, left_block in enumerate(left):
            out[i : i + right.shape[0]] += left_block @ right

        return MatrixLaurentPoly(out, min_degree=self.min_degree + other.min_degree)

    def is_su2_on_circle(self, num_samples: int = 512, tol: float = 1e-10) -> bool:
        """Whether ``F(e^{iθ})`` is approximately in ``SU(2)`` on uniform circle samples."""
        num_samples_int = int(num_samples)
        if num_samples_int < 1:
            raise ValueError("num_samples must be >= 1.")
        if tol < 0.0:
            raise ValueError("tol must be non-negative.")

        theta = np.linspace(0.0, 2.0 * np.pi, num_samples_int, endpoint=False)
        z = np.exp(1j * theta)
        fz = self.evaluate(z)
        fz_dag = np.conjugate(np.swapaxes(fz, -1, -2))
        gram = np.matmul(fz_dag, fz)
        eye = np.eye(2, dtype=np.complex128)
        unitary_err = float(np.max(np.abs(gram - eye)))
        det_err = float(np.max(np.abs(np.linalg.det(fz) - 1.0)))
        return unitary_err <= tol and det_err <= tol

coeffs property

Dense coefficient array (copy).

Same encapsulation pattern as :class:~qsp_proc.polynomials.laurent.LaurentPoly.

min_degree property

Minimum Laurent exponent with a stored block.

max_degree property

Maximum Laurent exponent with a stored block.

exponents property

Stored Laurent exponents as a dense integer range.

degree property

Half-width degree max(|min_degree|, |max_degree|).

__init__(coeffs, min_degree)

Initialize from dense (num_coeffs, 2, 2) blocks and minimum Laurent exponent.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
def __init__(self, coeffs: np.ndarray, min_degree: int) -> None:
    """Initialize from dense ``(num_coeffs, 2, 2)`` blocks and minimum Laurent exponent."""
    if not isinstance(min_degree, Integral):
        raise TypeError("min_degree must be an integer.")

    coeffs_arr = np.asarray(coeffs, dtype=np.complex128)
    if coeffs_arr.ndim != 3 or coeffs_arr.shape[1:] != (2, 2):
        raise ValueError("coeffs must have shape (num_coeffs, 2, 2).")
    if coeffs_arr.shape[0] <= 0:
        raise ValueError("coeffs must contain at least one 2x2 coefficient matrix.")

    self._coeffs = coeffs_arr
    self._min_degree = int(min_degree)

coefficient(k)

Return coefficient matrix C_k, or a 2x2 zero matrix if outside support.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
66
67
68
69
70
71
def coefficient(self, k: int) -> np.ndarray:
    """Return coefficient matrix ``C_k``, or a ``2x2`` zero matrix if outside support."""
    idx = int(k) - self._min_degree
    if 0 <= idx < self._coeffs.shape[0]:
        return self._coeffs[idx]
    return np.zeros((2, 2), dtype=np.complex128)

evaluate(z)

Evaluate at scalar or batched z with trailing output shape (2, 2).

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
73
74
75
76
77
78
79
80
81
82
83
def evaluate(self, z: np.ndarray | complex) -> np.ndarray:
    """Evaluate at scalar or batched ``z`` with trailing output shape ``(2, 2)``."""
    z_arr = np.asarray(z, dtype=np.complex128)
    powers = z_arr[..., None] ** self.exponents
    values = np.asarray(
        np.einsum("...k,kij->...ij", powers, self._coeffs, optimize=True),
        dtype=np.complex128,
    )
    if z_arr.ndim == 0:
        return values.reshape(2, 2)
    return values

__matmul__(other)

Multiply two matrix-valued Laurent polynomials via coefficient convolution.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def __matmul__(self, other: object) -> MatrixLaurentPoly:
    """Multiply two matrix-valued Laurent polynomials via coefficient convolution."""
    if not isinstance(other, MatrixLaurentPoly):
        return NotImplemented

    left = self._coeffs
    right = other._coeffs
    out_len = left.shape[0] + right.shape[0] - 1
    out = np.zeros((out_len, 2, 2), dtype=np.complex128)

    for i, left_block in enumerate(left):
        out[i : i + right.shape[0]] += left_block @ right

    return MatrixLaurentPoly(out, min_degree=self.min_degree + other.min_degree)

is_su2_on_circle(num_samples=512, tol=1e-10)

Whether F(e^{iθ}) is approximately in SU(2) on uniform circle samples.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def is_su2_on_circle(self, num_samples: int = 512, tol: float = 1e-10) -> bool:
    """Whether ``F(e^{iθ})`` is approximately in ``SU(2)`` on uniform circle samples."""
    num_samples_int = int(num_samples)
    if num_samples_int < 1:
        raise ValueError("num_samples must be >= 1.")
    if tol < 0.0:
        raise ValueError("tol must be non-negative.")

    theta = np.linspace(0.0, 2.0 * np.pi, num_samples_int, endpoint=False)
    z = np.exp(1j * theta)
    fz = self.evaluate(z)
    fz_dag = np.conjugate(np.swapaxes(fz, -1, -2))
    gram = np.matmul(fz_dag, fz)
    eye = np.eye(2, dtype=np.complex128)
    unitary_err = float(np.max(np.abs(gram - eye)))
    det_err = float(np.max(np.abs(np.linalg.det(fz) - 1.0)))
    return unitary_err <= tol and det_err <= tol

decompose(p, q, convention='gqsp', tol=1e-10)

Unified decomposition entry point for supported QSP conventions.

Parameters:

Name Type Description Default
p LaurentPoly

Primary Laurent polynomial P.

 required
q LaurentPoly

Complementary Laurent polynomial / Fejér factor from completion.

 required
convention str

Convention selector ("gqsp", "laurent", or "auto").

 'gqsp'
tol float

Numerical tolerance used by decomposition routines and verification.

 1e-10

Returns:

Type Description
DecompositionResult

Recursive carving decomposition result.

Raises:

Type Description
ValueError

If convention is unsupported.
Source code in src\qsp_proc\decomposition\carving.py 
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
def decompose(
    p: LaurentPoly,
    q: LaurentPoly,
    convention: str = "gqsp",
    tol: float = 1e-10,
) -> DecompositionResult:
    """Unified decomposition entry point for supported QSP conventions.

    Args:
        p: Primary Laurent polynomial ``P``.
        q: Complementary Laurent polynomial / Fejér factor from completion.
        convention: Convention selector (``"gqsp"``, ``"laurent"``, or ``"auto"``).
        tol: Numerical tolerance used by decomposition routines and verification.

    Returns:
        Recursive carving decomposition result.

    Raises:
        ValueError: If ``convention`` is unsupported.
    """
    conv = convention.lower()
    if conv not in {"gqsp", "laurent", "auto"}:
        raise ValueError(
            f"Unsupported convention '{convention}'. "
            "Expected one of {'gqsp', 'laurent', 'auto'}."
        )

    # ``auto`` defaults to the most general currently supported pathway.
    target_conv = "gqsp" if conv == "auto" else conv

    if target_conv == "gqsp":
        matrix_poly = build_matrix_poly_gqsp(p, q)
    else:
        from qsp_proc.polynomials.chebyshev import ChebyshevPoly

        first_kind, second_kind = p.to_appendix_a_components()
        a_laurent = LaurentPoly.from_chebyshev(first_kind)
        b_laurent = LaurentPoly.from_appendix_a_components(
            first_kind=ChebyshevPoly([0.0]),
            second_kind_coeffs=second_kind,
        )
        gamma = q
        matrix_poly = build_matrix_poly_laurent_qsp(a_laurent, b_laurent, gamma)

    # Decompose the convention-built matrix polynomial directly in the original
    # Laurent variable so reconstruction/verification map back to the same P, Q.
    result = recursive_carve(matrix_poly, tol=tol)
    result.convention = target_conv

    if conv == "auto":
        max_error = verify_decomposition(
            p=p,
            q=q,
            result=result,
            convention="gqsp",
            tol=tol,
        )
        if max_error > tol:
            LOGGER.warning(
                "Decomposition verification error %.3e exceeds tol %.3e in auto mode.",
                max_error,
                tol,
            )

    return result

reconstruct_from_decomposition(result, num_samples=256)

Reconstruct F(z) from E_0 and carved projectors.

Uses F(z) = E_0 @ prod_{j=1}^{m} E_{P_j}(z), where E_{P_j}(z) = (I - P_j) + z^{s_j} P_j, where exponents follow result.carve_sides (top carve -> s_j = +1, bottom carve -> s_j = -1). If carve_sides is unavailable (legacy results), parity fallback is used.

Parameters:

Name Type Description Default
result DecompositionResult

Decomposition output containing e0 and projector sequence.

 required
num_samples int

Reserved for API compatibility and optional future checks.

 256

Returns:

Type Description
MatrixLaurentPoly

Reconstructed matrix-valued Laurent polynomial.
Source code in src\qsp_proc\decomposition\carving.py 
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
def reconstruct_from_decomposition(
    result: DecompositionResult, num_samples: int = 256
) -> MatrixLaurentPoly:
    """Reconstruct ``F(z)`` from ``E_0`` and carved projectors.

    Uses
    ``F(z) = E_0 @ prod_{j=1}^{m} E_{P_j}(z)``, where
    ``E_{P_j}(z) = (I - P_j) + z^{s_j} P_j``, where exponents follow
    ``result.carve_sides`` (top carve -> ``s_j = +1``, bottom carve -> ``s_j = -1``).
    If ``carve_sides`` is unavailable (legacy results), parity fallback is used.

    Args:
        result: Decomposition output containing ``e0`` and projector sequence.
        num_samples: Reserved for API compatibility and optional future checks.

    Returns:
        Reconstructed matrix-valued Laurent polynomial.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")

    e0 = np.asarray(result.e0, dtype=np.complex128)
    if e0.shape != (2, 2):
        raise ValueError("result.e0 must have shape (2, 2).")

    reconstructed = MatrixLaurentPoly(e0[np.newaxis, ...], min_degree=0)
    identity = np.eye(2, dtype=np.complex128)

    projectors = result.projectors
    carve_sides = result.carve_sides
    if carve_sides is not None and len(carve_sides) != len(projectors):
        raise ValueError("result.carve_sides must have the same length as projectors.")

    if carve_sides is None:
        # Backward-compatibility fallback for older decomposition outputs.
        carve_sides = [idx % 2 == 0 for idx in range(len(projectors))]

    for idx in range(len(projectors) - 1, -1, -1):
        projector = projectors[idx]
        from_top = carve_sides[idx]
        proj = np.asarray(projector, dtype=np.complex128)
        if proj.shape != (2, 2):
            raise ValueError("Each projector must have shape (2, 2).")

        q_proj = identity - proj

        if from_top:
            # If carved with E_P^{-1}(z)=Q+z^{-1}P, multiply back by E_P(z)=Q+zP.
            factor_coeffs = np.stack([q_proj, proj], axis=0)
            factor = MatrixLaurentPoly(factor_coeffs, min_degree=0)
        else:
            # If carved with E_P^{-1}(z)=Q+zP, multiply back by E_P(z)=Q+z^{-1}P.
            factor_coeffs = np.stack([proj, q_proj], axis=0)
            factor = MatrixLaurentPoly(factor_coeffs, min_degree=-1)

        reconstructed = reconstructed @ factor

    return reconstructed

verify_decomposition(p, q, result, convention='gqsp', num_samples=512, tol=1e-08)

Verify a decomposition by sampled operator-norm reconstruction error.

This computes: max_z ||F(z) - F_recon(z)||_op over uniformly sampled points on U(1).

Parameters:

Name Type Description Default
p LaurentPoly

Primary Laurent polynomial P.

 required
q LaurentPoly

Complementary Laurent polynomial Q.

 required
result DecompositionResult

Decomposition result.

 required
convention str

QSP convention selector; currently only "gqsp" is supported.

 'gqsp'
num_samples int

Number of unit-circle samples.

 512
tol float

Tolerance used for optional debug logging threshold.

 1e-08

Returns:

Type Description
float

Maximum sampled operator-norm error.
Source code in src\qsp_proc\decomposition\carving.py 
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
476
477
478
479
480
def verify_decomposition(
    p: LaurentPoly,
    q: LaurentPoly,
    result: DecompositionResult,
    convention: str = "gqsp",
    num_samples: int = 512,
    tol: float = 1e-8,
) -> float:
    """Verify a decomposition by sampled operator-norm reconstruction error.

    This computes:
    ``max_z ||F(z) - F_recon(z)||_op`` over uniformly sampled points on ``U(1)``.

    Args:
        p: Primary Laurent polynomial ``P``.
        q: Complementary Laurent polynomial ``Q``.
        result: Decomposition result.
        convention: QSP convention selector; currently only ``"gqsp"`` is supported.
        num_samples: Number of unit-circle samples.
        tol: Tolerance used for optional debug logging threshold.

    Returns:
        Maximum sampled operator-norm error.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")
    if tol < 0.0:
        raise ValueError("tol must be non-negative.")

    conv = convention.lower()
    if conv != "gqsp":
        raise NotImplementedError(
            f"verify_decomposition currently supports only 'gqsp', got '{convention}'."
        )

    f_original = build_matrix_poly_gqsp(p, q)
    f_reconstructed = reconstruct_from_decomposition(result, num_samples=num_samples)

    theta = np.linspace(0.0, 2.0 * np.pi, num_samples, endpoint=False, dtype=np.float64)
    z = np.exp(1j * theta)

    diff = f_original(z) - f_reconstructed(z)
    op_norms = np.linalg.norm(diff, ord=2, axis=(-2, -1))
    max_error = float(np.max(op_norms))

    if max_error > tol:
        LOGGER.debug(
            "Decomposition verification error %.3e exceeds tol %.3e",
            max_error,
            tol,
        )

    return max_error

qsp_proc.decomposition.carving

QSP decomposition via recursive carving (Skelton 2025, Appendix D, Algorithm 2).

Core routines factor a matrix-valued Laurent polynomial by rank-1 projectors and extract a constant SU(2) residue. Matrix construction lives in :mod:qsp_proc.decomposition.builders; the polynomial type in :mod:qsp_proc.decomposition.matrix_laurent_poly.

extract_projector_from_rank1(C, tol=1e-10)

Recover rank-1 projector P satisfying C @ P = 0 via SVD null-space extraction.

Parameters:

Name Type Description Default
C ndarray

Input matrix with shape (2, 2) expected to be rank-1.

 required
tol float

Numerical tolerance for rank detection and projector validation.

 1e-10

Returns:

Type Description
ndarray

Projector matrix P with shape (2, 2).

Raises:

Type Description
ValueError

If C is not 2x2 or is numerically zero.
Source code in src\qsp_proc\decomposition\carving.py 
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
def extract_projector_from_rank1(C: np.ndarray, tol: float = 1e-10) -> np.ndarray:
    """Recover rank-1 projector ``P`` satisfying ``C @ P = 0`` via SVD null-space extraction.

    Args:
        C: Input matrix with shape ``(2, 2)`` expected to be rank-1.
        tol: Numerical tolerance for rank detection and projector validation.

    Returns:
        Projector matrix ``P`` with shape ``(2, 2)``.

    Raises:
        ValueError: If ``C`` is not ``2x2`` or is numerically zero.
    """
    C_arr = np.asarray(C, dtype=np.complex128)
    if C_arr.shape != (2, 2):
        raise ValueError("C must have shape (2, 2).")

    _, singular_values, vh = np.linalg.svd(C_arr, full_matrices=True)

    if singular_values[0] < tol and singular_values[1] < tol:
        raise ValueError("C is numerically zero; cannot extract a unique projector.")

    if singular_values[1] > tol:
        warnings.warn(
            "Input matrix is not numerically rank-1; using smallest-singular-vector null-space approximation.",
            RuntimeWarning,
            stacklevel=2,
        )

    # Right null-space vector for the smallest singular value.
    v = vh[1, :].conjugate().reshape(2, 1)
    v_norm = np.linalg.norm(v)
    if v_norm < tol:
        raise ValueError("Failed to extract a stable null-space vector from C.")
    v = v / v_norm

    P = v @ np.conjugate(v.T)

    # Validate projector properties.
    if not np.allclose(P @ P, P, atol=tol, rtol=0.0):
        raise ValueError("Extracted projector failed idempotency check (P^2 = P).")
    if not np.allclose(P, np.conjugate(P.T), atol=tol, rtol=0.0):
        raise ValueError("Extracted projector failed Hermiticity check (P^dagger = P).")
    if not np.isclose(np.trace(P), 1.0, atol=tol, rtol=0.0):
        raise ValueError("Extracted projector failed trace check (Tr(P) = 1).")

    return P

snap_to_projector(near_proj, tol=1e-10)

Snap a near-projector to an exact rank-1 projector via eigendecomposition.

Following Theorem 2 / Algorithm 2 (Skelton 2025), this selects the eigenvector associated with the eigenvalue closest to 1 and forms P = |e0><e0| to enforce exact projector structure.

Parameters:

Name Type Description Default
near_proj ndarray

Approximate Hermitian projector with shape (2, 2).

 required
tol float

Numerical tolerance for input validation and spectral sanity checks.

 1e-10

Returns:

Type Description
ndarray

Exact rank-1 projector matrix with shape (2, 2).

Raises:

Type Description
ValueError

If input shape is invalid, not approximately Hermitian, or projector validation fails.
Source code in src\qsp_proc\decomposition\carving.py 
 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
138
139
140
def snap_to_projector(near_proj: np.ndarray, tol: float = 1e-10) -> np.ndarray:
    """Snap a near-projector to an exact rank-1 projector via eigendecomposition.

    Following Theorem 2 / Algorithm 2 (Skelton 2025), this selects the eigenvector
    associated with the eigenvalue closest to ``1`` and forms
    ``P = |e0><e0|`` to enforce exact projector structure.

    Args:
        near_proj: Approximate Hermitian projector with shape ``(2, 2)``.
        tol: Numerical tolerance for input validation and spectral sanity checks.

    Returns:
        Exact rank-1 projector matrix with shape ``(2, 2)``.

    Raises:
        ValueError: If input shape is invalid, not approximately Hermitian, or
            projector validation fails.
    """
    near_proj_arr = np.asarray(near_proj, dtype=np.complex128)
    if near_proj_arr.shape != (2, 2):
        raise ValueError("near_proj must have shape (2, 2).")

    if not np.allclose(
        near_proj_arr, np.conjugate(near_proj_arr.T), atol=tol, rtol=0.0
    ):
        raise ValueError("near_proj must be approximately Hermitian.")

    eigenvalues, eigenvectors = np.linalg.eigh(near_proj_arr)

    idx_e0 = int(np.argmin(np.abs(eigenvalues - 1.0)))
    e0 = eigenvectors[:, idx_e0]
    e0_norm = np.linalg.norm(e0)
    if e0_norm <= tol:
        raise ValueError("Failed to extract a stable eigenvector near eigenvalue 1.")
    e0 = e0 / e0_norm

    projector = np.outer(e0, np.conjugate(e0))

    machine_tol = 128.0 * np.finfo(np.float64).eps
    if not np.allclose(projector @ projector, projector, atol=machine_tol, rtol=0.0):
        raise ValueError("Snapped projector failed idempotency check (P^2 = P).")
    if not np.isclose(np.trace(projector), 1.0, atol=machine_tol, rtol=0.0):
        raise ValueError("Snapped projector failed trace check (Tr(P) = 1).")

    return projector

carve_one_step(F, projector, from_top=True)

Apply one recursive carving step via right-multiplication by E_P^{-1}.

Implements Appendix D.1.3 (Eq. 88): - from_top=True: E_P^{-1}(z) = Q + z^{-1} P - from_top=False: E_P^{-1}(z) = Q + z P where Q = I - P.

Parameters:

Name Type Description Default
F MatrixLaurentPoly

Current matrix-valued Laurent polynomial F(z) = sum_k C_k z^k.

 required
projector ndarray

Rank-1 2x2 projector P.

 required
from_top bool

If True, carve the highest-degree side; else lowest-degree side.

 True

Returns:

Type Description
MatrixLaurentPoly

Updated matrix-valued Laurent polynomial after one carving step.

Raises:

Type Description
ValueError

If projector has invalid shape or if the carved-side leading coefficient is not numerically zero.
Source code in src\qsp_proc\decomposition\carving.py 
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
def carve_one_step(
    F: MatrixLaurentPoly, projector: np.ndarray, from_top: bool = True
) -> MatrixLaurentPoly:
    """Apply one recursive carving step via right-multiplication by ``E_P^{-1}``.

    Implements Appendix D.1.3 (Eq. 88):
    - ``from_top=True``: ``E_P^{-1}(z) = Q + z^{-1} P``
    - ``from_top=False``: ``E_P^{-1}(z) = Q + z P``
    where ``Q = I - P``.

    Args:
        F: Current matrix-valued Laurent polynomial ``F(z) = sum_k C_k z^k``.
        projector: Rank-1 ``2x2`` projector ``P``.
        from_top: If ``True``, carve the highest-degree side; else lowest-degree side.

    Returns:
        Updated matrix-valued Laurent polynomial after one carving step.

    Raises:
        ValueError: If ``projector`` has invalid shape or if the carved-side leading
            coefficient is not numerically zero.
    """
    projector_arr = np.asarray(projector, dtype=np.complex128)
    if projector_arr.shape != (2, 2):
        raise ValueError("projector must have shape (2, 2).")

    matrix_poly = F
    q_projector = np.eye(2, dtype=np.complex128) - projector_arr

    coeffs = matrix_poly.coeffs
    min_degree = int(matrix_poly.min_degree)
    max_degree = int(matrix_poly.max_degree)

    # Support after multiplying by Q (same degree) and shifted P term.
    new_min_degree = min_degree - 1 if from_top else min_degree
    new_max_degree = max_degree if from_top else max_degree + 1

    new_coeffs = np.zeros(
        (new_max_degree - new_min_degree + 1, 2, 2), dtype=np.complex128
    )

    for idx, c_k in enumerate(coeffs):
        k = min_degree + idx

        # C_k @ Q contributes at degree k.
        same_idx = k - new_min_degree
        new_coeffs[same_idx] += c_k @ q_projector

        # C_k @ P contributes at degree k-1 (top carve) or k+1 (bottom carve).
        shifted_degree = k - 1 if from_top else k + 1
        shifted_idx = shifted_degree - new_min_degree
        new_coeffs[shifted_idx] += c_k @ projector_arr

    tol = 1e-10
    if from_top:
        carved_coeff = new_coeffs[-1]  # degree = original max_degree
        if np.max(np.abs(carved_coeff)) > tol:
            raise ValueError(
                "Top carving failed: highest-degree coefficient is not numerically zero."
            )
        return MatrixLaurentPoly(new_coeffs[:-1], min_degree=new_min_degree)

    carved_coeff = new_coeffs[0]  # degree = original min_degree
    if np.max(np.abs(carved_coeff)) > tol:
        raise ValueError(
            "Bottom carving failed: lowest-degree coefficient is not numerically zero."
        )
    return MatrixLaurentPoly(new_coeffs[1:], min_degree=new_min_degree + 1)

recursive_carve(F, tol=1e-10)

Run the full recursive carving loop (Algorithm 2, Skelton 2025).

Starting from a padded matrix-valued Laurent polynomial, this repeatedly extracts a rank-1 projector from the active boundary coefficient, snaps it to an exact Hermitian projector, and carves one layer until a constant matrix remains.

Parameters:

Name Type Description Default
F MatrixLaurentPoly

Input matrix-valued Laurent polynomial to decompose.

 required
tol float

Numerical tolerance for projector-snapping and diagnostics.

 1e-10

Returns:

Type Description
DecompositionResult

Decomposition result containing snapped constant unitary E0 and extracted
DecompositionResult

projectors in carve order.
Source code in src\qsp_proc\decomposition\carving.py 
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
293
294
295
296
297
298
def recursive_carve(F: MatrixLaurentPoly, tol: float = 1e-10) -> DecompositionResult:
    """Run the full recursive carving loop (Algorithm 2, Skelton 2025).

    Starting from a padded matrix-valued Laurent polynomial, this repeatedly extracts a
    rank-1 projector from the active boundary coefficient, snaps it to an exact
    Hermitian projector, and carves one layer until a constant matrix remains.

    Args:
        F: Input matrix-valued Laurent polynomial to decompose.
        tol: Numerical tolerance for projector-snapping and diagnostics.

    Returns:
        Decomposition result containing snapped constant unitary ``E0`` and extracted
        projectors in carve order.
    """
    matrix_poly = F
    current = matrix_poly
    projectors: list[np.ndarray] = []
    carve_sides: list[bool] = []
    step = 0
    su2_check_every = 10
    identity = np.eye(2, dtype=np.complex128)

    while int(current.degree) > 0:
        current = _trim_zero_boundaries(current, tol=tol)
        if int(current.degree) <= 0:
            break

        step += 1
        from_top = int(current.max_degree) > 0
        boundary_degree = int(current.max_degree if from_top else current.min_degree)
        side = "top" if from_top else "bottom"

        c_boundary = current.coefficient(boundary_degree)
        # ``extract_projector_from_rank1`` returns the null-space projector ``Q`` with
        # ``C_boundary @ Q ≈ 0``. ``carve_one_step`` expects ``P`` in
        # ``E_P^{-1}(z) = Q + z^{±1} P``, so use the complementary projector.
        q_boundary = snap_to_projector(extract_projector_from_rank1(c_boundary))
        projector = identity - q_boundary

        discarded_coeff_norm = float(np.max(np.abs(c_boundary @ q_boundary)))

        LOGGER.info(
            "Carving step %d: degree=%d, side=%s, discarded-coeff-norm=%.3e",
            step,
            int(current.degree),
            side,
            discarded_coeff_norm,
        )
        if discarded_coeff_norm > tol:
            LOGGER.warning(
                "Discarded coefficient norm %.3e exceeds tol %.3e; "
                "decomposition may be inaccurate.",
                discarded_coeff_norm,
                tol,
            )

        current = carve_one_step(current, projector, from_top=from_top)
        projectors.append(projector)
        carve_sides.append(from_top)

        if step % su2_check_every == 0 and LOGGER.isEnabledFor(logging.DEBUG):
            theta = np.linspace(0.0, 2.0 * np.pi, 128, endpoint=False, dtype=np.float64)
            z = np.exp(1j * theta)
            values = current(z)
            values_dag = np.conjugate(np.swapaxes(values, -1, -2))
            gram = np.matmul(values_dag, values)
            eye = np.eye(2, dtype=np.complex128)
            su2_err = float(np.max(np.abs(gram - eye)))
            det_err = float(np.max(np.abs(np.linalg.det(values) - 1.0)))
            LOGGER.debug(
                "SU(2) check at step %d: gram_err=%.3e, det_err=%.3e",
                step,
                su2_err,
                det_err,
            )

    e0_raw = current.coefficient(0)
    e0_unitary, _ = polar(e0_raw)

    return DecompositionResult(
        e0=e0_unitary,
        projectors=projectors,
        convention=getattr(matrix_poly, "convention", "gqsp"),
        carve_sides=carve_sides,
    )

decompose(p, q, convention='gqsp', tol=1e-10)

Unified decomposition entry point for supported QSP conventions.

Parameters:

Name Type Description Default
p LaurentPoly

Primary Laurent polynomial P.

 required
q LaurentPoly

Complementary Laurent polynomial / Fejér factor from completion.

 required
convention str

Convention selector ("gqsp", "laurent", or "auto").

 'gqsp'
tol float

Numerical tolerance used by decomposition routines and verification.

 1e-10

Returns:

Type Description
DecompositionResult

Recursive carving decomposition result.

Raises:

Type Description
ValueError

If convention is unsupported.
Source code in src\qsp_proc\decomposition\carving.py 
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
def decompose(
    p: LaurentPoly,
    q: LaurentPoly,
    convention: str = "gqsp",
    tol: float = 1e-10,
) -> DecompositionResult:
    """Unified decomposition entry point for supported QSP conventions.

    Args:
        p: Primary Laurent polynomial ``P``.
        q: Complementary Laurent polynomial / Fejér factor from completion.
        convention: Convention selector (``"gqsp"``, ``"laurent"``, or ``"auto"``).
        tol: Numerical tolerance used by decomposition routines and verification.

    Returns:
        Recursive carving decomposition result.

    Raises:
        ValueError: If ``convention`` is unsupported.
    """
    conv = convention.lower()
    if conv not in {"gqsp", "laurent", "auto"}:
        raise ValueError(
            f"Unsupported convention '{convention}'. "
            "Expected one of {'gqsp', 'laurent', 'auto'}."
        )

    # ``auto`` defaults to the most general currently supported pathway.
    target_conv = "gqsp" if conv == "auto" else conv

    if target_conv == "gqsp":
        matrix_poly = build_matrix_poly_gqsp(p, q)
    else:
        from qsp_proc.polynomials.chebyshev import ChebyshevPoly

        first_kind, second_kind = p.to_appendix_a_components()
        a_laurent = LaurentPoly.from_chebyshev(first_kind)
        b_laurent = LaurentPoly.from_appendix_a_components(
            first_kind=ChebyshevPoly([0.0]),
            second_kind_coeffs=second_kind,
        )
        gamma = q
        matrix_poly = build_matrix_poly_laurent_qsp(a_laurent, b_laurent, gamma)

    # Decompose the convention-built matrix polynomial directly in the original
    # Laurent variable so reconstruction/verification map back to the same P, Q.
    result = recursive_carve(matrix_poly, tol=tol)
    result.convention = target_conv

    if conv == "auto":
        max_error = verify_decomposition(
            p=p,
            q=q,
            result=result,
            convention="gqsp",
            tol=tol,
        )
        if max_error > tol:
            LOGGER.warning(
                "Decomposition verification error %.3e exceeds tol %.3e in auto mode.",
                max_error,
                tol,
            )

    return result

reconstruct_from_decomposition(result, num_samples=256)

Reconstruct F(z) from E_0 and carved projectors.

Uses F(z) = E_0 @ prod_{j=1}^{m} E_{P_j}(z), where E_{P_j}(z) = (I - P_j) + z^{s_j} P_j, where exponents follow result.carve_sides (top carve -> s_j = +1, bottom carve -> s_j = -1). If carve_sides is unavailable (legacy results), parity fallback is used.

Parameters:

Name Type Description Default
result DecompositionResult

Decomposition output containing e0 and projector sequence.

 required
num_samples int

Reserved for API compatibility and optional future checks.

 256

Returns:

Type Description
MatrixLaurentPoly

Reconstructed matrix-valued Laurent polynomial.
Source code in src\qsp_proc\decomposition\carving.py 
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
def reconstruct_from_decomposition(
    result: DecompositionResult, num_samples: int = 256
) -> MatrixLaurentPoly:
    """Reconstruct ``F(z)`` from ``E_0`` and carved projectors.

    Uses
    ``F(z) = E_0 @ prod_{j=1}^{m} E_{P_j}(z)``, where
    ``E_{P_j}(z) = (I - P_j) + z^{s_j} P_j``, where exponents follow
    ``result.carve_sides`` (top carve -> ``s_j = +1``, bottom carve -> ``s_j = -1``).
    If ``carve_sides`` is unavailable (legacy results), parity fallback is used.

    Args:
        result: Decomposition output containing ``e0`` and projector sequence.
        num_samples: Reserved for API compatibility and optional future checks.

    Returns:
        Reconstructed matrix-valued Laurent polynomial.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")

    e0 = np.asarray(result.e0, dtype=np.complex128)
    if e0.shape != (2, 2):
        raise ValueError("result.e0 must have shape (2, 2).")

    reconstructed = MatrixLaurentPoly(e0[np.newaxis, ...], min_degree=0)
    identity = np.eye(2, dtype=np.complex128)

    projectors = result.projectors
    carve_sides = result.carve_sides
    if carve_sides is not None and len(carve_sides) != len(projectors):
        raise ValueError("result.carve_sides must have the same length as projectors.")

    if carve_sides is None:
        # Backward-compatibility fallback for older decomposition outputs.
        carve_sides = [idx % 2 == 0 for idx in range(len(projectors))]

    for idx in range(len(projectors) - 1, -1, -1):
        projector = projectors[idx]
        from_top = carve_sides[idx]
        proj = np.asarray(projector, dtype=np.complex128)
        if proj.shape != (2, 2):
            raise ValueError("Each projector must have shape (2, 2).")

        q_proj = identity - proj

        if from_top:
            # If carved with E_P^{-1}(z)=Q+z^{-1}P, multiply back by E_P(z)=Q+zP.
            factor_coeffs = np.stack([q_proj, proj], axis=0)
            factor = MatrixLaurentPoly(factor_coeffs, min_degree=0)
        else:
            # If carved with E_P^{-1}(z)=Q+zP, multiply back by E_P(z)=Q+z^{-1}P.
            factor_coeffs = np.stack([proj, q_proj], axis=0)
            factor = MatrixLaurentPoly(factor_coeffs, min_degree=-1)

        reconstructed = reconstructed @ factor

    return reconstructed

verify_decomposition(p, q, result, convention='gqsp', num_samples=512, tol=1e-08)

Verify a decomposition by sampled operator-norm reconstruction error.

This computes: max_z ||F(z) - F_recon(z)||_op over uniformly sampled points on U(1).

Parameters:

Name Type Description Default
p LaurentPoly

Primary Laurent polynomial P.

 required
q LaurentPoly

Complementary Laurent polynomial Q.

 required
result DecompositionResult

Decomposition result.

 required
convention str

QSP convention selector; currently only "gqsp" is supported.

 'gqsp'
num_samples int

Number of unit-circle samples.

 512
tol float

Tolerance used for optional debug logging threshold.

 1e-08

Returns:

Type Description
float

Maximum sampled operator-norm error.
Source code in src\qsp_proc\decomposition\carving.py 
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
476
477
478
479
480
def verify_decomposition(
    p: LaurentPoly,
    q: LaurentPoly,
    result: DecompositionResult,
    convention: str = "gqsp",
    num_samples: int = 512,
    tol: float = 1e-8,
) -> float:
    """Verify a decomposition by sampled operator-norm reconstruction error.

    This computes:
    ``max_z ||F(z) - F_recon(z)||_op`` over uniformly sampled points on ``U(1)``.

    Args:
        p: Primary Laurent polynomial ``P``.
        q: Complementary Laurent polynomial ``Q``.
        result: Decomposition result.
        convention: QSP convention selector; currently only ``"gqsp"`` is supported.
        num_samples: Number of unit-circle samples.
        tol: Tolerance used for optional debug logging threshold.

    Returns:
        Maximum sampled operator-norm error.
    """
    if num_samples <= 0:
        raise ValueError("num_samples must be positive.")
    if tol < 0.0:
        raise ValueError("tol must be non-negative.")

    conv = convention.lower()
    if conv != "gqsp":
        raise NotImplementedError(
            f"verify_decomposition currently supports only 'gqsp', got '{convention}'."
        )

    f_original = build_matrix_poly_gqsp(p, q)
    f_reconstructed = reconstruct_from_decomposition(result, num_samples=num_samples)

    theta = np.linspace(0.0, 2.0 * np.pi, num_samples, endpoint=False, dtype=np.float64)
    z = np.exp(1j * theta)

    diff = f_original(z) - f_reconstructed(z)
    op_norms = np.linalg.norm(diff, ord=2, axis=(-2, -1))
    max_error = float(np.max(op_norms))

    if max_error > tol:
        LOGGER.debug(
            "Decomposition verification error %.3e exceeds tol %.3e",
            max_error,
            tol,
        )

    return max_error

qsp_proc.decomposition.builders

Convention-specific construction of matrix-valued Laurent polynomials for decomposition.

build_matrix_poly_gqsp(p, q)

Build G-QSP matrix polynomial F = [[P, Q], [-Q_tilde, P_tilde]].

For P(z) = sum_k p_k z^k and Q(z) = sum_k q_k z^k, this constructs coefficient blocks C_k = [[p_k, q_k], [-conj(q_{-k}), conj(p_{-k})]].

Parameters:

Name Type Description Default
p LaurentPoly

Laurent polynomial P.

 required
q LaurentPoly

Laurent polynomial Q.

 required

Returns:

Type Description
MatrixLaurentPoly

Matrix-valued Laurent polynomial F.
Source code in src\qsp_proc\decomposition\builders.py 
29
30
31
32
33
34
35
36
37
38
39
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
def build_matrix_poly_gqsp(p: LaurentPoly, q: LaurentPoly) -> MatrixLaurentPoly:
    """Build G-QSP matrix polynomial ``F = [[P, Q], [-Q_tilde, P_tilde]]``.

    For
    ``P(z) = sum_k p_k z^k`` and ``Q(z) = sum_k q_k z^k``, this constructs
    coefficient blocks
    ``C_k = [[p_k, q_k], [-conj(q_{-k}), conj(p_{-k})]]``.

    Args:
        p: Laurent polynomial ``P``.
        q: Laurent polynomial ``Q``.

    Returns:
        Matrix-valued Laurent polynomial ``F``.
    """
    p_min, p_max = _support_bounds(p)
    q_min, q_max = _support_bounds(q)

    min_degree = min(p_min, q_min, -p_max, -q_max)
    max_degree = max(p_max, q_max, -p_min, -q_min)

    num_coeffs = max_degree - min_degree + 1
    coeff_blocks = np.zeros((num_coeffs, 2, 2), dtype=np.complex128)

    for k in range(min_degree, max_degree + 1):
        idx = k - min_degree
        p_k = _coeff_at(p, k)
        q_k = _coeff_at(q, k)
        p_minus_k = _coeff_at(p, -k)
        q_minus_k = _coeff_at(q, -k)

        coeff_blocks[idx, 0, 0] = p_k
        coeff_blocks[idx, 0, 1] = q_k
        coeff_blocks[idx, 1, 0] = -np.conjugate(q_minus_k)
        coeff_blocks[idx, 1, 1] = np.conjugate(p_minus_k)

    matrix_poly = MatrixLaurentPoly(coeff_blocks, min_degree=min_degree)

    if not matrix_poly.is_su2_on_circle():
        warnings.warn(
            "Constructed G-QSP matrix polynomial is not approximately SU(2) on U(1).",
            RuntimeWarning,
            stacklevel=2,
        )

    return matrix_poly

build_matrix_poly_laurent_qsp(a, b, gamma)

Build Laurent-QSP matrix polynomial from A, B and gamma.

Implements Appendix D.1.2 (Skelton 2025): F(z) = A(z)I + iB(z)sigma_X + iC(z)sigma_Y + iD(z)sigma_Z, where (C, D) = gamma_to_cd(gamma). Equivalently, F = [[A + iD, iB + C], [iB - C, A - iD]].

Parameters:

Name Type Description Default
a LaurentPoly

Laurent polynomial A.

 required
b LaurentPoly

Laurent polynomial B.

 required
gamma LaurentPoly

Completion output Laurent polynomial gamma.

 required

Returns:

Type Description
MatrixLaurentPoly

Matrix-valued Laurent polynomial F.
Source code in src\qsp_proc\decomposition\builders.py 
 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
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
def build_matrix_poly_laurent_qsp(
    a: LaurentPoly, b: LaurentPoly, gamma: LaurentPoly
) -> MatrixLaurentPoly:
    """Build Laurent-QSP matrix polynomial from ``A, B`` and ``gamma``.

    Implements Appendix D.1.2 (Skelton 2025):
    ``F(z) = A(z)I + iB(z)sigma_X + iC(z)sigma_Y + iD(z)sigma_Z``, where
    ``(C, D) = gamma_to_cd(gamma)``. Equivalently,
    ``F = [[A + iD, iB + C], [iB - C, A - iD]]``.

    Args:
        a: Laurent polynomial ``A``.
        b: Laurent polynomial ``B``.
        gamma: Completion output Laurent polynomial ``gamma``.

    Returns:
        Matrix-valued Laurent polynomial ``F``.
    """
    c, d = gamma_to_cd(gamma)
    a_min, a_max = _support_bounds(a)
    b_min, b_max = _support_bounds(b)
    c_min, c_max = _support_bounds(c)
    d_min, d_max = _support_bounds(d)

    min_degree = min(a_min, b_min, c_min, d_min)
    max_degree = max(a_max, b_max, c_max, d_max)

    coeff_blocks = np.zeros((max_degree - min_degree + 1, 2, 2), dtype=np.complex128)

    # Coefficient block for z^k:
    # [[a_k + i d_k, i b_k + c_k], [i b_k - c_k, a_k - i d_k]]
    for k in range(min_degree, max_degree + 1):
        idx = k - min_degree
        a_k = _coeff_at(a, k)
        b_k = _coeff_at(b, k)
        c_k = _coeff_at(c, k)
        d_k = _coeff_at(d, k)

        coeff_blocks[idx, 0, 0] = a_k + 1j * d_k
        coeff_blocks[idx, 0, 1] = 1j * b_k + c_k
        coeff_blocks[idx, 1, 0] = 1j * b_k - c_k
        coeff_blocks[idx, 1, 1] = a_k - 1j * d_k

    matrix_poly = MatrixLaurentPoly(coeff_blocks, min_degree=min_degree)

    if not matrix_poly.is_su2_on_circle():
        warnings.warn(
            "Constructed Laurent-QSP matrix polynomial is not approximately SU(2) on U(1).",
            RuntimeWarning,
            stacklevel=2,
        )

    return matrix_poly

gamma_to_cd(gamma)

Convert gamma into reciprocal/anti-reciprocal Laurent polynomials.

Implements Appendix D.1.1 (Skelton 2025): C(z) = (gamma(z) + gamma(1/z)) / 2, D(z) = (gamma(z) - gamma(1/z)) / (2i).

If gamma(z) = sum_k g_k z^k, then c_k = (g_k + g_-k)/2 and d_k = (g_k - g_-k)/(2i).

Parameters:

Name Type Description Default
gamma LaurentPoly

Input Laurent polynomial gamma.

 required

Returns:

Type Description
tuple[LaurentPoly, LaurentPoly]

Tuple (C, D) as Laurent polynomials.
Source code in src\qsp_proc\decomposition\builders.py 
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
def gamma_to_cd(gamma: LaurentPoly) -> tuple[LaurentPoly, LaurentPoly]:
    """Convert ``gamma`` into reciprocal/anti-reciprocal Laurent polynomials.

    Implements Appendix D.1.1 (Skelton 2025):
    ``C(z) = (gamma(z) + gamma(1/z)) / 2``,
    ``D(z) = (gamma(z) - gamma(1/z)) / (2i)``.

    If ``gamma(z) = sum_k g_k z^k``, then
    ``c_k = (g_k + g_-k)/2`` and ``d_k = (g_k - g_-k)/(2i)``.

    Args:
        gamma: Input Laurent polynomial ``gamma``.

    Returns:
        Tuple ``(C, D)`` as Laurent polynomials.
    """
    g_min, g_max = _support_bounds(gamma)

    cd_min = min(g_min, -g_max)
    cd_max = max(g_max, -g_min)

    c_coeffs = np.zeros(cd_max - cd_min + 1, dtype=np.complex128)
    d_coeffs = np.zeros(cd_max - cd_min + 1, dtype=np.complex128)

    for k in range(cd_min, cd_max + 1):
        idx = k - cd_min
        g_k = _coeff_at(gamma, k)
        g_minus_k = _coeff_at(gamma, -k)
        c_coeffs[idx] = 0.5 * (g_k + g_minus_k)
        d_coeffs[idx] = (g_k - g_minus_k) / (2.0j)

    return LaurentPoly(c_coeffs, min_degree=cd_min), LaurentPoly(d_coeffs, min_degree=cd_min)

pad_to_half_degree(F)

Apply Algorithm 2 zero-padding: F(z)=sum_k C_k z^k -> F'(w)=sum_k C_k w^{2k}.

This doubles the Laurent exponent grid by inserting zero 2x2 blocks at odd exponents. If n = degree(F), the output support is exactly [-2n, 2n] with coefficient-array shape (4n + 1, 2, 2).

Parameters:

Name Type Description Default
F MatrixLaurentPoly

Input matrix-valued Laurent polynomial.

 required

Returns:

Type Description
MatrixLaurentPoly

Re-parameterized matrix-valued Laurent polynomial F'.
Source code in src\qsp_proc\decomposition\builders.py 
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
def pad_to_half_degree(F: MatrixLaurentPoly) -> MatrixLaurentPoly:
    """Apply Algorithm 2 zero-padding: ``F(z)=sum_k C_k z^k -> F'(w)=sum_k C_k w^{2k}``.

    This doubles the Laurent exponent grid by inserting zero ``2x2`` blocks at
    odd exponents. If ``n = degree(F)``, the output support is exactly
    ``[-2n, 2n]`` with coefficient-array shape ``(4n + 1, 2, 2)``.

    Args:
        F: Input matrix-valued Laurent polynomial.

    Returns:
        Re-parameterized matrix-valued Laurent polynomial ``F'``.
    """
    matrix_poly = F
    n = int(matrix_poly.degree)
    padded_min_degree = -2 * n
    padded_coeffs = np.zeros((4 * n + 1, 2, 2), dtype=np.complex128)

    source_coeffs = matrix_poly.coeffs
    source_min_degree = matrix_poly.min_degree

    for idx, coeff_block in enumerate(source_coeffs):
        k = source_min_degree + idx
        target_degree = 2 * k
        target_idx = target_degree - padded_min_degree
        padded_coeffs[target_idx] = coeff_block

    return MatrixLaurentPoly(padded_coeffs, min_degree=padded_min_degree)

qsp_proc.decomposition.matrix_laurent_poly

Matrix-valued Laurent polynomial storage and carving output types.

DecompositionResult dataclass

Output of recursive carving: constant residue, projectors, and carve metadata.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
11
12
13
14
15
16
17
18
@dataclass
class DecompositionResult:
    """Output of recursive carving: constant residue, projectors, and carve metadata."""

    e0: np.ndarray
    projectors: list[np.ndarray]
    convention: str
    carve_sides: list[bool] | None = None

MatrixLaurentPoly

Two-by-two matrix-valued Laurent polynomial F(z) = sum_k C_k z^k.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 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
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
class MatrixLaurentPoly:
    """Two-by-two matrix-valued Laurent polynomial ``F(z) = sum_k C_k z^k``."""

    def __init__(self, coeffs: np.ndarray, min_degree: int) -> None:
        """Initialize from dense ``(num_coeffs, 2, 2)`` blocks and minimum Laurent exponent."""
        if not isinstance(min_degree, Integral):
            raise TypeError("min_degree must be an integer.")

        coeffs_arr = np.asarray(coeffs, dtype=np.complex128)
        if coeffs_arr.ndim != 3 or coeffs_arr.shape[1:] != (2, 2):
            raise ValueError("coeffs must have shape (num_coeffs, 2, 2).")
        if coeffs_arr.shape[0] <= 0:
            raise ValueError("coeffs must contain at least one 2x2 coefficient matrix.")

        self._coeffs = coeffs_arr
        self._min_degree = int(min_degree)

    @property
    def coeffs(self) -> np.ndarray:
        """Dense coefficient array (copy).

        Same encapsulation pattern as :class:`~qsp_proc.polynomials.laurent.LaurentPoly`.
        """
        return self._coeffs.copy()

    @property
    def min_degree(self) -> int:
        """Minimum Laurent exponent with a stored block."""
        return self._min_degree

    @property
    def max_degree(self) -> int:
        """Maximum Laurent exponent with a stored block."""
        return self._min_degree + self._coeffs.shape[0] - 1

    @property
    def exponents(self) -> np.ndarray:
        """Stored Laurent exponents as a dense integer range."""
        return np.arange(self.min_degree, self.max_degree + 1, dtype=np.int64)

    @property
    def degree(self) -> int:
        """Half-width degree ``max(|min_degree|, |max_degree|)``."""
        return max(abs(self.min_degree), abs(self.max_degree))

    def coefficient(self, k: int) -> np.ndarray:
        """Return coefficient matrix ``C_k``, or a ``2x2`` zero matrix if outside support."""
        idx = int(k) - self._min_degree
        if 0 <= idx < self._coeffs.shape[0]:
            return self._coeffs[idx]
        return np.zeros((2, 2), dtype=np.complex128)

    def evaluate(self, z: np.ndarray | complex) -> np.ndarray:
        """Evaluate at scalar or batched ``z`` with trailing output shape ``(2, 2)``."""
        z_arr = np.asarray(z, dtype=np.complex128)
        powers = z_arr[..., None] ** self.exponents
        values = np.asarray(
            np.einsum("...k,kij->...ij", powers, self._coeffs, optimize=True),
            dtype=np.complex128,
        )
        if z_arr.ndim == 0:
            return values.reshape(2, 2)
        return values

    __call__ = evaluate

    def __matmul__(self, other: object) -> MatrixLaurentPoly:
        """Multiply two matrix-valued Laurent polynomials via coefficient convolution."""
        if not isinstance(other, MatrixLaurentPoly):
            return NotImplemented

        left = self._coeffs
        right = other._coeffs
        out_len = left.shape[0] + right.shape[0] - 1
        out = np.zeros((out_len, 2, 2), dtype=np.complex128)

        for i, left_block in enumerate(left):
            out[i : i + right.shape[0]] += left_block @ right

        return MatrixLaurentPoly(out, min_degree=self.min_degree + other.min_degree)

    def is_su2_on_circle(self, num_samples: int = 512, tol: float = 1e-10) -> bool:
        """Whether ``F(e^{iθ})`` is approximately in ``SU(2)`` on uniform circle samples."""
        num_samples_int = int(num_samples)
        if num_samples_int < 1:
            raise ValueError("num_samples must be >= 1.")
        if tol < 0.0:
            raise ValueError("tol must be non-negative.")

        theta = np.linspace(0.0, 2.0 * np.pi, num_samples_int, endpoint=False)
        z = np.exp(1j * theta)
        fz = self.evaluate(z)
        fz_dag = np.conjugate(np.swapaxes(fz, -1, -2))
        gram = np.matmul(fz_dag, fz)
        eye = np.eye(2, dtype=np.complex128)
        unitary_err = float(np.max(np.abs(gram - eye)))
        det_err = float(np.max(np.abs(np.linalg.det(fz) - 1.0)))
        return unitary_err <= tol and det_err <= tol

coeffs property

Dense coefficient array (copy).

Same encapsulation pattern as :class:~qsp_proc.polynomials.laurent.LaurentPoly.

min_degree property

Minimum Laurent exponent with a stored block.

max_degree property

Maximum Laurent exponent with a stored block.

exponents property

Stored Laurent exponents as a dense integer range.

degree property

Half-width degree max(|min_degree|, |max_degree|).

__init__(coeffs, min_degree)

Initialize from dense (num_coeffs, 2, 2) blocks and minimum Laurent exponent.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
def __init__(self, coeffs: np.ndarray, min_degree: int) -> None:
    """Initialize from dense ``(num_coeffs, 2, 2)`` blocks and minimum Laurent exponent."""
    if not isinstance(min_degree, Integral):
        raise TypeError("min_degree must be an integer.")

    coeffs_arr = np.asarray(coeffs, dtype=np.complex128)
    if coeffs_arr.ndim != 3 or coeffs_arr.shape[1:] != (2, 2):
        raise ValueError("coeffs must have shape (num_coeffs, 2, 2).")
    if coeffs_arr.shape[0] <= 0:
        raise ValueError("coeffs must contain at least one 2x2 coefficient matrix.")

    self._coeffs = coeffs_arr
    self._min_degree = int(min_degree)

coefficient(k)

Return coefficient matrix C_k, or a 2x2 zero matrix if outside support.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
66
67
68
69
70
71
def coefficient(self, k: int) -> np.ndarray:
    """Return coefficient matrix ``C_k``, or a ``2x2`` zero matrix if outside support."""
    idx = int(k) - self._min_degree
    if 0 <= idx < self._coeffs.shape[0]:
        return self._coeffs[idx]
    return np.zeros((2, 2), dtype=np.complex128)

evaluate(z)

Evaluate at scalar or batched z with trailing output shape (2, 2).

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
73
74
75
76
77
78
79
80
81
82
83
def evaluate(self, z: np.ndarray | complex) -> np.ndarray:
    """Evaluate at scalar or batched ``z`` with trailing output shape ``(2, 2)``."""
    z_arr = np.asarray(z, dtype=np.complex128)
    powers = z_arr[..., None] ** self.exponents
    values = np.asarray(
        np.einsum("...k,kij->...ij", powers, self._coeffs, optimize=True),
        dtype=np.complex128,
    )
    if z_arr.ndim == 0:
        return values.reshape(2, 2)
    return values

__matmul__(other)

Multiply two matrix-valued Laurent polynomials via coefficient convolution.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def __matmul__(self, other: object) -> MatrixLaurentPoly:
    """Multiply two matrix-valued Laurent polynomials via coefficient convolution."""
    if not isinstance(other, MatrixLaurentPoly):
        return NotImplemented

    left = self._coeffs
    right = other._coeffs
    out_len = left.shape[0] + right.shape[0] - 1
    out = np.zeros((out_len, 2, 2), dtype=np.complex128)

    for i, left_block in enumerate(left):
        out[i : i + right.shape[0]] += left_block @ right

    return MatrixLaurentPoly(out, min_degree=self.min_degree + other.min_degree)

is_su2_on_circle(num_samples=512, tol=1e-10)

Whether F(e^{iθ}) is approximately in SU(2) on uniform circle samples.

Source code in src\qsp_proc\decomposition\matrix_laurent_poly.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def is_su2_on_circle(self, num_samples: int = 512, tol: float = 1e-10) -> bool:
    """Whether ``F(e^{iθ})`` is approximately in ``SU(2)`` on uniform circle samples."""
    num_samples_int = int(num_samples)
    if num_samples_int < 1:
        raise ValueError("num_samples must be >= 1.")
    if tol < 0.0:
        raise ValueError("tol must be non-negative.")

    theta = np.linspace(0.0, 2.0 * np.pi, num_samples_int, endpoint=False)
    z = np.exp(1j * theta)
    fz = self.evaluate(z)
    fz_dag = np.conjugate(np.swapaxes(fz, -1, -2))
    gram = np.matmul(fz_dag, fz)
    eye = np.eye(2, dtype=np.complex128)
    unitary_err = float(np.max(np.abs(gram - eye)))
    det_err = float(np.max(np.abs(np.linalg.det(fz) - 1.0)))
    return unitary_err <= tol and det_err <= tol