nlft_qsp.solvers.completion.weiss

Module for polynomial completion based on Weiss' algorithm.

Classes:

Name	Description
WeissConvergenceError	This error is thrown where Weiss' algorithm does not converge, more precisely,

Functions:

Name	Description
complete	Uses Weiss' algorithm to find a complementary polynomial to the given one. The polynomial will also be the unique outer, positive-mean polynomial with this property, according to arXiv:2407.05634.
ratio	Uses Weiss' algorithm to compute \(b/a\), where \(a\) is the unique outer, positive-mean polynomial such that \(|a|^2 + |b|^2 = 1\), up to working precision.
weiss_internal	Internal function for Weiss' algorithm. The user should call weiss.complete, or weiss.ratio.

WeissConvergenceError

Bases: Exception

This error is thrown where Weiss' algorithm does not converge, more precisely, when the error in Weiss' algorithm does not improve after WEISS_MAX_ATTEMPTS steps.

Source code in nlft_qsp/solvers/completion/weiss.py

class WeissConvergenceError(Exception):
    """This error is thrown where Weiss' algorithm does not converge, more precisely,
    when the error in Weiss' algorithm does not improve after `WEISS_MAX_ATTEMPTS` steps."""

    def __init__(self, *args):
        super().__init__(*args)

complete(b: Polynomial, eps: float = -1, verbose=False) -> Polynomial

Uses Weiss' algorithm to find a complementary polynomial to the given one. The polynomial will also be the unique outer, positive-mean polynomial with this property, according to arXiv:2407.05634.

Parameters:

Name	Type	Description	Default
b	Polynomial	The polynomial to complete.	required
eps	float	The desired tolerance. If not specified, it will be set to working precision.	-1
verbose	bool	verbosity during the procedure. Defaults to False.	False

Returns:

Type	Description
Polynomial	A polynomial \(a(z)\) satisfying \(|a|^2 + |b|^2 = 1\) on the unit circle (up to eps).

Source code in nlft_qsp/solvers/completion/weiss.py

def complete(b: Polynomial, eps:float=-1, verbose=False) -> Polynomial:
    """Uses Weiss' algorithm to find a complementary polynomial to the given one. The polynomial will also be the unique outer, positive-mean polynomial with this property, according to [arXiv:2407.05634](https://arxiv.org/abs/2407.05634).

    Args:
        b (Polynomial): The polynomial to complete.
        eps (float): The desired tolerance. If not specified, it will be set to working precision.
        verbose (bool, optional): verbosity during the procedure. Defaults to False.

    Returns:
        A polynomial $a(z)$ satisfying $|a|^2 + |b|^2 = 1$ on the unit circle (up to eps).
    """
    return weiss_internal(b, eps, verbose=verbose)

ratio(b: Polynomial, eps: float = -1, verbose=False) -> Polynomial

Uses Weiss' algorithm to compute \(b/a\), where \(a\) is the unique outer, positive-mean polynomial such that \(|a|^2 + |b|^2 = 1\), up to working precision.

Parameters:

Name	Type	Description	Default
b	Polynomial	The polynomial to complete.	required
eps	float	The desired tolerance. If not specified, it will be set to working precision.	-1
verbose	bool	verbosity during the procedure. Defaults to False.	False

Returns:

Type	Description
Polynomial	A polynomial \(a(z)\) satisfying \(|a(z)|^2 + |b(z)|^2 = 1\) on the unit circle (up to working precision), and a polynomial \(c\) that approximates \(b/a\).

Source code in nlft_qsp/solvers/completion/weiss.py

def ratio(b: Polynomial, eps:float=-1, verbose=False) -> Polynomial:
    """Uses Weiss' algorithm to compute $b/a$, where $a$ is the unique outer, positive-mean polynomial such that $|a|^2 + |b|^2 = 1$, up to working precision.

    Args:
        b (Polynomial): The polynomial to complete.
        eps (float): The desired tolerance. If not specified, it will be set to working precision.
        verbose (bool, optional): verbosity during the procedure. Defaults to False.

    Returns:
        A polynomial $a(z)$ satisfying $|a(z)|^2 + |b(z)|^2 = 1$ on the unit circle (up to working precision), and a polynomial $c$ that approximates $b/a$.
    """
    return weiss_internal(b, eps, compute_ratio=True, verbose=verbose)

weiss_internal(b: Polynomial, eps: float = -1, compute_ratio=False, verbose=False) -> Polynomial | tuple[Polynomial, Polynomial]

Internal function for Weiss' algorithm. The user should call weiss.complete, or weiss.ratio.

Parameters:

Name	Type	Description	Default
b	Polynomial	The starting polynomial to complete.	required
eps	float	The desired tolerance. If not specified, it will be set to working precision.	-1
compute_ratio	bool	If True, then also a Polynomial approximating \(b/a\) will be returned.	False
verbose	bool	verbosity during the procedure.	False

Returns:

Name	Type	Description
Polynomial	Polynomial | tuple[Polynomial, Polynomial]	A polynomial \(a(z)\) satisfying \(|a|^2 + |b|^2 = 1\) on the unit circle (up to working precision). If compute_ratio=True, a second polynomial \(c(z)\) that approximates \(b/a\) is returned.

Source code in nlft_qsp/solvers/completion/weiss.py

def weiss_internal(b: Polynomial, eps:float=-1, compute_ratio=False, verbose=False) -> Polynomial | tuple[Polynomial, Polynomial]:
    """Internal function for Weiss' algorithm. The user should call `weiss.complete`, or `weiss.ratio`.

    Args:
        b (Polynomial): The starting polynomial to complete.
        eps (float): The desired tolerance. If not specified, it will be set to working precision.
        compute_ratio (bool, optional): If True, then also a `Polynomial` approximating $b/a$ will be returned.
        verbose (bool, optional): verbosity during the procedure.

    Returns:
        Polynomial: A polynomial $a(z)$ satisfying $|a|^2 + |b|^2 = 1$ on the unit circle (up to working precision). If `compute_ratio=True`, a second polynomial $c(z)$ that approximates $b/a$ is returned.
    """
    d = b.effective_degree()
    if eps < 0:
        eps = bd.machine_threshold() * 10e-3

    eta = 1 - b.sup_norm(max(4000, 4*d))

    N = next_power_of_two(int(d/eta))//2 # Exponential search on N
    threshold = 1
    attempts = 0
    while threshold > eps:
        N *= 2
        if verbose:
            print(f"N = {N}")

        b_points = b.eval_at_roots_of_unity(N)

        R = laurent_approximation([np.log(1 - bz * np.conj(bz))/2 for bz in b_points])

        G = R.schwarz_transform()
        G_points = G.eval_at_roots_of_unity(N)

        a = laurent_approximation([np.exp(gz) for gz in G_points])
        a = a.truncate(-b.effective_degree(), 0) # a and b must have the same support

        new_thr = (a * a.conjugate() + b * b.conjugate() - 1).l2_norm()
        if verbose:
            print(f"N = {N:>7}, threshold = {new_thr}")

        if threshold <= new_thr:
            attempts += 1
            if attempts >= WEISS_MAX_ATTEMPTS:
                raise WeissConvergenceError()
        else:
            threshold = new_thr
            attempts = 0

    if compute_ratio:
        c = laurent_approximation([bz * np.exp(-gz) for bz, gz in zip(b_points, G_points)])
        return a, c.truncate(c.support_start, b.support().stop - 1)
    else:
        return a