A Chebyshev Spectral Method for Nonlinear Fourier Transform: Norming Constants

In this paper, we present a Chebyshev based spectral method for the computation of the Jost solutions corresponding to complex values of the spectral parameter in the Zakharov–Shabat scattering problem. The discrete framework is then used to devise a new algorithm based on a minimum total variation (MTV) principle for the computation of the norming constants which comprise the discrete part of the nonlinear Fourier spectrum. The method relies on the MTV principle to find the points where the expressions for norming constants are numerically well-conditioned.

09/09/2019

A Fast Chebyshev Spectral Method for Nonlinear Fourier Transform

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

Fast Nonlinear Fourier Transform using Chebyshev Polynomials

We explore the class of exponential integrators known as exponential tim...
12/05/2018

Nonlinear Fourier Transform of Truncated Multi-Soliton Pulses

Multi-soliton pulses, as special solutions of the Nonlinear Schroedinger...
01/31/2019

Statistics of the Nonlinear Discrete Spectrum of a Noisy Pulse

In the presence of additive Gaussian noise, the statistics of the nonlin...
03/12/2018

Minimum bias multiple taper spectral estimation

Two families of orthonormal tapers are proposed for multi-taper spectral...
11/08/2020

Fourier-based and Rational Graph Filters for Spectral Processing

Data are represented as graphs in a wide range of applications, such as ...
06/02/2020

Uncertainty Principle based optimization; new metaheuristics framework

To more flexibly balance between exploration and exploitation, a new met...

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, 8] (b) exponential time differencing (ETD) method based exponential integrators [9]. 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 [10] which achieves spectral accuracy at quasilinear complexity per sample of the signal. In this paper, we extend the recently proposed spectral method [11] for the computation continuous spectrum to compute the norming constants. It is well-known that the determination of the point where the expression which defines the norming constant is numerically well-conditioned is non-trivial problem. In the previous works [5, 6, 7] this point was taken to be origin, however, this choice can be shown to fail for carefully constructed examples. In order to remedy this problem, we propose a minimum total variation (MTV) principle to determine a set of points where the expression for the norming constants are well-conditioned. Note that total variation of the quantities in question are identically zero at the continuous level; therefore, it makes sense to seek the minima of TV for a sliding window of fixed size which traverses the sampling grid. The size of the sliding window can be adaptively reduced which adds to the effectiveness of the algorithm.

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) [1] scattering problem which can be stated as follows: For and ,

 vt=−iζσ3v+U(t)v, (1)

where is one of the Pauli matrices defined in the beginning of this article. The potential is defined by and . 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) with a prescribed behavior at or . The Jost solutions of the first kind, denoted by and , are the linearly independent solutions of (1) which have the following asymptotic behavior as : and . The Jost solutions of the second kind, denoted by and , are the linearly independent solutions of (1) which have the following asymptotic behavior as : and . The scattering coefficients are defined by

 (2)

for . The analytic continuation of the Jost solution with respect to is possible provided the potential is decays sufficiently fast or has a compact support. If the potential has a compact support, the Jost solutions have analytic continuation in the entire complex plane. Consequently, the scattering coefficients , , , are analytic functions of  [12].

In general, the nonlinear Fourier spectrum for the signal comprises a discrete and a continuous spectrum. The discrete spectrum consists of the so called eigenvalues , such that , and, the norming constants such that . For compactly supported potentials, . The continuous spectrum, also referred to as the reflection coefficient, is defined by for .

Ii The numerical scheme

Introducing the “local” scattering coefficients and such that , the ZS scattering problem can be written as and . Let the scattering potential be supported in . Such signals have been studied in more detail in [5, 13, 14, 15] also in order to understand the consequence of domain truncation for general signals. For , the ‘initial’ conditions for the Jost solution are: and . The scattering coefficients and are given by and . In the following, we describe a numerical scheme based on the Chebyshev polynomials to solve the coupled Volterra integral equations

 a(t;ζ)=1+∫t−1q(s)˘b(s;ζ)ds,˘b(t;ζ)=2iζ∫t−1˘b(s;ζ)ds+∫t−1r(s)a(s;ζ)ds, (3)

