We present results of a long-term team collaboration of mathematicians and biologists. We focus on building a mathematical framework for the shape space constituted by a collection of homologous bones or teeth from many species. The biological application is to quantitative morphological understanding of the evolutionary history of primates in particular, and mammals more generally. Similar to the practice of biologists, we leverage the power of the whole collection for results that are more robust than can be obtained by only pairwise comparisons, using tools from differential geometry and machine learning. This paper concentrates on the mathematical framework. We review methods for comparing anatomical surfaces, discuss the problem of registration and alignment, and address the computation of different distances. Next, we cover broader questions related to cross-dataset landmark selection, shape segmentation, and shape classification analysis. This paper summarizes the work of many team members other than the authors; in this paper that unites (for the first time) all their results in one joint context, space restrictions prevent a full description of the mathematical details, which are thoroughly covered in the original articles. Although our application is to the study of anatomical surfaces, we believe our approach has much wider applicability.

Mathematics is very human—we all enjoy using our brains, thinking, making connections, and checking to see if our conclusions are true. For me, mathematics is also a way to experience creativity, which is what makes me happy. In this article, I want to share this delight with you, through the story of wavelets. Wavelets are mathematical tools that my colleagues and I developed, which are very useful in analyzing images and other signals like medical scans or audio files. I hope that, by the end of this article, you will see that mathematics is much more interesting and diverse than memorizing formulas and performing calculations. It is a creative endeavor that advances society and can make our lives richer and more beautiful.

This paper proposes a novel kernel-based optimization scheme to handle tasks in the analysis, e.g., signal spectral estimation and single-channel source separation of 1D non-stationary oscillatory data. The key insight of our optimization scheme for reconstructing the time-frequency information is that when a nonparametric regression is applied on some input values, the output regressed points would lie near the oscillatory pattern of the oscillatory 1D signal only if these input values are a good approximation of the ground-truth phase function. In this work, Gaussian Process (GP) is chosen to conduct this nonparametric regression: the oscillatory pattern is encoded as the Pattern-inducing Points (PiPs) which act as the training data points in the GP regression; while the targeted phase function is fed in to compute the correlation kernels, acting as the testing input. Better approximated phase function generates more precise kernels, thus resulting in smaller optimization loss error when comparing the kernel-based regression output with the original signals. To the best of our knowledge, this is the first algorithm that can satisfactorily handle fully non-stationary oscillatory data, close and crossover frequencies, and general oscillatory patterns. Even in the example of a signal produced by slow variation in the parameters of a trigonometric expansion, we show that PiPs admits competitive or better performance in terms of accuracy and robustness than existing state-of-the-art algorithms.

In the desire to quantify the success of neural networks in deep learning and other applications, there is a great interest in understanding which functions are efficiently approximated by the outputs of neural networks. By now, there exists a variety of results which show that a wide range of functions can be approximated with sometimes surprising accuracy by these outputs. For example, it is known that the set of functions that can be approximated with exponential accuracy (in terms of the number of parameters used) includes, on one hand, very smooth functions such as polynomials and analytic functions and, on the other hand, very rough functions such as the Weierstrass function, which is nowhere differentiable. In this paper, we add to the latter class of rough functions by showing that it also includes refinable functions. Namely, we show that refinable functions are approximated by the outputs of deep ReLU neural networks with a fixed width and increasing depth with accuracy exponential in terms of their number of parameters. Our results apply to functions used in the standard construction of wavelets as well as to functions constructed via subdivision algorithms in Computer Aided Geometric Design.

