A Fast Chebyshev Spectral Method for Nonlinear Fourier Transform

by   Vishal Vaibhav, et al.

In this letter, we present a fast and well-conditioned spectral method based on the Chebyshev polynomials for computing the continuous part of the nonlinear Fourier spectrum. The algorithm achieves a complexity of O(N_iter.N N) per spectral node for N samples of the signal at the Chebyshev nodes where N_iter. is the number of iterations of the biconjugate gradient stabilized method.


page 1

page 2

page 3

page 4


Fast Nonlinear Fourier Transform using Chebyshev Polynomials

We explore the class of exponential integrators known as exponential tim...

A Chebyshev Spectral Method for Nonlinear Fourier Transform: Norming Constants

In this paper, we present a Chebyshev based spectral method for the comp...

A Fast and Accurate Nonlinear Spectral Method for Image Recognition and Registration

This article addresses the problem of two- and higher dimensional patter...

Robust, Nonparametric, Efficient Decomposition of Spectral Peaks under Distortion and Interference

We propose a decomposition method for the spectral peaks in an observed ...

Efficient Nonlinear Fourier Transform Algorithms of Order Four on Equispaced Grid

We explore two classes of exponential integrators in this letter to desi...

Successive Eigenvalue Removal for Multi-Soliton Spectral Amplitude Estimation

Optical nonlinear Fourier transform-based communication systems require ...

3rd-order Spectral Representation Method: Part I – Multi-dimensional random fields with fast Fourier transform implementation

This paper introduces a generalised 3rd-order Spectral Representation Me...

I Introduction

This paper considers the Zakharov and Shabat (ZS) [1] scattering problem which forms the basis for defining a nonlinear generalization of the conventional Fourier transform dubbed as the nonlinear Fourier transform (NFT). In an optical fiber communication system the nonlinear Fourier (NF) spectrum offers a novel way of encoding information in optical pulses where the nonlinear effects are adequately taken into account as opposed to being treated as a source of signal distortion [2, 3]. One of the challenges that has emerged in realizing these ideas is the development of a fast and well-conditioned NFT algorithm that can offer spectral accuracy at low complexity. Such an algorithm would prove extremely useful for system design and benchmarking. Currently, there are primarily two successful approaches proposed in the literature for computing the continuous NF spectrum which are capable of achieving algebraic orders convergence at quasilinear complexity: (a) the integrating factor (IF) based exponential integrators [4, 5, 6, 7] (b) exponential time differencing (ETD) method based exponential integrators [8]. Note that while the IF schemes uses fast polynomial arithmetic in the monomial basis, the ETD schemes use fast polynomial arithmetic in the Chebyshev basis. For the inverse transform, a sampling series based approach for computing the “radiative” part has been proposed in [9] which achieves spectral accuracy at quasilinear complexity per sample of the signal. In this paper, we propose a spectral method for direct NFT that achieves quasilinear complexity per sample of the continuous spectrum. The signal in this case must be sampled at the so-called Chebyshev–Gauss–Lobatto nodes.

Introducing the “local” scattering coefficients and such that , the ZS scattering problem can be written as and . Let the scattering potential (with where ) be supported in . Specializing to the real line , the initial conditions for the Jost solution are: and . The scattering coefficients and are given by and so that the reflection coefficients is . Let and . In the following, we describe a numerical scheme based on the Chebyshev polynomials to solve the coupled Volterra integral equations given by


derived from the ZS problem using the aforementioned initial conditions where the dependence on is suppressed for the sake of brevity of presentation. We refer the reader to [10] for the nonlinear Fourier analysis of time-limited signals.

Ii The numerical scheme

In this section, we describe the numerical scheme and carry out the theoretical analysis of its convergence and stability. Let us set the following notations: and (the set of whole numbers), () denotes the discrete Lebesgue spaces with norm and .

The first step is to obtain a discrete representation of an integral operator in the Chebyshev basis using the relations

Let and where , we have


Now setting

, and, introducing the infinite dimensional vectors

and , we have


