An “algebraic” reconstruction of piecewisesmooth functions from integral measurements
Abstract.
This paper presents some results on a wellknown problem in Algebraic Signal Sampling and in other areas of applied mathematics: reconstruction of piecewisesmooth functions from their integral measurements (like moments, Fourier coefficients, Radon transform, etc.). Our results concern reconstruction (from the moments or Fourier coefficients) of signals in two specific classes: linear combinations of shifts of a given function, and “piecewise finite functions” which satisfy on each continuity interval a linear differential equation with polynomial coefficients. In each case the problem is reduced to a solution of a certain type of nonlinear algebraic system of equations (“Pronytype system”). We recall some known methods for explicitly solving such systems in one variable, and provide extensions to some multidimensional cases. Finally, we investigate the local stability of solving the Pronytype systems.
Key words and phrases:
nonlinear moment inversion, Algebraic Sampling, Prony system2000 Mathematics Subject Classification:
94A12, 94A201. Introduction
It is well known that the error in the best approximation of a function by an th degree Fourier polynomial is of order . The same holds for algebraic polynomial approximation and for other basic approximation tools (see e.g. [21, Chapters IV, VI], and [32, Vol.I, Chapter 3, Theorem 13.6]). However, for with singularities, in particular, with discontinuities, the error is much larger: its order is . Considering the socalled Kolmogorow width of families of signals with moving discontinuities one can show that any linear approximation method provides the same order of error, if we do not fix a priori the discontinuities’ position (see [11], Theorem 2.10). Another manifestation of the same problem is the “Gibbs effect”  a relatively strong oscillation of the approximating function near the discontinuities. Practically important signals usually do have discontinuities, so the above feature of linear representation methods presents a serious problem in signal reconstruction. In particular, it visibly appears near the edges of images compressed by JPEG, as well as in the noise and low resolution of the CT and MRI images.
Recent nonlinear reconstruction methods, in particular, Compressed Sensing ([6, 7]) and Algebraic Sampling ([8, 19, 23, 10, 15]), address this problem in many cases. Both approaches appeal to an a priori information on the character of the signals to be reconstructed, assuming their “simplicity” in one or another sense. Compressed sensing assumes only a sparse representation in a certain (wavelets) basis, and thus it presents a rather general and “universal” approach. Algebraic Sampling usually requires more specific a priori assumptions on the structure of the signals, but it promises a better reconstruction accuracy. In fact, we believe that ultimately the Algebraic Sampling approach has a potential to reconstruct “simple signals with singularities” as good as smooth ones. The most difficult problem in this approach seems to be estimating the accuracy of solution of the nonlinear systems arising.
Our purpose in this paper is to further substantiate the Algebraic Sampling approach. On one hand, we present two algebraic reconstruction methods for generic classes of signals. The first one is reconstruction of combinations of shifts of a given function from the moments and Fourier coefficients (Section 2). The second one concerns piecewise finite moment inversion (Section 4). On the other hand, we consider some typical nonlinear systems arising in these reconstruction schemes. We describe the methods of their solution (Section 3), and provide some explicit bounds on their local stability (Section 5). We also present results of some numerical experiments in Section 6.
Our ultimate goal may be stated in terms of the following conjecture (which seems to be supported also by the results of [9, 18, 14, 30, 23]):
There is a nonlinear algebraic procedure reconstructing any signal in a class of piecewise functions (of one or several variables) from its first Fourier coefficients, with the overall accuracy of order . This includes the discontinuities’ positions, as well as the smooth pieces over the continuity domains.
Recently [4] we have shown that at least half the conjectured accuracy can be achieved. However, the question of maximal possible accuracy remains open. Our results presented in this paper can be considered as an additional step in this direction.
2. Linear combinations of shifts of a given function
Reconstruction of this class of signals from sampling has been described in [8, 19]. We study a rather similar problem of reconstruction from the moments. Our method is based on the following approach: we construct convolution kernels dual to the monomials. Applying these kernels, we get a Pronytype system of equations on the shifts and amplitudes.
Let us restate a general reconstruction problem, as it appears in our specific setting. We want to reconstruct signals of the form
(1) 
where is a known function of , while the number of the shifts and the form (1) of the signal are known a priori. The parameters are to be found from a finite number of “measurements”, i.e. of linear (usually integral) functionals like polynomial moments, Fourier moments, shifted kernels, evaluation over some grid etc.
In this paper we consider only linear combinations of shifts of one known function . Reconstruction of shifts of several functions based on “decoupling” via sampling at zeroes of the Fourier transforms of the shifted functions is presented in [29].
In what follows , is a scalar index, while and are multiindices. Partial ordering of multiindices is given by
Assume that a multisequence of functions is fixed. We consider the measurements provided by
(2) 
Our approach now works as follows: given and we now try to find an “convolution dual” system of functions in a form of certain “triangular” linear combinations
(3) 
More accurately, we try to find the coefficients in (3) in such a way that
(4) 
We shall call a sequence satisfying (3), (4)  convolution dual to . Below we shall find explicitly convolution dual systems to the usual and exponential monomials.
We consider a general problem of finding convolution dual sequences to a given sequence of measurements as an important step in the reconstruction problem. Notice that it can be generalized by dropping the requirement of a specific representation (3): . Instead we can require only that be expressible in terms of the measurements sequence . Also in (4) can be replaced by another a priori chosen sequence . This problem leads, in particular, to certain functional equations, satisfied by polynomials and exponential functions, as well as exponential polynomials and some kinds of elliptic functions (see [28]).
Now we have the following result:
Theorem 1.
Let a sequence be convolution dual to . Define by Then the parameters and in (1) satisfy the following system of equations (“generalized Prony system”):
(5) 
According to Theorem 1, in order to reduce the reconstruction problem with the measurements (2) and for signals of the a priori known form (1) to a solution of the generalized Prony system we have to find an  convolution dual system to the measurements kernels . In fact we need only the coefficients . Having these coefficients, we compute and get system (5).
Solvability of system (5) and robustness of its solution depend on the measurements kernels . Specific examples are considered below.
2.1. Reconstruction from moments
We are given a finite number of moments of a signal as in (1) in the form
(6) 
So here for each multiindex . We look for the dual functions satisfying the convolution equation
(7) 
for each multiindex . To solve this equation we apply Fourier transform to both sides of (7). Assuming that and that has the derivatives of all the orders at we find (see [28]) that there is a unique solution to (7) provided by
(8) 
where
The assumption is essential in the construction of the convolution dual system for the monomials as well as for other measurement kernels (as well as in the study of the shifts of a given function in general). The above calculation can be applied, with proper modifications, in more general situations (see [28]). On the other hand, the assumption of differentiability of at zero is not very restrictive, in particular, if we work with signals with finite support.
Returning to the moments reconstruction, we set the generalized polynomial moments as
(9) 
and obtain, as in Theorem 1, the following system of equations:
(10) 
This system is called “multidimensional Prony system”. It appears in numerous problems of theoretical and applied mathematics. In Section 3 below we recall one of the classical methods for its solution in onedimensional case, and describe, under certain assumptions, a method for its solution in several dimensions. Local stability of the solutions of the onedimensional Prony system is discussed in Section 5.
2.2. Fourier case
The signal has the same form as in section 2.1: . The measurement kernels we now choose as the harmonics . So our measurements are the Fourier coefficients
Let us assume now that satisfies the condition for all integer . Then the Fourier harmonics turn out to form essentially selfdual system: indeed, taking we get immediately that
According to our general scheme we put now and get a system of the form
(11) 
where . This is once more a multidimensional Prony system as in (10), with the nodes on the complex unit circle.
2.3. Further extensions
The approach above can be extended in the following directions:

Reconstruction of signals built from several functions or with the addition of dilations also can be investigated (a perturbation approach where the dilations are approximately 1 is studied in [26]).

Further study of “convolution duality” in [26] provides a certain extension of the class of signals and measurements allowing for a closed  form reconstruction.
3. Prony system
3.1. Onedimensional case
Let us start with the classical case of onedimensional Prony system:
(12) 
This system appears in many branches of mathematics (see [20]). There are several solution methods available, for instance direct nonlinear minimization (see e.g. [9]), the polynomial realization (original Prony method, [24]) or the statespace approach ([25]). Let us give a sketch of a method based on Padé approximation techniques which is rather close to the original Prony method. To simplify the presentation we shall assume that we know a priori that all the nodes are pairwise distinct. General case is treated similarly ([22], see also [31]).
Consider a “moments generating function” Assuming the equations (12) are satisfied for all and summing up the geometric progressions we get
(13) 
So is a rational function of degree tending to zero as . The unknowns and in (12) are the poles and the residues of , respectively.
Now in order to find explicitly from the first moments we use the Padé approximation approach (see [22]): write as with polynomials and of degrees and , respectively.
Multiplying by we have . Now equating the coefficients on both sides we get the following system of linear equations:
………………………..
…………………………
The rest of the equations in this system are obtained by further shifts of the indices of the moments, and so they form a Hankeltype matrix.
Now, being a rational function of degree , is uniquely defined by its first Taylor coefficients (the difference of two such functions cannot vanish at zero with the order higher than ). We conclude that the linear system consisting of the first equations as above is uniquely solvable up to a common scaling of and (of course, this fact follows also form a general Padé approximation theory  see [22]).
Now a solution procedure for the Prony system can be described as follows:
1. Solve a linear system of the first equations as above (with the coefficients  the known moments ) to find the moments generating function in the form .
2. Represent in a standard way as the sum of elementary fractions . (Equivalently, find poles and residues of ). Besides algebraic operations, this requires just finding the roots of the polynomial . Then form the unique solution of the Prony system (12).
3.2. Multidimensional case
Let us recall our multidimensional notations: , is a scalar index, while and are multiindices. Partial ordering of multiindices is given by
Let With the above notation, the multidimensional Prony system has exactly the same form as in the onedimensional case:
(14) 
Exactly as above we get that for the moments generating function is a rational function of degree of the form
(15) 
Representing as we get exactly in the same way as above an infinite system of linear equations for the coefficients of and , with a Hankeltype matrix formed by the moments . By the same consideration as above, after we take enough equations in this system the solution is unique up to a rescaling (see [1, 22] and references therein).
However, from this point the multidimensional situation becomes essentially more complicated. While in dimension one can be, essentially, any rational function of degree (naturally represented as the sum of elementary fractions), in several variables has a very special form given by (15). This fact can be easily understood via counting the degrees of freedom: the number of unknowns in (14) is while a rational function of degree in variables has much more coefficients than that, for . It would be desirable to use as many equations from (14) as the number of unknowns, but the method outlined above ignores a specific structure of and requires as many equations as in a Padé reconstruction for a general rational function of degree .
We treat this problem in [28] analyzing the structure of singularities of and on this base proposing a robust reconstruction algorithm. Let us give here one simple special case of this algorithm, which separates the variables in the problem.
3.2.1. Separation of variables in the multidimensional Prony system
Let us assume that we know a priori that the solution of multidimensional Prony system (14) is such that all the coordinates of the points are pairwise distinct. Moreover, we assume that for . Under these assumptions we proceed as follows:
Consider “partial moment generating functions” , defined by
(16) 
where is a multiindex with for and otherwise. We have the following simple fact:
Proposition 2.
is a onedimensional moments generating function of the form
(17) 
It coincides with the restriction of to the th coordinate axis in .
Proof.
Let us evaluate along the th coordinate axis, that is on the line with as above and . We get
which is the moments generating function (17). Now, to express through the multidimensional moments we notice that any monomial , with vanishes identically on the th coordinate axis. Hence
This shows that and completes the proof of the proposition. ∎
Now applying the method described in Section (3.1) above we find for each the coordinates and (repeatedly) the coefficients . It remains to arrange these coordinates into the points . This presents a certain combinatorial problem, since Prony system (14) is invariant under permutations of the index . Under the assumptions above we proceed as follows: for each we have obtained the (unordered) collection of the pairs . By assumptions for . Hence we can arrange in a unique way all the pairs into the string This gives us the desired solution of multidimensional Prony system (14).
Notice that the assumption for is essential here. Indeed, for and we have for on each of the coordinate axes. So the (unique up to permutations of the index ) solution of the Prony system cannot be reconstructed from these moments only.
Another remark is that the separation of variables as described above requires knowledge of moments ( on each of the coordinate axes). This is almost twice more than unknowns. This number can be significantly reduced in some cases. See [28] for further investigation in both of these directions.
Stability estimates for the solution of onedimensional Prony system (Section 5) can be easily extended to the case considered in the present subsection. We do not give here explicit statement of this result.
4. Reconstruction of piecewise finite functions from moments
In this section we present an overview of our previous findings on algebraic reconstruction of a certain general class of signals. See [3] for the complete details.
Let consist of “pieces” with jump points
Furthermore, let satisfy on each continuity interval some linear homogeneous differential equation with polynomial coefficients: where
(18) 
Each may be therefore written as a linear combination of functions which are a basis for the space :
(19) 
We term such functions “piecewise finite”. Many “simple” functions may be represented in this framework, such as polynomials, trigonometric polynomials and algebraic functions.
The sequence is given by the usual moments
Piecewise finite Reconstruction Problem.
Given and the moment sequence of a piecewise finite function , reconstruct all the parameters .
The above problem can be solved as follows.

