Discrete nonlinear Fourier transforms and their inverses

We study two discretisations of the nonlinear Fourier transform of AKNS-ZS type, ${\cal F}^E$ and ${\cal F}^D$. Transformation ${\cal F}^D$ is suitable for studying the distributions of the form $u = \sum_{n = 1}^N u_n \, \delta_{x_n}$, where $\delta _{x_n}$ are delta functions. The poles $x_n$ are not equidistant. The central result of the paper is the construction of recursive algorithms for inverses of these two transformations. The algorithm for $({\cal F}^D)^{- 1}$ is numerically more demanding than that for $({\cal F}^E)^{- 1}$. We describe an important symmetry property of ${\cal F}^D$. It enables the reduction of the nonlinear Fourier analysis of the constant mass distributions $u = \sum_{n = 1}^N u_c \, \delta _{x_n}$ for the numerically more efficient ${\cal F}^E$ and its inverse.


Introduction
The theory of nonlinear Fourier analysis stems from the inverse scattering method for solving nonlinear integrable partial differential equations. The inverse scattering method is a nonlinear version of the original Fourier's idea for solving the initial value problems for linear partial differential equations. The Fourier approach first uses the Fourier transform to recast the initial condition u(x, 0) into data whose time evolution is easily found. We find the transformed data for future time t and then apply the inverse Fourier transform to find the value of the solution u(x, t) of our problem at time t. Roughly speaking, the inverse scattering method has the same basic structure, except that the linear Fourier transform has to be replaced by a nonlinear one and the inverse linear transform by its nonlinear analogue. Therefore, finding the inverse nonlinear Fourier transform is an important problem. The inverse scattering method approach was discovered by Gardner, Greene, Kruskal and Miura in [6], [7], where they solved the Korteweg-deVries equation. There are many versions of the nonlinear Fourier transforms. One of the most popular is the transform, associated with the AKNS-ZS systems. This was introduced in the pioneering work of Ablowitz, Kaup, Newell, and Segur, and simultaneously by Zakharov and Shabat, see [1], [2], and [23].
In this study, we consider two discretisations of this transform and their inverses for functions, defined on a finite interval. The first discretisation F E is a simple Euler type construction. The second discretisation F D is more interesting. It can be used to transform the finite linear combinations of delta functions of the form u = N n=1 u n δ xn . Here, the points x n are not presumed to be equidistant. The arrays of values {u n } n=1,...,N and {x n } n=1,...,N are both the unknowns of the inverse of F D . Distributions of the form u n δ nn are of great interest in many fields of mathematics and their applications especially for signal processing. Crystalline measures and one-dimensional Fourier quasicrystals are special cases of infinite versions of such distributions. For an exciting work on these objects, see [8], [10], [9], [11], to mention but a few. The main tool of the experimental work in this field is finite sums of the form u = N n=1 u n δ xn . The transform F D can also be used as an approximation of the nonlinear transform of continuous functions when we apply it to the distributions of the form u = N n=1 u n ∆x n δ xn , where ∆x n = x n−1 − x n .
The main results of this study are Theorems 1 and 2. They provide recursive algorithms for the evaluation of the inverse transformations of F D and F E . Not surprisingly, the algorithm for inverse (F E ) −1 is numerically much cheaper than the one for (F D ) −1 . In [14] and [15] we provided perturbative constructions for the evaluations of two types of continuous AKNS-ZS nonlinear Fourier transforms. Our algorithms in this study are recursive, but they are not perturbative. In principle, they yield the exact solutions to the two inverse problems in finite time.
Transforms F E and F D send functions or distributions into functions z → A (z), where A (z) are matrices of the form All such functions are not elements of the images of F E or F D . As corollaries of theorems 1 and 2, we give the tests which decide when a matrix function , whose poles are equidistant. Then it is not difficult to see that one can essentially use the faster method for ( The question arises whether this can also be done in the case where the distributions in question are the constant mass distributions, that is, those of the form w = N n=1 u c δ xn . The answer to this question is negative. However, in F D , the roles of the variables u n and ∆x n = x n+1 − x n are almost symmetrical. Reversing the roles of u n and ∆x n yields the nonlinear Fourier transformation F D d which can be considered dual to F D , as explained in Section 6. Theorem 3 shows the following. Let B(z) be equal to F D d [w] for a constant mass distribution w with non-equidistant poles. Then can be computed using a faster algorithm for the evaluation of (F E ) −1 . We hope our results will contribute some new insight to the developing field of nonlinear Fourier analysis. In particular, the transformation F D could provide a helpful tool for the nonlinear Fourier analysis of distributions. Nonlinear Fourier analysis is becoming an essential part of the field of harmonic analysis and its applications in a broad sense; see the references [3], [12], [17], [18], [19], [20], to name but a few. The applicability of the linear Fourier transform is enormous, and the relation between the linear and nonlinear Fourier transforms is better understood in time, see the references [4], [13], [21], [22], and many other works. Undoubtedly, the role of nonlinear Fourier analysis will become ever more critical in studying various nonlinear problems.
The second section recalls the definition and basic properties of AKNS-ZS nonlinear Fourier transform for functions in finite intervals. In Section 3, we introduce discretisations F E and F D , respectively. The central part of the paper is sections 4 and 5, in which we construct the algorithms for evaluating (F D ) −1 and (F E ) −1 , respectively. In section 6, we describe the symmetric structure of F D and introduce the dual transformation F D d . We apply it to the study of the constant mass distributions of the form w = N n=1 u c δ xn . We conclude the paper with some brief remarks on the computational complexity of our algorithms. We also mention possible directions for further research.

Nonlinear Fourier transform
Let LSU(2) denote the space of functions of a real variable z, with values in the group of SU(2) matrices given by where the functions a(z) and b(z) are elements of the L 2 -space on a finite interval, say [0, 1]. Let x → u(x) be an L 2 [0, 1] complex-valued function. The nonlinear Fourier transform associated with the AKNS-ZS integrable systems on a finite interval is the map given by the following rule: Let Φ[u](x; z) be the solution of the family of linear initial value problems where the coeficient matrix is given by Then, The above initial value problem can be solved by a convergent function series, namely the Dyson series. We have Above ∆ d (x) denotes the ordered d dimensional simplex The linearisation of the nonlinear Fourier transform F [u] around u ≡ 0 is the usual linear Fourier transform. To see this, it is helpful to change the gauge using the gauge transformation matrix G(x, z) = diag(e −πixz , e πixz ). The transformed matrix L G of coefficients is given by The nonlinear Fourier transformation F G in the new gauge is given by

Discretizations of nonlinear Fourier transform
We shall describe two simple discretizations of the nonlinear Fourier transform.
Euler type of discretization: In the initial value problem (1) we can replace the differential equation using a difference equation. We replace the function u(x) with a function of the discrete variable n → u n , n = 0, . . . , N − 1, for some N ∈ Z, and obtain Φ(n + 1, z) − Φ(n, z) This gives and by recursion We restrict the values of z to {0, 1, . . . , N − 1}. Matrix multiplication yields the polynomial analogue of the Dyson series, namely, Let us denote We have Using (2) and (6), we collect all the exponential factors E N (−n, z) on the left. Then, (4) becomes Then the linear part of (7) is essentially the usual linear discrete Fourier transform. A similar approach was used in the continuous case in [5].
Transform of a sum of delta functions: Here, the spectral variable z takes value in the entire R. One of the most important ways to represent signals is to assign a suitable distribution of the form Let us introduce two discrete functions.
The values of u n are arbitrary complex numbers. The points x n are not assumed to be equidistant. We shall define {∆x n } n=0,...,N by Let now u ǫ (x) be the step function, given by Here, ǫ is a number that is smaller than the smallest ∆x n . So, the function u ǫ (x) is equal to zero at intervals [0, , it assumes large values of un ǫ . In terms of the theory of distributions, we have The nonlinear Fourier transform, defined by (1), of a step function u ǫ , is the product of the matrix exponentiations of constant pieces of u ǫ . We have where In the definition of E n , we have ǫ = ǫ for n = 1, . . . , N − 1 and ǫ = ǫ 2 for n = 0, N. A straightforward calculation gives R n (ǫ, z) = cos (A(u n , ǫz)) + iπǫz sinc(A(u n , ǫz)) u n sinc(A(u n , ǫz)) −u n sinc(A(u n , ǫz)) cos (A(u n , ǫz)) − iπǫz sinc(A(u n , ǫz)) , In the limit ǫ → 0, we have If we rewrite u n in polar coordinates, u n = e iφn r n , we get Definition 1 Let the distribution u be given by The above calculations provide proof of the following proposition.
. . , n n are arbitrary complex numbers. The nonlinear Fourier transform of u(x) is given by Proof: The proposition follows immediately from formulae (12) and (13), (14). ✷ Rewriting the formula for F D [u] in terms of the variables x n instead of ∆x n yields the following expression: The nonlinear Fourier transform of u = N n=1 u n δ xn can be expressed as the formula More explicitly, ✷ For easier writing, it will be handy to introduce a slightly modified, "reduced" transform At first glance, formula (15) appears to be different from the Euler-type discretisation F E of F . However, (17) enables us to find the analogue of Dyson's polynomial for F D [u], which makes the relation to the Euler discretization clearer. Let us introduce the matrix U n = 0 e iφn sin r n −e −iφn sin r n 0 .
We have cos r n e −2πixnz e iφn sin r n −e 2πixnz e −iφn sin r n cos r n = cos r n · I + E(−2x n , z) · U n . Denote cos r n .
Using the commutation rule the reduced (18) can be rewritten in the form The expressions (7) and (19) are indeed similar. In particular, we note that expanding V n (u) into the Taylor series around zero, with respect to u and keeping only the linear terms in the second summand = d k=1 tan (∆x n k u n k ) into Taylor series around zero with respect to the variables ∆x n u n and keep only the linear terms. Let us, further, set x n = n N for all n. Then, also ∆x n = 1 N , and the expression (19) becomes (7).

The inverse of F D
In this section, we describe an algorithm for evaluating the inverse of the nonlinear Fourier transform F D . The following theorem gives our result: Theorem 1 Suppose we are given an SU(2)-valued function can be computed from A (z) by the following recursion.
by the following construction. (iii) Set where matrices E and R are given by (13) and (14), respectively.
We repeat the sequence of operations (i) to (iv) until we reach A N −1 . Then the values, obtained in (iii), provide the wanted array In the proof of our theorem we will need the following lemma: Lemma 1 For every k = 0, . . . , ⌈ N 2 ⌉ and for every ordered array n 1 , n 2 , . . . n 2k−1 of odd number of elements from {n j } we have Proof: Consider the real numbers Because of the orderings of {n j } j=1,...,N and {x j } j=1,...,N , we have for every monotonically increasing array n 1 , n 2 , . . . n 2k−1 of integers, Indeed, on the right, the positive numbers x n j+1 − x n j are subtracted from x n 2k−1 . The number x N is the largest among all x n , so the lemma is proved. Let A (z) = F D r [u](z) for some u = N n=1 u n δ xn . We have n 2k−1 >n 2k >...>n 1 e −2πi(xn 2k−1 −xn 2k +...−xn 2 +xn 1 )·z S n 1 ,n 2 ,...,n 2k−1 (u). (21) Let w(z) = N j=1 w j e −2πix j z . Recall that the inverse linear Fourier transform of w is given by We can apply this formula to the function B[u](z), given by (21), and obtain where N (2k−1) is the set of all increasing integer arrays 0 < n 1 < n 2 < . . . < n 2k−1 ≤ N and x n 1 ,...,n 2k−1 = x n 2k−1 − x n 2k−2 + . . . − x n 2 + x n 1 . Therefore, for a generic choice of {x n } n=1,...,N , the linear inverse Fourier transform of B[u] is the sum of the 2 N −1 delta functions. Lemma 1 tells us that S N (u) δ x N is precisely the rightmost of these delta functions.

From the linear inverse transforms F −1 [B[u]] and F −1 [C (u) + A[u]]
, we can read the values x N , S N (u), and C (u). From the latter two we get We have obtained our first pair (x N , u N ) of the array U = {(x n , u n ); n = 1, . . . N}. To find another pair of U , we can proceed as follows: Recall the formula (18) Thus, the matrix function is the nonlinear Fourier transform F D r [u (1) ] of the distribution Now we can repeat the above procedure with A 1 = F D r [u (1) ] instead of A (z), and obtain the next pair (x N −1 , u N −1 ) of the array U . Repeating this procedure N times yields all the pairs in U . Hence our theorem is proved.
✷ An immediate consequence of the theorem is the following: We obtain the array U by dividing the second components of U by ∆x n which we obtain by subtracting the appropriate first components; ∆x n = x n+1 − x n .

The inverse of F E
In this section, z takes values in the set {0, 1, . . . , N − 1}. Let us denote and let us recall the definitions (5) of the matrices E N (n, z) and U n . We shall prove the following theorem: Theorem 2 Let A (z), z = 0, 1, . . . , N − 1 be an N-tuple of matrices of the form (22). Suppose we know that there exists a vector u = (u 0 , u 1 , . . . , u N −1 ) such that Then u can be computed by the following algorithm.
Then u = (u 0 , u 1 , . . . , u N −1 ), where un N are given by (ii), is the solution of the equation Proof: Recall that the discrete Fourier transform F E [u](z) of a discrete function u = (u 0 , u 1 , . . . , u N −1 ) can be given by its Dyson form where ∆ D d denotes the discrete ordered simplex ∆ D d = {(n 1 , n 2 , . . . , n d ) ∈ Z d ; 0 ≤ n 1 < n 2 < . . . < n d ≤ N − 1}. We have rewritten F E [u] more explicitly as It is easily seen from (3) that for every z, the transform F E [u](z) is a matrix of the form (22). (In general, the determinant of F E [u](z) is not equal to 1 and is therefore not an element of SU (2).) The entries α E (z) and β E (z) which determine the matrix are given by Note that in the above equations the expressions inside the large brackets, next to the factor e −2πi lz N are functions of variable l, because they are determined by the sets D n (l) Note also that the k = 1 term in (27) is the usual linear discrete Fourier transform F D [u](z) of u. Let us denote (n 1 ,...,n 2k−1 )∈D 2k−1 (l) u n 2k−1 u n 2k−2 · · · u n 2 u n 1 . Then Therefore, applying the linear inverse discrete transform to both sides of the above equation gives The right-hand sides of the above equations are known. Therefore, we have N equations for the N unknowns u 0 , u 1 , . . . , u N −1 . The above system is highly nonlinear. The terms S(l, k) are polynomials of degrees 2k−1 in the variables u 0 , u 1 , . . . , u N −1 . The number of terms of degree 2k −1 of such polynomial is equal to the number of points in the stratum D 2k−1 (l). Recall that the number of points in D 2k−1 (l) is the number of solutions of the equation where N − 1 ≥ n 2k−1 > n 2k−2 > . . . > n 1 ≥ 0 are integers. For l = N − 1, the above equation has only one solution: n 1 = N − 1. Therefore, the last equation of system (28) is: We have found the last component u N −1 of the unknown discrete function (u 0 , u 1 , . . . , u N −1 ).
Recall that F E [u] is given by (3); that is, Let for every z = 0, 1, . . . , N − 1, A rather straightforward calculation shows that is u with components cyclically permuted by one step. Now we repeat the procedure with A 1 in place of A : We calculate the linear discrete inverse Fourier transform ( and set Repeating this procedure N times yields all the wanted values (u 0 , u 1 , . . . , u N −1 ). ✷ for every z.

Remark 2
The starting point of the algorithm of Theorem 2 is the application of the inverse linear discrete Fourier transform to A (z). This yields the nonlinear system for l = 0, . . . , N − 1 and for the unknowns u 0 , u 1 , . . . u N −1 . The left-hand sides for l close to N 2 are huge. More precisely, the number of elements in each D 2k−1 (l) is equal to For N = 100, k = 10 and l = 25 we have ♯D 2k−1 (l) ∼ 2, 3 × 10 17 . It turns out that for fixed N and k, the numbers ♯D 2k−1 (l) are essentially (after normalization) distributed according to a discretisation of Beta distribution. For more details on this, see [16].

Dual transforms F D , F D d and the constant mass distributions
Throughout this section, we shall restrict ourselves to distributions u = N n=1 u n δ xn with real positive coefficients, u n > 0, and N n=1 u n < 1. As before, let 0 < x 1 < x 2 < . . . < x n < 1. The central point will be the fact that the roles of variables u n and ∆x n in F D are, to an extent, symmetrical. For easier reading, we denote: ξ n = ∆x n in this section.
We can define two transforms of the pair of functions (x, v) by and R(v) = R(v, 1 πz ) and E(x) = E(x, 1 πz ) = diag(e ix , e −ix ). Let v be an arbitrary function from (29) and let v r be its restriction on {1, 2, . . . , N}. Then we denote where F D [u](z) is the transform studied above and introduced in Definition 1.
is called the dual nonlinear Fourier transform, associated to F D [ξ, u](z). In variables (ξ, u) we shall denote it by F D d [ξ, u](z). Note that to obtain F D d from F D , we have to choose a value u 0 . The value u N +1 is then determined by the condition v(N + 1) = 1. A more explicit expression of the dual transform is which is analogous to the expression (16). Let now where ζ = −z and R(ξ n ) = cos ξ n −i sin ξ n −i sin ξ n cos ξ n = cos ξ n I + sin ξ n L.
Similar to page 13, we obtain from (33): and, taking into account that M − 1 = N, we have Here F D denotes the linear discrete Fourier transform. From here on, we follow the algorithm of Theorem 2 without changes to obtain all ξ n , and x n = n−1 k=0 ξ k . ✷

Conclusions
We described two algorithms for evaluating the inverses of the discretisation F E and F D of the nonlinear Fourier transform F . The two algorithms have similar structures. In both cases, we apply the inverse linear Fourier transform on the given matrix function A (z), where z is discrete in the case of F E and continuous in the case of F D . This yields highly nonlinear systems for the unknown discrete function {u n } n=0,...,N −1 , and also for {x n } n=0,...N −1 in the case of F D . However, in both cases, the "last" equation in these systems is trivial and yields the unknowns u N −1 or (x N −1 , u N −1 ). The structures of F E and F D allow us to find the other unknowns recursively. Our primary interest in these algorithms is their mathematical content and not numerical efficiency. We nevertheless comment shortly on the computational complexity of the algorithms. Of these two, the case of F D is numerically more demanding. We have to calculate good approximations of continuous linear Fourier transforms N times. These would, of course, be calculated using FFT, but the data size of these FFTs must be considerably larger than N.
The inverse of F E can be solved by N evaluations of the FFT with data size N: so the order of complexity is N 2 log N. Suppose we have a situation in which where H denotes the hyperfactorial function. This modification is an improvement. For N = 1000, the difference between the complexities of the modified and unmodified is roughly 3, 7 × 10 6 .
State of the art numerical methods concerning the nonlinear Fourier transform can be found in [12] and the references therein.
Our discretizations do not stem from some integrable discretisation of the integrable systems. However, in our opinion, they shed light on the structure of nonlinear Fourier analysis. Also, nonintegrable tools are often used in the study of integrable systems, see [3], [17] and many other references.
In a forthcoming paper, we intend to investigate what the results of this paper can tell us about the continuous nonlinear Fourier transformation and its applications in nonlinear dynamics. We also intend to explore the usefulness of F D and (F D ) −1 in the nonlinear Fourier analysis of distributions.

Acknowledgement
The research for this paper was supported in part by the research project, ARRS: J1-3005 from ARRS, Republic of Slovenia.