X-radiography (X-ray imaging) is a widely used imaging technique in art investigation. It can provide information about the condition of a painting as well as insights into an artist’s techniques and working methods, often revealing hidden information invisible to the naked eye. X-radiograpy of double-sided paintings results in a mixed X-ray image and this paper deals with the problem of separating this mixed image. Using the visible color images (RGB images) from each side of the painting, we propose a new Neural Network architecture, based upon ‘connected’ auto-encoders, designed to separate the mixed X-ray image into two simulated X-ray images corresponding to each side. This connected auto-encoders architecture is such that the encoders are based on convolutional learned iterative shrinkage thresholding algorithms (CLISTA) designed using algorithm unrolling techniques, whereas the decoders consist of simple linear convolutional layers; the encoders extract sparse codes from the visible image of the front and rear paintings and mixed X-ray image, whereas the decoders reproduce both the original RGB images and the mixed X-ray image. The learning algorithm operates in a totally self-supervised fashion without requiring a sample set that contains both the mixed X-ray images and the separated ones. The methodology was tested on images from the double-sided wing panels of the Ghent Altarpiece, painted in 1432 by the brothers Hubert and Jan van Eyck. These tests show that the proposed approach outperforms other state-of-the-art X-ray image separation methods for art investigation applications.

In this paper, we focus on X-ray images (X-radiographs) of paintings with concealed sub-surface designs (e.g., deriving from reuse of the painting support or revision of a composition by the artist), which therefore include contributions from both the surface painting and the concealed features. In particular, we propose a self-supervised deep learning-based image separation approach that can be applied to the X-ray images from such paintings to separate them into two hypothetical X-ray images. One of these reconstructed images is related to the X-ray image of the concealed painting, while the second one contains only information related to the X-ray image of the visible painting. The proposed separation network consists of two components: the analysis and the synthesis sub-networks. The analysis sub-network is based on learned coupled iterative shrinkage thresholding algorithms (LCISTA) designed using algorithm unrolling techniques, and the synthesis sub-network consists of several linear mappings. The learning algorithm operates in a totally self-supervised fashion without requiring a sample set that contains both the mixed X-ray images and the separated ones. The proposed method is demonstrated on a real painting with concealed content, Do na Isabel de Porcel by Francisco de Goya, to show its effectiveness.

Diffusion maps (DM) constitute a classic dimension reduction technique, for data lying on or close to a (relatively) low-dimensional manifold embedded in a much larger dimensional space. The DM procedure consists in constructing a spectral parametrization for the manifold from simulated random walks or diffusion paths on the data set. However, DM is hard to tune in practice. In particular, the task to set a diffusion time t when constructing the diffusion kernel matrix is critical. We address this problem by using the semigroup property of the diffusion operator. We propose a semigroup criterion for picking t. Experiments show that this principled approach is effective and robust.

Phase retrieval is known to always be unstable when using a frame or continuous frame for an infinite dimensional Hilbert space. We consider a generalization of phase retrieval to the setting of subspaces of L_2 which coincides with using a continuous frame for phase retrieval when the subspace is the range of the analysis operator of a continuous frame. We then prove that there do exist infinite dimensional subspaces of L_2 where phase retrieval is stable. That is, we give a method for constructing an infinite dimensional subspace Y⊆ L_2 such that there exists C≥ 1 so that min(f-g_L_2,f+g_L_2)≤ C |f|-|g| _L_2 for all f,g∈ Y. This construction also leads to new results on uniform stability of phase retrieval in finite dimensions. Our construction has a deterministic component and a random component. When using sub-Gaussian random variables we achieve phase retrieval with high probability and stability constant independent of the dimension n when using m on the order of n random vectors. Without sub-Gaussian or any other higher moment assumptions, we are able to achieve phase retrieval with high probability and stability constant independent of the dimension n when using m on the order of nlog(n) random vectors.

In this paper, we study the stability of phase retrieval problems via a family of locally stable phase retrieval frame measurements in Banach spaces, which we call “locally stable and conditionally connected” (LSCC) measurement schemes. For any signal f in the Banach space, we associate it with a weighted graph Gf, defined by the LSCC measurement scheme, and show that the phase retrievability of the signal f is determined by the connectivity of Gf. We quantify the phase retrieval stability of the signal by two common measures of graph connectivity: The Cheeger constant for real-valued signals, and algebraic connectivity for complex-valued signals. We then use our results to study the stability of two phase retrieval models. In the first model, we study a finite-dimensional phase retrieval problem from locally supported measurements such as the windowed Fourier transform. We show that signals “without large holes” are phase retrievable, and that for such signals in Rd the phase retrieval stability constant grows proportionally to d1/2, while in Cd it grows proportionally to d. The second model we consider is an infinite-dimensional phase retrieval problem in a shift-invariant space. In infinite-dimension spaces, even phase retrievable signals can have the Cheeger constant being zero, and hence have an infinite stability constant. We give an example of signals with monotone polynomial decay which has the Cheeger constant being zero, and an example with exponential decay which has a strictly positive Cheeger constant.

