Efficient Nonlinear Fourier Transform Algorithms of Order Four on Equispaced Grid

by   Vishal Vaibhav, et al.

We explore two classes of exponential integrators in this letter to design nonlinear Fourier transform (NFT) algorithms with a desired accuracy-complexity trade-off and a convergence order of 4 on an equispaced grid. The integrating factor based method in the class of Runge-Kutta methods yield algorithms with complexity O(N^2N) (where N is the number of samples of the signal) which have superior accuracy-complexity trade-off than any of the fast methods known currently. The integrators based on Magnus series expansion, namely, standard and commutator-free Magnus methods yield algorithms of complexity O(N^2) that have superior error behavior even for moderately small step-sizes and higher signal strengths.



page 2


Fast Nonlinear Fourier Transform using Chebyshev Polynomials

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

A Fast Chebyshev Spectral Method for Nonlinear Fourier Transform

In this letter, we present a fast and well-conditioned spectral method b...

Fast Nonlinear Fourier Transform Algorithms Using Higher Order Exponential Integrators

The nonlinear Fourier transform (NFT) has recently gained significant at...

Efficient computation of the zeros of the Bargmann transform under additive white noise

We study the computation of the zero set of the Bargmann transform of a ...

Aliasing error of the exp(β√(1-z^2)) kernel in the nonuniform fast Fourier transform