which are equivalent to the ZS problem with the aforementioned initial conditions. The numerical scheme computes the approximations to the Jost solutions in terms of the Chebyshev polynomials which can then used to compute the scattering coefficients. In the following, for the sake of brevity of presentation, we assume that the eigenvalues are known. The method of computing the eigenvalues using a Chebyshev based spectral would be presented in a future publication. Here we focus entirely on the computation of the norming constants. In principle, the norming constants can be computed simply by evaluating the numerical approximation to at the eigenvalue to obtain , however, becomes negligibly small which must now be multiplied with which grows exponentially. These intermediate quantities can easily suffer from lack of precision leading to inaccurate determination of .

Let . Following [11], the integral operator defined by in the Chebyshev basis is given by

 d(t)=[C0−14C1−∞∑n=2(−1)nCnn2−1]T0(t)+[C0−12C2]T1(t)+∞∑n=212n[Cn−1−Cn+1]Tn(t). (4)

In the matrix form, has the representation

 K=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1−14−13+18−115…10−12140−14160−16⋱⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (5)

The next step in the discretization of (3) involves expanding the signal in the Chebyshev basis. Let and where . A truncated expansion upto terms can be accomplished by sampling the potentials at the Chebyshev–Gauss–Lobatto (CGL) nodes given by and carrying out discrete Chebyshev transform which can be implemented using an FFT of size  [16].

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 (3) are the products and

which must be represented as linear operations on the unknown coefficient vectors

and . Again following [11], let and ; then, , and

 2Gl=l−1∑k=0Rl−kAk+2R0Al+∞∑k=1RkAk+l+∞∑k=0Rk+lAk,2Hl=l−1∑k=0Ql−kBk+2Q0Bl+∞∑k=1QkBl+k+∞∑k=0Qk+lBk. (6)

for . These relations define the operator , which comprises a Töplitz and an almost Hankel matrix given by

 (7)

Similarly, the representation of the operator can be defined. Setting , the discrete version of (3) can be stated as

 (I−ΛΛ∗I−2iζK)(AB)=S(AB)=(E00), (8)

where

is the identity matrix,

and . Noting that and setting , we have where we have used the fact that . The numerical scheme can be obtained as follows: truncation of the Chebyshev expansion of to terms, truncation of , and to -dimensional vectors, and, truncation of and to matrix where . If a direct sparse solver is used, the complexity of the numerical scheme would be lower than .

Remark II.1.

Let us remark that within the iterative approach and using the structured nature of the matrices involved it is possible to lower the complexity of the linear solver as in [11]. Our preliminary investigation indicate that the formulation

 [I+(I−2iζK)−1Γ]B=(I−2iζK)−1KR, (9)

works better for iterative solvers. However, the number of iterations needed for the stabilized biconjugate gradient method [17] can be as large as . Let us also mention that fast versions of the direct solvers for banded systems with certain number of filled first rows have been proposed [18]. We defer these aspects to a future publication.

The recipe discussed above accomplishes the computation of the Jost function of second kind , let us now show that by solving the scattering problem in the setting described above for it is possible to compute the Jost solution of the first kind . If denotes the Jost solution of the second kind for , we have  [5]. For convenience, let so that, for , the ‘initial’ conditions are: and . The scattering coefficients and are given by and . Again, we do not attempt to evaluate at the eigenvalues . Finally, let us note that the Chebyshev coefficients of are given by which follows from the symmetry property of the Chebyshev polynomials: .

Having discussed the computation of the Jost solutions, we turn to the computation of the norming. To this end, let . For convenience, we set and introduce the functions

 fδ(τ;ζk)=ln∣∣∣a(τ;ζk)˘c(τ;ζk)∣∣∣+2Imζkτ,gδ(τ;ζk)=ln∣∣ ∣∣˘b(τ;ζk)d(τ;ζk)∣∣ ∣∣+2Imζkτ,fθ(τ;ζk)=arg[a(τ;ζk)˘c(τ;ζk)e−2iReζkτ]gθ(τ;ζk)=arg[˘b(τ;ζk)d(τ;ζk)e−2iReζkτ], (10)

