What happens when generative machine learning models are pretrained on web-scale datasets containing data generated by earlier models? Some prior work warns of “model collapse” as the web is overwhelmed by synthetic data; other work suggests the problem can be contained (i.e. collapse can be avoided) by managing how available data are used in pretraining. In this paper, we report experiments on three ways of using data (training-workflows), across three generative model task-settings (multivariate Gaussian estimation, kernel density estimation, and language-model fine-tuning) to further confirm the possibility of containment: (a) we confirm that the training-workflow of {\it replacing} all real data by successive generations of purely synthetic data indeed suffers model collapse in all task-settings studied; (b) we consider the training-workflow of {\it accumulating} synthetic data alongside real data and training on all data combined and confirming that, although the proportion of real data eventually becomes zero, models remain stable and their test losses do not diverge under this training-workflow; (c) we consider a training-workflow where real and synthetic data accumulate together but successive generations of pretraining are constrained to use fixed-size data subsets each generation. In this workflow, we observe slow and gradual rather than explosive degradation of test loss performance across generations. Our insights are particularly important when forecasting whether future frontier generative models will collapse or thrive, and our results open avenues for empirically and mathematically studying the context-dependent value of synthetic data.

Science progresses by iteratively advancing and correcting humanity's understanding of the world. In machine learning (ML) research, rapid advancements have led to an explosion of publications, but have also led to misleading, incorrect, flawed or perhaps even fraudulent studies being accepted and sometimes highlighted at ML conferences due to the fallibility of peer review. While such mistakes are understandable, ML conferences do not offer robust processes to help the field systematically correct when such errors are made.This position paper argues that ML conferences should establish a dedicated "Refutations and Critiques" (R&C) Track. This R&C Track would provide a high-profile, reputable platform to support vital research that critically challenges prior research, thereby fostering a dynamic self-correcting research ecosystem.We discuss key considerations including track design, review principles, potential pitfalls, and provide an illustrative example submission concerning a recent ICLR 2025 Oral.We conclude that ML conferences should create official, reputable mechanisms to help ML research self-correct.

Many applications seek to recover low-rank approximations of noisy tensor data. We consider several practical and effective matricization strategies which construct specific matrices from such tensors and then apply spectral methods; the strategies include tensor unfolding, partial tracing, power iteration, and recursive unfolding. We settle the behaviors of unfolding and partial tracing, identifying sharp thresholds in signal-to-noise ratio above which the signal is partially recovered. In particular, we extend previous results to a much larger class of tensor shapes where axis lengths may be different. For power iteration and recursive unfolding, we prove that under conditions where previous algorithms partially recovery the signal, these methods achieve (asymptotically) exact recovery. Our analysis deploys random matrix theory to obtain sharp thresholds which elude perturbation and concentration bounds. Specifically, we rely upon recent disproportionate random matrix results, which describe sequences of matrices with diverging aspect ratio.

We propose a method to perform two-sample tests on discrete distributions, i.e., an unsupervised discriminator testing whether two frequency tables were each sampled from a single unspecified probability distribution. Our proposal takes feature-by-feature P-values based on a binomial allocation model, and combines the P-values using Higher Criticism. Performance on real-world data (e.g. authorship attribution challenges) shows this to be an effective unsupervised untrained discriminator even in violations of the binomial allocation model. The method has interesting theoretical properties, in the ‘rare/weak departures’ setting where, if two distributions are actually different, they differ only in relatively few features and only by relatively subtle amounts. We perform a phase diagram analysis in which the phase space quantifies how rare and how weak such departures are. Although our proposal does not require any formal specification of an alternative hypothesis, nor does it require any specification of a baseline or null hypothesis, in the limit where word counts are high, the method delivers the optimal phase diagram in the rare/weak setting: it is asymptotically fully powerful inside the region of phase space where a formally specified test would have been fully powerful. In the limit where counts are low, we derive the phase diagram as well, although the optimality of the resulting diagram is not discussed here.

A matrix network is a family of matrices, with their relations modeled as a weighted graph. We consider the task of completing a partially observed matrix network. The observation comes from a novel sampling scheme where a fraction of matrices might be completely unobserved. How can we recover the entire matrix network from incomplete observations? This mathematical problem arises in many applications including medical imaging and social networks. To recover the matrix network, we propose a structural assumption that the matrices are low-rank after the graph Fourier transform on the network. We formulate a convex optimization problem and prove an exact recovery guarantee for the optimization problem. Furthermore, we numerically characterize the exact recovery regime for varying rank and sampling rate and discover a new phase transition phenomenon. Then we give an iterative imputation algorithm to efficiently solve optimization problem and complete large scale matrix networks. We demonstrate the algorithm with a variety of applications such as MRI and Facebook user network.

Recovering high-resolution images from limited sensory data typically leads to a serious ill-posed inverse problem, demanding inversion algorithms that effectively capture the prior information. Learning a good inverse mapping from training data faces severe challenges, including: (i) scarcity of training data; (ii) need for plausible reconstructions that are physically feasible; (iii) need for fast reconstruction, especially in real-time applications. We develop a successful system solving all these challenges, using as basic architecture the repetitive application of alternating proximal and data fidelity constraints. We learn a proximal map that works well with real images based on residual networks with recurrent blocks. Extensive experiments are carried out under different settings: (a) reconstructing abdominal MRI of pediatric patients from highly undersampled k-space data and (b) super-resolving natural face images. Our key findings include: 1. a recurrent ResNet with a single residual block (10-fold repetition) yields an effective proximal which accurately reveals MR image details. 2. Our architecture significantly outperforms conventional non-recurrent deep ResNets by 2dB SNR; it is also trained much more rapidly. 3. It outperforms state-of-the-art compressed-sensing Wavelet-based methods by 4dB SNR, with 100x speedups in reconstruction time.

