Mark Iwen

Mark Iwen, Aditya Viswanathan, Yang Wang

We develop a fast phase retrieval method which can utilize a large class of local phaseless correlation-based measurements in order to recover a given signal ${\bf x} \in \mathbb{C}^d$ (up to an unknown global phase) in near-linear $\mathcal{O} \left( d \log^4 d \right)$-time. Accompanying theoretical analysis proves that the proposed algorithm is guaranteed to deterministically recover all signals ${\bf x}$ satisfying a natural flatness (i.e., non-sparsity) condition for a particular choice of deterministic correlation-based measurements. A randomized version of these same measurements is then shown to provide nonuniform probabilistic recovery guarantees for arbitrary signals ${\bf x} \in \mathbb{C}^d$. Numerical experiments demonstrate the method's speed, accuracy, and robustness in practice -- all code is made publicly available. Finally, we conclude by developing an extension of the proposed method to the sparse phase retrieval problem; specifically, we demonstrate a sublinear-time compressive phase retrieval algorithm which is guaranteed to recover a given $s$-sparse vector ${\bf x} \in \mathbb{C}^d$ with high probability in just $\mathcal{O}(s \log^5 s \cdot \log d)$-time using only $\mathcal{O}(s \log^4 s \cdot \log d)$ magnitude measurements. In doing so we demonstrate the existence of compressive phase retrieval algorithms with near-optimal linear-in-sparsity runtime complexities.

Mark A. Iwen, Brian Preskitt, Rayan Saab et al.

We improve a phase retrieval approach that uses correlation-based measurements with compactly supported measurement masks [27]. The improved algorithm admits deterministic measurement constructions together with a robust, fast recovery algorithm that consists of solving a system of linear equations in a lifted space, followed by finding an eigenvector (e.g., via an inverse power iteration). Theoretical reconstruction error guarantees from [27] are improved as a result for the new and more robust reconstruction approach proposed herein. Numerical experiments demonstrate robustness and computational efficiency that outperforms competing approaches on large problems. Finally, we show that this approach also trivially extends to phase retrieval problems based on windowed Fourier measurements.

Mark Iwen, Fadil Santosa, Rachel Ward

We investigate the problem of decoding a bar code from a signal measured with a hand-held laser-based scanner. Rather than formulating the inverse problem as one of binary image reconstruction, we instead incorporate the symbology of the bar code into the reconstruction algorithm directly, and search for a sparse representation of the UPC bar code with respect to this known dictionary. Our approach significantly reduces the degrees of freedom in the problem, allowing for accurate reconstruction that is robust to noise and unknown parameters in the scanning device. We propose a greedy reconstruction algorithm and provide robust reconstruction guarantees. Numerical examples illustrate the insensitivity of our symbology-based reconstruction to both imprecise model parameters and noise on the scanned measurements.

M. A. Iwen

A compressed sensing method consists of a rectangular measurement matrix, $M \in \mathbbm{R}^{m \times N}$ with $m \ll N$, together with an associated recovery algorithm, $\mathcal{A}: \mathbbm{R}^m \rightarrow \mathbbm{R}^N$. Compressed sensing methods aim to construct a high quality approximation to any given input vector ${\bf x} \in \mathbbm{R}^N$ using only $M {\bf x} \in \mathbbm{R}^m$ as input. In particular, we focus herein on instance optimal nonlinear approximation error bounds for $M$ and $\mathcal{A}$ of the form $ \| {\bf x} - \mathcal{A} (M {\bf x}) \|_p \leq \| {\bf x} - {\bf x}^{\rm opt}_k \|_p + C k^{1/p - 1/q} \| {\bf x} - {\bf x}^{\rm opt}_k \|_q$ for ${\bf x} \in \mathbbm{R}^N$, where ${\bf x}^{\rm opt}_k$ is the best possible $k$-term approximation to ${\bf x}$. In this paper we develop a compressed sensing method whose associated recovery algorithm, $\mathcal{A}$, runs in $O((k \log k) \log N)$-time, matching a lower bound up to a $O(\log k)$ factor. This runtime is obtained by using a new class of sparse binary compressed sensing matrices of near optimal size in combination with sublinear-time recovery techniques motivated by sketching algorithms for high-volume data streams. The new class of matrices is constructed by randomly subsampling rows from well-chosen incoherent matrix constructions which already have a sub-linear number of rows. As a consequence, fewer random bits than previously required are needed in order to select the rows utilized by the fast reconstruction algorithms considered herein.

