Schedule for: 16w5038 - Sparse Interpolation, Rational Approximation and Exponential Analysis

A welcome drink will be served at the hotel.

We consider the interpolation of an $n$-variate exponential sum $$ F(x_1, \ldots, x_n) = \sum_{j=1}^t c_j e^{f_{j,1} x_1 + f_{j,2} x_2 + \cdots + f_{j,n} x_n}. $$ In the univariate case, $n=1$, there is an entire branch of algorithms, which can be traced back to Prony's method dated in the 18th century, devoted to the recovery of the $2t$ unknowns, $c_1, \ldots, c_t, f_{1}, \ldots, f_{t}$, in $$ F(x) = \sum_{j=1}^t c_j e^{f_{j} x}. $$ In the multivariate case, $n>1$, it remains an active research topic to identify and separate distinct multivariate parameters from results computed by a Prony-like method from samples along projections. On top of the above, if the $f_{j,k}$ are allowed to be complex, the evaluations of the imaginary parts of distinct $f_{j,k}$ can also collide. This aliasing phenomenon can occur in either the univariate or the multivariate case. Our method interpolates $F(x_1, \ldots, x_n)$ from $(n+1) \cdot t$ evaluations. Since the total number of parameters $c_j$ and $f_{j,k}$ is exactly $(n+1) \cdot t$, we interpolate $F(x_1, \ldots, x_n)$ from the minimum possible number of evaluations. The method can also be used to recover the correct frequencies from aliased results. Essentially, we offer a scheme that can be embedded in any Prony-like algorithm, such as the least squares Prony version, ESPRIT, the matrix pencil approach, etc., thus can be viewed as a new tool offering additional possibilities in exponential analysis. This is joint work with A. Cuyt (Antwerp)