We consider the recovery of low-rank matrices from noisy data by shrinkage of singular values, in which a single, univariate nonlinearity is applied to each of the empirical singular values. We adopt an asymptotic framework, in which the matrix size is much larger than the rank of the signal matrix to be recovered, and the signal-to-noise ratio of the low-rank piece stays constant. For a variety of loss functions, including Mean Square Error (MSE) - (square Frobenius norm), the nuclear norm loss and the operator norm loss, we show that in this framework there is a well-defined asymptotic loss that we evaluate precisely in each case. In fact, each of the loss functions we study admits a unique admissible shrinkage nonlinearity dominating all other nonlinearities. We provide a general method for evaluating these optimal nonlinearities, and demonstrate our framework by working out simple, explicit formulas for the optimal nonlinearities in the Frobenius, nuclear and operator norm cases. For example, for a square low-rank n-by-n matrix observed in white noise with level σ, the optimal nonlinearity for MSE loss simply shrinks each data singular value y to √(y2-4nσ2) (or to 0 if y2, which compares favorably with optimally tuned hard thresholding and optimally tuned soft thresholding, providing guarantees of 3nrσ2and 6nrσ2, respectively. Our general method also allows one to evaluate optimal shrinkers numerically to arbitrary precision. As an example, we compute optimal shrinkers for the Schatten-p norm loss, for any p > 0.

The increasing availability of access to large-scale computing clusters - for example via the cloud - is changing the way scientific research can be conducted, enabling experiments of a scale and scope that would have been inconceivable several years ago. An ambitious data scientist today can carry out projects involving several million CPU hours. In the near future, we anticipate a typical Ph.D. in computational science may be expected or even required to offer findings based on at least 1 million CPU hours of computations. The massive scale of these soon-to-be-upon-us computational experiments demands that we change how we organize our experimental practices. Traditionally, and still the dominant paradigm today, the end-to-end process of experiment design and execution involves a significant amount of manual intervention and situational tweaking, cutting and pasting, and the use of disparate disconnected tools, much of which is undocumented and easily lost. This makes it difficult to detect and understand possible failure points in the computational workflow, making it virtually impossible to correct, let alone simply rerun the experiment. This is an amazing state of affairs, considering the ubiquity of error in scientific computation and in research generally. Following such unstructured and undocumented research practices limits the ability of the researcher to exploit cluster and cloud-based paradigms, as each increase in scale under the dominant paradigm is likely to lead to ever more errors and misunderstandings. A better paradigm will integrate the design of large experiments seamlessly with job management, output harvesting, data analysis, reporting, and publication of code and data. In particular such a paradigm would submerge the details of all the processing, harvesting, and management while exposing transparently the description of the discovery process itself, including details such as the parameter space exploration. Reproducing any job would be a push-button affair, and creating a new experiment from a previous one might involve only changes of a line or two of code followed again by push-button execution and reporting. Even though such experiments would be operating at a much greater scale than today, under such a paradigm they would be easier to conduct, obtain a lower error rate, and offer a much greater opportunity for `outsiders' to understand the results. In this article, we discuss the challenges of massive computational experimentation and present a taxonomy of some of the desiderata which such paradigms should offer. We then present ClusterJob (CJ), an efficient computing environment that we and other researchers have used to conduct and share million-CPU-hour experiments in a painless and reproducible way.

We consider recovery of low-rank matrices from noisy data by hard thresholding of singular values, in which empirical singular values below a threshold λ are set to 0. We study the asymptotic mean squared error (AMSE) in a framework, where the matrix size is large compared with the rank of the matrix to be recovered, and the signal-to-noise ratio of the low-rank piece stays constant. The AMSE-optimal choice of hard threshold, in the case of n-by-n matrix in white noise of level σ, is simply (4/√3)√nσ ≈ 2.309√nσ when σ is known, or simply 2.858 · ymedwhen σ is unknown, where ymedis the median empirical singular value. For nonsquare, m by n matrices with m ≠ n the thresholding coefficients 4/√3 and 2.858 are replaced with different provided constants that depend on m/n. Asymptotically, this thresholding rule adapts to unknown rank and unknown noise level in an optimal manner: it is always better than hard thresholding at any other value, and is always better than ideal truncated singular value decomposition (TSVD), which truncates at the true rank of the low-rank matrix we are trying to recover. Hard thresholding at the recommended value to recover an n-by-n matrix of rank r guarantees an AMSE at most 3 nrσ2. In comparison, the guarantees provided by TSVD, optimally tuned singular value soft thresholding and the best guarantee achievable by any shrinkage of the data singular values are 5 nrσ2, 6 nrσ2, and 2 nrσ2, respectively. The recommended value for hard threshold also offers, among hard thresholds, the best possible AMSE guarantees for recovering matrices with bounded nuclear norm. Empirical evidence suggests that performance improvement over TSVD and other popular shrinkage rules can be substantial, for different noise distributions, even in relatively small n.

Recent work on Approximate Message Passing algorithms in compressed sensing focuses on `ideal' algorithms which at each iteration face a subproblem of recovering an unknown sparse signal in Gaussian white noise. The noise level in each subproblem changes from iteration to iteration in a way that depends on the underlying signal (which we don't know!). For such algorithms to be used in practice, it seems we need an estimator that achieves the MMSE when the noise level is unknown. In this paper we solve this problem using convex optimization, Stein Unbiased Risk Estimates and Huber Splines.

