Skip to content

nlft_qsp.poly

Definitions of complex sequences, polynomials and Chebyshev expansions, as well as all operations between them.

Classes:

Name Description
ChebyshevTExpansion

Linear combination of Chebyshev polynomials of the first kind.
ComplexL0Sequence

Represents a sequence of complex numbers or complex matrices index by \(\Z\), whose support is finite.
Polynomial

Represents a general Laurent polynomial of one complex variable.

ChebyshevTExpansion

Bases: ComplexL0Sequence

Linear combination of Chebyshev polynomials of the first kind.

Parameters:

Name Type Description Default
coeffs list[complex_type] | Polynomial

Either the coefficients of the linear combination, or the symmetric Laurent polynomial \(P(z)\) which is equal up to the change of variable \(x = \frac{z + z^{-1}}{2}\).

 required
Note

The change of variable between Laurent polynomials and Chebyshev expansions is given by the relation \(T_k(\cos \theta) = \cos k\theta = \frac{z^k + z^{-k}}{2}\). Thus the expansion is

\[\sum_{k = 0}^n c_k T_k(x) = \sum_{k = 0}^n c_k \frac{z^k + z^{-k}}{2}\]

Methods:

Name Description
__call__

Evaluates the Chebyshev expansion at the given number.
__str__

Converts the expansion to a human-readable string representation.
from_laurent_polynomial

Returns the Chebyshev expansion \(T\) satisfying \(T(x) = \frac{P(z) + P^*(z)}{2}\).
from_polynomial

Returns the Chebyshev expansion \(T\) satisfying \(T(x) = P(x)\).
to_laurent

Returns the Laurent polynomial \(P(e^{i\theta}) = T(\cos \theta)\) where \(T(x)\) is represented by self.
to_polynomial

Returns the polynomial \(P\) satisfying \(P(x) = T(x)\).
Source code in nlft_qsp/poly.py 
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
@serializable(
    type_tag="@qspx/chebyshev_t",
    fields={
        "coeffs": "coeffs"
    }
)
class ChebyshevTExpansion(ComplexL0Sequence):
    r"""Linear combination of Chebyshev polynomials of the first kind.

    Args:
        coeffs: Either the coefficients of the linear combination, or the symmetric Laurent polynomial $P(z)$ which is equal up to the change of variable $x = \frac{z + z^{-1}}{2}$.

    Note:
        The change of variable between Laurent polynomials and Chebyshev expansions is given by the relation $T_k(\cos \theta) = \cos k\theta = \frac{z^k + z^{-k}}{2}$. Thus the expansion is

        $$\sum_{k = 0}^n c_k T_k(x) = \sum_{k = 0}^n c_k \frac{z^k + z^{-k}}{2}$$
    """
    def __init__(self, coeffs: list[complex_type] | Polynomial):
        if isinstance(coeffs, list) or isinstance(coeffs, np.ndarray):
            super().__init__(coeffs, support_start=0)
        elif isinstance(coeffs, Polynomial):
            if not coeffs.is_symmetric():
                raise ValueError("The given Laurent polynomial is not symmetric.")

            coeffs = [2*coeffs[k] for k in range(coeffs.support().stop)]
            coeffs[0] /= 2
            super().__init__(coeffs, support_start=0)
        else:
            raise ValueError("Only a coefficient vector or symmetric Laurent polynomials are allowed.")

    def degree(self) -> int:
        return self.coeffs.shape[0] - 1

    def __call__(self, x: float_type) -> complex_type:
        """Evaluates the Chebyshev expansion at the given number.

        Args:
            x (float_type): The point at which to evaluate the expansion.

        Returns:
            complex: The evaluated result.
        """
        theta = np.arccos(x)

        return sum(self[k] * np.cos(k * theta) for k in self.support())

    def __str__(self) -> str:
        """Converts the expansion to a human-readable string representation.

        Returns:
            str: The string representation of the expansion.
        """
        return ' + '.join(f"{c} T_{self.support_start + k}(x)" for k, c in enumerate(self.coeffs))

    def to_laurent(self) -> Polynomial:
        r"""Returns the Laurent polynomial $P(e^{i\theta}) = T(\cos \theta)$ where $T(x)$ is represented by `self`."""
        P = Polynomial(np.concat([np.flipud(self.coeffs), self.coeffs[1:]]), support_start=-self.coeffs.shape[0]+1)
        P[0] *= 2
        return P/2

    @classmethod
    def from_polynomial(cls, P: Polynomial) -> "ChebyshevTExpansion":
        r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = P(x)$."""
        return ChebyshevTExpansion(np.polynomial.chebyshev.poly2cheb(P.coeffs))

    @classmethod
    def from_laurent_polynomial(cls, P: Polynomial):
        r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = \frac{P(z) + P^*(z)}{2}$.

        Raises:
            ValueError: if $P$ is not symmetric.

        Note: $P$ must be symmetric."""
        if not P.is_symmetric():
            raise ValueError("The given Laurent polynomial is not symmetric.")

        coeffs = [2*P[k] for k in range(P.support().stop)]
        coeffs[0] /= 2

        return ChebyshevTExpansion(coeffs)

    def to_polynomial(self) -> Polynomial:
        """Returns the polynomial $P$ satisfying $P(x) = T(x)$."""
        return Polynomial(np.polynomial.chebyshev.cheb2poly(self.coeffs))

__call__(x: float_type) -> complex_type

Evaluates the Chebyshev expansion at the given number.

Parameters:

Name Type Description Default
x float_type

The point at which to evaluate the expansion.

 required

Returns:

Name Type Description
complex complex_type

The evaluated result.
Source code in nlft_qsp/poly.py 
541
542
543
544
545
546
547
548
549
550
551
552
def __call__(self, x: float_type) -> complex_type:
    """Evaluates the Chebyshev expansion at the given number.

    Args:
        x (float_type): The point at which to evaluate the expansion.

    Returns:
        complex: The evaluated result.
    """
    theta = np.arccos(x)

    return sum(self[k] * np.cos(k * theta) for k in self.support())

__str__() -> str

Converts the expansion to a human-readable string representation.

Returns:

Name Type Description
str str

The string representation of the expansion.
Source code in nlft_qsp/poly.py 
554
555
556
557
558
559
560
def __str__(self) -> str:
    """Converts the expansion to a human-readable string representation.

    Returns:
        str: The string representation of the expansion.
    """
    return ' + '.join(f"{c} T_{self.support_start + k}(x)" for k, c in enumerate(self.coeffs))

from_laurent_polynomial(P: Polynomial) classmethod

Returns the Chebyshev expansion \(T\) satisfying \(T(x) = \frac{P(z) + P^*(z)}{2}\).

