nlft_qsp.solvers.riemann_hilbert

Module for inverse NLFT based on Riemann-Hilbert factorizations.

Functions:

Name	Description
factorize	Computes the (right) Riemann-Hilbert factorization of \((a, b)\), i.e., the pair of polynomials \((a_+, b_+)\) such that \((a, b) = (a_-, b_-) (a_+, b_+)\), where \(b_+\) has support in \([k, n]\).
inlft	Compute the inverse nonlinear Fourier transform using the Riemann-Hilbert algorithm. See arXiv:2407.05634 for an explanation of the method.
system_matrix	Returns a block matrix of the form
toeplitz	Returns the \((n-k+1) \times (n-k+1)\) Toeplitz matrix constructed with the coefficients of \(c\) (\(n\) is the index of the last coefficient of \(c\)). Specifically, the first column of the matrix will be \((c_n, c_{n-1}, ..., c_k)^T\), i.e., the reversed order of its coefficients in \([k, n]\).

factorize(c: Polynomial, k: int, normalize: bool = False) -> tuple[Polynomial, Polynomial]

Computes the (right) Riemann-Hilbert factorization of \((a, b)\), i.e., the pair of polynomials \((a_+, b_+)\) such that \((a, b) = (a_-, b_-) (a_+, b_+)\), where \(b_+\) has support in \([k, n]\).

Note

The two polynomials are given up to a common multiplicative constant.

Parameters:

Name	Type	Description	Default
c	Polynomial	The polynomial \(c\) that approximates the ratio \(b/a\).	required
k	int	The displacement, i.e., the point in the sequence where the cut should be made.	required
normalize	bool	whether to normalize the polynomials so that \(`|A_+|^2 + |B_+|^2 = 1\).	False

Returns:

Type	Description
tuple[Polynomial, Polynomial]	A pair of polynomials \((A_+, B_+)\) equal (up to a \(\sqrt{A_+(\infty)}\) factor, if normalize=False) to the right Riemann-Hilbert factorization of \((a, b)\), whose ratio is approximated by \(c\).

Source code in nlft_qsp/solvers/riemann_hilbert.py

def factorize(c: Polynomial, k: int, normalize: bool = False) -> tuple[Polynomial, Polynomial]:
    r"""Computes the (right) Riemann-Hilbert factorization of $(a, b)$, i.e., the pair of polynomials $(a_+, b_+)$ such that $(a, b) = (a_-, b_-) (a_+, b_+)$, where $b_+$ has support in $[k, n]$.

    Note:
        The two polynomials are given up to a common multiplicative constant.

    Args:
        c (Polynomial): The polynomial $c$ that approximates the ratio $b/a$.
        k (int): The displacement, i.e., the point in the sequence where the cut should be made.
        normalize (bool): whether to normalize the polynomials so that $`|A_+|^2 + |B_+|^2 = 1$.

    Returns:
        A pair of polynomials $(A_+, B_+)$ equal (up to a $\sqrt{A_+(\infty)}$ factor, if normalize=False) to the right Riemann-Hilbert factorization of $(a, b)$, whose ratio is approximated by $c$.
    """
    n = c.support().stop - 1
    d = n - k

    A = system_matrix(c, k)
    x = np.linalg.solve(A, [0] * (2*d+1) + [1])

    Ap = Polynomial(x[d+1:2*d+2], support_start=-d)
    Bp = Polynomial(x[0:d+1])

    if normalize:
        a_inf = np.sqrt(Ap[0])
        Bp /= a_inf
        Ap /= a_inf

    return Ap, Bp

inlft(b: Polynomial, c: Polynomial) -> NonLinearFourierSequence

Compute the inverse nonlinear Fourier transform using the Riemann-Hilbert algorithm. See arXiv:2407.05634 for an explanation of the method.

Parameters:

Name	Type	Description	Default
b	Polynomial	The starting polynomial, such that \((a, b)\) is the NLFT we want to compute the sequence for.	required
c	Polynomial	A polynomial approximating the ratio \(b/a\). The end of its support must coincide with the one of \(b\).	required

Returns:

Type	Description
NonLinearFourierSequence	A sequence whose NLFT is equal to \((a, b)\) (up to working precision).

Source code in nlft_qsp/solvers/riemann_hilbert.py

def inlft(b: Polynomial, c: Polynomial) -> NonLinearFourierSequence:
    """Compute the inverse nonlinear Fourier transform using the Riemann-Hilbert algorithm. See [arXiv:2407.05634](https://arxiv.org/abs/2407.05634) for an explanation of the method.

    Args:
        b (Polynomial): The starting polynomial, such that $(a, b)$ is the NLFT we want to compute the sequence for.
        c (Polynomial): A polynomial approximating the ratio $b/a$. The end of its support must coincide with the one of $b$.

    Returns:
        A sequence whose NLFT is equal to $(a, b)$ (up to working precision).
    """
    if c.support().stop != b.support().stop:
        raise ValueError("The supports of b and c must end at the same point.")

    F = []
    for k in b.support():
        Ap, Bp = factorize(c, k)

        F.append(Bp[0]/Ap[0])

    return NonLinearFourierSequence(F, b.support_start)

system_matrix(c: Polynomial, k: int)

Returns a block matrix of the form

\[\begin{pmatrix} I & -T^T \\ \overline{T} & I \end{pmatrix}\]

where \(T\) is the Toeplitz matrix of \(c\) and and \(I\) is the identity matrix. Solving the associated linear system yields the Riemann-Hilbert factorization of \((a, b)\), whose ratio is approximated by \(c\). The Toeplitz matrix will be computed using toeplitz(c, k), see the documentation for toeplitz for more information.

Parameters:

Name	Type	Description	Default
c	Polynomial	The Laurent polynomial \(c\) that approximates \(b/a\). It should have support in \((-\infty, d]\).	required
k	int	Parameter passed to the toeplitz function. The function will compute the Riemann-Hilbert factorization from the \(k\)-th element of the NLFT. This parameter should not go over the end of the support of \(c\).	required

Source code in nlft_qsp/solvers/riemann_hilbert.py

def system_matrix(c: Polynomial, k: int):
    r"""Returns a block matrix of the form 

    $$\begin{pmatrix} I & -T^T \\ \overline{T} & I \end{pmatrix}$$

    where $T$ is the Toeplitz matrix of $c$ and and $I$ is the identity matrix.
    Solving the associated linear system yields the Riemann-Hilbert factorization
    of $(a, b)$, whose ratio is approximated by $c$. The Toeplitz matrix will be
    computed using `toeplitz(c, k)`, see the documentation for `toeplitz` for more
    information.

    Args:
        c (Polynomial): The Laurent polynomial $c$ that approximates $b/a$. It should have support in $(-\infty, d]$.
        k (int): Parameter passed to the `toeplitz` function. The function will compute the Riemann-Hilbert factorization from the $k$-th element of the NLFT. This parameter should not go over the end of the support of $c$.
    """
    n = c.support().stop - 1
    d = n - k


    T = toeplitz(c, k)
    M = np.zeros(shape=(2*d+2, 2*d+2), dtype=complex_type)
    for i in range(d+1):
        for j in range(d+1):
            M[i, (d+1)+j] = -T[j, i]
            M[(d+1)+i, j] = np.conj(T[i, j])

    for i in range(2*d+2):
        M[i, i] = 1

    return M

toeplitz(c: Polynomial, k: int) -> np.ndarray

Returns the \((n-k+1) \times (n-k+1)\) Toeplitz matrix constructed with the coefficients of \(c\) (\(n\) is the index of the last coefficient of \(c\)). Specifically, the first column of the matrix will be \((c_n, c_{n-1}, ..., c_k)^T\), i.e., the reversed order of its coefficients in \([k, n]\).

Note

Assuming that \(c\) is supported on \((-\infty, n]\), this matrix will act as the linear operator \(f \mapsto z^k P (c z^{-k} f)\), where \(P\) is the Cauchy projection annihilating any negative-degree term.

Parameters:

Name	Type	Description	Default
c	Polynomial	The polynomial \(c\) whose coefficients will be arranged in the matrix.	required

Source code in nlft_qsp/solvers/riemann_hilbert.py

def toeplitz(c: Polynomial, k: int) -> np.ndarray:
    r"""Returns the $(n-k+1) \times (n-k+1)$ Toeplitz matrix constructed with the coefficients of $c$ ($n$ is the index of the last coefficient of $c$). Specifically, the first column of the matrix will be $(c_n, c_{n-1}, ..., c_k)^T$, i.e., the reversed order of its coefficients in $[k, n]$.

    Note:
        Assuming that $c$ is supported on $(-\infty, n]$, this matrix will act as the linear operator $f \mapsto z^k P (c z^{-k} f)$, where $P$ is the Cauchy projection annihilating any negative-degree term.

    Args:
        c (Polynomial): The polynomial $c$ whose coefficients will be arranged in the matrix.
    """
    n = c.support().stop - 1

    return np.array([[c[n + i - j] for i in range(0, n-k+1)] for j in range(0, n-k+1)], dtype=complex_type)