Alp Ozdemir, Ali Zare, Mark A. Iwen et al.

The widespread use of multisensor technology and the emergence of big datasets have created the need to develop tools to reduce, approximate, and classify large and multimodal data such as higher-order tensors. While early approaches focused on matrix and vector based methods to represent these higher-order data, more recently it has been shown that tensor decomposition methods are better equipped to capture couplings across their different modes. For these reasons, tensor decomposition methods have found applications in many different signal processing problems including dimensionality reduction, signal separation, linear regression, feature extraction, and classification. However, most of the existing tensor decomposition methods are based on the principle of finding a low-rank approximation in a linear subspace structure, where the definition of the rank may change depending on the particular decomposition. Since many datasets are not necessarily low-rank in a linear subspace, this often results in high approximation errors or low compression rates. In this paper, we introduce a new adaptive, multi-scale tensor decomposition method for higher order data inspired by hybrid linear modeling and subspace clustering techniques. In particular, we develop a multi-scale higher-order singular value decomposition (MS-HoSVD) approach where a given tensor is first permuted and then partitioned into several sub-tensors each of which can be represented as a low-rank tensor with increased representational efficiency. The proposed approach is evaluated for dimensionality reduction and classification for several different real-life tensor signals with promising results.

Mark A. Iwen, Sami Merhi, Michael Perlmutter

In this short note, we consider the worst case noise robustness of any phase retrieval algorithm which aims to reconstruct all nonvanishing vectors $\mathbf{x} \in \mathbb{C}^d$ (up to a single global phase multiple) from the magnitudes of an arbitrary collection of local correlation measurements. Examples of such measurements include both spectrogram measurements of $\mathbf{x}$ using locally supported windows and masked Fourier transform intensity measurements of $\mathbf{x}$ using bandlimited masks. As a result, the robustness results considered herein apply to a wide range of both ptychographic and Fourier ptychographic imaging scenarios. In particular, the main results imply that the accurate recovery of high-resolution images of extremely large samples using highly localized probes is likely to require an extremely large number of measurements in order to be robust to worst case measurement noise, independent of the recovery algorithm employed. Furthermore, recent pushes to achieve high-speed and high-resolution ptychographic imaging of integrated circuits for process verification and failure analysis will likely need to carefully balance probe design (e.g., their effective time-frequency support) against the total number of measurements acquired in order for their imaging techniques to be stable to measurement noise, no matter what reconstruction algorithms are applied.

Jade Larriva-Latt, Angela Morrison, Alison Radgowski et al.

Certain applications such as Magnetic Resonance Imaging (MRI) require the reconstruction of functions from Fourier spectral data. When the underlying functions are piecewise-smooth, standard Fourier approximation methods suffer from the Gibbs phenomenon - with associated oscillatory artifacts in the vicinity of edges and an overall reduced order of convergence in the approximation. This paper proposes an edge-augmented Fourier reconstruction procedure which uses only the first few Fourier coefficients of an underlying piecewise-smooth function to accurately estimate jump information and then incorporate it into a Fourier partial sum approximation. We provide both theoretical and empirical results showing the improved accuracy of the proposed method, as well as comparisons demonstrating superior performance over existing state-of-the-art sparse optimization-based methods. Extensions of the proposed techniques to functions of several variables are also addressed preliminarily. All code used to generate the results in this report are made publicly available.

Santhosh Karnik, Rongrong Wang, Mark Iwen