We address the structure identification and the uniform approximation of sums of ridge functions $f(x)=\sum _{i=1}^m g_i(\langle a_i,x\rangle )$ on ${\mathbb{R}}^d$ , representing a general form of a shallow feed-forward neural network, from a small number of query samples. Higher order differentiation, as used in our constructive approximations, of sums of ridge functions or of their compositions, as in deeper neural network, yields a natural connection between neural network weight identification and tensor product decomposition identification. In the case of the shallowest feed-forward neural network, second-order differentiation and tensors of order two (i.e., matrices) suffice as we prove in this paper. We use two sampling schemes to perform approximate differentiation—active sampling, where the sampling points are universal, actively and randomly designed, and passive sampling, where sampling points were preselected at random from a distribution with known density. Based on multiple gathered approximated first- and second-order differentials, our general approximation strategy is developed as a sequence of algorithms to perform individual sub-tasks. We first perform an active subspace search by approximating the span of the weight vectors $a_1,\dots ,a_m$ . Then we use a straightforward substitution, which reduces the dimensionality of the problem from $d$ to $m$ . The core of the construction is then the stable and efficient approximation of weights expressed in terms of rank- $1$ matrices $a_i \otimes a_i$ , realized by formulating their individual identification as a suitable nonlinear program. We prove the successful identification by this program of weight vectors being close to orthonormal and we also show how we can constructively reduce to this case by a whitening procedure, without loss of any generality. We finally discuss the implementation and the performance of the proposed algorithmic pipeline with extensive numerical experiments, which illustrate and confirm the theoretical results.

X-ray images are widely used in the study of paintings. When a painting has hidden sub-surface features (e.g., reuse of the canvas or revision of a composition by the artist), the resulting X-ray images can be hard to interpret as they include contributions from both the surface painting and the hidden design. In this paper we propose a self-supervised deep learning-based image separation approach that can be applied to the X-ray images from such paintings (‘mixed X-ray images’) to separate them into two hypothetical X-ray images, one containing information related to the visible painting only and the other containing the hidden features. The proposed approach involves two steps: (1) separation of the mixed X-ray image into two images, guided by the combined use of a reconstruction and an exclusion loss; (2) even allocation of the error map into the two individual, separated X-ray images, yielding separation results that have an appearance that is more familiar in relation to Xray images. The proposed method was demonstrated on a real painting with hidden content, Doña Isabel de Porcel by Francisco de Goya, to show its effectiveness.

To help understand the underlying mechanisms of neural networks (NNs), several groups have studied the number of linear regions ℓ of piecewise linear (PwL) functions, generated by deep neural networks (DNN). In particular, they showed that ℓ can grow exponentially with the number of network parameters p, a property often used to explain the advantages of deep over shallow NNs. Nonetheless, a dimension argument shows that DNNs cannot generate all PwL functions with ℓ linear regions when ℓ > p.It is thus natural to seek to characterize specific families of functions with ℓ > p linear regions that can be constructed by DNNs. Iterated Function Systems (IFS) recursively construct a sequence of PwL functions Fk with a number of linear regions which is exponential in k. We show that Fk can be generated by a NN using only O(k) parameters. IFS are used extensively to generate natural-looking landscape textures in artificial images as well as for compression of natural images. The surprisingly good performance of this compression suggests that human visual system may lock in on self-similarities. The combination of this phenomenon with the capacity of DNNs to efficiently approximate IFS may contribute to the success of DNNs, particularly striking for image processing tasks.