It turns out that the moment sequence of every piecewise finite function satisfies a linear recurrence relation, such that the coefficients of annihilating every piece of and the discontinuity locations can be recovered from the moments by solving a linear systems of equations plus a nonlinear step of polynomial roots finding. As we shall explain below, in many cases this step is equivalent to solving a certain generalized form of the previously mentioned Prony system (12).

The function itself can be subsequently reconstructed by numerically calculating the basis for the space and solving an additional linear system of equations to recover the coefficients .
The above algorithm has been tested on reconstruction of piecewise polynomials, piecewise sinusoids and rational functions  see [3] for details. The results reported there were relatively accurate for low noise levels. In this paper we continue to explore the numerical stability of the method  see Sections 5 and 6.
Now let us show how the Prony system arises in the piecewise finite reconstruction method. Consider the distribution . It is of the form (see [3, Theorem 2.12])
(20) 
where are the discontinuity locations, are the associated “jump magnitudes” which depend on the values , and is the Dirac function. In particular, when is just a piecewise constant function, we have and so
(21) 
where . Let us now apply the moment transform to the equations (21) and (20). We get, correspondingly,
(22)  
(23) 
The system (22) is of course identical to the previously considered (12). However, the following question arises: how are the numbers in the lefthand side of (22) and (23) related to the known quantities ? It turns out that the numbers are certain linear combinations of these known moments, with coefficients which are determined by in a welldefined way. The conclusion is as follows: if the operator which annihilates every piece of is known (but the other parameters are not^{1}^{1}1This assumption is realistic, for instance when reconstructing piecewisepolynomials or piecewisesinusoids.), then one can recover the discontinuity locations of by solving the Pronylike system (23). We term the system (23) “confluent Prony system”. It can be solved in a similar manner to the standard Prony system. A unique solution exists whenever all the ’s are distinct and the highest coefficients are nonzero. We provide the details in [2].
5. Accuracy of solving Pronylike systems
In the preceding sections we have shown that several algebraic methods for nonlinear reconstruction can be reduced to solving certain types of systems of equations, the most basic one of which is the Prony system (12), (22). A crucial factor for the eventual applicability of the reconstruction methods is the sensitivity of solving these systems to measurement errors.
In this section we consider two such systems  (22) and (23) and provide some theoretical results regarding the local sensitivity of their solution to noise. These results will hopefully enable further “global” analysis.
We consider the following question: if the measurements in the lefthand side of (22) are known with error at most , how well can one recover the unknowns ? We regard this stability problem to be absolutely central in Algebraic Sampling. To our best knowledge, no general treatment of this problem exists, therefore we consider our results below to be a step in this direction.
For simplicity, let us assume that the number of equations in (22) equals the number of unknowns, which in this case is . Let us further consider the “measurement map” given by (22) (we call it the “Prony map”):
Then, one possible answer (by no means a complete one) to the question posed above can be given in terms of the local Lipschitz constant of the “solution map” , whenever this inverse is defined. We then obtain the following result.
Theorem 3.
Let be the exact unperturbed moments of the model (21). Assume that all the ’s are distinct and also for . Now let be perturbations of the above moments such that . Then, for sufficiently small , the perturbed Prony system has a unique solution which satisfies:
where is an explicit constant depending only on the geometry of .
Proof.
Consider the Jacobian determinant of the map . It is easily factorized as follows:
where is a block
The matrix is a special case of the socalled confluent Vandermonde matrix, wellknown in numerical analysis ([5, 12, 13, 16]). In particular, the paper [12, Theorem 3] contains the following estimate for the norm of :
where
Now since the ’s are distinct and , the Jacobian is nonsingular and so in a sufficiently small neighbourhood of the exact solution, the map is approximately linear. By the inverse function theorem
and so taking completes the proof. ∎
By a similar technique with slightly more involved calculations we obtain the following result for the system (23).
Theorem 4.
Let be the exact unperturbed moments of the model (20), where . Assume that all the ’s are distinct and also for . Now let be perturbations of the above moments such that . Then, for sufficiently small , the perturbed confluent Prony system has a unique solution which satisfies:
where is a constant depending only on the nodes and the multiplicities .
Proof outline.
As before, we obtain a factorization of the Jacobian determinant as a product of a confluent Vandermonde matrix (defined in a similar manner but with each column having “confluencies”) and the block diagonal matrix where
Inverting this expression and taking completes the proof. ∎
6. Numerical experiments
We have tested the piecewise finite reconstruction method on a simple case of a piecewiseconstant signal. The implementation details are identical to those used in [3, Appendix]. As can be seen from Figure 1, the reconstruction is accurate even in the presence of mediumlevel noise.
We have already mentioned that for a piecewise constant function , the distribution is of the form (21). Therefore, the piecewise finite reconstruction problem for essentially reduces to solving the Prony system (22), and so the local estimates of Theorem 3 should apply in the case of a very small noise level. This prediction is partially confirmed by the results of our second experiment, presented in Figure 2. In particular, it can be seen that indeed , while the accuracy of other jump points does not depend on for . While this is certainly an encouraging result, more investigation is clearly needed in order to fully understand the dependence of the error on all the parameters of the problem. Such an investigation should, in our opinion, concentrate on the global structure of the mapping .
References
 [1] G. A. Baker, P. GravesMorris, Padé Approximants. Cambridge U.P., 1996.
 [2] D.Batenkov, Y.Yomdin, Accuracy of nonlinear reconstruction for discontinuous models, in preparation.
 [3] D.Batenkov, Moment inversion problem for piecewise finite functions, Inverse Problems, 25(10):105001, October 2009.
 [4] Batenkov, D., Yomdin, Y. Algebraic Fourier reconstruction of piecewise smooth functions, to appear in Mathematics of Computation.
 [5] Å. Björck and T. Elfving. Algorithms for confluent Vandermonde systems. Numerische Mathematik, 21(2):130–137, 1973.
 [6] E. J. Candes̀. Compressive sampling. Proceedings of the International Congress of Mathematicians, Madrid, Spain, 2006. Vol. III, 1433–1452, Eur. Math. Soc., Zurich, 2006.
 [7] D. Donoho, Compressed sensing. IEEE Trans. Inform. Theory 52 (2006), no. 4, 1289–1306.
 [8] P.L. Dragotti, M. Vetterli and T. Blu, Sampling Moments and Reconstructing Signals of Finite Rate of Innovation: Shannon Meets StrangFix, IEEE Transactions on Signal Processing, Vol. 55, Nr. 5, Part 1, pp. 17411757, 2007.
 [9] K. Eckhoff, Accurate reconstructions of functions of finite regularity from truncated Fourier series expansions, Math. Comp. 64 (1995), no. 210, 671–690.
 [10] M. Elad, P. Milanfar, G. H. Golub, Shape from moments—an estimation theory perspective, IEEE Trans. Signal Process. 52 (2004), no. 7, 1814–1829.
 [11] B. Ettinger, N. Sarig. Y. Yomdin, Linear versus nonlinear acqusition of stepfunctions, J. of Geom. Analysis, 18 (2008), 2, 369399.
 [12] W. Gautschi. On inverses of Vandermonde and confluent Vandermonde matrices. Numerische Mathematik, 4(1):117–123, 1962.
 [13] W. Gautschi. On inverses of Vandermonde and confluent Vandermonde matrices III. Numerische Mathematik, 29(4):445–450, 1978.
 [14] A. Gelb, E. Tadmor, Detection of edges in spectral data II. Nonlinear enhancement, SIAM J. Numer. Anal. 38 (2000), 13891408.
 [15] B. Gustafsson, Ch. He, P. Milanfar, M. Putinar, Reconstructing planar domains from their moments. Inverse Problems 16 (2000), no. 4, 1053–1070.
 [16] D. Kalman. The generalized Vandermonde matrix. Mathematics Magazine, pages 15–21, 1984.
 [17] V. Kisunko, Cauchy type integrals and a moment problem. C.R. Math. Acad. Sci. Soc. R. Can. 29 (2007), no. 4, 115–122.
 [18] G. Kvernadze, T. Hagstrom, H. Shapiro, Locating discontinuities of a bounded function by the partial sums of its Fourier series., J. Sci. Comput. 14 (1999), no. 4, 301–327.
 [19] I. Maravic and M. Vetterli, Exact Sampling Results for Some Classes of Parametric NonBandlimited 2D Signals, IEEE Transactions on Signal Processing, Vol. 52, Nr. 1, pp. 175189, 2004.
 [20] Y.I. Lyubich. The SylvesterRamanujan System of Equations and The Complex Power Moment Problem. The Ramanujan Journal, 8(1):23–45, 2004.
 [21] Natanson, I., Constructive Function Theory (in Russian), Gostekhizdat, 1949.
 [22] E. M. Nikishin, V. N. Sorokin, Rational Approximations and Orthogonality, Translations of Mathematical Monographs, Vol 92, AMS, 1991.
 [23] P. Prandoni, M. Vetterli, Approximation and compression of piecewise smooth functions, R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci. 357 (1999), no. 1760, 2573–2591.
 [24] R. Prony. Essai experimental et analytique. J. Ec. Polytech.(Paris), 2:24–76, 1795.
 [25] B.D. Rao and K.S. Arun. Model based processing of signals: A state space approach. Proceedings of the IEEE, 80(2):283–309, 1992.
 [26] N. Sarig, Signal Acquisition from Measurements via Convolution Duals, preprint, Weizmann Institute, Dec. 2010.
 [27] N. Sarig, Y. Yomdin, Signal Acquisition from Measurements via NonLinear Models, C. R. Math. Rep. Acad. Sci. Canada Vol. 29 (4) (2007), 97114.
 [28] N. Sarig and Y. Yomdin, Multidimensional Prony systems for uniform and nonuniform sampling, in preparation.
 [29] N. Sarig and Y. Yomdin, Decoupling of a Reconstruction Problem for Shifts of Several Signals via Non Uniform Sampling, in preparation.
 [30] E. Tadmor, High resolution methods for time dependent problems with piecewise smooth solutions. Proceedings of the International Congress of Mathematicians, Vol. III (Beijing, 2002), 747–757, Higher Ed. Press, Beijing, 2002.
 [31] Y. Yomdin, Singularities in Algebraic Data Acquisition, Real and Complex Singularities, M. Manoel, M. C. Romero Fuster, C. T. C. Wall, Editors, London Mathematical Society Lecture Note Series, No. 380, 2010, 378394.
 [32] Zygmund, A., Trigonometric Series. Vols. I, II, Cambridge University Press, New York, 1959.