Consider the compressed sensing problem of estimating an unknown k-sparse n-vector from a set of m noisy linear equations. Recent work focused on the noise sensitivity of particular algorithms - the scaling of the reconstruction error with added noise. In this paper, we study the minimax noise sensitivity - the minimum is over all possible recovery algorithms and the maximum is over all vectors obeying a sparsity constraint. This fundamental quantity characterizes the difficulty of recovery when nothing is known about the vector other than the fact that it has at most k nonzero entries. Assuming random sensing matrices (i.i.d. Gaussian), we obtain non-asymptotic bounds which show that the minimax noise sensitivity is finite if m ≥ k + 3 and infinite if m ≤ k + 1. We also study the large system behavior where δ = m/n ∈ (0,1) denotes the undersampling fraction and k/n = ε ∈ (0,1) denotes the sparsity fraction. There is a phase transition separating successful and unsuccessful recovery: the minimax noise sensitivity is bounded for any δ > ε and is unbounded for any δ <; ε. One consequence of our results is that the Bayes optimal phase transitions of Wu and Verdu can be obtained uniformly over the class of all sparse vectors.

Compressed sensing posits that, within limits, one can undersample a sparse signal and yet reconstruct it accurately. Knowing the precise limits to such undersampling is important both for theory and practice. We present a formula that characterizes the allowed undersampling of generalized sparse objects. The formula applies to approximate message passing (AMP) algorithms for compressed sensing, which are here generalized to employ denoising operators besides the traditional scalar soft thresholding denoiser. This paper gives several examples including scalar denoisers not derived from convex penalization-the firm shrinkage nonlinearity and the minimax nonlinearity-and also nonscalar denoisers-block thresholding, monotone regression, and total variation minimization. Let the variables ε = k/N and δ = n/N denote the generalized sparsity and undersampling fractions for sampling the k-generalized-sparse N-vector x0according to y=Ax0. Here, A is an n×N measurement matrix whose entries are iid standard Gaussian. The formula states that the phase transition curve δ = δ(ε) separating successful from unsuccessful reconstruction of x0by AMP is given by δ = M(ε|Denoiser) where M(ε|Denoiser) denotes the per-coordinate minimax mean squared error (MSE) of the specified, optimally tuned denoiser in the directly observed problem y = x + z. In short, the phase transition of a noiseless undersampling problem is identical to the minimax MSE in a denoising problem. We prove that this formula follows from state evolution and present numerical results validating it in a wide range of settings. The above formula generates numerous new insights, both in the scalar and in the nonscalar cases.

We study the compressed sensing reconstruction problem for a broad class of random, band-diagonal sensing matrices. This construction is inspired by the idea of spatial coupling in coding theory. As demonstrated heuristically and numerically by Krzakala [30], message passing algorithms can effectively solve the reconstruction problem for spatially coupled measurements with undersampling rates close to the fraction of nonzero coordinates. We use an approximate message passing (AMP) algorithm and analyze it through the state evolution method. We give a rigorous proof that this approach is successful as soon as the undersampling rate δ exceeds the (upper) Rényi information dimension of the signal, d̅(pX). More precisely, for a sequence of signals of diverging dimension n whose empirical distribution converges to pX, reconstruction is with high probability successful from d̅(pX) n+o(n) measurements taken according to a band diagonal matrix. For sparse signals, i.e., sequences of dimension n and k(n) nonzero entries, this implies reconstruction from k(n)+o(n) measurements. For “discrete” signals, i.e., signals whose coordinates take a fixed finite set of values, this implies reconstruction from o(n) measurements. The result is robust with respect to noise, does not apply uniquely to random signals, but requires the knowledge of the empirical distribution of the signal pX.