so that

 δk=fδ(τ;ζk)=gδ(τ;ζk),θk=fθ(τ;ζk)=gθ(τ;ζk), (11)

for . The functions defined above are constant with respect to for any given eigenvalue , however, at the discrete level they may vary. The choice of where the expressions above must be evaluated depends on the numerical conditioning of the quantities involved. Let a sub-grid be said to be admissible with respect to if its total variation over , given by

 V[fδ;J]=∑n|fδ(τn+1;ζk)−fδ(τn;ζk)|, (12)

satisfies the condition . In the numerical implementation, we introduce a sliding window given by where is fixed. If an appropriate tolerance cannot be guessed a priori, we simply choose that corresponds to . The point can then be reported as the optimal point for the computation of and . We label this algorithm as the minimum total variation (MTV) algorithm.

Iii Numerical Tests

For numerical tests, we consider the chirped secant-hyperbolic potential [19] given by for where

 f(s)=A0exp[−2iμA0log(coshs)]/cosh(s),s∈R. (13)

Here, controls how well the signal is supported in . Let and and set . Further, set and let ; then, the eigenvalues are given by and the corresponding norming constants are given by

 bk=ωe−2iμA0(log2)+iπkk−1∏j=1(ωA0−jω∗A0−j)e2iζkt0, (14)

for . We set so that and

 ζk=i(K+1/2−k)W,k=1,2,…,K. (15)

The signal for the choice of parameters , , and is shown in Fig. 2. The two expression for and provided in (11) are plotted in Fig. 3 which correspond to the numerically computed Jost solution for shown in Fig. 4 with and . It is evident from Fig. 3 that the choice of is non-trivial and the choice of is certainly not admissible. The MTV algorithm applied to the functions and with sliding-window size finds the optimal points to be , respectively. The errors and for each of the choices of are of the order .

A second example to demonstrate the effectiveness of the MTV algorithm is that of an -soliton solution whose discrete spectrum in terms of the triplets is listed in Table I. The signal computed using the classical Darboux transformation [20, 5, 14, 15] for the choice is shown in Fig. 5 and numerically computed Jost solutions in Fig. 6. The two expression for and provided in (11) are plotted in Fig. 7 which correspond to the numerically computed Jost solution for with and . It is evident from Fig. 7 that the choice of is again non-trivial and the choice of is once again not admissible. The MTV algorithm applied to the functions and with sliding-window size finds the optimal points to be , respectively. The errors and for each of the choices of are of the order .

The third example to demonstrate the effectiveness of the MTV algorithm is derived from the last example by multiplying a linear phase factor to the signal so that the eigenvalues acquire a shift of which we choose to set (see Table I). We let the numerical parameters to be the same as in the last example. The signal shown in Fig. 8, the variation of are shown in Fig. 9 which correspond to the numerically computed Jost solution (see Fig. 10) for . The MTV algorithm applied to the functions and finds the optimal points to be , respectively. The errors and for each of the choices of are of the order and , respectively.

The final numerical tests are meant to verify the spectral convergence of the numerical scheme. To this end, we quantify the error by

 eθ=1KK∑k=1|θk−θnum.k|,eδ=1KK∑k=1|δk−δnum.k|. (16)

We set and consider the set of values .

The first profile we would like to use for the convergence analysis is chirped hyperbolic profile. Let us set the parameters as , and . The choice of these parameters makes optimal for the computation of the norming constants. The results of the convergence analysis is shown in Fig. 11 which confirms the spectral convergence of the numerical scheme. The plateauing of the error seen in these plots are on account of the lack of compact support of .