We propose a deterministic stable sparse FFT algorithm for the inverse discrete Fourier transform of length $N$. If the vector ${\bf x} \in {\bf R}^{N}$ is non-negative and $M$-sparse, we need at most $\min \{M \log(N), N \}$ Fourier values and (for $M^{2} \le N$) only ${\cal O}(M^{2} N)$ arithmatical operations to recover ${\bf x}$. The algorithm works iteratively and does not incorporate any a priori knowledge on the sparsity $M$ of ${\bf x}$. Each iteration step only involves the solution of a linear system of size at most $M$. Choosing the Fourier samples for reconstruction adaptively at each iteration level, we can develop a strategy to ensure that the coefficient matrix in the linear system is well-conditioned. These results have been obtained jointly with Katrin Wannenwetsch (University of G\"ottingen), Annie Cuyt and Wen-Shin Lee (University of Antwerp).

We derive a new method for optimal $\ell^2$-approximation of discrete signals on $\mathbb{N}_0$ whose entries can be represented as an exponential sum of finite length. Our approach employs Prony's method in a first step to recover the exponential sum that is determined by the signal. In the second step we use the theory of Adamjan, Arov and Krein (AAK-theory) to derive an algorithm for computing a shorter exponential sum that approximates the original signal in the $\ell^{2}$-norm well. AAK-theory originally determines best approximations of bounded periodic functions in Hardy-subspaces. We reformulate these ideas for our purposes and present the theory using only basic tools from linear algebra and Fourier analysis. The new algorithm is tested numerically in different examples. These results have been obtained jointly with Gerlind Plonka (Institute of Numerical and Applied Mathematics, University of G\"ottingen).

Motivated by a problem from physics, in~1795 Prony proposed a method to reconstruct a linear combination of exponentials given only a finite set of samples of sufficient cardinality. Prony's method translates this into the problem of solving a univariate polynomial equation. Variants and generalizations, in particular to the multivariate case, have been studied recently by several authors. We present a generalization of Prony's method to the multivariate case that translates the problem into a system of multivariate polynomial equations. In particular, we consider some of its algebraic properties. Joint work with Stefan Kunis, Thomas Peter, H. Michael M\"oller and Tim R\"omer

While there exists a growing literature dealing with the interpolation of sparse polynomials, fewer references address problems related to the evaluation of such polynomials. A notable exception is a result by Canny, Kaltofen and Lakshman, which shows that this can be done in amortized polylogarithmic time, provided we use very specific evaluation points. In this talk, I will discuss some questions around the evaluation of sparse polynomials for general evaluation points. This is joint work with Dorian Nogneng and Larry Li.

Ankur Moitra in his paper at STOC 2015 has given an in-depth analysis of how oversampling improves the conditioning of the arising Prony systems for sparse interpolation and signal recovery from numeric data. Moitra assumes that oversampling is done for a number of samples beyond the actual sparsity of the polynomial/signal. We give an algorithm that can be used to compute the sparsity and estimate the minimal number of samples needed in numerical sparse interpolation. The early termination strategy of polynomial interpolation has been incorporated in the algorithm: by oversampling at a small number of extra sample points we can diagnose that the sparsity has not been reached. Our algorithm still has to make a guess, the number $\zeta$ of oversamples, and we show by example that if $\zeta$ is guessed too small, premature termination can occur, but our criterion is numerically more accurate than that by Kaltofen, Lee and Yang (Proc.\ SNC 2011, ACM) but not as efficiently computable. For heuristic justification one has available the multivariate early termination theorem by Kaltofen and Lee (JSC vol.\ 36(3--4) 2003) for exact arithmetic, and the numeric Schwartz-Zippel Lemma by Kaltofen, Yang and Zhi (Proc.\ SNC 2007, ACM). A main contribution here is a modified proof of the Theorem by Kaltofen and Lee that permits starting the sequence at the point $(1,\ldots,1)$, for scalar fields of characteristic $\ne 2$ (in characteristic~$2$ counter-examples are given). Joint work with Z. Hao (KLMM, Beijing) and E. Kaltofen (NC state)

Several problems with applications in signal processing, functional approximation involve series representations, which are sums of polynomial-exponential functions. Their decomposition often allows to recover the structure of the underlying signal or data. Using the duality between polynomials and formal power series, we will describe the correspondence between polynomial-exponential series and Artinian Gorenstein algebras. We will give a generalisation of Kronecker theorem to the multivariate case, showing that the symbol of a Hankel operator of finite rank is a polynomial-exponential series and by connecting the rank of the Hankel operator with the decomposition of the symbol. We will describe algorithms for computing the polynomial-exponential decomposition of a series from its first coefficients, exploiting eigenvector methods for solving polynomial equations and a Gram-Schmidt orthogonalization process. A key ingredient of the approach is the flat extension criteria, which leads to a rank condition for a Carath{\'e}odory-Fej{\'e}r decomposition of multivariate Hankel matrices. The approach will be illustrated in different contexts: sparse interpolation of polylog functions from values, decomposition of polynomial-exponential functions from values, tensor decomposition, decomposition of symbols of convolution (or cross-correlation) operators of finite rank, reconstruction of measures as weighted sums of Dirac measures from moments. Numerical computation will be given to show the behaviour of the method.

A quadrature is an approximation of the definite integral of a function by a weighted sum of function values at specified points, or nodes, within the domain of integration. Gaussian quadratures are constructed to yield exact results for any polynomials of degree $2r-1$ or less by a suitable choice of $r$ nodes and weights. Cubature is a generalization of quadrature in higher dimension. Constructing a cubature amounts to find a linear form\\[0pt] \hspace*{3cm}$\displaystyle \Lambda : \R[x] \to \R, p \mapsto \sum_{j=1}^r a_j \, p(\xi_j) $ \\[0pt] from the knowledge of its restriction to $\R[x]_{\leq d}$. The unknowns are the weights $a_j$ and the nodes $\xi_j$. An approach based on moment matrices was proposed in \cite{Fialkow05,Lasserre12,Abril-Bucero16c}. We give a basis-free version in terms of the Hankel operator $\mathcal{H}$ associated to $\Lambda$. The existence of a cubature of degree $d$ with $r$ nodes boils down to conditions of ranks and positive semidefiniteness on $\mathcal{H}$. We then recognize the nodes as the solutions of a generalized eigenvalue problem. Standard domains of integration are symmetric under the action of a finite group. It is natural to look for cubatures that respect this symmetry \cite{Cools97,Gatermann89}. They are exact for all anti-symmetric functions beyond the degree of the cubature. Introducing adapted bases obtained from representation theory, the symmetry constraint allows to block diagonalize the Hankel operator $\mathcal{H}$. The size of the blocks is explicitly related to the orbit types of the nodes. From the computational point of view, we then deal with smaller-sized matrices both for securing the existence of the cubature and computing the nodes. This is joint work with Mathieu Collowald.

Tour of Monte Alban for morning

Given the scalar sequence $\{f_m\}^\infty_{m=0}$ that satisfies $$f_m=\sum^k_{j=1}a_j\zeta_j^m,\quad m=0,1,\ldots,$$ where $a_j, \zeta_j\in \mathbb{C}$, $j=1,\ldots,k,$ are independent of $m$, and $\zeta_j$ are distinct, the algorithm of Prony concerns the determination of the $a_j$ and the $\zeta_j$ from the $f_m$. In this talk, we discuss ways of extending Prony's algorithm to sequences of vectors $\{\mathbf{f}_m\}^\infty_{m=0}$ in $\mathbb{C}^N$ that satisfy $$\mathbf{f}_m=\sum^k_{j=1}\mathbf{a}_j\zeta_j^m, \quad m=0,1,\ldots,$$ where $\mathbf{a}_j\in\mathbb{C}^N$ and $\zeta_j\in\mathbb{C}$, $j=1,\ldots,k.$ We consider different approaches that enable us to determine the $\mathbf{a}_j$ and $\zeta_j$ whether the $\mathbf{a}_j$ are linearly independent or not. In so doing, we concentrate on extensions that take into account the possibility of the components of the $\mathbf{a}_j$ being coupled. The methods suggested here can be extended to vector sequences in infinite dimensional spaces in a straightforward manner.

We construct a frame of wave-atom-like multi-scale Gaussian beams following the dyadic parabolic decomposition of phase space. With this frame we generate a parametrix for the wave equation and study the approximate solution of the Cauchy initial value problem. We provide a natural estimate in terms of the scale parameter which is of order $-\sfrac{1}{2}$. We then introduce the notion of well-spread families of multi-scale Gaussian beam parameters. These families enable the generation of approximate solutions to boundary value problems with the natural error estimate. We show an application of this result in reverse-time-migration type imaging with wavefield reconstruction of seismic array data. joint work with Michele Berra and Jos\'{e}-Luis Romero

In this talk, integral representation formula (IRF) for poroelastic materials and two-phase elastic composite materials will be presented. Poroelastic materials refer to two-phase composites of solid and fluid. Wave propagation with wave-length much longer than the scale of the micro-geometry can be described by the Biot-Johnson-Kaplik-Dashen (Biot-JKD) equations, which can be considered as the homogenized wave equation for composite materials. All the coefficients in the equations depend on the microgeometry. Complex function theory, esp. Herglotz functions, provides an elegant and effective way for quantifying this dependence. This link will be one of the main topics of the talk. Another main topic in the presentation will be the numerical algorithm that utilizes the analytic structure implied by the IRF of the dynamic permeability/tortuosity functions of poroelastic materials to approximate the memory term in the wave equations, which describes the drag force due to solid-fluid interaction at the microscopic level, by sums of Prony type.

In this talk we consider the problem of working with rational functions in a numeric environment. A particular problem when modelling with such functions is the existence of Froissart doublets, where a zero is close to a pole. We discuss three different parameters which allow one to monitor the absence of Froissart doublets for a given general rational function. These include the euclidean condition number of an underlying Sylvester-type matrix, a parameter for determining coprimeness of two numerical polynomials and bounds on the spherical derivative. We show that our parameters sharpen those found in previous papers (1) and (2) (1) B. Beckermann and G. Labahn, When are two numerical polynomials relatively prime? {\em Journal of Symbolic Computation} {\bf 26} (1998) 677-689. (2) B. Beckermann and A. Matos, Algebraic properties of robust Pad\'e approximants, {\em Journal of Approximation Theory }{\bf 190} (2015) 91-115. This is joint work with B. Beckermann (Lille) and G. Labahn (Waterloo)

In 1853, Sylvester introduced double sum expressions for two finite sets of indeterminates and subresultants for univariate polynomials, showing the relationship between both notions in several but not all cases. Here we show how Sylvester's double sums can be interpreted in terms of symmetric multivariate Lagrange interpolation, allowing to recover in a natural way the full description of all cases. We will also report on preliminary results on extensions to symmetric multivariate Hermite interpolation, and applications. Joint work with T. Krick and M. Valdettaro (Universidad de Buenos Aires)

This work arose out of joint work with Jonathan M. Borwein (1951--2016), undertaken while he was the Distinguished Visiting Scholar in Residence at Western (April-August 2016). That work was a survey of articles in the American Mathematical Monthly on the~$\Gamma$ function, of which there were nearly fifty. The survey identified some gaps: visualization, computation, and perhaps most surprisingly \textsl{nothing} on the functional inverse of~$\Gamma$, which we denote~$y=\invG(x)$ if~$x=\Gamma(y)$. This is multivalued with infinitely many branches. In this talk, I show that one can invert Stirling's original approximation (which is more accurate than the popular approximation by that name) to give a surprisingly accurate asymptotic formula for~$\invG$.

We consider multivariate trigonometric polynomials \begin{align*} p(x)&=\sum_{k\in I}\hat{p}_{k}\mathrm{e}^{2\pi\mathrm{i} k\cdot x},\\ p&\in\Pi_I:=\operatorname{span}\{\mathrm{e}^{2\pi\mathrm{i}k\circ}\colon k\in I\}, \end{align*} with arbitrary frequency index set $I\subset\mathbb{Z}^d$, i.e., $I$ has finite cardinality and there is not required a specific structure of $I$. Our aim are spatial discretizations such that all trigonometric polynomials belonging to $\Pi_I$ can be uniquely reconstructed from their sampling values at the discretization nodes. In the recent years, we used single rank-1 lattices as spatial discretizations and developed a workable component--by--component strategy in order to determine such rank-1 lattices for given frequency index sets $I$. The structure of rank-1 lattices allows for the efficient computation of the corresponding discrete Fourier transforms. Moreover, the computation is stable in the sense of well-conditioned Fourier matrices. However, detailed investigations on the structure of rank-1 lattices reveal a unsuitable huge oversampling in specific situations. More precise, the number $M$ of necessarily needed rank-1 lattice nodes may be nearly $\frac{1}{4}|I|^2$ in the worst case. Motivated by the sparse grids idea of Smolyak, we would like to combine multiple rank-1 lattices to a spatial discretization similar to the combination of small full grids to sparse grids. We collect different interpretations of our approach and present the most urgently and---at least for the author---outstanding issues. Many similarities to the interpolation of a straight-line program as a sparse algebraic polynomial suggest that some of these issues are easily solvable by adapting results of the SLP community.

TBA

In sparse polynomial interpolation, we can improve the conditioning of the problem by splitting the overall problem into several subproblems of smaller sizes. In this talk, we extend this idea to a more general scenario where the exponents are not necessarily discrete. We investigate the size-reduction strategy for this setting, which also leads us to a specific case of digital filters. This is joint work with Annie Cuyt and Wen-shin Lee (Antwerp)

The study of quaternions, quaternion polynomials and quaternion matrices has a long history back to 1940's, and has recently gained more and more interest in many areas such as quantum mechanics and signal processing. In this talk, we discuss the computing approximate greatest common right divisors of quaternion polynomials. We explore some properties of Sylvester matrices in quaternion case. The difficult thing is that there are no ``good" definitions of the determinate for matrices over division rings. Base on the singular value decomposition for quaternion matrices, we present algorithms to compute approximate GCRDs for quaternion polynomials.