We study the compressed sensing reconstruction problem for a broad class of random, band-diagonal sensing matrices. This construction is inspired by the idea of spatial coupling in coding theory. As demonstrated heuristically and numerically by Krzakala et al. [11], message passing algorithms can effectively solve the reconstruction problem for spatially coupled measurements with undersampling rates close to the fraction of non-zero coordinates. We use an approximate message passing (AMP) algorithm and analyze it through the state evolution method. We give a rigorous proof that this approach is successful as soon as the undersampling rate δ exceeds the (upper) Rényi information dimension of the signal, d̅(pX). More precisely, for a sequence of signals of diverging dimension n whose empirical distribution converges to pX, reconstruction is with high probability successful from d̅(pX)n+o(n) measurements taken according to a band diagonal matrix. For sparse signals, i.e. sequences of dimension n and k(n) non-zero entries, this implies reconstruction from k(n) + o(n) measurements. For `discrete' signals, i.e. signals whose coordinates take a fixed finite set of values, this implies reconstruction from o(n) measurements. The result is robust with respect to noise, does not apply uniquely to random signals, but requires the knowledge of the empirical distribution of the signal pX.

Many papers studying compressed sensing consider the noisy underdetermined system of linear equations: y = Ax0+ z, with n × N measurement matrix A, n2I). Both y and A are known, both x0and z are unknown, and we seek an approximation to x0; we let δ = n/N ϵ (0,1) denote the undersampling fraction. In the popular strict sparsity model of compressed sensing, such papers further assume that x0has at most a specified fraction ε of nonzeros. In this paper, we relax the assumption of strict sparsity by assuming the vector x0is close in mean p-th power to a sparse signal. We study how this relaxation affects the performance of ℓ1-penalized ℓ2minimization, in which the reconstruction x1,λsolves min min ||y - Ax||22/2 + λ||x||1. We study asymptotic mean-squared error (AMSE), the large-system limit of the MSE of x1, λ. Using recently developed tools based on Approximate Message Passing (AMP), we develop expressions for minimax AMSE Mϵ,p*(δ, ξ, σ) - max over all approximately sparse signals, min over penalizations λ, where ξ measures the deviation from strict sparsity. There is of course a phase transition curve δ*= δ*(ε); only above this curve, δ >; δ*(ε), can we have exact recovery even in the noiseless-data strict-sparsity setting. It turns out that the minimax AMSE can be characterized succinctly by a coefficient sensp*(ε, δ) which we refer to as the sparsity-relaxation sensitivity. We give explicit expressions for sensp*(ε, δ), compute them, and interpret them. Our approach yields precise formulas in place of loose order bounds based on restricted isometry property and instance optimality results. Our formulas reveal that sensitivity is finite everywhere exact recovery is possible under strict sparsity, and that sensitivity to added random noise in the measurements y is smaller than the sensitivity to adding a comparable amount of noise to the estimand x0. Our methods can also treat the mean q-th power loss. The methods themselves are based on minimax decision theory and seem of independent interest.

Finding the sparsest solution to underdetermined systems of linear equationsy= Φxis NP-hard in general. We show here that for systems with “typical”/“random” Φ, a good approximation to the sparsest solution is obtained by applying a fixed number of standard operations from linear algebra. Our proposal, Stagewise Orthogonal Matching Pursuit (StOMP), successively transforms the signal into a negligible residual. Starting with initial residualr0=y, at thes-th stage it forms the “matched filter” ΦTrs-1, identifies all coordinates with amplitudes exceeding a specially chosen threshold, solves a least-squares problem using the selected coordinates, and subtracts the least-squares fit, producing a new residual. After a fixed number of stages (e.g., 10), it stops. In contrast to Orthogonal Matching Pursuit (OMP), many coefficients can enter the model at each stage in StOMP while only one enters per stage in OMP; and StOMP takes a fixed number of stages (e.g., 10), while OMP can take many (e.g.,n). We give both theoretical and empirical support for the large-system effectiveness of StOMP. We give numerical examples showing that StOMP rapidly and reliably finds sparse solutions in compressed sensing, decoding of error-correcting codes, and overcomplete representation.

We consider the compressed sensing problem where the object x0∈ ℝNis to be recovered from incomplete measurements y = Ax0+z. Here the sensing matrix A is an n×N random matrix with Gaussian entries and n1-penalized least-squares reconstruction (aka LASSO, Basis Pursuit). Suppose that xοis sparse in the sense of having ℓρnorm bounded by ξ · N1/ρfor some fixed 0; 0. In both the noisy (ziiid N(0, σ2)) and noiseless (z = 0) cases, we evaluate the worst-case asymptotic mean square error (AMSE) for optimally tuned ℓ1penalized least-squares, and we exhibit the least-favorable object xο(hardest sparse signal to recover) and the maximin penalization. Our explicit formulas yield precise relations. For example, in the noiseless case z = 0, for vectors xοof ℓρnorm bounded by 1, we show that mm max∥ x̑λ- x0∥2= n1-2/p(2log(N/n))2/p-1{1+oN(l)} where n, N → ∞, n/N → 0 slowly. The complete formulas applies to a general scaling limit n/N → δ and sparsity parameter ξ, and unexpectedly involve quantities from statistical decision theory. This reflects a deep connection between ℓ1-penalized ℓ2minimization and scalar soft thresholding.

Consider the noisy underdetermined system of linear equations: y = Ax0+ z, with A an n × N measurement matrix, n2I) a Gaussian white noise. Both y and A are known, both x0and z are unknown, and we seek an approximation to x0. When x0has few nonzeros, useful approximations are often obtained by ℓ1-penalized ℓ2minimization, in which the reconstruction x̂1,λsolves min{||y - Ax||22/2 + λ||x||1}. Consider the reconstruction mean-squared error MSE = E|| x̂1,λ- x0||22/N, and define the ratio MSE/σ2as the noise sensitivity. Consider matrices A with i.i.d. Gaussian entries and a large-system limit in which n, N → ∞ with n/N → δ and k/n → ρ. We develop exact expressions for the asymptotic MSE of x̂1,λ, and evaluate its worst-case noise sensitivity over all types of k-sparse signals. The phase space 0 ≤ 8, ρ ≤ 1 is partitioned by the curve ρ = ρMSE(δ) into two regions. Formal noise sensitivity is bounded throughout the region ρ = ρMSE(δ) and is unbounded throughout the region ρ = ρMSE(δ). The phase boundary ρ = ρMSE(δ) is identical to the previously known phase transition curve for equivalence of ℓ1- ℓ0minimization in the k-sparse noiseless case. Hence, a single phase boundary describes the fundamental phase transitions both for the noise less and noisy cases. Extensive computational experiments validate these predictions, including the existence of game-theoretical structures underlying it (saddlepoints in the payoff, least-favorable signals and maximin penalization). Underlying our formalism is an approximate message passing soft thresholding algorithm (AMP) introduced earlier by the authors. Other papers by the authors detail expressions for the formal MSE of AMP and its close connection to ℓ1-penalized reconstruction. The focus of the present paper is on computing the minimax formal MSE within the class of sparse signals x0.

This special issue presents many exciting and surprising developments in signal and image processing owing to sparse representations. While the technology is new, much of the intellectual heritage of these recent developments can be traced back generations. This foreword sketches the ¿prehistory of sparsity,¿ a fascinating range of early discoveries ranging from pure mathematics to oil exploration to biological vision.

Undersampling theorems state that we may gather far fewer samples than the usual sampling theorem while exactly reconstructing the object of interest-provided the object in question obeys a sparsity condition, the samples measure appropriate linear combinations of signal values, and we reconstruct with a particular nonlinear procedure. While there are many ways to crudely demonstrate such undersampling phenomena, we know of only one mathematically rigorous approach which precisely quantifies the true sparsity-undersampling tradeoff curve of standard algorithms and standard compressed sensing matrices. That approach, based on combinatorial geometry, predicts the exact location in sparsity-undersampling domain where standard algorithms exhibitphase transitionsin performance. We review the phase transition approach here and describe the broad range of cases where it applies. We also mention exceptions and state challenge problems for future research. Sample result: one can efficiently reconstruct a k-sparse signal of length N from n measurements, provided n ?? 2k ?? log(N/n), for (k,n,N) large.k ?? N.AMS 2000 subject classifications. Primary: 41A46, 52A22, 52B05, 62E20, 68P30, 94A20; Secondary: 15A52, 60F10, 68P25, 90C25, 94B20.

In“Counting faces of randomly projected polytopes when the projection radically lowers dimension”the authors proved an asymptoticsampling theorem for sparse signals, showing that$n$random measurements permit to reconstruct an$N$-vector having$k$nonzeros provided$$ n > 2 \cdot k \cdot \log(N/n) (1+o(1))$$reconstruction uses$\ell_1$minimization. They also proved anasymptotic rate theorem, showing existence of real error-correcting codes for messages of length$N$which can correct all possible$k$-element error patterns using just$n$generalized checksum bits, where$$ n > 2e\cdot k \log(N/n) (1 + o(1))$$decoding uses$\ell_1$minimization. Both results require an asymptotic framework, with$N$growing large. For applications, on the other hand, we are concerned with specific triples$k, n, N$. We exhibit triples$(k,n,N)$for which Compressed Sensing Matrices and Real Error-Correcting Codes surely exist and can be obtained with high probability by random sampling. These derive from exponential bounds on the probability of drawing ‘bad’ matrices. The bounds give conditions effective at finite-$N$, and converging to the known sharp asymptotic conditions for large$N$. Compared to other finite-$N$bounds known to us, they are much stronger, and much more explicit. Our bounds derive from asymptotics in“Counting faces of randomly projected polytopes when the projection radically lowers dimension”counting the expected number of$k$-dimensional faces of the randomly projected simplex$T^{N-1}$and cross-polytope$C^N$. We develop here finite-$N$bounds on the expected discrepancy between the number of$k$-faces of the projected polytope$AQ$and its generator$Q$, for$Q=T^{N-1}$and$C^N$. Our bounds also imply existence of interesting geometric objects. Thus, we exhibit triples$(k,n,N)$for which polytopes with$2N$vertices can be centrally$k$-neighborly.

The minimum lscr1-norm solution to an underdetermined system of linear equationsy=Axis often, remarkably, also the sparsest solution to that system. This sparsity-seeking property is of interest in signal processing and information transmission. However, general-purpose optimizers are much too slow for lscr1minimization in many large-scale applications.In this paper, the Homotopy method, originally proposed by Osborne et al. and Efron et al., is applied to the underdetermined lscr1-minimization problem min parxpar1subject toy=Ax. Homotopy is shown to run much more rapidly than general-purpose LP solvers when sufficient sparsity is present. Indeed, the method often has the followingk-step solution property: if the underlying solution has onlyknonzeros, the Homotopy method reaches that solution in onlykiterative steps. Thisk-step solution property is demonstrated for several ensembles of matrices, including incoherent matrices, uniform spherical matrices, and partial orthogonal matrices. These results imply that Homotopy may be used to rapidly decode error-correcting codes in a stylized communication system with a computational budget constraint. The approach also sheds light on the evident parallelism in results on lscr1minimization and orthogonal matching pursuit (OMP), and aids in explaining the inherent relations between Homotopy, least angle regression (LARS), OMP, and polytope faces pursuit.

In a recent paper, a method called morphological component analysis (MCA) has been proposed to separate the texture from the natural part in images. MCA relies on an iterative thresholding algorithm, using a threshold which decreases linearly towards zero along the iterations. This paper shows how the MCA convergence can be drastically improved using the mutual incoherence of the dictionaries associated to the different components. This modified MCA algorithm is then compared to basis pursuit, and experiments show that MCA and BP solutions are similar in terms of sparsity, as measured by the l1 norm, but MCA is much faster and gives us the possibility of handling large scale data sets.

Many applications in signal processing lead to the optimization problems min parxpar1subject to y = Ax, and min parxpar1subject to pary - Axpar les epsi, where A is a given d times n matrix, d1minimization by using LARS, Lasso, and homotopy methods (Efron et al., Tibshirani, Osborne et al.). While these methods were first proposed for use in statistical model selection, we show that under certain conditions these methods find the sparsest solution rapidly, as opposed to conventional general purpose optimizers which are prohibitively slow. We define a phase transition diagram which shows how algorithms behave for random problems, as the ratio of unknowns to equations and the ratio of the sparsity to equations varies. We find that whenever the number k of nonzeros in the sparsest solution is less than d/2log(n) then LARS/homotopy obtains the sparsest solution in k steps each of complexity O(d2)

The classical multivariate linear regression problem assumes p variables X1, X2,... ,Xpand a response vector y, each with n observations, and a linear relationship between the two: y = Xbeta + z, where z ~ N(0, sigma2). We point out that when p > n, there is a breakdown point for standard model selection schemes, such that model selection only works well below a certain critical complexity level depending on n/p. We apply this notion to some standard model selection algorithms (Forward Stepwise, LASSO, LARS) in the case where pGtn. We find that 1) the breakdown point is well-de ned for random X-models and low noise, 2) increasing noise shifts the breakdown point to lower levels of sparsity, and reduces the model recovery ability of the algorithm in a systematic way, and 3) below breakdown, the size of coef cient errors follows the theoretical error distribution for the classical linear model.

Consider a d times n matrix A, with d < n. The problem of solving for x in y = Ax is underdetermined, and has infinitely many solutions (if there are any). Given y, the minimum Kolmogorov complexity solution (MKCS) of the input x is defined to be an input z (out of many) with minimum Kolmogorov-complexity that satisfies y = Az. One expects that if the actual input is simple enough, then MKCS will recover the input exactly. This paper presents a preliminary study of the existence and value of the complexity level up to which such a complexity-based recovery is possible. It is shown that for the set of all d times n binary matrices (with entries 0 or 1 and d < n), MKCS exactly recovers the input for an overwhelming fraction of the matrices provided the Kolmogorov complexity of the input is O(d). A weak converse that is loose by a log n factor is also established for this case. Finally, we investigate the difficulty of finding a matrix that has the property of recovering inputs with complexity of O(d) using MKCS

Suppose x is an unknown vector in Ropfm(a digital image or signal); we plan to measure n general linear functionals of x and then reconstruct. If x is known to be compressible by transform coding with a known transform, and we reconstruct via the nonlinear procedure defined here, the number of measurements n can be dramatically smaller than the size m. Thus, certain natural classes of images with m pixels need only n=O(m1/4log5/2(m)) nonadaptive nonpixel samples for faithful recovery, as opposed to the usual m pixel samples. More specifically, suppose x has a sparse representation in some orthonormal basis (e.g., wavelet, Fourier) or tight frame (e.g., curvelet, Gabor)-so the coefficients belong to an lscrpball for 02error O(N1/2-1p/). It is possible to design n=O(Nlog(m)) nonadaptive measurements allowing reconstruction with accuracy comparable to that attainable with direct knowledge of the N most important coefficients. Moreover, a good approximation to those N important coefficients is extracted from the n measurements by solving a linear program-Basis Pursuit in signal processing. The nonadaptive measurements have the character of "random" linear combinations of basis/frame elements. Our results use the notions of optimal recovery, of n-widths, and information-based complexity. We estimate the Gel'fand n-widths of lscrpballs in high-dimensional Euclidean space in the case 0

Overcomplete representations are attracting interest in signal processing theory, particularly due to their potential to generate sparse representations of signals. However, in general, the problem of finding sparse representations must be unstable in the presence of noise. This paper establishes the possibility of stable recovery under a combination of sufficient sparsity and favorable structure of the overcomplete system. Considering an ideal underlying signal that has a sufficiently sparse representation, it is assumed that only a noisy version of it can be observed. Assuming further that the overcomplete system is incoherent, it is shown that the optimally sparse approximation to the noisy data differs from the optimally sparse decomposition of the ideal noiseless signal by at most a constant multiple of the noise level. As this optimal-sparsity method requires heavy (combinatorial) computational effort, approximation algorithms are considered. It is shown that similar stability is also available using the basis and the matching pursuit algorithms. Furthermore, it is shown that these methods result in sparse approximation of the noisy data that contains only terms also appearing in the unique sparsest representation of the ideal noiseless sparse signal.

We briefly recall previous literature about the interaction between sparsity and /spl lscr/ /sup 1/ minimization. We then discuss /spl lscr/ /sup 1/ minimization in geometry separation, in compressed sensing, and in compressed sensing of separable signals.

The images generated by varying the underlying articulation parameters of an object (pose, attitude, light source position, and so on) can be viewed as points on a low-dimensional image parameter articulation manifold (IPAM) in a high-dimensional ambient space. In this paper, we develop theory and methods for the inverse problem of estimating, from a given image on or near an IPAM, the underlying parameters that produced it. Our approach is centered on the observation that, while typical image manifolds are not differentiable, they have an intrinsic multiscale geometric structure. In fact, each IPAM has a family of approximate tangent spaces, each one good at a certain resolution. Putting this structural aspect to work, we develop a new algorithm for high-accuracy parameter estimation based on a coarse-to-fine Newton iteration through the family of approximate tangent spaces. We test the algorithm in several idealized registration and pose estimation problems.

A new class of related algorithms for deblocking block-transform compressed images and video sequences is proposed in this paper. The algorithms apply weighted sums on pixel quartets, which are symmetrically aligned with respect to block boundaries. The basic weights, which are aimed at very low bit-rate images, are obtained from a two-dimensional function which obeys predefined constraints. Using these weights on images compressed at higher bit rates produces a deblocked image which contains blurred "false" edges near real edges. We refer to this phenomenon as the ghosting effect. In order to prevent its occurrences, the weights of pixels, which belong to nonmonotone areas, are modified by dividing each pixel's weight by a predefined factor called a grade. This scheme is referred to as weight adaptation by grading (WABG). Better deblocking of monotone areas is achieved by applying three iterations of the WABG scheme on such areas followed by a fourth iteration which is applied on the rest of the image. We refer to this scheme as deblocking frames of variable size (DFOVS). DFOVS automatically adapts itself to the activity of each block. This new class of algorithms produces very good subjective results and PSNR results which are competitive relative to available state-of-the-art methods.

The separation of image content into semantic parts plays a vital role in applications such as compression, enhancement, restoration, and more. In recent years, several pioneering works suggested such a separation be based on variational formulation and others using independent component analysis and sparsity. This paper presents a novel method for separating images into texture and piecewise smooth (cartoon) parts, exploiting both the variational and the sparsity mechanisms. The method combines the basis pursuit denoising (BPDN) algorithm and the total-variation (TV) regularization scheme. The basic idea presented in this paper is the use of two appropriate dictionaries, one for the representation of textures and the other for the natural scene parts assumed to be piecewise smooth. Both dictionaries are chosen such that they lead to sparse representations over one type of image-content (either texture or piecewise smooth). The use of the BPDN with the two amalgamed dictionaries leads to the desired separation, along with noise removal as a by-product. As the need to choose proper dictionaries is generally hard, a TV regularization is employed to better direct the separation process and reduce ringing artifacts. We present a highly efficient numerical scheme to solve the combined optimization problem posed by our model and to show several experimental results that validate the algorithm's performance.

We construct detectors for "geometric" objects in noisy data. Examples include a detector for presence of a line segment of unknown length, position, and orientation in two-dimensional image data with additive white Gaussian noise. We focus on the following two issues. i) The optimal detection threshold-i.e., the signal strength below which no method of detection can be successful for large dataset size n. ii) The optimal computational complexity of a near-optimal detector, i.e., the complexity required to detect signals slightly exceeding the detection threshold. We describe a general approach to such problems which covers several classes of geometrically defined signals; for example, with one-dimensional data, signals having elevated mean on an interval, and, in d-dimensional data, signals with elevated mean on a rectangle, a ball, or an ellipsoid. In all these problems, we show that a naive or straightforward approach leads to detector thresholds and algorithms which are asymptotically far away from optimal. At the same time, a multiscale geometric analysis of these classes of objects allows us to derive asymptotically optimal detection thresholds and fast algorithms for near-optimal detectors.

Traffic matrices are required inputs for many IP network management tasks, such as capacity planning, traffic engineering, and network reliability analysis. However, it is difficult to measure these matrices directly in large operational IP networks, so there has been recent interest in inferring traffic matrices from link measurements and other more easily measured data. Typically, this inference problem is ill-posed, as it involves significantly more unknowns than data. Experience in many scientific and engineering fields has shown that it is essential to approach such ill-posed problems via "regularization". This paper presents a new approach to traffic matrix estimation using a regularization based on "entropy penalization". Our solution chooses the traffic matrix consistent with the measured data that is information-theoretically closest to a model in which source/destination pairs are stochastically independent. It applies to both point-to-point and point-to-multipoint traffic matrix estimation. We use fast algorithms based on modern convex optimization theory to solve for our traffic matrices. We evaluate our algorithm with real backbone traffic and routing data, and demonstrate that it is fast, accurate, robust, and flexible.

Beamlets are combined with zero-tree coding algorithm to create a new coding method, particularly suitable for lines and curves. Beamlets are a multiscale collection of line segments at a range of scales, location, having a variety of lengths and orientations. The new coding scheme - named JBEAM is more efficient - in terms of bit rates - in coding binary curvy images in simulations.

We interpret non-negative matrix factorization geometrically, as the problem of finding a simplicial cone which contains a cloud of data points and which is contained in the positive orthant. We show that under certain conditions, basically requiring that some of the data are spread across the faces of the positive orthant, there is a unique such simpli- cial cone. We give examples of synthetic image articulation databases which obey these conditions; these require separated support and facto- rial sampling. For such databases there is a generative model in terms of ‘parts’ and NMF correctly identifies the ‘parts’. We show that our theoretical results are predictive of the performance of published NMF code, by running the published algorithms on one of our synthetic image articulation databases.

Traffic matrices are required inputs for many IP network management tasks: for instance, capacity planning, traffic engineering and network reliability analysis. However, it is difficult to measure these matrices directly, and so there has been recent interest in inferring traffic matrices from link measurements and other more easily measured data. Typically, this inference problem is ill-posed, as it involves significantly more unknowns than data. Experience in many scientific and engineering fields has shown that it is essential to approach such ill-posed problems via "regularization". This paper presents a new approach to traffic matrix estimation using a regularization based on "entropy penalization". Our solution chooses the traffic matrix consistent with the measured data that is information-theoretically closest to a model in which source/destination pairs are stochastically independent. We use fast algorithms based on modern convex optimization theory to solve for our traffic matrices. We evaluate the algorithm with real backbone traffic and routing data, and demonstrate that it is fast, accurate, robust, and flexible.

We present in this paper a new method for contrast enhancement based on the curvelet transform. The curvelet transform represents edges better than wavelets, and is therefore well-suited for multiscale edge enhancement. We compare this approach with enhancement based on the wavelet transform, and the Multiscale Retinex. In a range of examples, we use edge detection and segmentation, among other processing applications, to provide for quantitative comparative evaluation. Our findings are that curvelet based enhancement out-performs other enhancement methods on noisy images, but on noiseless or near noiseless images curvelet based enhancement is not remarkably better than wavelet based enhancement.

We consider the problem of recovering a binary image consisting of many filaments or linear fragments in the presence of severe binary noise. Our approach exploits beamlets—a dyadically organized, multiscale system of line segments—and associated fast algorithms for beamlet analysis. It considers models based on beamlet-decorated recursive dyadic partitions, and models the image as a Bernoulli random process with spatially variant success probability, which is “high” within the beamlet complexity-penalized model fitting. Simulation results demonstrate the effectiveness of the method.

A new algorithm for the removal of blocking artifacts in block-DCT compressed images and video sequences is proposed in this paper. The algorithm uses deblocking frames of variable size (DFOVS). A deblocking frame is a square of pixels which overlaps image blocks. Deblocking is achieved by applying weighted summation on pixels quartets which reside in deblocking frames. The pixels in a quartet are symmetrically aligned with respect to block boundaries. The weights are determined according to a special 2-D function and a predefined factor, we refer to as a grade. A pixel's grade is determined according to the amount of details in its neighborhood. Deblocking of monotone areas is achieved by iteratively applying deblocking frames of decreasing sizes on such areas. This new algorithm produces very good subjective results and PSNR results which are competitive relative to available state-of-the-art methods.

We describe approximate digital implementations of two new mathematical transforms, namely, the ridgelet transform and the curvelet transform. Our implementations offer exact reconstruction, stability against perturbations, ease of implementation, and low computational complexity. A central tool is Fourier-domain computation of an approximate digital Radon transform. We introduce a very simple interpolation in the Fourier space which takes Cartesian samples and yields samples on a rectopolar grid, which is a pseudo-polar sampling set based on a concentric squares geometry. Despite the crudeness of our interpolation, the visual performance is surprisingly good. Our ridgelet transform applies to the Radon transform a special overcomplete wavelet pyramid whose wavelets have compact support in the frequency domain. Our curvelet transform uses our ridgelet transform as a component step, and implements curvelet subbands using a filter bank of a; trous wavelet filters. Our philosophy throughout is that transforms should be overcomplete, rather than critically sampled. We apply these digital transforms to the denoising of some standard images embedded in white noise. In the tests reported here, simple thresholding of the curvelet coefficients is very competitive with "state of the art" techniques based on wavelets, including thresholding of decimated or undecimated wavelet transforms and also including tree-based Bayesian posterior mean methods. Moreover, the curvelet reconstructions exhibit higher perceptual quality than wavelet-based reconstructions, offering visually sharper images and, in particular, higher quality recovery of edges and of faint linear and curvilinear features. Existing theory for curvelet and ridgelet transforms suggests that these new approaches can outperform wavelet methods in certain image reconstruction problems. The empirical results reported here are in encouraging agreement.

Suppose a discrete-time signal S(t), 0/spl les/t<N, is a superposition of atoms taken from a combined time-frequency dictionary made of spike sequences 1/sub {t=/spl tau/}/ and sinusoids exp{2/spl pi/iwt/N}//spl radic/N. Can one recover, from knowledge of S alone, the precise collection of atoms going to make up S? Because every discrete-time signal can be represented as a superposition of spikes alone, or as a superposition of sinusoids alone, there is no unique way of writing S as a sum of spikes and sinusoids in general. We prove that if S is representable as a highly sparse superposition of atoms from this time-frequency dictionary, then there is only one such highly sparse representation of S, and it can be obtained by solving the convex optimization problem of minimizing the l/sup 1/ norm of the coefficients among all decompositions. Here "highly sparse" means that N/sub t/+N/sub w/

We propose a best basis algorithm for signal enhancement in white Gaussian noise. The best basis search is performed in families of orthonormal bases constructed with wavelet packets or local cosine bases. We base our search for the "best" basis on a criterion of minimal reconstruction error of the underlying signal. This approach is intuitively appealing, because the enhanced or estimated signal has an associated measure of performance, namely, the resulting mean-square error. Previous approaches in this framework have focused on obtaining the most "compact" signal representations, which consequently contribute to effective denoising. These approaches, however, do not possess the inherent measure of performance which our algorithm provides. We first propose an estimator of the mean-square error, based on a heuristic argument and subsequently compare the reconstruction performance based upon it to that based on the Stein (1981) unbiased risk estimator. We compare the two proposed estimators by providing both qualitative and quantitative analyses of the bias term. Having two estimators of the mean-square error, we incorporate these cost functions into the search for the "best" basis, and subsequently provide a substantiating example to demonstrate their performance.

We apply basis pursuit, an atomic decomposition technique, for spectrum estimation. Compared with several modern time series methods, our approach can greatly reduce the problem of power leakage; it is able to superresolve; moreover, it works well with noisy and unevenly sampled signals. We present experiments on bizarrely spaced radial velocity data from one of the newly-discovered extrasolar planetary systems.

Automatic modulation classification (or recognition) is an intrinsically interesting problem with a variety of regulatory and military applications. We developed a method which is simple, fast, efficient and robust. The feature being used is the counts of signals falling into different parts of the signal plane. Compared with the likelihood method and the high order correlation method, it is much easier to be implemented, and the execution is much faster. When the channel model is correct, our method is efficient, in the sense that it will achieve the "optimal" classification rate. When unknown contamination is present, our method can automatically overcome it to certain degree. At SNRs of 10 and 15 dB, examples of classifying two modulation types-QAM4 and PSK6-are given. Simulations demonstrate its ability to deal with unknown noise.

In this paper we review some recent interactions between harmonic analysis and data compression. The story goes back of course to Shannon's R(D) theory in the case of Gaussian stationary processes, which says that transforming into a Fourier basis followed by block coding gives an optimal lossy compression technique; practical developments like transform-based image compression have been inspired by this result. In this paper we also discuss connections perhaps less familiar to the information theory community, growing out of the field of harmonic analysis. Recent harmonic analysis constructions, such as wavelet transforms and Gabor transforms, are essentially optimal transforms for transform coding in certain settings. Some of these transforms are under consideration for future compression standards. We discuss some of the lessons of harmonic analysis in this century. Typically, the problems and achievements of this field have involved goals that were not obviously related to practical data compression, and have used a language not immediately accessible to outsiders. Nevertheless, through an extensive generalization of what Shannon called the "sampling theorem", harmonic analysis has succeeded in developing new forms of functional representation which turn out to have significant data compression interpretations. We explain why harmonic analysis has interacted with data compression, and we describe some interesting recent ideas in the field that may affect data compression in the future.

We propose a best basis algorithm for signal enhancement in white Gaussian noise. We base our search of best basis on a criterion of minimal reconstruction error of the underlying signal. We subsequently compare our simple error criterion to the Stein (1981) unbiased risk estimator, and provide a substantiating example to demonstrate its performance. A review is also given of noise removal by thresholding and of wavepacket orthonormal bases.

Donoho and Johnstone (1994) proposed a method for reconstructing an unknown function f on [0,1] from noisy data d/sub i/=f(t/sub i/)+/spl sigma/z/sub i/, i=0, ..., n-1,t/sub i/=i/n, where the z/sub i/ are independent and identically distributed standard Gaussian random variables. The reconstruction f/spl circ/*/sub n/ is defined in the wavelet domain by translating all the empirical wavelet coefficients of d toward 0 by an amount /spl sigma//spl middot//spl radic/(2log (n)/n). The authors prove two results about this type of estimator. [Smooth]: with high probability f/spl circ/*/sub n/ is at least as smooth as f, in any of a wide variety of smoothness measures. [Adapt]: the estimator comes nearly as close in mean square to f as any measurable estimator can come, uniformly over balls in each of two broad scales of smoothness classes. These two properties are unprecedented in several ways. The present proof of these results develops new facts about abstract statistical inference and its connection with an optimal recovery model.>

If the measure of the support of a function f is small, its symmetric decreasing rearrangement f* is more nearly bandlimited to low frequencies than f, while their norms are equal. An immediate corollary is that the time-limited zero-order prolate spheroidal wavefunction is the extremal function for a new optimization problem involving time- and bandlimiting. The result has an application in exploration seismology.>