The remarkable successes of neural networks in a huge variety of inverse problems have fueled their adoption in disciplines ranging from medical imaging to seismic analysis over the past decade. However, the high dimensionality of such inverse problems has simultaneously left current theory, which predicts that networks should scale exponentially in the dimension of the problem, unable to explain why the seemingly small networks used in these settings work as well as they do in practice. To reduce this gap between theory and practice, we provide a general method for bounding the complexity required for a neural network to approximate a Hölder (or uniformly) continuous function defined on a high-dimensional set with a low-complexity structure. The approach is based on the observation that the existence of a Johnson-Lindenstrauss embedding $A\in\mathbb{R}^{d\times D}$ of a given high-dimensional set $S\subset\mathbb{R}^D$ into a low dimensional cube $[-M,M]^d$ implies that for any Hölder (or uniformly) continuous function $f:S\to\mathbb{R}^p$, there exists a Hölder (or uniformly) continuous function $g:[-M,M]^d\to\mathbb{R}^p$ such that $g(Ax)=f(x)$ for all $x\in S$. Hence, if one has a neural network which approximates $g:[-M,M]^d\to\mathbb{R}^p$, then a layer can be added that implements the JL embedding $A$ to obtain a neural network that approximates $f:S\to\mathbb{R}^p$. By pairing JL embedding results along with results on approximation of Hölder (or uniformly) continuous functions by neural networks, one then obtains results which bound the complexity required for a neural network to approximate Hölder (or uniformly) continuous functions on high dimensional sets. The end result is a general theoretical framework which can then be used to better explain the observed empirical successes of smaller networks in a wider variety of inverse problems than current theory allows.

Mark A. Iwen, Mark Philip Roach

Let $\mathcal{M}$ be a compact $d$-dimensional submanifold of $\mathbb{R}^N$ with reach $τ$ and volume $V_{\mathcal M}$. Fix $ε\in (0,1)$. In this paper we prove that a nonlinear function $f: \mathbb{R}^N \rightarrow \mathbb{R}^{m}$ exists with $m \leq C \left(d / ε^2 \right) \log \left(\frac{\sqrt[d]{V_{\mathcal M}}}τ \right)$ such that $$(1 - ε) \| {\bf x} - {\bf y} \|_2 \leq \left\| f({\bf x}) - f({\bf y}) \right\|_2 \leq (1 + ε) \| {\bf x} - {\bf y} \|_2$$ holds for all ${\bf x} \in \mathcal{M}$ and ${\bf y} \in \mathbb{R}^N$. In effect, $f$ not only serves as a bi-Lipschitz function from $\mathcal{M}$ into $\mathbb{R}^{m}$ with bi-Lipschitz constants close to one, but also approximately preserves all distances from points not in $\mathcal{M}$ to all points in $\mathcal{M}$ in its image. Furthermore, the proof is constructive and yields an algorithm which works well in practice. In particular, it is empirically demonstrated herein that such nonlinear functions allow for more accurate compressive nearest neighbor classification than standard linear Johnson-Lindenstrauss embeddings do in practice.

Edem Boahen, Simone Brugiapaglia, Hung-Hsu Chou et al.

Motivated by applications such as sparse PCA, in this paper we present provably-accurate one-pass algorithms for the sparse approximation of the top eigenvectors of extremely massive matrices based on a single compact linear sketch. The resulting compressive-sensing-based approaches can approximate the leading eigenvectors of huge approximately low-rank matrices that are too large to store in memory based on a single pass over its entries while utilizing a total memory footprint on the order of the much smaller desired sparse eigenvector approximations. Finally, the compressive sensing recovery algorithm itself (which takes the gathered compressive matrix measurements as input, and then outputs sparse approximations of its top eigenvectors) can also be formulated to run in a time which principally depends on the size of the sought sparse approximations, making its runtime sublinear in the size of the large matrix whose eigenvectors one aims to approximate. Preliminary experiments on huge matrices having $\sim 10^{16}$ entries illustrate the developed theory and demonstrate the practical potential of the proposed approach.

M. A. Iwen

We study the problem of estimating the best B term Fourier representation for a given frequency-sparse signal (i.e., vector) $\textbf{A}$ of length $N \gg B$. More explicitly, we investigate how to deterministically identify B of the largest magnitude frequencies of $\hat{\textbf{A}}$, and estimate their coefficients, in polynomial$(B,\log N)$ time. Randomized sub-linear time algorithms which have a small (controllable) probability of failure for each processed signal exist for solving this problem. However, for failure intolerant applications such as those involving mission-critical hardware designed to process many signals over a long lifetime, deterministic algorithms with no probability of failure are highly desirable. In this paper we build on the deterministic Compressed Sensing results of Cormode and Muthukrishnan (CM) \cite{CMDetCS3,CMDetCS1,CMDetCS2} in order to develop the first known deterministic sub-linear time sparse Fourier Transform algorithm suitable for failure intolerant applications. Furthermore, in the process of developing our new Fourier algorithm, we present a simplified deterministic Compressed Sensing algorithm which improves on CM's algebraic compressibility results while simultaneously maintaining their results concerning exponential decay.