Raises:

Type Description
ValueError

if \(P\) is not symmetric.

Note: \(P\) must be symmetric.

Source code in nlft_qsp/poly.py 
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
@classmethod
def from_laurent_polynomial(cls, P: Polynomial):
    r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = \frac{P(z) + P^*(z)}{2}$.

    Raises:
        ValueError: if $P$ is not symmetric.

    Note: $P$ must be symmetric."""
    if not P.is_symmetric():
        raise ValueError("The given Laurent polynomial is not symmetric.")

    coeffs = [2*P[k] for k in range(P.support().stop)]
    coeffs[0] /= 2

    return ChebyshevTExpansion(coeffs)

from_polynomial(P: Polynomial) -> ChebyshevTExpansion classmethod

Returns the Chebyshev expansion \(T\) satisfying \(T(x) = P(x)\).

Source code in nlft_qsp/poly.py 
568
569
570
571
@classmethod
def from_polynomial(cls, P: Polynomial) -> "ChebyshevTExpansion":
    r"""Returns the Chebyshev expansion $T$ satisfying $T(x) = P(x)$."""
    return ChebyshevTExpansion(np.polynomial.chebyshev.poly2cheb(P.coeffs))

to_laurent() -> Polynomial

Returns the Laurent polynomial \(P(e^{i\theta}) = T(\cos \theta)\) where \(T(x)\) is represented by self.

Source code in nlft_qsp/poly.py 
562
563
564
565
566
def to_laurent(self) -> Polynomial:
    r"""Returns the Laurent polynomial $P(e^{i\theta}) = T(\cos \theta)$ where $T(x)$ is represented by `self`."""
    P = Polynomial(np.concat([np.flipud(self.coeffs), self.coeffs[1:]]), support_start=-self.coeffs.shape[0]+1)
    P[0] *= 2
    return P/2

to_polynomial() -> Polynomial

Returns the polynomial \(P\) satisfying \(P(x) = T(x)\).

Source code in nlft_qsp/poly.py 
589
590
591
def to_polynomial(self) -> Polynomial:
    """Returns the polynomial $P$ satisfying $P(x) = T(x)$."""
    return Polynomial(np.polynomial.chebyshev.cheb2poly(self.coeffs))

ComplexL0Sequence

Represents a sequence of complex numbers or complex matrices index by \(\Z\), whose support is finite.

Attributes:

Name Type Description
coeffs list[complex_type]

List of complex coefficients.
support_start int

Index of the first element of the sequence.
shape tuple[int]

The shape of the coefficients. This is always () for polynomials, but is used for matrix polynomials.

Methods:

Name Description
__getitem__

Returns the \(k\)-th element of the sequence, i.e., \(F_k\).
__init__

Initializes a complex sequence.
__setitem__

Sets the coefficient of \(z^k\) to be \(c\), allocating space if needed.
is_imaginary

Whether the sequence has only imaginary elements.
is_real

Whether the sequence has only real elements.
is_symmetric

Whether the sequence satisfies \(F_k = F_{-k}\).
l2_norm

Computes the \(\ell_2\) norm.
l2_squared_norm

Computes the squared \(\ell_2\) norm.
support

Returns the range in \(\Z\) where the sequence is non-zero.
Source code in nlft_qsp/poly.py 
 16
 17
 18
 19
 20
 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
119
120
121
122
123
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
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
class ComplexL0Sequence:
    r"""Represents a sequence of complex numbers or complex matrices index by $\Z$, whose support is finite.

    Attributes:
        coeffs (list[complex_type]): List of complex coefficients.
        support_start (int): Index of the first element of the sequence.
        shape (tuple[int]): The shape of the coefficients. This is always () for polynomials, but is used for matrix polynomials.
    """

    def __init__(self, coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None):
        """Initializes a complex sequence.

        Args:
            coeffs: List of complex numbers as coefficients. `coeffs[k, :]` will be the coefficient of $z^(s + k)$ where $s$ is `support_start`.
            shape: The shape of the sequence. If it is none, then the shape of coeffs is taken. This is only used for matrix sequences, and is ignored otherwise.
            support_start (optional): Index of the first element of the sequence. Defaults to $0$.

        Note:
            Please provide exactly one of coeffs or shape.
        """
        if coeffs is not None and shape is not None or coeffs is None and shape is None:
            raise ValueError("Please provide exactly one between coeffs or shape.")

        self.support_start = support_start

        if shape is not None:
            self.coeffs = np.zeros((1, *shape), dtype=complex_type)
            self.shape = shape
            return

        if isinstance(coeffs, list):
            if any(not isinstance(ck, Number) for ck in coeffs):
                raise ValueError(f"The list of coefficients must be a single list or np.ndarray of complex values. Got {coeffs}. If you want to create a matrix sequence, please provide a np.ndarray of shape (num_coeffs, *shape) or specify the shape argument.")

            self.coeffs = np.array(coeffs, dtype=complex_type)
            self.shape = ()
            return

        if isinstance(coeffs, np.ndarray):
            if coeffs.ndim == 0:
                raise ValueError(f"There must be at least one coefficient. Got {coeffs}.")

            self.coeffs = np.array(coeffs, dtype=complex_type)
            self.shape = coeffs.shape[1:]
            return

        raise ValueError(f"Coefficients must be provided as a list or np.ndarray. Got {coeffs}.")

    def support(self) -> range:
        r"""Returns the range in $\Z$ where the sequence is non-zero.

        Note:
            This simply checks the allocated array of coefficients but does not check leading or trailing zeros.
        """
        return range(self.support_start, self.support_start + self.coeffs.shape[0])

    def __getitem__(self, k: int | tuple) -> complex_type | np.ndarray:
        r"""Returns the $k$-th element of the sequence, i.e., $F_k$.

        Args:
            k (int): The index of the sequence.

        Returns:
            complex_type | np.ndarray: The coefficient $F_k$, or $0$ if $k$ is out of the support.

        Note: in case of matrix sequence, numpy slices can be used on the coefficients.
        """
        if not isinstance(k, tuple):
            k = (k,)

        supp = self.support()
        if isinstance(k[0], slice):
            sl1 = max(k[0].start, supp.start) - self.support_start if k[0].start else 0
            sl2 = min(k[0].stop, supp.stop) - self.support_start if k[0].stop else self.coeffs.shape[0]

            return Polynomial(self.coeffs[(slice(sl1, sl2, k[0].step), *k[1:])], support_start=supp.start + sl1)

        if k[0] in supp:
            return self.coeffs[(k[0] - self.support_start, *k[1:])]

        if self.shape != ():
            return np.zeros(self.shape, dtype=complex_type)[k[1:]]

        return 0

    def __setitem__(self, k: int, c: complex_type | np.ndarray):
        r"""Sets the coefficient of $z^k$ to be $c$, allocating space if needed.

        Args:
            k (int): The exponent of $z$.
            c (complex): The coefficient to set.
        """
        if isinstance(k, int):
            k = (k,)

        if self.support_start + self.coeffs.shape[0] <= k[0]:
            pad_width = [(0, k[0] - self.support_start - self.coeffs.shape[0] + 1)] + [(0, 0)] * (self.coeffs.ndim - 1)
            self.coeffs = np.pad(self.coeffs, pad_width)
        elif self.support_start > k[0]:
            pad_width = [(self.support_start - k[0], 0)] + [(0, 0)] * (self.coeffs.ndim - 1)
            self.coeffs = np.pad(self.coeffs, pad_width)
            self.support_start = k[0]

        self.coeffs[(k[0] - self.support_start, *k[1:])] = c

    def l2_norm(self) -> float_type:
        r"""Computes the $\ell_2$ norm."""
        return np.sqrt(self.l2_squared_norm())

    def l2_squared_norm(self) -> float_type:
        r"""Computes the squared $\ell_2$ norm.

        Note:
            For matrix sequences, this is the sum of the squared Frobenius norms of the coefficients.
        """
        if self.shape == ():
            return sum(c * np.conj(c) for c in self.coeffs)

        return sum(np.linalg.norm(c)**2 for c in self.coeffs)

    def is_real(self) -> bool:
        """Whether the sequence has only real elements."""
        if self.shape == ():
            return all(np.abs(np.imag(F)) <= bd.machine_threshold() for F in self.coeffs)

        return all(np.abs(np.imag(F)) <= bd.machine_threshold() for F in self.coeffs.flatten())

    def is_imaginary(self) -> bool:
        """Whether the sequence has only imaginary elements."""
        if self.shape == ():
            return all(np.abs(np.real(F)) <= bd.machine_threshold() for F in self.coeffs)

        return all(np.abs(np.real(F)) <= bd.machine_threshold() for F in self.coeffs.flatten())

    def is_symmetric(self) -> bool:
        """Whether the sequence satisfies $F_k = F_{-k}$."""
        if self.shape == ():
            return all(np.abs(self[k] - self[-k]) <= bd.machine_threshold() for k in self.support())

        return all(np.linalg.norm(self[k] - self[-k]) <= bd.machine_threshold() for k in self.support())

    def __add__(self, other):
        if isinstance(other, Number):
            q = self.duplicate()
            q[0] += other

            return q
        elif isinstance(other, np.ndarray):
            if other.shape != self.shape:
                raise ValueError(f"The shape of the given matrix {other.shape} does not match the expected shape {self.shape}.")

            q = self.duplicate()
            q[0] += other

            return q
        elif not isinstance(other, ComplexL0Sequence):
            raise TypeError("Sequence addition admits only other sequences or scalars.")

        self_end = self.support_start + self.coeffs.shape[0]
        other_end = other.support_start + other.coeffs.shape[0]

        sum_start = min(self.support_start, other.support_start)
        sum_end = max(self_end, other_end)

        sum_coeffs = []
        for k in range(sum_start, sum_end):
            res = 0

            if self.support_start <= k and k < self_end:
                res += self.coeffs[k - self.support_start]

            if other.support_start <= k and k < other_end:
                res += other.coeffs[k - other.support_start]

            sum_coeffs.append(res)

        return Polynomial(np.array(sum_coeffs, dtype=complex_type), sum_start) # TODO This should return the lowest class between self and other.

    def __radd__(self, other):
        return self + other

    def __neg__(self):
        return type(self)(-self.coeffs, self.support_start)

    def __sub__(self, other):
        return self + (-other)

    def __rsub__(self, other):
        return self + (-other)

__getitem__(k: int | tuple) -> complex_type | np.ndarray

Returns the \(k\)-th element of the sequence, i.e., \(F_k\).

Parameters:

Name Type Description Default
k int

The index of the sequence.

 required

Returns:

Type Description
complex_type | ndarray

complex_type | np.ndarray: The coefficient \(F_k\), or \(0\) if \(k\) is out of the support.

Note: in case of matrix sequence, numpy slices can be used on the coefficients.

Source code in nlft_qsp/poly.py 
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
def __getitem__(self, k: int | tuple) -> complex_type | np.ndarray:
    r"""Returns the $k$-th element of the sequence, i.e., $F_k$.

    Args:
        k (int): The index of the sequence.

    Returns:
        complex_type | np.ndarray: The coefficient $F_k$, or $0$ if $k$ is out of the support.

    Note: in case of matrix sequence, numpy slices can be used on the coefficients.
    """
    if not isinstance(k, tuple):
        k = (k,)

    supp = self.support()
    if isinstance(k[0], slice):
        sl1 = max(k[0].start, supp.start) - self.support_start if k[0].start else 0
        sl2 = min(k[0].stop, supp.stop) - self.support_start if k[0].stop else self.coeffs.shape[0]

        return Polynomial(self.coeffs[(slice(sl1, sl2, k[0].step), *k[1:])], support_start=supp.start + sl1)

    if k[0] in supp:
        return self.coeffs[(k[0] - self.support_start, *k[1:])]

    if self.shape != ():
        return np.zeros(self.shape, dtype=complex_type)[k[1:]]

    return 0

__init__(coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None)

Initializes a complex sequence.

Parameters:

Name Type Description Default
coeffs list[complex_type] | ndarray

List of complex numbers as coefficients. coeffs[k, :] will be the coefficient of \(z^(s + k)\) where \(s\) is support_start.

 None
shape tuple[int]

The shape of the sequence. If it is none, then the shape of coeffs is taken. This is only used for matrix sequences, and is ignored otherwise.

 None
support_start optional

Index of the first element of the sequence. Defaults to \(0\).

 0
Note

Please provide exactly one of coeffs or shape.

Source code in nlft_qsp/poly.py 
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
def __init__(self, coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None):
    """Initializes a complex sequence.

    Args:
        coeffs: List of complex numbers as coefficients. `coeffs[k, :]` will be the coefficient of $z^(s + k)$ where $s$ is `support_start`.
        shape: The shape of the sequence. If it is none, then the shape of coeffs is taken. This is only used for matrix sequences, and is ignored otherwise.
        support_start (optional): Index of the first element of the sequence. Defaults to $0$.

    Note:
        Please provide exactly one of coeffs or shape.
    """
    if coeffs is not None and shape is not None or coeffs is None and shape is None:
        raise ValueError("Please provide exactly one between coeffs or shape.")

    self.support_start = support_start

    if shape is not None:
        self.coeffs = np.zeros((1, *shape), dtype=complex_type)
        self.shape = shape
        return

    if isinstance(coeffs, list):
        if any(not isinstance(ck, Number) for ck in coeffs):
            raise ValueError(f"The list of coefficients must be a single list or np.ndarray of complex values. Got {coeffs}. If you want to create a matrix sequence, please provide a np.ndarray of shape (num_coeffs, *shape) or specify the shape argument.")

        self.coeffs = np.array(coeffs, dtype=complex_type)
        self.shape = ()
        return

    if isinstance(coeffs, np.ndarray):
        if coeffs.ndim == 0:
            raise ValueError(f"There must be at least one coefficient. Got {coeffs}.")

        self.coeffs = np.array(coeffs, dtype=complex_type)
        self.shape = coeffs.shape[1:]
        return

    raise ValueError(f"Coefficients must be provided as a list or np.ndarray. Got {coeffs}.")

__setitem__(k: int, c: complex_type | np.ndarray)

Sets the coefficient of \(z^k\) to be \(c\), allocating space if needed.

Parameters:

Name Type Description Default
k int

The exponent of \(z\).

 required
c complex

The coefficient to set.

 required
Source code in nlft_qsp/poly.py 
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def __setitem__(self, k: int, c: complex_type | np.ndarray):
    r"""Sets the coefficient of $z^k$ to be $c$, allocating space if needed.

    Args:
        k (int): The exponent of $z$.
        c (complex): The coefficient to set.
    """
    if isinstance(k, int):
        k = (k,)

    if self.support_start + self.coeffs.shape[0] <= k[0]:
        pad_width = [(0, k[0] - self.support_start - self.coeffs.shape[0] + 1)] + [(0, 0)] * (self.coeffs.ndim - 1)
        self.coeffs = np.pad(self.coeffs, pad_width)
    elif self.support_start > k[0]:
        pad_width = [(self.support_start - k[0], 0)] + [(0, 0)] * (self.coeffs.ndim - 1)
        self.coeffs = np.pad(self.coeffs, pad_width)
        self.support_start = k[0]

    self.coeffs[(k[0] - self.support_start, *k[1:])] = c

is_imaginary() -> bool

Whether the sequence has only imaginary elements.

Source code in nlft_qsp/poly.py 
143
144
145
146
147
148
def is_imaginary(self) -> bool:
    """Whether the sequence has only imaginary elements."""
    if self.shape == ():
        return all(np.abs(np.real(F)) <= bd.machine_threshold() for F in self.coeffs)

    return all(np.abs(np.real(F)) <= bd.machine_threshold() for F in self.coeffs.flatten())

is_real() -> bool

Whether the sequence has only real elements.

Source code in nlft_qsp/poly.py 
136
137
138
139
140
141
def is_real(self) -> bool:
    """Whether the sequence has only real elements."""
    if self.shape == ():
        return all(np.abs(np.imag(F)) <= bd.machine_threshold() for F in self.coeffs)

    return all(np.abs(np.imag(F)) <= bd.machine_threshold() for F in self.coeffs.flatten())

is_symmetric() -> bool

Whether the sequence satisfies \(F_k = F_{-k}\).

Source code in nlft_qsp/poly.py 
150
151
152
153
154
155
def is_symmetric(self) -> bool:
    """Whether the sequence satisfies $F_k = F_{-k}$."""
    if self.shape == ():
        return all(np.abs(self[k] - self[-k]) <= bd.machine_threshold() for k in self.support())

    return all(np.linalg.norm(self[k] - self[-k]) <= bd.machine_threshold() for k in self.support())

l2_norm() -> float_type

Computes the \(\ell_2\) norm.

Source code in nlft_qsp/poly.py 
121
122
123
def l2_norm(self) -> float_type:
    r"""Computes the $\ell_2$ norm."""
    return np.sqrt(self.l2_squared_norm())

l2_squared_norm() -> float_type

Computes the squared \(\ell_2\) norm.

Note

For matrix sequences, this is the sum of the squared Frobenius norms of the coefficients.

Source code in nlft_qsp/poly.py 
125
126
127
128
129
130
131
132
133
134
def l2_squared_norm(self) -> float_type:
    r"""Computes the squared $\ell_2$ norm.

    Note:
        For matrix sequences, this is the sum of the squared Frobenius norms of the coefficients.
    """
    if self.shape == ():
        return sum(c * np.conj(c) for c in self.coeffs)

    return sum(np.linalg.norm(c)**2 for c in self.coeffs)

support() -> range

Returns the range in \(\Z\) where the sequence is non-zero.

Note

This simply checks the allocated array of coefficients but does not check leading or trailing zeros.

Source code in nlft_qsp/poly.py 
64
65
66
67
68
69
70
def support(self) -> range:
    r"""Returns the range in $\Z$ where the sequence is non-zero.

    Note:
        This simply checks the allocated array of coefficients but does not check leading or trailing zeros.
    """
    return range(self.support_start, self.support_start + self.coeffs.shape[0])

Polynomial

Bases: ComplexL0Sequence

Represents a general Laurent polynomial of one complex variable.

Attributes:

Name Type Description
coeffs list[complex_type]

List of complex coefficients.
shape tuple

The shape of the coefficients. This is always () for scalar polynomials, but is used for matrix polynomials.
support_start int

Minimum degree that appears in the polynomial.

Methods:

Name Description
__call__

Evaluates the polynomial using Horner's method.
__init__

Initializes a Polynomial instance.
__str__

Converts the polynomial to a human-readable string representation.
analytic_part

Discards all the negative degrees, keeping only the non-negative ones.
anti_analytic_part

Discards all the positive degrees, keeping only the non-positive ones.
block_matrix

Returns a matrix polynomial containing the given polynomials/constants as blocks.
conjugate

Returns the conjugate polynomial on the unit circle. If \(p(z) = \sum_k p_k z^k\), then its conjugate is defined as \(p^*(z) = \sum_k p_k^* z^{-k}\). This is also known as the Schwarz reflection across \(\mathbb{T}\).
diagonal_block_matrix

Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.
duplicate

Creates a duplicate of the current polynomial.
effective_degree

Returns the size of the support of the polynomial minus 1 (max degree - min degree).
eval_at_roots_of_unity

Evaluates the polynomial at the \(N\)-th roots of unity using the inverse fast Fourier transform.
hilbert_transform

Returns the polynomial \(Q\) such that \(P + Q\) yields an analytic polynomial (\(P\) being self).
only_negative_degrees

DEPRECATED: use anti_analytic_part() instead.
only_positive_degrees

DEPRECATED: use analytic_part() instead.
schwarz_transform

Returns the anti-analytic polynomial whose real part gives the current polynomial.
sharp

Same as conjugate(), but support_start is left unchanged.
shift

Creates a new polynomial equal to the current one, multiplied by z^k.
sup_norm

Estimates the supremum norm of the polynomial over the unit circle
truncate

Keeps only the coefficients in \([m, n]\), discarding the others.
Source code in nlft_qsp/poly.py 
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
293
294
295
296
297
298
299
300
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
366
367
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
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
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
@serializable(
    type_tag="@qspx/polynomial",
    fields={
        "support_start": "support_start",
        "coeffs": "coeffs"
    }
)
class Polynomial(ComplexL0Sequence):
    """Represents a general Laurent polynomial of one complex variable.

    Attributes:
        coeffs (list[complex_type]): List of complex coefficients.
        shape (tuple): The shape of the coefficients. This is always `()` for scalar polynomials, but is used for matrix polynomials.
        support_start (int): Minimum degree that appears in the polynomial.
    """

    def __init__(self, coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None):
        """Initializes a Polynomial instance.

        Args:
            coeffs: List of complex numbers as coefficients.
            shape (optional): Shape of the coefficient array.
            support_start (optional): Minimum degree in the polynomial. Defaults to 0.

        Note:
            Please provide exactly one of coeffs or shape.
        """
        super().__init__(coeffs, support_start, shape)

    def duplicate(self) -> "Polynomial":
        """Creates a duplicate of the current polynomial.

        Returns:
            Polynomial: A new Polynomial instance with the same coefficients and support.
        """
        return type(self)(self.coeffs, self.support_start)

    def shift(self, k: int):
        """Creates a new polynomial equal to the current one, multiplied by `z^k`."""
        return type(self)(self.coeffs, self.support_start + k)

    def effective_degree(self) -> int:
        """Returns the size of the support of the polynomial minus 1 (max degree - min degree).

        Note:
            This does not check for leading or trailing zeros in the coefficient array.

        Returns:
            int: The effective degree of the polynomial.
        """
        return self.coeffs.shape[0] - 1

    def conjugate(self) -> "Polynomial":
        r"""Returns the conjugate polynomial on the unit circle. If $p(z) = \sum_k p_k z^k$, then its conjugate is defined as $p^*(z) = \sum_k p_k^* z^{-k}$. This is also known as the Schwarz reflection across $\mathbb{T}$.

        Note:
            For a matrix polynomial, each coefficient is conjugate-transposed.
        """
        if self.shape == ():
            conj_coeffs = np.flip(np.conj(self.coeffs), axis=0)
        else:
            axes = [0] + list(range(1, self.coeffs.ndim))[::-1] # transpose = reverse the order of the axes except the first one, which is the index of the coefficient
            conj_coeffs = np.flip(np.transpose(np.conj(self.coeffs), axes=axes), axis=0)

        return Polynomial(conj_coeffs, -(self.support_start + self.coeffs.shape[0] - 1))

    def sharp(self) -> "Polynomial":
        r"""Same as `conjugate()`, but `support_start` is left unchanged.

        Returns:
            Polynomial: The sharp-conjugate polynomial.
        """
        p = self.conjugate()
        p.support_start += self.effective_degree() + 1
        return p

    def schwarz_transform(self) -> "Polynomial":
        r"""Returns the anti-analytic polynomial whose real part gives the current polynomial.

        Note:
            This is equivalent to adding $i \mathcal{H}[p]`, where $\mathcal{H}[p]$ is the Hilbert transform of $p$.
        """
        schwarz_coeffs = []
        for k in self.support():
            if k < 0:
                schwarz_coeffs.append(2*self[k])
            elif k == 0:
                schwarz_coeffs.append(self[k])

        return Polynomial(np.array(schwarz_coeffs, dtype=complex_type), self.support_start)

    def hilbert_transform(self) -> "Polynomial":
        r"""Returns the polynomial $Q$ such that $P + Q$ yields an analytic polynomial ($P$ being `self`).

        Note:
            This is actually $i \mathcal{H}[P]$, i.e., the Hilbert transform as returned is already multiplied by $i$.
        """
        hilbert_coeffs = []
        for k in self.support():
            if k < 0:
                hilbert_coeffs.append(-self[k])
            elif k > 0:
                hilbert_coeffs.append(self[k])
            else:
                hilbert_coeffs.append(0)

        return Polynomial(np.array(hilbert_coeffs, dtype=complex_type), self.support_start)

    def __mul__(self, other):
        if isinstance(other, Number):
            return Polynomial(self.coeffs * other, self.support_start)
        elif isinstance(other, np.ndarray):
            if other.shape != self.shape:
                raise ValueError(f"The shape of the given matrix {other.shape} does not match the expected shape {self.shape}.")

            return Polynomial(self.coeffs * other, self.support_start)
        elif not isinstance(other, Polynomial):
            raise TypeError("Polynomial multiplication admits only other polynomials or constants with compatible shape.")

        # Pad so they end up with the same length
        target_len = self.coeffs.shape[0] + other.coeffs.shape[0] - 1
        pad_a = [(0, target_len - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
        pad_b = [(0, target_len - other.coeffs.shape[0])] + [(0, 0)] * (other.coeffs.ndim - 1)

        coeffs_a = np.fft.fft(np.pad(self.coeffs, pad_width=pad_a), axis=0)
        coeffs_b = np.fft.fft(np.pad(other.coeffs, pad_width=pad_b), axis=0)

        # Multiply in the Fourier domain
        coeffs_c = [a * b for a, b in zip(coeffs_a, coeffs_b)]

        # Inverse FFT to get the result
        new_coeffs = np.fft.ifft(coeffs_c, axis=0)
        support_start = self.support_start + other.support_start  # Lowest degree of the new poly

        return Polynomial(new_coeffs, support_start)

    def __rmul__(self, other):
        return self * other

    def __truediv__(self, other):
        if isinstance(other, Number):
            return Polynomial(self.coeffs / other, self.support_start)

        raise TypeError("Polynomial division is only possible with scalars.")

    def __matmul__(self, other):
        if isinstance(other, np.ndarray):
            if self.coeffs.shape[-1] != other.shape[-2]:
                raise ValueError(f"Incompatible shapes: {self.shape} vs {other.shape}.")

            return Polynomial(self.coeffs @ other, self.support_start)
        elif not isinstance(other, Polynomial):
            raise TypeError("Polynomial multiplication admits only other polynomials or constants with compatible shape.")

        if self.coeffs.shape[-1] != other.shape[-2]:
            raise ValueError(f"Incompatible shapes: {self.shape} vs {other.shape}.")

        # Pad so they end up with the same length
        target_len = self.coeffs.shape[0] + other.coeffs.shape[0] - 1
        pad_a = [(0, target_len - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
        pad_b = [(0, target_len - other.coeffs.shape[0])] + [(0, 0)] * (other.coeffs.ndim - 1)

        coeffs_a = np.fft.fft(np.pad(self.coeffs, pad_width=pad_a), axis=0)
        coeffs_b = np.fft.fft(np.pad(other.coeffs, pad_width=pad_b), axis=0)

        # Multiply in the Fourier domain
        coeffs_c = [a @ b for a, b in zip(coeffs_a, coeffs_b)]

        # Inverse FFT to get the result
        new_coeffs = np.fft.ifft(coeffs_c, axis=0)
        support_start = self.support_start + other.support_start  # Lowest degree of the new poly

        return Polynomial(new_coeffs, support_start)

    def __rmatmul__(self, other):
        if isinstance(other, np.ndarray):
            if other.shape[-1] != self.coeffs.shape[-2]:
                raise ValueError(f"The shape of the given matrix {other.shape} does not match the expected shape {self.shape}.")

            return Polynomial(other @ self.coeffs, self.support_start)

        raise TypeError("Polynomial multiplication admits only other polynomials or constants with compatible shape.")

    def __call__(self, z) -> complex_type | np.ndarray:
        """Evaluates the polynomial using Horner's method.

        Args:
            z (complex): The point at which to evaluate the polynomial.

        Returns:
            complex: The evaluated result.
        """
        res = self.coeffs[-1]
        for k in reversed(range(self.coeffs.shape[0] - 1)):
            res = res * z + self.coeffs[k]
        return res * (z ** self.support_start)

    def eval_at_roots_of_unity(self, N: int) -> list[complex_type]: # TODO remove power-of-two assumption
        r"""Evaluates the polynomial at the $N$-th roots of unity using the inverse fast Fourier transform.

        Args:
            N (int): A power of two specifying the number of roots. If $N$ is not a power of two, then the next power of two is implicitly taken.

        Returns:
            list[complex]: List of evaluations at the N-th roots of unity.
            The $k$-th element of the list will be $P(e^{2\pi i k/N})$, $P$ being `self`.
        """
        N = next_power_of_two(N)
        M = next_power_of_two(max(N, self.coeffs.shape[0]))

        pad_width = [(0, M - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
        coeffs = np.pad(self.coeffs, pad_width=pad_width)
        coeffs = np.roll(coeffs, self.support_start, axis=0)
        # This has the effect of having everything multiplied by z^s

        evals = np.fft.ifft(coeffs, norm='forward', axis=0) # M evaluations at the M-th roots of unity
        return evals[::M//N]

    def sup_norm(self, N=1024) -> float_type:
        """Estimates the supremum norm of the polynomial over the unit circle

        Args:
            N (int): the number of samples to compute the maximum from. If $N$ is not a power of two, then the next power of two is taken.

        Returns:
            float_type: An estimate for the supremum norm of the polynomial over the unit circle.
        """
        if self.shape == ():
            return np.max([np.abs(sample) for sample in self.eval_at_roots_of_unity(N)])

        return np.max([np.linalg.norm(sample, ord=2) for sample in self.eval_at_roots_of_unity(N)])

    def truncate(self, m: int, n: int) -> "Polynomial":
        """Keeps only the coefficients in $[m, n]$, discarding the others.

        Returns:
            Polynomial: A new, truncated polynomial.
        """
        if self.shape == ():
            return Polynomial([self[k] for k in range(m, n+1)], m)

        return Polynomial(np.array([self[k] for k in range(m, n+1)], dtype=complex_type), m)

    def only_positive_degrees(self) -> "Polynomial":
        """DEPRECATED: use `analytic_part()` instead."""
        return self.analytic_part()

    def analytic_part(self) -> "Polynomial":
        """Discards all the negative degrees, keeping only the non-negative ones.

        Returns:
            Polynomial: A new polynomial containing only the positive-degree coefficients."""
        return self.truncate(0, self.support_start + self.coeffs.shape[0] - 1)

    def only_negative_degrees(self) -> "Polynomial":
        """DEPRECATED: use `anti_analytic_part()` instead."""
        return self.anti_analytic_part()

    def anti_analytic_part(self) -> "Polynomial":
        """Discards all the positive degrees, keeping only the non-positive ones.

        Returns:
            Polynomial: A new polynomial containing only the non-positive-degree coefficients."""
        return self.truncate(self.support_start, 0)

    def __str__(self) -> str:
        """Converts the polynomial to a human-readable string representation.

        Returns:
            str: The string representation of the polynomial.
        """
        return ' + '.join(f"{c} z^{self.support_start + k}" for k, c in enumerate(self.coeffs))

    @classmethod
    def block_matrix(cls, blocks: list[list["Polynomial"]]):
        """Returns a matrix polynomial containing the given polynomials/constants as blocks.

        Raises: ValueError if the block shapes do not match."""
        joint_support = range(min(min(p.support().start for p in l) for l in blocks),
                              max(max(p.support().stop for p in l) for l in blocks))

        coeffs = []
        for k in joint_support:
            coeffs.append(np.block([[p[k] for p in l] for l in blocks]))

        return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

    @classmethod
    def diagonal_block_matrix(cls, blocks: list["Polynomial"]):
        """Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.

        Raises: ValueError if the block shapes do not match."""
        joint_support = range(min(p.support().start for p in blocks), max(p.support().stop for p in blocks))

        coeffs = []
        for k in joint_support:
            coeffs.append(sp.linalg.block_diag(*[p[k] for p in blocks]))

        return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

__call__(z) -> complex_type | np.ndarray

Evaluates the polynomial using Horner's method.

Parameters:

Name Type Description Default
z complex

The point at which to evaluate the polynomial.

 required

Returns:

Name Type Description
complex complex_type | ndarray

The evaluated result.
Source code in nlft_qsp/poly.py 
390
391
392
393
394
395
396
397
398
399
400
401
402
def __call__(self, z) -> complex_type | np.ndarray:
    """Evaluates the polynomial using Horner's method.

    Args:
        z (complex): The point at which to evaluate the polynomial.

    Returns:
        complex: The evaluated result.
    """
    res = self.coeffs[-1]
    for k in reversed(range(self.coeffs.shape[0] - 1)):
        res = res * z + self.coeffs[k]
    return res * (z ** self.support_start)

__init__(coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None)

Initializes a Polynomial instance.

Parameters:

Name Type Description Default
coeffs list[complex_type] | ndarray

List of complex numbers as coefficients.

 None
shape optional

Shape of the coefficient array.

 None
support_start optional

Minimum degree in the polynomial. Defaults to 0.

 0
Note

Please provide exactly one of coeffs or shape.

Source code in nlft_qsp/poly.py 
223
224
225
226
227
228
229
230
231
232
233
234
def __init__(self, coeffs: list[complex_type] | np.ndarray = None, support_start: int = 0, shape: tuple[int] = None):
    """Initializes a Polynomial instance.

    Args:
        coeffs: List of complex numbers as coefficients.
        shape (optional): Shape of the coefficient array.
        support_start (optional): Minimum degree in the polynomial. Defaults to 0.

    Note:
        Please provide exactly one of coeffs or shape.
    """
    super().__init__(coeffs, support_start, shape)

__str__() -> str

Converts the polynomial to a human-readable string representation.

Returns:

Name Type Description
str str

The string representation of the polynomial.
Source code in nlft_qsp/poly.py 
472
473
474
475
476
477
478
def __str__(self) -> str:
    """Converts the polynomial to a human-readable string representation.

    Returns:
        str: The string representation of the polynomial.
    """
    return ' + '.join(f"{c} z^{self.support_start + k}" for k, c in enumerate(self.coeffs))

analytic_part() -> Polynomial

Discards all the negative degrees, keeping only the non-negative ones.

Returns:

Name Type Description
Polynomial Polynomial

A new polynomial containing only the positive-degree coefficients.
Source code in nlft_qsp/poly.py 
454
455
456
457
458
459
def analytic_part(self) -> "Polynomial":
    """Discards all the negative degrees, keeping only the non-negative ones.

    Returns:
        Polynomial: A new polynomial containing only the positive-degree coefficients."""
    return self.truncate(0, self.support_start + self.coeffs.shape[0] - 1)

anti_analytic_part() -> Polynomial

Discards all the positive degrees, keeping only the non-positive ones.

Returns:

Name Type Description
Polynomial Polynomial

A new polynomial containing only the non-positive-degree coefficients.
Source code in nlft_qsp/poly.py 
465
466
467
468
469
470
def anti_analytic_part(self) -> "Polynomial":
    """Discards all the positive degrees, keeping only the non-positive ones.

    Returns:
        Polynomial: A new polynomial containing only the non-positive-degree coefficients."""
    return self.truncate(self.support_start, 0)

block_matrix(blocks: list[list[Polynomial]]) classmethod

Returns a matrix polynomial containing the given polynomials/constants as blocks.

Raises: ValueError if the block shapes do not match.

Source code in nlft_qsp/poly.py 
480
481
482
483
484
485
486
487
488
489
490
491
492
@classmethod
def block_matrix(cls, blocks: list[list["Polynomial"]]):
    """Returns a matrix polynomial containing the given polynomials/constants as blocks.

    Raises: ValueError if the block shapes do not match."""
    joint_support = range(min(min(p.support().start for p in l) for l in blocks),
                          max(max(p.support().stop for p in l) for l in blocks))

    coeffs = []
    for k in joint_support:
        coeffs.append(np.block([[p[k] for p in l] for l in blocks]))

    return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

conjugate() -> Polynomial

Returns the conjugate polynomial on the unit circle. If \(p(z) = \sum_k p_k z^k\), then its conjugate is defined as \(p^*(z) = \sum_k p_k^* z^{-k}\). This is also known as the Schwarz reflection across \(\mathbb{T}\).

Note

For a matrix polynomial, each coefficient is conjugate-transposed.

Source code in nlft_qsp/poly.py 
259
260
261
262
263
264
265
266
267
268
269
270
271
def conjugate(self) -> "Polynomial":
    r"""Returns the conjugate polynomial on the unit circle. If $p(z) = \sum_k p_k z^k$, then its conjugate is defined as $p^*(z) = \sum_k p_k^* z^{-k}$. This is also known as the Schwarz reflection across $\mathbb{T}$.

    Note:
        For a matrix polynomial, each coefficient is conjugate-transposed.
    """
    if self.shape == ():
        conj_coeffs = np.flip(np.conj(self.coeffs), axis=0)
    else:
        axes = [0] + list(range(1, self.coeffs.ndim))[::-1] # transpose = reverse the order of the axes except the first one, which is the index of the coefficient
        conj_coeffs = np.flip(np.transpose(np.conj(self.coeffs), axes=axes), axis=0)

    return Polynomial(conj_coeffs, -(self.support_start + self.coeffs.shape[0] - 1))

diagonal_block_matrix(blocks: list[Polynomial]) classmethod

Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.

Raises: ValueError if the block shapes do not match.

Source code in nlft_qsp/poly.py 
494
495
496
497
498
499
500
501
502
503
504
505
@classmethod
def diagonal_block_matrix(cls, blocks: list["Polynomial"]):
    """Returns a matrix polynomial containing the given polynomials/constants as blocks along the diagonal.

    Raises: ValueError if the block shapes do not match."""
    joint_support = range(min(p.support().start for p in blocks), max(p.support().stop for p in blocks))

    coeffs = []
    for k in joint_support:
        coeffs.append(sp.linalg.block_diag(*[p[k] for p in blocks]))

    return Polynomial(np.array(coeffs, dtype=complex_type), support_start=joint_support.start)

duplicate() -> Polynomial

Creates a duplicate of the current polynomial.

Returns:

Name Type Description
Polynomial Polynomial

A new Polynomial instance with the same coefficients and support.
Source code in nlft_qsp/poly.py 
236
237
238
239
240
241
242
def duplicate(self) -> "Polynomial":
    """Creates a duplicate of the current polynomial.

    Returns:
        Polynomial: A new Polynomial instance with the same coefficients and support.
    """
    return type(self)(self.coeffs, self.support_start)

effective_degree() -> int

Returns the size of the support of the polynomial minus 1 (max degree - min degree).

Note

This does not check for leading or trailing zeros in the coefficient array.

Returns:

Name Type Description
int int

The effective degree of the polynomial.
Source code in nlft_qsp/poly.py 
248
249
250
251
252
253
254
255
256
257
def effective_degree(self) -> int:
    """Returns the size of the support of the polynomial minus 1 (max degree - min degree).

    Note:
        This does not check for leading or trailing zeros in the coefficient array.

    Returns:
        int: The effective degree of the polynomial.
    """
    return self.coeffs.shape[0] - 1

eval_at_roots_of_unity(N: int) -> list[complex_type]

Evaluates the polynomial at the \(N\)-th roots of unity using the inverse fast Fourier transform.

Parameters:

Name Type Description Default
N int

A power of two specifying the number of roots. If \(N\) is not a power of two, then the next power of two is implicitly taken.

 required

Returns:

Type Description
list[complex_type]

list[complex]: List of evaluations at the N-th roots of unity.
list[complex_type]

The \(k\)-th element of the list will be \(P(e^{2\pi i k/N})\), \(P\) being self.
Source code in nlft_qsp/poly.py 
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
def eval_at_roots_of_unity(self, N: int) -> list[complex_type]: # TODO remove power-of-two assumption
    r"""Evaluates the polynomial at the $N$-th roots of unity using the inverse fast Fourier transform.

    Args:
        N (int): A power of two specifying the number of roots. If $N$ is not a power of two, then the next power of two is implicitly taken.

    Returns:
        list[complex]: List of evaluations at the N-th roots of unity.
        The $k$-th element of the list will be $P(e^{2\pi i k/N})$, $P$ being `self`.
    """
    N = next_power_of_two(N)
    M = next_power_of_two(max(N, self.coeffs.shape[0]))

    pad_width = [(0, M - self.coeffs.shape[0])] + [(0, 0)] * (self.coeffs.ndim - 1)
    coeffs = np.pad(self.coeffs, pad_width=pad_width)
    coeffs = np.roll(coeffs, self.support_start, axis=0)
    # This has the effect of having everything multiplied by z^s

    evals = np.fft.ifft(coeffs, norm='forward', axis=0) # M evaluations at the M-th roots of unity
    return evals[::M//N]

hilbert_transform() -> Polynomial

Returns the polynomial \(Q\) such that \(P + Q\) yields an analytic polynomial (\(P\) being self).

Note

This is actually \(i \mathcal{H}[P]\), i.e., the Hilbert transform as returned is already multiplied by \(i\).

Source code in nlft_qsp/poly.py 
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
def hilbert_transform(self) -> "Polynomial":
    r"""Returns the polynomial $Q$ such that $P + Q$ yields an analytic polynomial ($P$ being `self`).

    Note:
        This is actually $i \mathcal{H}[P]$, i.e., the Hilbert transform as returned is already multiplied by $i$.
    """
    hilbert_coeffs = []
    for k in self.support():
        if k < 0:
            hilbert_coeffs.append(-self[k])
        elif k > 0:
            hilbert_coeffs.append(self[k])
        else:
            hilbert_coeffs.append(0)

    return Polynomial(np.array(hilbert_coeffs, dtype=complex_type), self.support_start)

only_negative_degrees() -> Polynomial

DEPRECATED: use anti_analytic_part() instead.

Source code in nlft_qsp/poly.py 
461
462
463
def only_negative_degrees(self) -> "Polynomial":
    """DEPRECATED: use `anti_analytic_part()` instead."""
    return self.anti_analytic_part()

only_positive_degrees() -> Polynomial

DEPRECATED: use analytic_part() instead.

Source code in nlft_qsp/poly.py 
450
451
452
def only_positive_degrees(self) -> "Polynomial":
    """DEPRECATED: use `analytic_part()` instead."""
    return self.analytic_part()

schwarz_transform() -> Polynomial

Returns the anti-analytic polynomial whose real part gives the current polynomial.

Note

This is equivalent to adding $i \mathcal{H}[p]`, where \(\mathcal{H}[p]\) is the Hilbert transform of \(p\).