The most popular algorithm for the nonuniform fast Fourier transform (NU...

Multi-layered Cepstrum for Instantaneous Frequency Estimation

We propose the multi-layered cepstrum (MLC) method to estimate multiple ...

Characteristic Analysis of 1024-Point Quantized Radix-2 FFT/IFFT Processor

The precise analysis and accurate measurement of harmonic provides a rel...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

In a series of papers [1, 2, 3], it was shown recently that exponential linear multistep methods (LMMs) provide a natural setting for higher-order convergent fast nonlinear Fourier transform (NFT). This followed from a simple observation that the transfer matrices obtained are amenable to FFT-based fast polynomial arithmetic. Note that in the earliest works on fast NFTs [4], the Ablowitz-Ladik method can be interpreted as the exponential Euler method. In this paper, we use the exponential Runge-Kutta methods to obtain a family of fast NFTs (provided that the nodes are equispaced). In particular, we present two fast NFTs with fourth order accuracy based on fourth order Runge-Kutta methods. The structure of the transfer matrix reveals that such methods are superior to those based on LMMs in terms of complexity while the numerical tests reveal that they also have a superior accuracy-complexity trade-off. The first algorithm is based on the classical fourth order (explicit) Runge-Kutta method which has been studied by several authors in the context of NFTs [5, 6]. The second method uses the three-stage Lobatto IIIA (implicit) Runge-Kutta method.

Fig. 1: The figure shows the convergence analysis of various methods for the secant-hyperbolic profile. It is evident that the CF method with two exponentials (CF) outperforms fourth order Magnus method (M). Note that the complexity of methods M, M and CF is while the rest are of .
Fig. 2: The figure shows the accuracy-complexity trade-off for various methods for the secant-hyperbolic profile. The legends are same as that of Fig. 2. The plots demonstrate that the Runge-Kutta methods (ERK, IRK) turn out to be far superior than all the other ‘fast’ methods with regard to the accuracy-complexity trade-off. It is interesting to note that the ‘slow’ methods become quite competitive with increasing signal strength.
Fig. 1: The figure shows the convergence analysis of various methods for the secant-hyperbolic profile. It is evident that the CF method with two exponentials (CF) outperforms fourth order Magnus method (M). Note that the complexity of methods M, M and CF is while the rest are of .

For moderately small step-sizes, most fast methods yield poor accuracy specially corresponding to the large values of the spectral parameter. An error analysis of such integrating factor based methods [5] shows that the error terms contain positive powers of the spectral parameter which necessitates the use of smaller step-sizes in order to keep the error low. On the contrary, integrators based on Magnus series expansion are known to have error terms that contain negative powers of the spectral parameter. In this letter, we follow the recipe presented by Blanes et al. [7] to develop a fourth order Magnus method and a fourth order commutator-free (CF) Magnus method [8], both of which take samples of the potential on an equispaced grid in order to compute the NFT. Let us emphasize that our CF method is different from those considered in [9]

where the samples of the potential are needed on the Gauss-Legendre nodes (the authors generate the samples by interpolation on an equispaced grid, locally). The fourth-order (standard) Magnus method happens to be faster than the corresponding CF method on account of the fact that there is an additional matrix exponential introduced in the CF method in order to avoid the use of commutators. The accuracy-complexity trade-off, however, is similar for the two methods. Finally, we also present a fast variant of the CF method of order four (formally) by employing the fourth order splitting on the lines of 

[10, 6]. Despite the well-known limitation imposed on the order and stability of such techniques as demonstrated by Sheng [11], we do not find any reduction of order within the double precision arithmetic. It is noteworthy that the authors in [6] found the aforementioned splitting worsen in accuracy after a certain step-size. It is also not clear from their analysis if this scheme is convergent in their setting.

We begin our discussion with a brief review of the scattering theory closely following the formalism presented in [12]. The nonlinear Fourier transform of any signal is defined via the Zakharov-Shabat (ZS) scattering problem which can be stated as follows: For and ,


where . The potential is defined by and with (). Here, is known as the spectral parameter and is the complex-valued signal. The solution of the scattering problem (1), henceforth referred to as the ZS problem, consists in finding the so called scattering coefficients which are defined through special solutions of (1) known as the Jost solutions which are linearly independent solutions of (1) such that they have a plane-wave like behavior at or . The Jost solution of the second kind, denoted by , has the asymptotic behavior as . The asymptotic behavior as determines the scattering coefficients and for . In this letter, we primarily focus on the continuous spectrum, also referred to as the reflection coefficient, which is defined by for .

Ii The Numerical Scheme

Ii-a Runge-Kutta Method

In this section, we will develop the integrating factor based exponential Runge-Kutta (RK) method for the numerical solution of the ZS problem. Following [1, 3], we begin with the transformation so that (1) becomes with whose entries are and . Let the step size be and the quantities be ordered so that the nodes within the step can be stated as . For the potential sampled at these nodes, we use the convention , and . In order for the resulting discrete system to be amenable to FFT-based fast polynomial arithmetic, it is sufficient to have each of the

’s belong to the set of uniformly distributed nodes in

. A -stage RK method is characterized by the nodes , and, the weights and . Introducing the intermediate stage quantities for , we have


This system of equations can be solved in any computer algebra system to obtained the transfer matrix connecting the vectors

to .

Setting , the Lobatto IIIA method (labelled as IRK) of order  [13] simplifies to


The fourth order classical RK method (labelled as ERK) simplifies to the form (3) with and the entries of the transfer matrix are given by


where .

Fig. 3: The figure shows the error behavior of the NFT algorithms for as a function of for the secant-hyperbolic profile. The legends are same as that of Fig. 2. It is straightforward to conclude that for smaller number of samples, the ‘slower’ methods (M, M and CF) are far superior to the ‘fast’ methods for , the principal domain of a given method, with comparable run-times.

Ii-A1 Scattering coefficients

Let the computational domain be and set . Let the number of steps be and . Also, let and . The grid is defined by with where is the number of samples. Let the potential be supported in and we assume for convenience. In order to represent the Jost solutions, we introduce the polynomial vector


and the polynomial . Consider the Jost solution . For the Lobatto IIIA method, the Jost solution can be stated as , with


where is the constant part of and


From , it follows that the discrete scattering coefficients are given by


where . For the classical fourth order RK method, ; therefore, and . The principal branch for the discrete scattering coefficients here works out to be . This follows from the principle branch of the individual transfer matrices. The nodes lead to which is not in the standard form for FFT algorithms to be used. Therefore, we would like to work with nodes so that

and pad the input vector with zeros. Following as in 

[3], the complexity of computing the scattering coefficient works out to be .

Ii-B Standard and commutator-free Magnus methods

Ii-B1 Magnus method

Let us assume that the solution of the ODE (1) can be written as for , then has a series representation known as the Magnus series [7]. Truncating this series to achieve the desired order of accuracy yields a family of numerical schemes known as Magnus method. For the ZS problem, Magnus integrators preserve the Lie group structure of the Jost solution and its accuracy does not worsen with increasing . In designing a fourth order Magnus integrator (labelled as M), we follow the method due to Blanes et al. [7, 8]: Defining


the method proceeds by expanding in terms of the quantities using the Magnus series. For the fourth order method, setting , we have [7, 8] . Evaluating and upto fourth order accuracy using the three-point Gauss quadrature involving Legendre-Gauss-Lobatto (LGL) nodes, , the numerical scheme for the ZS problem can be stated as where,


where , and . For the purpose of comparison, we would also like to consider the Magnus method with one point Gauss quadrature (labelled as M) which is of order  [1].

Ii-B2 Commutator-free Magnus method

Blanes and Moan [8] have constructed fourth-order commutator-free (CF) methods that are based on Magnus method. Using the quantities defined above, the CF method (labelled as CF) can be stated as


Evaluating and upto fourth-order accuracy using the three-point Gauss quadrature involving LGL nodes, we have where ,


and . Given that there are two matrix exponentials involved, the CF Magnus method has higher complexity than that of the standard Magnus method.

Ii-B3 A fast variant

The CF Magnus method further allows us to obtain a fast NFT algorithm (labelled as SCF) via splitting of the matrix exponential. Consider a formally fourth order splitting [10]:


This splitting is stable and convergent [10]; however, its global order of convergence is  [11]. Introducing and , and, putting , the transfer matrix relation can be written as


where the entries of the matrix are


where . The discrete scattering coefficients can be written as and . The principal branch for the discrete scattering coefficients here works out to be . This again follows from the principle branch of the individual transfer matrices. As before, the nodes lead to which is not in the standard form for FFT algorithms to be used. Therefore, we would like to work with nodes so that and pad the input vector with zeros.

Iii Numerical Tests and Conclusion

For the numerical experiments, we employ the well-known secant-hyperbolic potential given by , () for which the scattering coefficients are given in [1]. We set the computational domain to be and let . Let be the principal branch; then, the error in computing is quantified by


where the integrals are computed using the trapezoidal rule. Similar consideration applies to . For the purpose of testing, we include the implicit Adams method presented in [3] which are labelled as IA with . The method IA is identical to the trapezoidal rule, therefore, we use the label TR. The convergence analysis is carried out in Fig. 2 and the trade-off between accuracy and complexity is presented in Fig. 2. In terms of accuracy, the CF outperforms every other method with M being a close second. However, the accuracy-complexity trade-off is similar for the two methods. The ‘fast’ methods evidently lower complexity at the cost of accuracy. The RK methods (ERK and IRK) outperform all the other ’fast’ methods in terms of accuracy-complexity trade-off (see Fig. 2); however, with increasing signal strength, the ‘slow’ methods becoming equally competitive. In fact, at moderately small step-sizes, the ‘slow’ methods far outperform the ‘fast’ methods (see Fig. 3) with increasing signal strength.


  • [1] V. Vaibhav, “Fast inverse nonlinear Fourier transformation using exponential one-step methods: Darboux transformation,” Phys. Rev. E, vol. 96, p. 063302, 2017.
  • [2] ——, “Fast inverse nonlinear Fourier transform,” Phys. Rev. E, vol. 98, p. 013304, 2018.
  • [3] ——, “Higher order convergent fast nonlinear Fourier transform,” IEEE Photonics Technol. Lett., vol. 30, no. 8, pp. 700–703, 2018.
  • [4] S. Wahls and H. V. Poor, “Fast numerical nonlinear Fourier transforms,” IEEE Trans. Inf. Theory, vol. 61, no. 12, pp. 6957–6974, 2015.
  • [5] S. Burtsev, R. Camassa, and I. Timofeyev, “Numerical algorithms for the direct spectral transform with applications to nonlinear Schrödinger type systems,” J. Comput. Phys., vol. 147, no. 1, pp. 166–186, 1998.
  • [6] S. Chimmalgi, P. J. Prins, and S. Wahls, “Fast nonlinear Fourier transform algorithms using higher order exponential integrators,” 2018, arXiv:1812.00703[eess.SP]. [Online]. Available: https://arxiv.org/abs/1812.00703
  • [7] S. Blanes, F. Casas, and J. Ros, “Improved high order integrators based on the magnus expansion,” BIT Numer. Math., vol. 40, no. 3, pp. 434–450, 2000.
  • [8] S. Blanes and P. C. Moan, “Fourth- and sixth-order commutator-free Magnus integrators for linear and non-linear dynamical systems,” Appl. Numer. Math., vol. 56, no. 12, pp. 1519–1537, 2006.
  • [9] S. Chimmalgi, P. J. Prins, and S. Wahls, “Nonlinear Fourier transform algorithm using a higher order exponential integrator,” in Signal Processing in Photonic Communications, 2018, p. SpM4G.5.
  • [10] S. Descombes, “Convergence of a splitting method of high order for reaction-diffusion systems,” Mathematics of Computation, vol. 70, no. 236, pp. 1481–1501, 2001. [Online]. Available: http://www.jstor.org/stable/2698737
  • [11]

    Q. Sheng, “Solving linear partial differential equations by exponential splitting,”

    IMA J. Numer. Anal., vol. 9, no. 2, pp. 199–212, 1989.
  • [12] M. J. Ablowitz, D. J. Kaup, A. C. Newell, and H. Segur, “The inverse scattering transform - Fourier analysis for nonlinear problems,” Stud. Appl. Math., vol. 53, no. 4, pp. 249–315, 1974.
  • [13] J. C. Butcher,

    Numerical Methods for Ordinary Differential Equations

    , 2nd ed.   San Francisco: John Wiley & Sons, Ltd, 2003.