Santhosh Karnik, Anna Veselovska, Mark Iwen et al.

We provide a rigorous analysis of implicit regularization in an overparametrized tensor factorization problem beyond the lazy training regime. For matrix factorization problems, this phenomenon has been studied in a number of works. A particular challenge has been to design universal initialization strategies which provably lead to implicit regularization in gradient-descent methods. At the same time, it has been argued by Cohen et. al. 2016 that more general classes of neural networks can be captured by considering tensor factorizations. However, in the tensor case, implicit regularization has only been rigorously established for gradient flow or in the lazy training regime. In this paper, we prove the first tensor result of its kind for gradient descent rather than gradient flow. We focus on the tubal tensor product and the associated notion of low tubal rank, encouraged by the relevance of this model for image data. We establish that gradient descent in an overparametrized tensor factorization model with a small random initialization exhibits an implicit bias towards solutions of low tubal rank. Our theoretical findings are illustrated in an extensive set of numerical simulations show-casing the dynamics predicted by our theory as well as the crucial role of using a small random initialization.

Michael Perlmutter, Jieqian He, Mark Iwen et al.

The scattering transform is a wavelet-based model of Convolutional Neural Networks originally introduced by S. Mallat. Mallat's analysis shows that this network has desirable stability and invariance guarantees and therefore helps explain the observation that the filters learned by early layers of a Convolutional Neural Network typically resemble wavelets. Our aim is to understand what sort of filters should be used in the later layers of the network. Towards this end, we propose a two-layer hybrid scattering transform. In our first layer, we convolve the input signal with a wavelet filter transform to promote sparsity, and, in the second layer, we convolve with a Gabor filter to leverage the sparsity created by the first layer. We show that these measurements characterize information about signals with isolated singularities. We also show that the Gabor measurements used in the second layer can be used to synthesize sparse signals such as those produced by the first layer.

Mark A. Iwen, Benjamin Schmidt, Arman Tavakoli

Let $\mathcal{M}$ be a smooth $d$-dimensional submanifold of $\mathbb{R}^N$ with boundary that's equipped with the Euclidean (chordal) metric, and choose $m \leq N$. In this paper we consider the probability that a random matrix $A \in \mathbb{R}^{m \times N}$ will serve as a bi-Lipschitz function $A: \mathcal{M} \rightarrow \mathbb{R}^m$ with bi-Lipschitz constants close to one for three different types of distributions on the $m \times N$ matrices $A$, including two whose realizations are guaranteed to have fast matrix-vector multiplies. In doing so we generalize prior randomized metric space embedding results of this type for submanifolds of $\mathbb{R}^N$ by allowing for the presence of boundary while also retaining, and in some cases improving, prior lower bounds on the achievable embedding dimensions $m$ for which one can expect small distortion with high probability. In particular, motivated by recent modewise embedding constructions for tensor data, herein we present a new class of highly structured distributions on matrices which outperform prior structured matrix distributions for embedding sufficiently low-dimensional submanifolds of $\mathbb{R}^N$ (with $d \lesssim \sqrt{N}$) with respect to both achievable embedding dimension, and computationally efficient realizations. As a consequence we are able to present, for example, a general new class of Johnson-Lindenstrauss embedding matrices for $\mathcal{O}(\log^c N)$-dimensional submanifolds of $\mathbb{R}^N$ which enjoy $\mathcal{O}(N \log (\log N))$-time matrix vector multiplications.

M. A. Iwen, D. Needell, E. Rebrova et al.