The second profile we would like to use for the convergence analysis are multisoliton solutions. Let for . Then the eigenvalues are defined as . The norming constants are chosen as . The discrete spectrum corresponding to is depicted in Fig. 12. The potential can be computed with machine precision using the classical Darboux transformation algorithm [5]. The scaling parameter is set to be . Here, is known to be optimal for the computation of the norming constants. The results of the convergence analysis is shown in Fig. 13 which confirms the spectral convergence of the numerical scheme. The plateauing of the error seen in these plots are again on account of the lack of compact support of .

Iv Conclusion

In this paper, we presented a Chebyshev spectral method for the solution of the Zakharov–Shabat scattering problem for complex values of the spectral parameter. Within this discrete framework, we also proposed a robust algorithm for computing the norming constants. This algorithm is based on a minimum total variation principle and therefore the algorithm is abbreviated as the MTV algorithm for norming constants. Future work in this direction will focus on developing fast solvers for the linear system within the direct (relying on the structured nature of the system matrix) as well as iterative (relying on the fast matrix–vector multiplication for structured matrices) methods with and without preconditioning.

References

• [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, J. E. Prilepsky, S. T. Le, S. Wahls, L. L. Frumin, M. Kamalian, and S. A. Derevyanko, “Nonlinear Fourier transform for optical data processing and transmission: advances and perspectives,” Optica, vol. 4, no. 3, pp. 307–322, 2017.
• [4] V. Vaibhav and S. Wahls, “Introducing the fast inverse NFT,” in Optical Fiber Communication Conference.   Los Angeles, CA, USA: Optical Society of America, 2017, p. Tu3D.2.
• [5] V. Vaibhav, “Fast inverse nonlinear Fourier transformation using exponential one-step methods: Darboux transformation,” Phys. Rev. E, vol. 96, p. 063302, 2017.
• [6] V. Vaibhav, “Fast inverse nonlinear Fourier transform,” Phys. Rev. E, vol. 98, p. 013304, 2018.
• [7] V. Vaibhav, “Higher order convergent fast nonlinear Fourier transform,” IEEE Photonics Technol. Lett., vol. 30, no. 8, pp. 700–703, 2018.
• [8] V.Vaibhav, “Efficient nonlinear Fourier transform algorithms of order four on equispaced grid,” IEEE Photonics Technol. Lett., vol. 31, no. 15, pp. 1269–1272, 2019.
• [9] V. Vaibhav, “Fast nonlinear Fourier transform using Chebyshev polynomials,” arXiv, 2019, arXiv:1908.09811[physics.comp-ph]. [Online]. Available: https://arxiv.org/abs/1908.09811
• [10] V. Vaibhav, “Nonlinearly bandlimited signals,” J. Phys. A: Math. Theor., vol. 52, no. 10, p. 105202, 2019.
• [11] V. Vaibhav, “A fast Chebyshev spectral method for nonlinear Fourier transform,” arXiv, 2019, arXiv:1909.03710[physics.comp-ph]. [Online]. Available: https://arxiv.org/abs/1909.03710
• [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] V. Vaibhav, “Nonlinear Fourier transform of time-limited and one-sided signals,” J. Phys. A: Math. Theor., vol. 51, no. 42, p. 425201, 2018.
• [14] V. Vaibhav, “Exact solution of the Zakharov–Shabat scattering problem for doubly-truncated multisoliton potentials,” Commun. Nonlinear Sci. Numer. Simul., vol. 61, pp. 22–36, 2018.
• [15] V. Vaibhav, “Darboux transformation: new identities,” Physica Scripta, vol. 94, no. 11, p. 115504, 2019.
• [16] P. Giorgi, “On polynomial multiplication in Chebyshev basis,” IEEE Trans. Comput., vol. 61, no. 6, pp. 780–789, 2011.
• [17] 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.
• [18] S. Olver and A. Townsend, “A fast and well-conditioned spectral method,” SIAM Review, vol. 55, no. 3, pp. 462–489, 2013.
• [19] 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.
• [20] V. Vaibhav and W. Wahls, “Multipoint newton-type nonlinear Fourier transform for detecting multi-solitons,” in Optical Fiber Communication Conference.   Anaheim, CA, USA: Optical Society of America, 2016, p. W2A.34.