Source code in nlft_qsp/poly.py 
283
284
285
286
287
288
289
290
291
292
293
294
295
296
def schwarz_transform(self) -> "Polynomial":
    r"""Returns the anti-analytic polynomial whose real part gives the current polynomial.

    Note:
        This is equivalent to adding $i \mathcal{H}[p]`, where $\mathcal{H}[p]$ is the Hilbert transform of $p$.
    """
    schwarz_coeffs = []
    for k in self.support():
        if k < 0:
            schwarz_coeffs.append(2*self[k])
        elif k == 0:
            schwarz_coeffs.append(self[k])

    return Polynomial(np.array(schwarz_coeffs, dtype=complex_type), self.support_start)

sharp() -> Polynomial

Same as conjugate(), but support_start is left unchanged.

Returns:

Name Type Description
Polynomial Polynomial

The sharp-conjugate polynomial.
Source code in nlft_qsp/poly.py 
273
274
275
276
277
278
279
280
281
def sharp(self) -> "Polynomial":
    r"""Same as `conjugate()`, but `support_start` is left unchanged.

    Returns:
        Polynomial: The sharp-conjugate polynomial.
    """
    p = self.conjugate()
    p.support_start += self.effective_degree() + 1
    return p

shift(k: int)

Creates a new polynomial equal to the current one, multiplied by z^k.