We formulate a novel characterization of a family of invertible maps between two-dimensional domains. Our work follows two classic results: The Rad\'o-Kneser-Choquet (RKC) theorem, which establishes the invertibility of harmonic maps into a convex planer domain; and Tutte's embedding theorem for planar graphs - RKC's discrete counterpart - which proves the invertibility of piecewise linear maps of triangulated domains satisfying a discrete-harmonic principle, into a convex planar polygon. In both theorems, the convexity of the target domain is essential for ensuring invertibility. We extend these characterizations, in both the continuous and discrete cases, by replacing convexity with a less restrictive condition. In the continuous case, Alessandrini and Nesi provide a characterization of invertible harmonic maps into non-convex domains with a smooth boundary by adding additional conditions on orientation preservation along the boundary. We extend their results by defining a condition on the normal derivatives along the boundary, which we call the cone condition; this condition is tractable and geometrically intuitive, encoding a weak notion of local invertibility. The cone condition enables us to extend Alessandrini and Nesi to the case of harmonic maps into non-convex domains with a piecewise-smooth boundary. In the discrete case, we use an analog of the cone condition to characterize invertible discrete-harmonic piecewise-linear maps of triangulations. This gives an analog of our continuous results and characterizes invertible discrete-harmonic maps in terms of the orientation of triangles incident on the boundary.

X-radiography is a widely used imaging technique in art investigation, whether to investigate the condition of a painting or provide insights into artists' techniques and working methods. In this paper, we propose a new architecture based on the use of `connected' auto-encoders in order to separate mixed X-ray images acquired from double-sided paintings, where in addition to the mixed X-ray image one can also exploit the two RGB images associated with the front and back of the painting. This proposed architecture uses convolutional auto-encoders that extract features from the RGB images that can be employed to (1) reproduce both of the original RGB images, (2) reconstruct the associated separated X-ray images, and (3) regenerate the mixed X-ray image. It operates in a totally self-supervised fashion without the need for examples containing both the mixed X-ray images and the separated ones. Based on images from the double-sided wing panels from the famous Ghent Altarpiece, painted in 1432 by the brothers Hubert and Jan Van Eyck, the proposed algorithm has been experimentally verified to outperform state-of-the-art X-ray separation methods in art investigation applications.

The approximation of both geodesic distances and shortest paths on point cloud sampled from an embedded submanifold ℳ of Euclidean space has been a long-standing challenge in computational geometry. Given a sampling resolution parameter h, state-of-the-art discrete methods yield O(h) provable approximations. In this paper, we investigate the convergence of such approximations made by Manifold Moving Least-Squares (Manifold-MLS), a method that constructs an approximating manifold ℳ^h using information from a given point cloud that was developed by Sober & Levin in 2019. In this paper, we show that provided that ℳ∈ C^k and closed (i.e. ℳ is a compact manifold without boundary) the Riemannian metric of ℳ^h approximates the Riemannian metric of ℳ,. Explicitly, given points p_1, p_2 ∈ℳ with geodesic distance ρ_ℳ(p_1, p_2), we show that their corresponding points p_1^h, p_2^h ∈ℳ^h have a geodesic distance of ρ_ℳ^h(p_1^h,p_2^h) = ρ_ℳ(p_1, p_2)(1 + O(h^k-1)) (i.e., the Manifold-MLS is nearly an isometry). We then use this result, as well as the fact that ℳ^h can be sampled with any desired resolution, to devise a naive algorithm that yields approximate geodesic distances with a rate of convergence O(h^k-1). We show the potential and the robustness to noise of the proposed method on some numerical simulations.

Massive data sets have their own architecture. Each data source has an inherent structure, which we should attempt to detect in order to utilize it for applications, such as denoising, clustering, anomaly detection, knowledge extraction, or classification. Harmonic analysis revolves around creating new structures for decomposition, rearrangement and reconstruction of operators and functions—in other words inventing and exploring new architectures for information and inference. Two previous very successful workshops on applied harmonic analysis and sparse approximation have taken place in 2012 and in 2015. This workshop was the an evolution and continuation of these workshops and intended to bring together world leading experts in applied harmonic analysis, data analysis, optimization, statistics, and machine learning to report on recent developments, and to foster new developments and collaborations.