The next step in the discretization of (1) involves expanding the scattering potentials in the Chebyshev basis. Let and where . A truncated expansion upto terms can be accomplished by sampling the potentials at the Chebyshev–Gauss–Lobatto nodes given by and carrying out discrete Chebyshev transform which can be implemented using an FFT of size  [11]. Now, our final goal is to obtain an expansion of the local scattering coefficients in the Chebyshev basis: To this end, let and where and are to be determined (for fixed value of ). The last ingredient needed in the discretization of (1) are the products and which must be represented as a linear operation on the unknown coefficient vectors and . To this end, let and ; then, it follows that , and


for . Following Olver and Townsend [12], these relations define the operator , which comprises a Töplitz and an almost Hankel matrix given by


Similarly, the representation of the operator follows from the same convention where we recall . Setting , the discrete version of (1) can be stated as



is the identity matrix,

and . Noting that and setting , we have . Before we attempt to solve this equation, let us determine the conditions under which a solution exists. To this end, let us show that the infinite matrices and (or, equivalently, ) form bounded linear operators on : Let and consider , then


The dominant term, which is clearly

, satisfies the estimate

so that which yields . Next, assuming that (which corresponds to the fact that is absolutely continuous on  [13]), we have and for so that which yields .

Next, we would like to show that the inverse of is bounded on . To this end, let denote the -th column of ; then, for any arbitrary , we have . Therefore, it suffices to show that is bounded for all where the resolvent of defined by so that . Consider ; by definition


Let . Following [10] and taking into account the definition of the Chebyshev coefficients together with the property , we have

It is possible to conclude that for continuing as above yielding the estimate


From here, it also follows that .

Fig. 1: The figure depicts the sparsity structure of the truncated version of the block matrix defined in (6). Here where is the number of Chebyshev polynomials used for the potential function . The matrix is banded with exactly two filled rows each corresponding to the filled top row of .
Fig. 2: The figure shows the convergence analysis for the proposed algorithms using a chirped hyperbolic potential.

The numerical algorithm is now obtained by truncating the Chebyshev expansion for and, equivalently, to terms while truncating to an matrix where so that the discrete system in reads as where and .

The sparsity structure of the original linear system in (6) is depicted in Fig. 1. Let denote the resolvent of ; then, where ( being the –term Chebyshev approximation to ).

In order to study the convergence behavior of the numerical algorithm, we consider the error where, with slight abuse of notation, we have assumed that is infinite dimensional by setting the additional entries to . We also extend the matrix to an infinite dimensional matrix in the same fashion. The error can be estimated from the relation which yields


The convergence of the numerical scheme, therefore, depends on which can be estimated as follows:


If the first derivatives of are absolutely continuous, then  [13]. Using the estimates obtained so far, an estimate for the condition number of the matrix in the –norm works out to be which guarantees that the numerical scheme remains well-conditioned for all .

Turning to the solution of the linear system , there are two options, namely, the direct method (possibly a sparse solver which takes into account the sparsity of ) or an iterative solver such as the BiCGSTAB [14], a stable variant of the biconjugate gradient method. For each value of the spectral parameter, the complexity of computing the scattering coefficient using the direct sparse solver is less than while the same for the sparse iterative solver is less than . A fast variant of the iterative solver for can be easily developed by observing that the matrix–vector product can be computed with complexity using the fast polynomial multiplication in Chebyshev basis [11] (which in turn relies on the FFT algorithm) yielding a complexity of for each value of the spectral parameter. We label this fast iterative solver as FBiCGSTAB. Let us note that the FBiCGSTAB does not need to actually create the matrix which makes it efficient in terms of usage of memory. The tolerance for the iterative solvers is set to be and the maximum iteration threshold to .

Iii Numerical Tests

In this section, we conduct numerical tests in order to confirm the convergence behavior and the complexity of the algorithms proposed. For this purpose we consider the chirped secant-hyperbolic potential [15] given by for where for . We restrict to the case . The parameter controls how well the profile is supported in . Let so that the discrete spectrum is empty. It can be shown that for ; therefore, the interval can be considered as the effective width of the coefficient. We set in the following. The numerical error is quantified by where comprises samples over some finite sequence and denotes the numerically computed values.