In this paper new general modewise Johnson-Lindenstrauss (JL) subspace embeddings are proposed that are both considerably faster to generate and easier to store than traditional JL embeddings when working with extremely large vectors and/or tensors. Corresponding embedding results are then proven for two different types of low-dimensional (tensor) subspaces. The first of these new subspace embedding results produces improved space complexity bounds for embeddings of rank-$r$ tensors whose CP decompositions are contained in the span of a fixed (but unknown) set of $r$ rank-one basis tensors. In the traditional vector setting this first result yields new and very general near-optimal oblivious subspace embedding constructions that require fewer random bits to generate than standard JL embeddings when embedding subspaces of $\mathbb{C}^N$ spanned by basis vectors with special Kronecker structure. The second result proven herein provides new fast JL embeddings of arbitrary $r$-dimensional subspaces $\mathcal{S} \subset \mathbb{C}^N$ which also require fewer random bits (and so are easier to store - i.e., require less space) than standard fast JL embedding methods in order to achieve small $ε$-distortions. These new oblivious subspace embedding results work by $(i)$ effectively folding any given vector in $\mathcal{S}$ into a (not necessarily low-rank) tensor, and then $(ii)$ embedding the resulting tensor into $\mathbb{C}^m$ for $m \leq C r \log^c(N) / ε^2$. Applications related to compression and fast compressed least squares solution methods are also considered, including those used for fitting low-rank CP decompositions, and the proposed JL embedding results are shown to work well numerically in both settings.

Nada Sissouno, Florian Boßmann, Frank Filbir et al.

Measurements achieved with ptychographic imaging are a special case of diffraction measurements. They are generated by illuminating small parts of a sample with, e.g., a focused X-ray beam. By shifting the sample, a set of far-field diffraction patterns of the whole sample are then obtained. From a mathematical point of view those measurements are the squared modulus of the windowed Fourier transform of the sample. Thus, we have a phase retrieval problem for local Fourier measurements. A direct solver for this problem was introduced by Iwen, Viswanathan and Wang in 2016 and improved by Iwen, Preskitt, Saab and Viswanathan in 2018. Motivated by the applied perspective of ptychographic imaging, we present a generalization of this method and compare the different versions in numerical experiments. The new method proposed herein turns out to be more stable, particularly in the case of missing data.

Sami Merhi, Ruochuan Zhang, Mark A. Iwen et al.

In this paper we consider Sparse Fourier Transform (SFT) algorithms for approximately computing the best $s$-term approximation of the Discrete Fourier Transform (DFT) $\mathbf{\hat{f}} \in \mathbb{C}^N$ of any given input vector $\mathbf{f} \in \mathbb{C}^N$ in just $\left( s \log N\right)^{\mathcal{O}(1)}$-time using only a similarly small number of entries of $\mathbf{f}$. In particular, we present a deterministic SFT algorithm which is guaranteed to always recover a near best $s$-term approximation of the DFT of any given input vector $\mathbf{f} \in \mathbb{C}^N$ in $\mathcal{O} \left( s^2 \log ^{\frac{11}{2}} (N) \right)$-time. Unlike previous deterministic results of this kind, our deterministic result holds for both arbitrary vectors $\mathbf{f} \in \mathbb{C}^N$ and vector lengths $N$. In addition to these deterministic SFT results, we also develop several new publicly available randomized SFT implementations for approximately computing $\mathbf{\hat{f}}$ from $\mathbf{f}$ using the same general techniques. The best of these new implementations is shown to outperform existing discrete sparse Fourier transform methods with respect to both runtime and noise robustness for large vector lengths $N$.

Sami Merhi, Aditya Viswanathan, Mark Iwen

A novel phase retrieval method, motivated by ptychographic imaging, is proposed for the approximate recovery of a compactly supported specimen function $f:\mathbb{R}\rightarrow\mathbb{C}$ from its continuous short time Fourier transform (STFT) spectrogram measurements. The method, partially inspired by the well known PhaseLift algorithm, is based on a lifted formulation of the infinite dimensional problem which is then later truncated for the sake of computation. Numerical experiments demonstrate the promise of the proposed approach.

Mark Iwen, Aditya Viswanathan, Yang Wang

In this short note we propose a simple two-stage sparse phase retrieval strategy that uses a near-optimal number of measurements, and is both computationally efficient and robust to measurement noise. In addition, the proposed strategy is fairly general, allowing for a large number of new measurement constructions and recovery algorithms to be designed with minimal effort.