Source code in nlft_qsp/poly.py 
244
245
246
def shift(self, k: int):
    """Creates a new polynomial equal to the current one, multiplied by `z^k`."""
    return type(self)(self.coeffs, self.support_start + k)

sup_norm(N=1024) -> float_type

Estimates the supremum norm of the polynomial over the unit circle

Parameters:

Name Type Description Default
N int

the number of samples to compute the maximum from. If \(N\) is not a power of two, then the next power of two is taken.

 1024

Returns:

Name Type Description
float_type float_type

An estimate for the supremum norm of the polynomial over the unit circle.
Source code in nlft_qsp/poly.py 
425
426
427
428
429
430
431
432
433
434
435
436
437
def sup_norm(self, N=1024) -> float_type:
    """Estimates the supremum norm of the polynomial over the unit circle

    Args:
        N (int): the number of samples to compute the maximum from. If $N$ is not a power of two, then the next power of two is taken.

    Returns:
        float_type: An estimate for the supremum norm of the polynomial over the unit circle.
    """
    if self.shape == ():
        return np.max([np.abs(sample) for sample in self.eval_at_roots_of_unity(N)])

    return np.max([np.linalg.norm(sample, ord=2) for sample in self.eval_at_roots_of_unity(N)])

truncate(m: int, n: int) -> Polynomial

Keeps only the coefficients in \([m, n]\), discarding the others.

Returns:

Name Type Description
Polynomial Polynomial

A new, truncated polynomial.
Source code in nlft_qsp/poly.py 
439
440
441
442
443
444
445
446
447
448
def truncate(self, m: int, n: int) -> "Polynomial":
    """Keeps only the coefficients in $[m, n]$, discarding the others.

    Returns:
        Polynomial: A new, truncated polynomial.
    """
    if self.shape == ():
        return Polynomial([self[k] for k in range(m, n+1)], m)

    return Polynomial(np.array([self[k] for k in range(m, n+1)], dtype=complex_type), m)