For a first test of convergence, we restrict ourselves to moderate values: , and with and . The results of the convergence analysis for fixed are shown in Fig. 2 which confirms the spectral convergence of the numerical scheme. Note that the direct solver has been dropped for the case . The seed solution for the iterative algorithms were null vectors, however, the fast iterative algorithm uses one final iteration of BiCGSTAB with the output of the former as the seed solution in order to improve the accuracy. This step can be omitted for large when forming becomes costly. The result of the complexity analysis is plotted in Fig. 3 which shows that iterative solvers outperform the direct solver. Note that the iterative solvers become increasingly more attractive when the samples of the continuous spectrum are needed on a spectral grid so that the solution of the previous grid point can be used as the seed solution for the current.

Fig. 3: Complexity analysis of the proposed algorithms for , and, .

In the final test, we focus only on the fast method. Let and () while keeping . In order to avoid early plateauing of error, we set because the chirped secant-hyperbolic profile does not strictly have a compact support. For the spectral parameter , we choose a uniform grid over of points. The result of the convergence analysis is plotted in Fig. 4 (top) which once again confirms the spectral convergence of the numerical scheme. The complexity of the numerical scheme per spectral node can be estimated from Fig. 4 (bottom) which confirms the quasilinear complexity.

Let us conclude with the numerical evidence in Fig. 5 for the claim that the condition number of remains bounded.

Fig. 4: Convergence (top) and complexity analysis for the fast method.
Fig. 5: The figure shows the condition number for and .


  • [1] V. E. Zakharov and A. B. Shabat, “Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media,” Sov. Phys. JETP, vol. 34, pp. 62–69, 1972.
  • [2] M. I. Yousefi and F. R. Kschischang, “Information transmission using the nonlinear Fourier transform, Part I,” IEEE Trans. Inf. Theory, vol. 60, no. 7, pp. 4312–4369, 2014.
  • [3] S. K. Turitsyn et al., “Nonlinear Fourier transform for optical data processing and transmission: advances and perspectives,” Optica, vol. 4, no. 3, pp. 307–322, 2017.
  • [4] V. Vaibhav, “Fast inverse nonlinear Fourier transformation using exponential one-step methods: Darboux transformation,” Phys. Rev. E, vol. 96, p. 063302, 2017.
  • [5] ——, “Fast inverse nonlinear Fourier transform,” Phys. Rev. E, vol. 98, p. 013304, 2018.
  • [6] ——, “Higher order convergent fast nonlinear Fourier transform,” IEEE Photonics Technol. Lett., vol. 30, no. 8, pp. 700–703, 2018.
  • [7] ——, “Efficient nonlinear Fourier transform algorithms of order four on equispaced grid,” IEEE Photonics Technol. Lett., vol. 31, no. 15, pp. 1269–1272, 2019.
  • [8] ——, “Fast nonlinear Fourier transform using Chebyshev polynomials,” 2019, arXiv:1908.09811[physics.comp-ph]. [Online]. Available: https://arxiv.org/abs/1908.09811
  • [9] ——, “Nonlinearly bandlimited signals,” J. Phys. A: Math. Theor., vol. 52, no. 10, p. 105202, 2019.
  • [10] ——, “Nonlinear Fourier transform of time-limited and one-sided signals,” J. Phys. A: Math. Theor., vol. 51, no. 42, p. 425201, 2018.
  • [11] P. Giorgi, “On polynomial multiplication in Chebyshev basis,” IEEE Trans. Comput., vol. 61, no. 6, pp. 780–789, 2011.
  • [12] S. Olver and A. Townsend, “A fast and well-conditioned spectral method,” SIAM Review, vol. 55, no. 3, pp. 462–489, 2013.
  • [13] L. N. Trefethen, “Is Gauss quadrature better than Clenshaw–Curtis?” SIAM Review, vol. 50, no. 1, pp. 67–87, 2008.
  • [14] H. A. Van der Vorst, “Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems,” SIAM J. Sci. Stat. Comput., vol. 13, no. 2, pp. 631–644, 1992.
  • [15] A. Tovbis et al., “On semiclassical (zero dispersion limit) solutions of the focusing nonlinear Schrödinger equation,” Commun. Pure Appl. Math., vol. 57, no. 7, pp. 877–985, 2004.