Shape characterizers are metrics that quantify aspects of the overall geometry of a three‐dimensional (3D) digital surface. When computed for biological objects, the values of a shape characterizer are largely independent of homology interpretations and often contain a strong ecological and functional signal. Thus, shape characterizers are useful for understanding evolutionary processes. Dirichlet normal energy (DNE) is a widely used shape characterizer in morphological studies. Recent studies found that DNE is sensitive to various procedures for preparing 3D mesh from raw scan data, raising concerns regarding comparability and objectivity when utilizing DNE in morphological research. We provide a r obustly i mplemented a lgorithm for computing the Dirichlet energy of the normal (ariaDNE) on 3D meshes. We show through simulation that the effects of preparation‐related mesh surface attributes, such as triangle count, mesh representation, noise, smoothing and boundary triangles, are much more limited on ariaDNE than DNE. Furthermore, ariaDNE retains the potential of DNE for biological studies, illustrated by its effectiveness in differentiating species by dietary preferences. Use of ariaDNE can dramatically enhance the assessment of the ecological aspects of morphological variation by its stability under different 3D model acquisition methods and preparation procedure. Towards this goal, we provide scripts for computing ariaDNE and ariaDNE values for specimens used in previously published DNE analyses.

The problem of phase retrieval is to determine a signal $$f\in \mathcal {H}$$, with $$ \mathcal {H}$$ a Hilbert space, from intensity measurements $$|F(\omega )|$$, where $$F(\omega ):=\langle f, \varphi _\omega \rangle $$ are measurements of f with respect to a measurement system $$(\varphi _\omega )_{\omega \in \Omega }\subset \mathcal {H}$$. Although phase retrieval is always stable in the finite-dimensional setting whenever it is possible (i.e. injectivity implies stability for the inverse problem), the situation is drastically different if $$\mathcal {H}$$ is infinite-dimensional: in that case phase retrieval is never uniformly stable (Alaifari and Grohs in SIAM J Math Anal 49(3):1895–1911, 2017 ; Cahill et al. in Trans Am Math Soc Ser B 3(3):63–76, 2016 ); moreover, the stability deteriorates severely in the dimension of the problem (Cahill et al. 2016 ). On the other hand, all empirically observed instabilities are of a certain type: they occur whenever the function | F | of intensity measurements is concentrated on disjoint sets $$D_j\subset \Omega $$, i.e. when $$F= \sum _{j=1}^k F_j$$ where each $$F_j$$ is concentrated on $$D_j$$ (and $$k \ge 2$$). Motivated by these considerations, we propose a new paradigm for stable phase retrieval by considering the problem of reconstructing F up to a phase factor that is not global, but that can be different for each of the subsets $$D_j$$, i.e. recovering F up to the equivalence $$\begin{aligned} F \sim \sum _{j=1}^k e^{\mathrm {i}\alpha _j} F_j. \end{aligned}$$We present concrete applications (for example in audio processing) where this new notion of stability is natural and meaningful and show that in this setting stable phase retrieval can actually be achieved, for instance, if the measurement system is a Gabor frame or a frame of Cauchy wavelets.

As a means of improving analysis of biological shapes, we propose an algorithm for sampling a Riemannian manifold by sequentially selecting points with maximum uncertainty under a Gaussian process model. This greedy strategy is known to be near-optimal in the experimental design literature, and it appears to outperform the use of user-placed landmarks in representing the geometry of biological objects in our application. In the noiseless regime, we establish an upper bound for the mean squared prediction error (MSPE) in terms of the number of samples and geometric quantities of the manifold, demonstrating that the MSPE for our proposed sequential design decays at a rate comparable to the oracle rate achievable by any sequential or nonsequential optimal design; to the best of our knowledge this is the first result of this type for sequential experimental design. The key is to link the greedy algorithm to reduced basis methods in the context of model reduction for partial differential equations (PDEs). We expect this approach will find additional applications in other fields of research.