Matrix Factorization Techniques in Machine Learning, Signal Processing, and Statistics

Abstract

Compressed sensing is an alternative to Shannon/Nyquist sampling for acquiring sparse or compressible signals. Sparse coding represents a signal as a sparse linear combination of atoms, which are elementary signals derived from a predefined dictionary. Compressed sensing, sparse approximation, and dictionary learning are topics similar to sparse coding. Matrix completion is the process of recovering a data matrix from a subset of its entries, and it extends the principles of compressed sensing and sparse approximation. The nonnegative matrix factorization is a low-rank matrix factorization technique for nonnegative data. All of these low-rank matrix factorization techniques are unsupervised learning techniques, and can be used for data analysis tasks, such as dimension reduction, feature extraction, blind source separation, data compression, and knowledge discovery. In this paper, we survey a few emerging matrix factorization techniques that are receiving wide attention in machine learning, signal processing, and statistics. The treated topics are compressed sensing, dictionary learning, sparse representation, matrix completion and matrix recovery, nonnegative matrix factorization, the Nyström method, and CUR matrix decomposition in the machine learning framework. Some related topics, such as matrix factorization using metaheuristics or neurodynamics, are also introduced. A few topics are suggested for future investigation in this article.

1. Introduction

Matrix factorization is widely used for inferring the structure in multivariate data. Given a noisy measurement of the product of two matrices, the matrix factorization problem aims to estimate the original matrices. It represents an observed data matrix $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$ as

$$\mathbf{Y}={\mathbf{L}}^T\mathbf{F}+\mathbf{E},$$

(1)

where $\mathbf{L}\in {\mathbb{R}}^{k\times m}$, $\mathbf{F}\in {\mathbb{R}}^{k\times n}$, and the residual matrix is denoted as $\mathbf{E}\in {\mathbb{R}}^{m\times n}$, which is usually assumed to have normally distributed entries. For factor analysis, $\mathbf{L}$ is referred to as the loadings, and $\mathbf{F}$ is referred to as the factors.

Matrix factorization is a bilinear inverse problem in individual matrices. Model (1) has many applications. For the matrix completion problem, the estimations of $\mathbf{L}$ and $\mathbf{F}$ from partially observed $\mathbf{Y}$ provide a natural and simple way to estimate the missing entries; this is a typical scenario in matrix completion tasks. Another wide range of applications involve inferring and summarizing the structures in multivariate data in $\mathbf{Y}$, where each row of $\mathbf{Y}$ is approximated by a linear combination of the rows of $\mathbf{F}$, which are referred to as factors, corresponding to factor analysis, the principal component analysis (PCA), or dictionary learning.

The dictionary $\mathbf{A}$, corresponding to ${\mathbf{L}}^T$ in (1), is subject to some constraints. A prominent example is PCA [1], where $\mathbf{A}$ has orthogonal columns, representing the subspace where the signal in the given class is contained. Another example is sparse coding, where $\mathbf{A}$ typically consists of normalized columns that form an overcomplete basis of the signal space, and the signal $\mathbf{x}\in {\mathbb{R}}^k$, corresponding to $\mathbf{F}$ degenerating into a column vector (i.e., $n=1$), is assumed to be sparse. Matrix factorization is more commonly, represented by

$$\mathbf{Y}=\mathbf{A}\mathbf{X}+\mathbf{E},$$

(2)

where $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$ is the data matrix, $\mathbf{A}\in {\mathbb{R}}^{m\times k}$ is the dictionary matrix, and $\mathbf{X}\in {\mathbb{R}}^{k\times n}$ is the code matrix. Each column of $\mathbf{Y}$ is approximated by a linear combination of the columns of $\mathbf{A}$, where the coefficients are given by the corresponding column of matrix $\mathbf{X}$. This is illustrated in Figure 1.

Matrix factorization provides a low-rank approximation of a matrix. It arises in many machine learning and signal processing applications [2], such as singular value decomposition (SVD), factor analysis, PCA, blind source separation, independent component analysis (ICA), blind matrix calibration, dictionary learning, low-rank matrix completion, nonnegative matrix factorization (NMF), C-means clustering [3], unsupervised representation learning, and so on. A penalty or prior distribution is usually used to achieve sparse representations (e.g., sparse factor analysis and sparse PCA). As a related concept, matrix decomposition aims to recover, from a matrix, a low-rank matrix and a sparse matrix. These techniques fall under the category of unsupervised learning in the machine learning framework [4].

Representation learning is a concept behind many machine learning applications, including the deep learning framework. In the context of representation learning, dictionary learning is sometimes referred to as sparse coding. Sparse coding has been proposed as a theory for modeling the visual cortex and as an unsupervised algorithm for learning representations. These methods are widely used for feature extraction and knowledge discovery. Low-rank matrix approximation is a ubiquitous problem in data processing.

Compressed sensing is a powerful framework used for acquiring sparse signals. Learning a sparsifying dictionary or transforming from compressive measurements [5,6] requires fewer equations to determine the unknowns, compared to matrix factorization. Dictionary learning aims to find a good sparse representation for a dataset. It is a matrix factorization problem with a sparsity constraint. Compressive blind source separation [7] is another matrix factorization problem with a sparsity constraint. These problems can be solved by recovering a sparse and low-rank matrix from its linear measurements. NMF can be made equivalent to some clustering problems by reframing them slightly [8]. NMF and spectral clustering are two popular clustering techniques. However, NMF cannot deal with nonlinear data, and spectral clustering relies on post-processing.

Crowdsourcing is a scalable approach to collecting data from humans. Through crowdsourcing platforms such as Amazon Mechanical Turk, a large number of data tasks are assigned to workers for binary or multiclass labeling. The goal is to estimate the unknown ground truth from the input of the various workers. The unknown labels can be estimated through aggregation, such as majority voting or weighted majority voting.

Sparsity is a natural property that consists of many real signals. In the human brain, very few neurons (1–4%) are active at any time [9]. Many observed response properties in the primary visual cortex (V1), such as the classical receptive field structure [10] and nonclassical response modulations [11], are accounted for by using the sparse coding model of the primary visual cortex. In the sparse coding model, most natural images are encoded by very few learned dictionary elements, and high-dimensional visual inputs can be represented by a small number of active cortical neurons [10].

Sparse coding involves learning an overcomplete set of basis vectors, where each data point is represented as a sparse combination of the basis vectors. On the contrary, sparse recovery aims to reconstruct sparse signals from an underdetermined set of compressed linear measurements.

The recovery of sparse signals corrupted by additive noise covers a wide range of applications, such as image inpainting, super-resolution, signal separation, and recovery of signals that are impaired by clipping, impulse noise, or narrowband interference.

For linear inverse or compressed sensing problems, a series of representer theorems give the generic form of the solution, which depends on whether $L_1$- or $L_2$-norm regularization is under consideration [12]. $L_1$-norm solutions are proven to be intrinsically sparse, and the use of the $L_1$-norm regularization is much more favorable for incorporating prior knowledge compared to the $L_2$-norm scenario. $L_1$-norm has long been used for pruning neural network architecture [13].

Compressed sensing, also known as compressive sampling, is a recent sampling method [14,15]. If a signal is sufficiently sparse, it can be exactly reconstructed from very few random measurements. Compressed sensing serves as an alternative to classical Shannon/Nyquist sampling for sparse or compressible signals, allowing for the perfect recovery of sparse signals using only a small number of random measurements. Compressible signals can be well approximated by sparse signals. A non-sparse signal can be compressed by using compressive covariance sensing [16], and its second-order statistics can be recovered from the compressed signal without sparsity constraint.

Low-rank representation is usually used to recover data from corruption or outliers. PCA is the best low-rank representation in terms of $L_2$ errors. By minimizing the $L_2$-norm error, data are projected onto the fixed-rank low-dimensional space. Robust PCA [17] and GoDec [18,19] decompose data into low-rank components and sparse components that capture corruptions.

Matrix completion [20] recovers a matrix from a subset of its entries. It had been prevalent in computer vision, statistics, collaborative filtering, and manifold learning in the last decade. It is an ill-posed problem. A common constraint is applied to the underlying matrix. The task is formulated as a low-rank matrix approximation problem. The method is related to compressed sensing. The values of matrix entries may be discrete or quantized, such as in the Netflix problem and recommender systems.

Many real-life data or physical signals are represented by nonnegative numbers. When analyzing mixtures of such data, nonnegativity constraints on the individual components are applied. Nonnegative PCA, nonnegative ICA, and NMF [21] are techniques used for the analysis of such data, where nonnegative data are represented as nonnegative linear combinations of nonnegative bases. Nonnegativity is inspired by neuronal properties, such as the firing rate representation and signed synaptic weight [21].

The inferior temporal cortex is a critical region in the primate visual cortex for object recognition. Object representation in this region has two prominent features. An object is represented by a combination of the activities of columnar clusters of neurons, where each cluster represents component features or parts of objects [22].

NMF [21,23] factorizes a matrix as a product of two matrices, whose elements are all nonnegative. Nonnegativity prevents mutual cancellation between basis functions and, thus, generates a parts-based representation, in agreement with human thinking. NMF can be used for tasks such as blind source separation (BSS) of images and nonnegative signals [24], spectra recovery [25], feature extraction, and clustering [26].

Canonical correlation analysis, SVD, PCA, and ICA are classical matrix factorization methods, derived by the low-rank approximation to a matrix by minimizing the squared error. Latent semantic indexing [27] is an application that uses SVD for automatic indexing and retrieval. We do not describe them here. In this paper, we provide a survey on the recent matrix factorization techniques that are prevalent in machine learning and signal processing. In Section 2, we describe compressed sensing. Section 3 introduces sparse coding and dictionary learning. Section 4 extends sparse coding to matrix completion. In Section 5, low-rank representation is reviewed. In Section 6, NMF is introduced. Section 7 introduces techniques for symmetric positive semidefinite matrix approximation, including the Nyström method. Section 8 describes the CX Decomposition and CUR decomposition. Finally, a summary is given in Section 9.

2. Compressed Sensing

Compressed sensing seeks to recover sparse or compressible signals from undersampled linear measurements [15,28]. A sparse or compressible high-dimensional signal can be projected onto a low-dimensional space when applying a random observation matrix. Compressibility of data and acquisition of incoherent measurements are the two fundamental properties underlying compressed sensing.

Given a signal $\overrightarrow{x}$, if $\overrightarrow{\alpha }={\mathsf{\Phi }}^T\overrightarrow{x}$ is sparsely distributed for any dictionary $\mathsf{\Phi }$, $\overrightarrow{x}$ is said to be compressible.

2.1. Signal Model

In compressed sensing, a signal with N samples, $\overrightarrow{x}\in {\mathbb{R}}^N$, is derived from a set of linear measurements

$$\overrightarrow{y}=\mathbf{A}\overrightarrow{x}+\overrightarrow{n},$$

(3)

where $\mathbf{A}\in {\mathbb{R}}^{M\times N}$ denotes a random sampling, sensing, or measurement matrix, $\overrightarrow{y}\in {\mathbb{R}}^M$ is a measurement vector with M measurements, $M<N$, and $\overrightarrow{n}\in {\mathbb{R}}^M$ is noise.

The problem (3) is underdetermined. When the norm of each column of $\mathbf{A}$ is unity, all of the columns form an incomplete basis with $M\ll N$.

In order to preserve the information in sparse or compressible signals $\overrightarrow{x}$, and to ensure the stable recovery of such signals, $\mathbf{A}$ has to satisfy the so-called restricted isometry property (RIP) [14].

$\mathbf{A}$ should also satisfy the incoherence property. $M=O(klnN)$ measurements are sufficient for perfectly reconstructing a k-sparse vector $\overrightarrow{x}$, if $\mathbf{A}$ is perfectly incoherent (e.g., uniformly random Gaussian measurements).

Compressed sensing achieves the stable recovery of compressible, noisy signals by solving the $L_0$-norm regularized inverse problem, or the corresponding computationally tractable $L_1$-norm problem

$$\underset{\overrightarrow{x}}{min}\parallel \overrightarrow{x}{\parallel }_1\mathrm{or}\parallel \overrightarrow{x}{\parallel }_0\mathrm{subject} \mathrm{to}{\parallel \mathbf{A}\overrightarrow{x}-\overrightarrow{y}\parallel }_2^2\le {\epsilon }^2,$$

(4)

where $L_0$-norm $\parallel \overrightarrow{x}{\parallel }_0$ counts the nonzero entries in $\overrightarrow{x}$, $L_1$-norm is the convex envelope of $L_0$-norm, $\parallel \overrightarrow{x}{\parallel }_1={\sum }_{i=1}^N\left|x_i\right|$, and $\epsilon$ is a tolerance.

The $L_0/L_1$-regularized least squares (LS) approach is used to deal with linear inverse problems under sparsity constraints. Linear programming (LP) has the best sparsity–undersampling trade-off, at a cost of high computation complexity. The approximate message-passing algorithm [29] is an iterative thresholding algorithm corresponding to the LP procedure, but is dramatically faster.

The east absolute selection and shrinkage operator (LASSO) [30] and approximate message passing [29] are low-complexity reconstruction procedures. By minimizing a weighted sum of the residual norm and a regularization term $\parallel \overrightarrow{x}{\parallel }_1$, LASSO can reconstruct sparse solutions and recover the sparsity pattern exactly as the number of observations increases, asymptotically with probability one.

Standard compressed sensing guarantees robust signal recovery from $O(klog\frac{N}{k})$ measurements with provable performance guarantees [31]. Based on a model-based compressed sensing theory, two recovery algorithms incorporating wavelet trees and block sparsity are proven to offer robust recovery from $O\left(k\right)$ measurements [31].

In [32], sensing vectors are selected independently at random from a probability distribution $\mathcal{F}$. If $\mathcal{F}$ obeys an incoherence property and an isotropy property, approximately sparse signals can be faithfully recovered from a minimum number of noisy measurements. The recovery does not require the RIP to hold near the sparsity level, and it also does not need a random model for the signal. A k-sparse signal can be faithfully recovered from about $klogN$ noisy Fourier coefficients.

When the measurements are obtained using a matrix with i.i.d. Gaussian entries, the weighted $L_1$-norm minimization recovers the sparse signal with overwhelmingly high probability [33]. For any stationary process satisfying certain mixing conditions, if the sampling rate is greater than the information dimension of the source process, the minimum entropy pursuit (MEP) optimization approach for universal compressed sensing can reliably recover the source vector almost losslessly, without any prior information about its distribution [34].

2.2. RIP, ERC, and MIP

Conditions used in the compressed sensing literature include the RIP [28], exact recovery condition (ERC) [35], and mutual inheritance property (MIP).

2.2.1. RIP

The RIP of sampling matrices is a sufficient condition for the reliable reconstruction of sparse signals [15]. RIP matrices can be constructed from binary vectors [36]. Both the algorithmic and constructive aspects were pursued to connect to error correction codes [37,38].

A k-sparse vector $\overrightarrow{x}$ has, at most, k nonzero entries, $\parallel \overrightarrow{x}{\parallel }_0\le k$. For any k-sparse vector $\overrightarrow{x}\in {\mathbb{R}}^N$, if there exists a constant $\delta \in (0,1)$, such that [14,39]

$$(1-\delta )\parallel \overrightarrow{x}{\parallel }_2^2\le \parallel \mathbf{A}\overrightarrow{x}{\parallel }_2^2\le (1+\delta ){\parallel \overrightarrow{x}\parallel }_2^2,$$

(5)

a sensing matrix $\mathbf{A}\in {\mathbb{R}}^{M\times N}$ is said to satisfy the k-restricted isometry property (k-RIP).

The minimum of all $\delta \in (0,1)$, denoted as ${\delta }_k$, is referred to as the kth-order restricted isometry constant (RIC), or simply RIC, of $\mathbf{A}$. A smaller RIC corresponds to a transformation closer to an isometry. The RIC has monotonicity properties, ${\delta }_k\le {\delta }_{2k}$. Via a RIP analysis, one can find the maximum sparsity order k that guarantees the recovery of all sparse vectors.

k-RIP ensures that all $M\times k$ submatrices of $\mathbf{A}$ approximately satisfy the isometry property and, hence, distance preserving. When M is close to k, maximal signal compression is achieved. RIP can be used to measure the orthogonality of column vectors of a dictionary. The problem of determining whether $\mathbf{A}$ has RIP for any accuracy is proven to be NP-hard [40]. Constructing RIP matrices deterministically is a hard problem [41], but RIP matrices can be generated with high probability by simple random methods [42,43].

The RIP analysis shows that Gaussian measurement matrices are information-theoretically optimal since the required number of measurements for sparse recovery is minimal [28,44]. A k-RIP matrix must have at least $M=\mathsf{\Omega }(klog\frac{N}{k})$ rows. That is, $M=\mathsf{\Omega }(klog\frac{N}{k})$ samples are required for any recovery algorithm to approximate the signal $\overrightarrow{x}$ with an accuracy expressed by the $L_1$- or $L_2$-norm [45,46]. This bound is applicable to non-k-sparse signals. Random Gaussian or Bernoulli matrices provide, with high probability, the best-known accuracy for recovery that matches this lower bound [28,47,48].

When most $M\times k$ submatrices define a near-isometric map of ${\mathbb{R}}^k$ into ${\mathbb{R}}^M$, a matrix has a statistical RIP of order k [49]. Statistical RIP is particularly useful for the sparse signal recovery of deterministic sensing matrices. For many existing deterministic sampling matrices, $M=O\left(k\right)$ rows guarantee k-statistical RIP [49]. With the conditions of statistical RIP and statistical incoherence, stable sparse recovery by a basis pursuit can be proved [49].

The k-RIP concept is extended to the case of the $L_p$-norm [50]. For every $1\le p<\infty$, almost tight bounds on the minimum number of rows, M, which are necessary for the RIP to hold, as well as almost tight bounds for the column sparsity of RIP matrices, are obtained [50].

2.2.2. ERC

The ERC is a necessary and sufficient condition for exact support recovery in a worst-case analysis [35]. If a subset of atoms satisfies the ERC, then it can be recovered from any linear combination of the atoms in, at most, k steps. The ERC necessarily holds when the latter conditions are fulfilled since the ERC is a worst-case necessary condition for exact recovery.

When the ERC is met, the orthogonal least squares (OLS) algorithm is guaranteed to exactly recover unknown support in, at most, k iterations, where k denotes the support cardinality [51]. The authors of [51] provide a closer look at the analysis of both orthogonal matching pursuit (OMP) and OLS when the ERC is not fulfilled. The existence of dictionaries for which some subsets are never recovered by OMP is proved. This phenomenon also appears with basis pursuit, where support recovery depends on the sign patterns, but it does not occur for OLS. None of the OMP, OLS, and basis pursuit algorithms is uniformly better than the others, but for correlated dictionaries, OLS may achieve guaranteed exact recovery in fewer iterations compared to OMP.

2.2.3. MIP

The conditioning of the dictionary characterizes how different its atoms are [52]. The performance of a sparse recovery algorithm is affected by the conditioning. The conditioning of a matrix is commonly measured by its condition number, which characterizes the sensitivity of the solution of a system of linear equations to noise.

Mutual coherence and the RIP constant can be seen as two measures of the conditioning of the dictionary. A large mutual coherence corresponds to two similar atoms, implying a bad conditioning in the dictionary and, hence, difficulties in finding the sparse solution.

One common assumption in studying the statistical performance of the estimators is the MIP that requires the mutual incoherence $\mu \left(\mathbf{A}\right)$ to be small [53,54,55]. Mutual coherence [55,56] is defined as the maximum value among all the correlation coefficients of normalized columns of a dictionary $\mathbf{A}$,

$$\mu \left(\mathbf{A}\right)=\underset{i\ne j}{max}|<{\overrightarrow{a}}_i,{\overrightarrow{a}}_j>|,$$

(6)

where ${\overrightarrow{a}}_i$ and ${\overrightarrow{a}}_j$ are two columns of $\mathbf{A}$.

If the $L_0$-norm problem has a k-sparse solution ${\overrightarrow{x}}_0$, for which $k<\frac{1}{2}(1+\mu {\left(\mathbf{A}\right)}^{-1})$, then it is the unique solution for both the $L_1$-norm and $L_0$-norm minimization problems [55,56]. This condition can be replaced by ERC, but ERC is not easy to check because it depends on unknown support. The condition $k<\frac{1}{2}(1+{\mu }^{-1})$ is a sufficient condition for ERC [35] to hold, and is easy to check.

The MIP implies RIP and ERC but the converse is not true. For OMP, support recovery was considered in the noiseless case [35], where the MIP condition $\mu <\frac{1}{2k-1}$ is a sufficient condition for exactly recovering a k-sparse signal $\overrightarrow{x}$ in the noiseless case. This condition is, in fact, sharp [57]. Under the MIP condition $\mu <\frac{1}{2k-1}$ and a condition on the minimum magnitude of the nonzero coordinates of $\overrightarrow{x}$, the support of $\overrightarrow{x}$ can be recovered exactly by the OMP algorithm in the bounded noise cases and with high probability in the Gaussian case [58].

2.3. Sparse Recovery

A high-dimensional k-sparse signal $\overrightarrow{x}\in {\mathbb{R}}^N$ can be expressed as a linear combination of M atoms ($2k\le M\le N$), defined by the sensing matrix $\mathbf{A}={\left[{\overrightarrow{a}}_1,{\overrightarrow{a}}_2,\dots ,{\overrightarrow{a}}_M\right]}^T\in {\mathbb{R}}^{M\times N}$, as given by (3).

Nonlinear reconstruction algorithms are derived from the $L_0$-norm regularized problem,

$$\underset{\overrightarrow{x}\in {\mathbb{R}}^N}{min}\mathsf{\Psi }\left(\overrightarrow{x}\right)=\frac{1}{2}\parallel \mathbf{A}\overrightarrow{x}-\overrightarrow{y}{\parallel }^2\mathrm{subject} \mathrm{to}{\parallel \overrightarrow{x}\parallel }_0\le k.$$

(7)

When $\overrightarrow{y}$ is a k-sparse signal, the problem is known as the k-exact-sparse problem.

The $L_0$-norm regularized problem (4) is non-convex and NP-hard. It has many local minima; when there is zero noise, a unique global minimum is the k-sparse vector ${\overrightarrow{x}}^{*}$ [59]. Many suboptimal algorithms approximate its solution. They recover the true value of ${\overrightarrow{x}}^{*}$ when $\overrightarrow{x}$ is sufficiently sparse and the columns of $\mathbf{A}$ are incoherent.

The $L_0$-norm problem (7) is usually formulated as

$$\underset{\overrightarrow{x}\in {\mathbb{R}}^N}{min}{\parallel \overrightarrow{x}\parallel }_0\mathrm{subject} \mathrm{to}\mathbf{A}\overrightarrow{x}=\overrightarrow{y}.$$

(8)

The $L_0$-norm problem (8) is NP-hard [60]. It is usually relaxed, and is reconstructed via the $L_1$-norm minimization problem [15,61]:

$$\underset{\overrightarrow{x}\in {\mathbb{R}}^N}{min}{\parallel \overrightarrow{x}\parallel }_1\mathrm{subject} \mathrm{to}\mathbf{A}\overrightarrow{x}=\overrightarrow{y}.$$

(9)

This is the basis pursuit method [62]. By imposing appropriate constraints on $\mathbf{A}$, the basis pursuit generates the exact recovery of $\overrightarrow{x}$.

Sparsity is a basic type of regularization. Sparse approximation finds a k-sparse signal ${\overrightarrow{x}}^{*}$ to approximate $\overrightarrow{y}=\mathbf{A}\overrightarrow{x}$ while $k\ll N$. Compressed sensing is a type of sparse approximation problem.

Compared with $L_1$-norm regularization, sparsity is better achieved with $L_0$-norm penalties based on experiments [63]. Under the RIP condition, the solutions by $L_1$-norm and $L_0$-norm regularization are equal [44]. However, $L_1$-norm regularization over-penalizes large coefficients, yielding biased estimation [64,65].

$L_p$-norm ($0<p<1$) is non-convex. Alternatively, $L_p$-norm ($0<p<1$) is an approximation to the $L_0$-norm,

$$\underset{\overrightarrow{x}\in R^N}{min}{\parallel \overrightarrow{x}\parallel }_p^p\mathrm{subject} \mathrm{to}\mathbf{A}\overrightarrow{x}=\overrightarrow{y}.$$

(10)

Regarding the random Gaussian matrix $\mathbf{A}$, the recovering ability of the $L_p$-norm minimization ($p\in [0,1)$) was investigated in [66]. When $\alpha =\frac{M}{N}\to 1$, the sharp threshold of the sparsity ratio differentiates the success and failure via $L_p$-norm minimization. $L_p$-norm minimization succeeds below the threshold. For strong recovery, the threshold decreases strictly from $0.5$ to $0.239$ as p rises from 0 to 1, whereas for weak recovery, the threshold is $2/3$ for $p\in [0,1)$. The threshold is 1 for $L_1$-norm minimization. $L_p$-norm minimization can return a denser solution compared to $L_1$-norm minimization. For any $\alpha \in (0,1)$, thresholds of the sparsity ratio for strong recovery and weak recovery are, respectively, provided in [66]. For strong recovery, $L_p$-norm minimization has a higher threshold with smaller p; for sectional recovery, the threshold is the same for all p; for weak recovery, $L_1$-norm minimization can outperform $L_p$-norm minimization. $L_1$-norm minimization generally outperforms $L_p$-norm minimization for sparse recovery.

Suboptimal signal recovery methods are categorized into greedy pursuit, thresholding, and convex relaxation methods. Greedy pursuit methods, such as matching pursuit [67], OMP [68], OLS [69], subspace pursuit [70], and compressive sampling matching pursuit (CoSaMP) [71], tackle the $L_0$-norm problem directly. Iterative hard thresholding (IHT) is a thresholding method for the $L_0$-norm problem [72]. Convex relaxation methods, such as gradient projection [73,74], accelerated proximal gradient [75], iterative reweighted method [76], homotopy method [77], and least angle regression [78], solve the $L_1$-norm problem, and there are also relaxation methods for the $L_p$-norm problem [79].

Basis pursuit [62] is a greedy sparse approximation technique used for solving the $L_1$-norm problem. For a given set of basis vectors, the method solves the optimization problem by greedily searching for vectors to add or remove. OMP [68], also known as the fully corrective forward greedy selection or simply the forward selection [60], is a simple and effective greedy algorithm for sparse recovery/approximation. The atom selection steps in matching pursuit and the Frank–Wolfe algorithm are very similar [80]. If k is small enough, then at each iteration, both matching pursuit and OMP select an atom indexed by the support, thus ensuring recovery properties [35,81]. Under the condition $k<\frac{1}{2}(1+{\mu }^{-1})$, it has been proven that matching pursuit shows an exponential rate of convergence, and that OMP reaches convergence after exactly k iterations [35]. Under this same condition, it was proved that the Frank–Wolfe algorithm converges exponentially beyond a certain iteration, even though the function is not strongly convex [80].

CoSaMP and subspace pursuit improve upon OMP by selecting multiple coordinates and incorporating a pruning step at each iteration, and sparse signal recovery was based on RIP [28]. This selection strategy results in an optimal sample complexity. At each iteration, CoSaMP first prunes the gradient, then solves an LS program restricted on a small support set; finally, hard thresholding is implemented to form a k-sparse iterate for future updates. As a greedy method, the hard thresholding pursuit adds and prunes indices on a list [82].

Compared to batch solvers, $L_1$-norm-based stochastic algorithms struggle to preserve the sparse structure of the solution [83].

If the optimal solution is sufficiently sparse, the problems (9) and (8) have approximately the same solution [44,48,62,68]. The problem (9) can be effectively solved by LP methods.

For the $L_p$-norm problem, projected gradient [84], iterative reweighted method [85], and approximate operator [86,87] can be used. These methods converge to a global minimum for proper initial points, which is easier to choose for larger p [84].

For compressed sensing, a sparse signal is represented in a finite discrete dictionary. The true parameters, however, may be from a continuous dictionary. Thus, the loss of sparsity arises from spectral leakage along the Dirichlet kernel [88]. From a small random subset of its N time-domain samples, spectrally compressed sensing can recover a spectrally sparse signal. It is assumed that a signal can be represented as the sum of k complex multidimensional sinusoids.

Assuming that the signal frequencies are on a grid, it is guaranteed that the spectrally sparse signal can be faithfully recovered from $O(klogN)$ random time-domain samples by using compressed sensing algorithms derived from $L_1$-norm minimization [44,61], even in the presence of bounded noise [47] or sparse outliers [89]. Additionally, a total-variation norm minimization method can recover a sparse signal from low-frequency samples [90].

In [91], a greedy method for the sparse recovery of linear measurements corrupted by highly impulsive noise is designed to solve a minimum dispersion optimization problem by adopting the family of symmetric alpha-stable distributions.

Sparse recovery algorithms with invariance properties are less affected when the sensing matrix (i.e., the dictionary) is ill-conditioned [52]. There implicitly exists an equivalent well-conditioned problem. Some sparse recovery algorithms, such as smoothed $L_0$-norm [92], basis pursuit, FOCUSS [93], and hard thresholding algorithms, are invariant, while others, such as matching pursuit and the spectral projected gradient for $L_1$-norm minimization (SPGL1) [94], are not.

2.4. Iterative Hard Thresholding

A simple greedy technique known as IHT [72,95] can generate the steepest descent steps that are feasible for the $L_1$-norm problem (9). This is achieved by utilizing hard thresholding to project steps along the negative gradient direction of $\mathsf{\Psi }$ onto the $L_0$-norm constraint. In the hard thresholding method, all but the k largest magnitude elements of a vector are set to zero. IHT utilizes a proximal-point technique at each iteration, and it implements gradient projection with a constant step size. In normalized IHT [96], a self-adaptive step size is adopted to guarantee stability and performance.

In case the spectral norm of $\mathbf{A}$ is less than one, IHT converges to a fixed point of (9) [72]. IHT guarantees stable recovery, provided that $\mathbf{A}$ has RIP [95]. There are also other RIP-based recovery conditions for IHT [82] and normalized IHT [96]. The theory for compressed sensing assumes that samples are taken from linear measurements. Under similar conditions, IHT accurately recovers sparse or structured signals from a few nonlinear observations [97]. Global linear convergence for IHT is guaranteed from theory works based on matrix completion, which is based on standard properties of incoherence and uniform sampling.

For any $\mathbf{A}$, sufficient conditions for the convergence of IHT to a fixed point, as well as necessary conditions for the existence of fixed points, are given in [98]. Sparse signal recovery analysis can be performed using these conditions. The analysis has been extended to normalized IHT. A theoretical analysis of normalized IHT, when both $\mathbf{A}$ and $\overrightarrow{y}$ are quantized, is given in [99], and it is proved that a low-precision normalized IHT can provide recovery guarantees under mild conditions.

IHT [72] and iterative soft-thresholding [100] algorithms are for the $L_0$- and $L_1$-norm problems, respectively. Half thresholding [86] is given for $p=\frac{1}{2}$. These iterative thresholding algorithms are efficient for high-dimensional problems, and it is also relatively easy to specify the regularization parameter.

In the proximal gradient homotopy method, the solutions of the regularized problem are computed and traced along a continuous homotopy path, and the selection of a regularization parameter is not needed. IHT, when combined with the homotopy technique, avoids the requirement of choosing a regularization parameter [101].

The hard thresholding pursuit [82] is an iterative greedy selection procedure used for the sparse recovery of the $L_0$-norm problem. As a combination of CoSaMP and IHT, it outperforms both methods in terms of the RIP. The exact recovery of sparse signals is theoretically justified under conditions of restricted strong condition number bounding [102]. In [103], the hard thresholding pursuit is generalized to sparsity-constrained convex optimization. The algorithm includes iterations of a gradient descent step and a hard thresholding step.

Iterative thresholding algorithms have much lower computational complexity per iteration and lower storage requirements than interior point methods. The recursions are modifications of the gradient method used to solve a linear system but consist of a (hard or soft) shrinkage operator to promote the sparsity of the estimate at each iteration. Hard thresholding algorithms are always orders of magnitude faster than convex programs [104].

A theoretical analysis of hard thresholding algorithms is given in [105]. A tight bound used for characterizing hard thresholding algorithms is derived, and RIP and sparsity parameters are related. Parsimonious solutions are guaranteed by following a stochastic hard thresholding procedure. Global linear convergence is proved under certain mild assumptions [105].

IHT is a first-order greedy selection method that minimizes the primal formulation. The original non-convex problem can be equivalently or approximately solved in a concave dual formulation under certain conditions. The dual IHT algorithm [106] is a super-gradient ascent method used to solve the non-smooth dual problem. Dual IHT is superior to IHT in model estimation accuracy and computational efficiency.

The GradMP algorithm [107] generalizes the idea of CoSaMP [71]. Stochastic IHT and stochastic GradMP [108] are two stochastic variants of greedy algorithms. The expected linear convergence towards the solution within a specified tolerance is proven, providing methods that often outperform their deterministic counterparts. In stochastic GradMP [108], at each iteration, only the gradient of a function is evaluated, a subspace is searched based on the previous estimate, and then a new solution is obtained via solving a convex low-dimensional sub-optimization problem.

Sparse convex optimization solves the optimization of a convex function subject to a sparsity constraint $k\le k^{*}\gamma$, where $k^{*}$ is a target sparsity and $\gamma \ge 1$ is an approximation factor. The adaptively regularized hard thresholding (ARHT) algorithm [109] brings the bound down to $\gamma =O\left(\kappa \right)$, with $\kappa$ being the restricted condition number, which is tight for a general class of algorithms, including IHT, OMP, and LASSO. ARHT is comparable to the most efficient greedy algorithms in terms of runtime, as it requires a single function minimization per iteration. ARHT provides a strong trade-off between the RIP condition and the solution sparsity.

Newton step-based IHT and Newton step-based hard thresholding pursuit adopt the Newton-like search direction instead of the steepest descent direction [110]. Sufficient guarantees for these algorithms are established in terms of the RIP of a sensing matrix.

A class of distributed iterative thresholding algorithms for $L_0/L_1$-regularized LS optimization problems is presented in [111]. By introducing a suitably distributed and regularized LS functional, it has been demonstrated that the algorithms reach their minima based on the dynamical systems theory [111].

2.5. Orthogonal Matching Pursuit

OMP [68], also known as forward stepwise regression, is a fast greedy algorithm. At each iteration, the algorithm selects the column from matrix $\mathbf{A}$ that is maximally correlated with the current residual and adds it to the set of selected columns. That is, at each iteration, one new element of the dictionary is added and one orthogonal projection is made. The residuals are updated by projecting the observations $\overrightarrow{y}$ onto the linear subspace that is spanned by the columns already selected, and then the algorithm iterates. The stopping rule depends on the noise structure.

OMP can recover a k-sparse signal $\overrightarrow{x}\in {\mathbb{R}}^N$ from incomplete measurements $\overrightarrow{y}\in {\mathbb{R}}^M$ obeying (3), with $k\ll M\ll N$. OMP recovers the true signal with high probability for random matrices, including Gaussian, but it may fail for some deterministic sensing matrices [35,112].

OMP, such as OLS, constructs the dictionary in an incremental way. By adding one index to the list at a time, the support of the underlying sparse signal is identified, and the sparse coefficients over the enlarged support are estimated. In OLS, a candidate that leads to the most significant decrease in residual power is selected. In comparison, in OMP, a column that is the most strongly correlated to the residual is chosen. OLS outperforms OMP in terms of convergence, at a cost of higher computational complexity [51].

Some greedy methods, such as stagewise OMP [113], regularized OMP [114], and generalized OMP [115] (also referred to as the orthogonal super greedy algorithm [116]) add multiple indices per iteration. At each iteration, candidates are identified according to correlations between the residual vector and columns of $\mathbf{A}$. The multipath matching pursuit [117] extends OMP by recovering sparse signals with a tree-searching strategy. At each iteration, multiple candidate paths are traced and extended, and the candidate that minimizes the residual power is chosen. Multiple OLS [118] extends OLS by selecting multiple indices at each iteration, leading to convergence in fewer iterations. Stable sparse recovery can be guaranteed, provided that the signal-to-noise ratio (SNR) grows linearly with the sparsity level k of the input signals.

When MIP or SNR satisfies certain conditions, OLS and multiple OLS methods reliably recover k-sparse signals in, at most, k iterations, while block OLS succeeds in, at most, $k/d$ iterations, where d is the block length [119]. The theoretical analysis for block OLS utilizes the block-MIP to deal with block sparsity. An MIP-based theoretical analysis shows OLS and multiple OLS methods in [119].

Signal space matching pursuit [120] sequentially adjusts the support of jointly sparse vectors to minimize the subspace distance to the residual space. The method accurately reconstructs any row k-sparse matrix of rank r in the full row rank scenario when the maximum number of linearly independent columns in $\mathbf{A}$ is not less than $k+1$. The selection rule reduces to that of multiple OLS when $r=1$, and to that of OLS when $r=L=1$, where L is the number of indices chosen in each iteration.

OMP with replacement [121] is known as partial hard thresholding with parameter $r=1$ [122], which is a generalization of IHT. It is essentially a variant of OLS that includes replacement steps. Block OMP [123] recovers block sparse signals whose nonzero entries occur in a few blocks by assuming a uniform block size and known block boundaries. In [124], block OMP is implemented by using coarse-fine block localization of a nonzero cluster. OMP is used to reconstruct a class of structured sparse signals modeled by trigonometric polynomials in [125].

Regarding sparse approximation, a single sufficient condition, where both basis pursuit and OMP could reliably recover a sparse signal, was developed in [35]. OMP can faithfully recover a k-sparse signal $\overrightarrow{x}\in {\mathbb{R}}^N$ given $O(klnN)$ random linear measurements [112]. OMP yields exact support recovery under certain RIP assumptions [126], and several improvements to the condition were proposed [116,127].

2.6. LASSO

LASSO [30,128], originally proposed for estimation in the linear model $\overrightarrow{y}=\mathbf{A}\overrightarrow{x}+\overrightarrow{n}$, has become a popular supervised learning technique for the recovery of sparse signals from high-dimensional measurements and an unsupervised learning technique for the feature selection of high-dimensional samples. The $L_1$-norm regularizer in LASSO tends to generate sparse regression coefficients.

In the context of supervised learning, the LASSO formulation is equivalent to the SVM formulation [129]. For unsupervised learning, the LASSO regression has been applied in biclustering tasks [130].

LASSO minimizes the sum of squared errors, subject to a bound on the sum of the modulus of the regression coefficients. It can be formulated as an $L_1$-norm regularized LS problem. The convex optimization problem, (4) or (9), can be represented by an LS problem subject to $L_1$-norm penalty, which has the same formulation as LASSO [30]

$$\overrightarrow{x}={argmin}_{\tilde{\overrightarrow{x}}\in {\mathbb{R}}^N}\{\parallel \mathbf{A}\tilde{\overrightarrow{x}}-\overrightarrow{y}{\parallel }_2^2+\lambda \parallel \tilde{\overrightarrow{x}}{\parallel }_1\},$$

(11)

where $\overrightarrow{x}\in {\mathbb{R}}^N$ is a regression coefficient vector, and $\lambda >0$ is a regularization parameter. There are many methods used for solving (11), such as stochastic gradient descent and stochastic coordinate descent [83,131]. Software packages for LASSO are publicly available.

The LASSO estimator satisfies the well-known prediction bound

$$\lambda \ge \frac{2\parallel {\mathbf{A}}^T\overrightarrow{n}{\parallel }_{\infty }}{N}⟹\frac{1}{N}\parallel \mathbf{A}(\overrightarrow{x}-{\widehat{\overrightarrow{x}}}_{\lambda }){\parallel }_2^2\le 2\lambda {\parallel \overrightarrow{x}\parallel }_1.$$

(12)

We call $\frac{2\parallel {\mathbf{A}}^T\overrightarrow{n}{\parallel }_{\infty }}{N}$ LASSO’s effective noise. Such bounds are referred to as oracle inequalities. The effective noise plays an important role in finite-sample bounds for LASSO, the calibration of LASSO’s tuning parameter, and inference on the coefficient vector $\overrightarrow{x}$. A bootstrap-based estimator of the quantiles of the effective noise was developed in [132]. The estimator is fully data-driven, i.e., it does not need any additional tuning parameters. The estimator is equipped with finite-sample guarantees and is applied to the calibration of tuning parameters for LASSO as well as to high-dimensional inference $\overrightarrow{x}$.

To estimates the regression vector $\overrightarrow{x}$ in the generic linear model with $\overrightarrow{n}=\mathcal{N}(0,{\sigma }^2\mathbf{I})$, when the variance ${\sigma }^2$ is unknown, two LASSO-type methods that jointly estimate $\overrightarrow{x}$ and the variance are minimizers of the $L_1$-norm-penalized LS functional, where the relaxation parameter is tuned according to two strategies [133].

LASSO implicitly performs model selection and shares many connections with forward stepwise regression. The least angle regression [78] performs stepwise variable selection. At each iteration, the variable that is correlated with all of the residuals obtained thus far the most is put in the set of active variables, and the current update is in a direction that is equiangular with all other active variables. Unlike OMP, which maintains a variable permanently, the least angle regression continually modifies the coefficient of the most correlated variable until that variable is no longer the one that is most correlated with the recent residual. The entire LASSO regularization path is generated with a computational cost that is similar to that of standard LS via QR decomposition.

In order to handle nonlinearity, instance-wise nonlinear LASSO [134] applies a nonlinear function on an instance $\overrightarrow{x}$ to give a sparse solution, in terms of instances. On the other hand, the feature-wise nonlinear LASSO, also referred to as the feature vector machine [135], imposes a nonlinear transformation ‘feature-wisely’ to obtain sparsity in terms of features. Both methods use the kernel trick.

For large-scale LASSO regression problems, the Frank–Wolfe method [136] uses the randomized iteration, and it is superior to the coordinate descent method. It achieves a convergence rate of $O(1/k)$ (in terms of the expected value). The solutions are significantly more sparse compared with the competing methods while retaining the same accuracy.

Sparsity-inducing algorithms, such as LASSO, are not algorithmically stable [137]. To put it differently, each iteration of the leave-one-out cross-validation of the LASSO estimator may—each time—produce disparate results. The tuning parameter and the model have to be estimated separately by using the data twice. The LASSO estimator can be risk-consistent when a tuning parameter is chosen through cross-validation under certain restrictions [138]. For LASSO, the robust optimization formulation is related to kernel density estimation [139]. According to the no-free-lunch theorem, sparsity and algorithmic stability are contradictory requirements, thus LASSO is not stable [139]. Compared with the unbounded asymptotic variance of the LASSO estimator, some robust LASSO estimators have stabilized asymptotic variances in the presence of large variance noise [140].

Group LASSO [141] selects variables at the group level. The penalty is in an intermediate mode between the $L_1$-norm and $L_2$-norm penalties. Group LASSO minimizes the square loss plus a penalty term proportional to the sum of the Euclidean norms of groups of coefficients. Group square-root LASSO [142] minimizes the square root of the residual sum of squares plus the same penalty term for group LASSO. It is independent of the variance of the error terms. Square-root LASSO, with or without groups, achieves the same correct pattern recovery and prediction accuracy under similar conditions, but with a simplified tuning strategy, compared to the LASSO or group-LASSO methods. Group square-root LASSO, with proven convergence properties, scales well with the dimension of the problem.

LASSO belongs to a family of regularized linear regression methods, which also includes ridge regression [143] and elastic net [144]. LASSO is an $L_1$-regularized LS method. Ridge regression substitutes the $L_1$-norm by the squared $L_2$-norm ridge regularization on the coefficients. LASSO not only reduces the variance of coefficient estimates but also selects variables by setting those coefficients below a threshold to zero. elastic net regularization uses a linearly mixed penalty of $L_1$- and $L_2$-norms [144]. By regressing each dependent variable separately on each covariate, marginal regression is roughly two orders of magnitude faster than LASSO for sparse and high-dimensional regression problems [145].

2.7. Other Sparse Algorithms

PCA is a classic method, and we do not describe it in this paper. Sparse PCA is targeted to find a sparse basis in order to make the result easy to interpret. A trade-off needs to be made between statistical fidelity and interpretability.

The orthogonality of loadings is considered in the simplified component technique-LASSO (SCoTLASS) [146], loading rotation [147], simple thresholding [148], and augmented Lagrangian sparse PCA [149]. SCoTLASS [146] optimizes the objective function of PCA subject to a sparsity constraint on each loading. The loading rotation [147] rotates the PCA loadings using various criteria in order to find a simple structure. Simple thresholding [148] obtains sparse loadings by setting PCA loadings to zero. Augmented Lagrangian sparse PCA [149] solves an augmented Lagrangian optimization problem, where the explained variance, orthogonality, and correlation between principal components are simultaneously considered.

Examples of deflation methods are the greedy methods [150], SCoTLASS, rSVD [151,152], GPower [153], and TPower [154]. Greedy search and branch-and-bound methods can solve small problems exactly, but with the complexity of $O\left(N^4\right)$ [150]. PathSPCA [151] is an approximate alternative to the solution of [150], leading to a reduced complexity of $O\left(N^3\right)$. rSVD [152] solves a sequence of rank-1 matrix approximations, subject to a sparsity penalty, to obtain sparse loadings. GPower [153] maximizes a convex objective and solves it by the power method. TPower [154] and a related power method, referred to as iterative thresholding sparse PCA [155], are targeted at the recovery of the sparse principal subspace.

In [156], sparse PCA is formulated as a regression-type optimization so as to use LASSO or elastic net techniques. Direct sparse PCA [157] relaxes the problem into a semidefinite convex problem, which has a computational complexity of $O\left(N^4{(logN)}^{1/2}\right)$ for N variables. A variable elimination method [158] reduces the complexity to $O\left(N^3\right)$. A methodology for uncertainty quantification is proposed in [159] based on an M-estimator with the LASSO penalty. It achieves minimax optimal rates and is used to construct a de-biased sparse PCA estimator. The estimator has a Gaussian limiting distribution and can be used for hypothesis testing or support recovery of the first eigenvector. It outperforms PCA in moderately high-dimensional regimes.

The sparse LMS algorithm [160] penalizes the quadratic cost function of the LMS algorithm by two sparsity constraints. Recursive $L_1$-regularized LS [161] estimates a sparse tap-weight vector for adaptive filtering by using an EM-type algorithm. The method outperforms the RLS algorithm in terms of both MSE and computational complexity.

Sparse SVD [162] is based on iterative thresholding of singular vectors, and is robust to tuning parameters. The penalized matrix decomposition [163] penalizes the likelihood with the $L_1$-norm penalty on factors and/or loadings. softImpute [164] fits a regularized low-rank matrix using a nuclear norm penalty.

The sparse factor analysis [165] and nonparametric Bayesian sparse factor analysis [166] are Bayesian approaches with different prior specifications. These Bayesian methods are self-tuning. A general empirical Bayes approach to matrix factorization [167] estimates the sparsity by estimating prior distributions from the observed data, and uses a variational approximation to effectively solve a simpler so-called normal means problem.

2.8. Restricted Isometry Property for Signal Recovery Methods

For the linear regression problem, when ${\delta }_{2k}<\sqrt{2}-1\approx 0.41$, the $L_0$-norm and $L_1$-norm problems are equivalent [39]. The LASSO algorithm can recover a solution [15,28,39,47]. In [168], the condition improved to ${\delta }_{2k}<0.493$. Many of the later results either provided related guarantees for LASSO while improving the RIP upper bound [168,169,170], reaching a bound of ${\delta }_{2k}<0.6248$, or obtained similar results by using greedy algorithms under more strict RIP conditions, but typically converging faster than LASSO [71,82,95,114,121,171].

For linear regression, CoSaMP [71] achieves a bound that is similar to that in [39], but the implementation is more efficient. Their method is valid for the more restricted RIP upper bound of ${\delta }_{2k}<0.025$, or ${\delta }_{4k}<0.4782$, as improved by [172]. IHT achieves a bound similar to that of CoSaMP [95], with the condition ${\delta }_{3k}<0.067$, which is improved to ${\delta }_{2k}<\frac{1}{3}$ by [121] and to ${\delta }_{3k}<0.5774$ by [82].

We consider sufficient conditions for perfect signal recovery using OMP. In the noiseless case, OMP can exactly identify the support of a k-sparse signal in k iterations, provided that $\mathbf{A}$ satisfies the $k+1$th-order RIP with ${\delta }_{k+1}<\frac{1}{3\sqrt{k}}$ [126]. Since a smaller RIC leads to better reconstruction, sparse signal recovery with interference-nulling achieves better performance than what is predicted in [173]. The sufficient condition is relaxed to ${\delta }_{k+1}<\frac{1}{\sqrt{k}+1}$ [115,127,174].

Sufficient conditions specified with the RIC bound ${\delta }_{k+1}<\frac{1}{\sqrt{k}+1}$ and certain requirements on the minimal magnitude of signal entries guarantee exact support identification under measurement noise [175]. In the noisy case, a relaxed upper bound ${\delta }_{k+1}<\frac{\sqrt{4k+1}-1}{2k}$ as well, as relaxed requirements on the minimal magnitude of the signal entries, guarantee perfect support recovery by using OMP [176]. In the noiseless case, the relaxed bound guarantees exact support recovery in k iterations.

If $\mathbf{A}$ satisfies RIP with ${\delta }_{k+1}<\frac{1}{\sqrt{k+1}}$, then OMP faithfully recovers a k-sparse signal $\overrightarrow{x}$ in k iterations under constraints on the minimum magnitude of nonzero entries of $\overrightarrow{x}$ [177,178]. This sufficient condition on ${\delta }_{k+1}$ is sharp. If $\mathbf{A}$ satisfies RIP with ${\delta }_{k+1}<\frac{1}{\sqrt{k+1}}$, then OLS also exactly recovers $\overrightarrow{x}$ in k iterations [179].

In [180], an OMP-like algorithm is analyzed based on RIP. Based on the technique in [180], sparse approximation by greedy algorithms is studied in [181]. With high probability, the exact recovery of random k-sparse signals within $k(1+\epsilon )$ iterations of OMP is proved. Thus, OMP is almost optimal for the exact recovery in a probabilistic sense [181].

It has been proven that for a k-sparse $\overrightarrow{x}$ and matrices with a rank of at most k, if the RIC of $\mathbf{A}$, ${\delta }_{tk}\le \sqrt{\frac{t-1}{t}}$, where $t\ge \frac{4}{3}$, then the $L_1$-norm problem can recover $\overrightarrow{x}$ exactly in the noiseless case and stably in the noisy case [182]. ${\delta }_{tk}<\frac{t}{4-t}$ was connected to be a sharp condition for $0<t<\frac{4}{3}$, and it was only partially proved in [182]. The conjecture on the RIP constant ${\delta }_{tk}<\frac{t}{4-t}$$(0<t<\frac{4}{3})$ is completely proven in [183]. Thus, in the noiseless case, a complete characterization of sharp RIP constants ${\delta }_{tk}$ for all $t>0$ is obtained, ensuring the exact recovery of all k-sparse signals and matrices with a rank of at most k through $L_1$-norm minimization and nuclear norm minimization, respectively. Noisy cases and approximately sparse cases are also considered.

Multiple OLS ($L>1$) [118] recovers k-sparse signals faithfully in, at most, k iterations, if $\mathbf{A}$ obeys the RIP with ${\delta }_{Lk}<\frac{\sqrt{L}}{\sqrt{k}+2\sqrt{L}}$. OLS ($L=1$) guarantees exact recovery under ${\delta }_{k+1}<\frac{1}{\sqrt{k}+2}$. This bound is tight since even a slight relaxation disables OLS from guaranteeing exact recovery.

Multipath matching pursuit faithfully recovers all k-sparse signals, provided that $\mathbf{A}$ satisfies the $k+L$-order RIP with ${\delta }_{k+L}<\frac{\sqrt{L}}{\sqrt{k}+2\sqrt{L}}$, In the case of L child paths per candidate [117]. This bound is further improved to ${\delta }_{k+L}<\sqrt{\frac{L}{k+L}}$ [184].

For the subspace pursuit, RIP-based exact recovery guarantees in both noiseless and noisy cases are given in [70,185]. For block OMP, block RIP is used to derive some sufficient conditions for the exact or stable recovery of block sparse signals in [186]. In [124], the convergence of the coarse-fine block OMP is analyzed by defining a pseudoblock-interleaved block RIP and then imposing upper bounds on the corresponding RIC.

Signal space matching pursuit guarantees exact reconstruction in at most $k-r+\lceil \frac{r}{L}\rceil$ iterations, if $\mathbf{A}$ satisfies the RIP of order $L(k-r)+r+1$ with ${\delta }_{L(k-r)+r+1}<$$max\left\{\frac{\sqrt{r}}{\sqrt{k+\frac{r}{4}}+\sqrt{\frac{r}{4}}},\frac{\sqrt{L}}{\sqrt{k}+1.15\sqrt{L}}\right\}$ [120]. The RIC requirement becomes less restrictive as r increases, and is less restrictive than those for OLS and multiple OLS. In case of $r=1$ and more than k iterations, the performance guarantee can be improved to ${\delta }_{\lfloor 7.8k\rfloor }\le 0.155$. Under a suitable RIP condition, the reconstruction error is upper bounded by a constant multiple of the noise power [120].

The $L_p$-norm problem is investigated based on RIP [187,188]. For $0<p\le 1$, any k-sparse signal can be recovered if ${\delta }_{2k}<\delta \left(p\right)$, with $\delta \left(p\right)>0$ being a constant decided by p [188]. Sufficient conditions for the exact recovery were derived in terms of the RIC for $L_p$-norm minimization [169], CoSaMP [71], and regularized OMP [114].

In addition to RIP, notions such as the restricted orthogonality constant (ROC) [28] and null space property [189] have been used for the analysis of sparse recovery. The $L_1$-norm and $L_p$-norm problems have been studied by using the null space property [189,190,191,192]. A null space constant of less than 1 is a sufficient and necessary condition for guaranteeing the $L_p$-norm problem to exactly recover any k-sparse signal [189]. Based on the null space property, if the $L_p$-norm problem ($p=p_1$) can recover any k-sparse signal, then the $L_p$-norm problem ($p\le p_1$) can also work [190]. $L_p$-minimization with a sufficiently small p is equivalent to $L_0$-minimization for sparse recovery [191]. The $L_p$-norm problem ($p<p^{*}$), with $p^{*}$ from an upper bound on the null space constant, is guaranteed for exact recovery [192].

2.9. Related Topics

One-bit compressed sensing [193] adopts the compressed sensing model, but only the sign of each measurement is retained. k-sparse signals in ${\mathbb{R}}^N$ can be estimated (up to normalization) from $\mathsf{\Omega }(klog\frac{N}{k})$ one-bit measurements. Recovery algorithms can be based on nonlinear programming [194], linear programming [195], convex programming [196], and modifications of IHT [197]. A uniform $L_2$-reconstruction error of at most $\gamma >0$ can be achieved with $M\ge \frac{1}{\gamma }klog\frac{N}{k}$ one-bit measurements [196,197]. The optimal quantization scheme is obtained with respect to the mean square error (MSE) of the LASSO reconstruction [198]. In [199], the decay of the error is optimized as a function of the oversampling factor $\lambda =\frac{M}{klog\frac{N}{k}}$. The error in reconstructed signals from one-bit measurements is bounded below by $\mathsf{\Omega }\left(\frac{1}{\lambda }\right)$. The adaptive thresholding used for quantization can lower the error rate to $e^{-\mathsf{\Omega }\left(\lambda \right)}$, which improves upon other adaptive thresholding methods, such as the sigma-delta quantization. A general recursive strategy achieves this exponential decay, realized by two specific polynomial-time algorithms, one based on convex programming and one on hard thresholding.

For a deterministic finite alphabet vector $\overrightarrow{x}$, two convex optimization methods, namely, the regularization-based method and transform method, have been introduced for the recovery of finite alphabet signals via $L_1$-norm minimization [200]. When the alphabet sizes $p=2$ and $(M,N)$ grow proportionally, the conditions for high-probability signal recovery are the same for both methods.

Without prior knowledge of the sparsity basis in both the sampling and recovery processes, blind compressed sensing is ill-posed in general [5]. Some constraints on the sparsity basis can be added to guarantee a unique solution. The methods can achieve results similar to those of standard compressed sensing, as long as the signals are sparse enough.

3. Dictionary Learning

Sparse coding, also referred to as dictionary learning, represents a dense signal using only a few elements from an overcomplete dictionary [10]. Dictionary learning is targeted to recover the elementary signals (atoms, exemplars, words), collectively referred to as a dictionary, which efficiently represents a set of homogeneous signals. This is generally performed by imposing certain sparseness constraints on the representative coefficients. Dictionary learning is usually used to find a sparse, patch-level representation of an image [201]. It is useful in image de-noising.

3.1. Problem Formulation

Sparse approximation has a formulation similar to that of compressed sensing but with a different objective. A target signal $\overrightarrow{y}\in {\mathbb{R}}^M$ is represented by a linear combination of atoms in an overcomplete dictionary $\mathbf{A}\in {\mathbb{R}}^{M\times N}$,

$$\overrightarrow{y}=\mathbf{A}\overrightarrow{x},$$

(13)

where the basis matrix $\mathbf{A}=\left[{\overrightarrow{a}}_1,{\overrightarrow{a}}_2,\dots ,{\overrightarrow{a}}_N\right]$ ($N\gg M$), ${\overrightarrow{a}}_i\in {\mathbb{R}}^M$ represents a word with the unit norm, $\parallel {\overrightarrow{a}}_i{\parallel }_2=1$, and $\overrightarrow{x}\in {\mathbb{R}}^N$ is a representation of $\overrightarrow{y}$. Any vector $\overrightarrow{y}$ can be represented as a linear combination of words in the overcomplete dictionary.

The recovery of both $\mathbf{A}$ and $\overrightarrow{x}$ given $\overrightarrow{y}$ is an underdetermined problem. For random and sparse $\overrightarrow{x}$, both $\mathbf{A}$ and $\overrightarrow{x}$ can be recovered from $\overrightarrow{y}$ with a high probability for sufficiently large M [202]. A polynomial-time algorithm, referred to as the exact recovery of sparsely-used dictionaries (ER-SpUD), consists of an ER-SpUD step and a greedy step. ER-SpUD is proved to probably recover the dictionary and coefficient matrices for the sufficiently sparse coefficient matrix [202]. The method is valid for $M\ge C{N}^2{log}^{2}N$, $C>0$ is a constant, and it was conjectured that $M\ge CNlogN$ suffices from an information-theoretical view [202]. This bound improves to $M\ge CN{log}^{4}N$ [203]. In [204], an improved Er-SpUD algorithm faithfully recovers $\mathbf{A}$ and $\overrightarrow{x}$ with high probability when $M\ge CNlogN$.

The set of linear Equations (13) has no unique solution. A sufficiently sparse $\overrightarrow{x}$ can be uniquely obtained by solving the $L_0$-norm minimization problem given by (8) [15]. Under weak conditions of $\mathbf{A}$, the $L_0$-norm minimization problem (8) has a solution equal to that of the $L_1$-norm minimization problem (9) [44]. Solutions to the $L_0$-norm and $L_1$-norm problems have been discussed and solved in the previous section. As such, the problem of recovering sparse signals from compressed measurements is the same as that of constructing sparse approximation.

3.2. Dictionary Learning Methods

Some examples of dictionary learning methods are sparse coding [205], nonnegative sparse coding [206,207], $L_p$-sparse coding [208], K-SVD [201], hierarchical sparse coding [209], fused-LASSO-based dictionary learning [210], and elastic net-based dictionary learning [211].

The leading dictionary learning methods are convex optimization algorithms, such as the alternating direction method of multipliers (ADMM), which was used for solving (8), and greedy algorithms, such as matching pursuit [67] and OMP [68]), which were also used for solving (9).

ADMM alternatively minimizes the coefficients and atoms separately. The K-SVD method [201] is a popular ADMM algorithm used for solving $L_0$-norm-based problems. It sequentially updates atoms in the dictionary $\mathbf{A}$ by SVD and finds sparse coefficients $\overrightarrow{x}$ by OMP by alternating iterations. However, the computational cost of OMP is nontrivial, and K-SVD is not always convergent.

Proximal alternating methods [212,213] can be used for a class of non-convex optimization problems, leading to global convergence. Accelerated plain dictionary learning [214] is a multi-block alternating scheme for $L_0$-norm sparse coding, with global convergence. A multi-block hybrid proximal alternating scheme [214] combines ideas from multi-block coordinate descent, proximal alternating methods, and the K-SVD method.

Motivated by the K-SVD algorithm, dictionary learning over positive definite matrices is solved by the alternating minimization approach [215]. Coordinate descent is implemented, and it is much faster than generic interior point methods.

For the recovery of sparse signals represented by a general dictionary that is corrupted by additive noise, the derived deterministic recovery guarantees depend on the signal and noise sparsity levels, on the coherence parameters of the involved dictionaries, and on the amount of prior knowledge about the signal and noise support sets [216]. When both signal and noise are sparse but in different domains, signal recovery is a non-convex and NP-hard problem. In [217], the problem is solved either by replacing $L_0$-norm with $L_1$-norm and then applying ADMM, or replacing $L_0$-norm with a smoothed $L_0$-norm and then applying the gradient projection method.

Sparse coding with latent variables described by discrete prior distributions was investigated in [218]. The sparse latent variables can take a value from a finite set of values and the prior probability of any value is learned from the data. Discrete sparse coding algorithms can scale efficiently with datasets.

Assume that the signals are generated as i.i.d.-random linear combinations of the K atoms from a complete reference dictionary ${\mathbf{D}}^{*}\in {\mathbb{R}}^{K\times K}$, where the linear combination coefficients are from either a Bernoulli-type model or an exact sparse model. A necessary and sufficient norm condition for ${\mathbf{D}}^{*}$ to be the unique sharp local minimum of the expected $L_1$-norm objective function is obtained, thus establishing the global property of $L_1$-norm dictionary learning [219]. The algorithm is based on block coordinate descent, guaranteeing a monotonic decrease in the objective function.

By casting dictionary learning as a classical (or frequentist) estimation problem, lower bounds on the worst-case MSE are derived by applying different generative models for the observed signals [220]. A lower bound on the worst-case MSE in terms of the SNR was obtained. The lower bounds are used to derive the required number of observations, such that dictionary learning is feasible.

Many signals cannot be sparsely represented using an orthonormal basis, but have sparse representations in a redundant dictionary $\mathbf{D}$. Standard compressive sensing methods can be extended to handle this case, provided that the dictionary is sufficiently incoherent or well conditioned, but fail in the case of a truly redundant or overcomplete dictionary [221]. The projected Landweber algorithm [222] extends IHT [95], and signal-space CoSaMP [223] extends CoSaMP, in order to operate in the signal space. They are oriented to recover the signal rather than its dictionary coefficients. $\mathbf{D}$-RIP [221] is a condition on the sensing matrix analogous to RIP. Both works assume that $\mathbf{A}$ satisfies the $\mathbf{D}$-RIP, which is a less-restrictive condition to satisfy than requiring $\mathbf{AD}$ to satisfy RIP. Implementing both algorithms requires the ability to compute projections of vectors in the signal space onto a sparse representation in the model family.

A non-convex generalization of the online matrix factorization algorithm for the i.i.d. data stream [224] is demonstrated to converge almost surely in the network dictionary learning algorithm [225]. A network dictionary learning algorithm [225] combines the online NMF and an MCMC algorithm for sampling motifs from networks. It extracts network dictionary patches from a given network. The convergence guarantee of the network dictionary learning algorithm is given.

The uniform spread of information is enforced over representative coefficients in robust encoding for digital communications. The $L_{\infty }$-norm penalty is used to naturally express anti-sparse regularization. A fully Bayesian formulation of anti-sparse coding is derived by using a prior known as the democratic prior to enhance anti-sparsity in a Gaussian linear model [226].

4. Matrix Completion

Matrix completion is a special case of the more general matrix recovery problem, which aims to reconstruct a matrix from generic and often random linear measurements.

Let $\mathbf{Y}$ represent the measured data that are corrupted by errors $\mathbf{E}$. The task is to recover a low-rank matrix $\mathbf{X}\in {\mathbb{R}}^{m_1\times n_1}$ from $\mathbf{Y}=F\left(\mathbf{X}\right)+\mathbf{E}\in {\mathbb{R}}^{m\times n}$ with the linear operator $F:{\mathbb{R}}^{m_1\times n_1}\to {\mathbb{R}}^{m\times n}$,

$$\underset{\mathbf{X},\mathbf{E}}{min}\mathrm{rank}\left(\mathbf{X}\right)+\lambda \parallel \mathbf{E}\parallel \mathrm{subject} \mathrm{to}\mathbf{Y}=F\left(\mathbf{X}\right)+\mathbf{E},$$

(14)

where $\lambda$ is a regularization parameter, and $\parallel ·\parallel$ is the $L_0$-norm [17] or $L_{2.0}$-norm [227] for sparsity.

When $F(.)$ is an identity operator, the model (14) is used for the low-rank and sparse matrix decomposition [17]. When $F\left(\mathbf{X}\right)=\mathbf{A}\mathbf{X}$ and $\mathbf{A}$ is a dictionary, it pertains to the low-rank representation [227,228]. When $F(.)$ is a sampling operator, it pertains to the low-rank matrix completion [229]. The matrix completion is a special case of the matrix recovery problem, which aims to recover a matrix from generic, random linear measurements. The problem (14) is NP-hard, owing to the discrete and non-convex nature of the rank function and $L_0$-norm (or $L_{2,0}$-norm) [229].

The rank function counts the number of nonzero singular values. A popular convex relaxation of the rank functional is the nuclear norm, also referred to as the trace norm, which is the sum of all singular values of a matrix [230]. The nuclear norm is the convex envelope of the rank function, and is the tightest convex lower bound of the rank function of a matrix [231]. Since the nuclear norm is a convex function, it can be efficiently optimized by semidefinite programming. Similar to the nuclear norm for the rank function, a convex relaxation of $L_0$-norm (or $L_{2,0}$-norm) is $L_1$-norm (or $L_{2,1}$-norm). The $L_{\infty }$-norm is a convex relaxation of the rank function for matrix completion under a uniform sampling distribution [232]. The rank function can also be relaxed by the Schatten p-norm.

Given an incomplete low-rank data matrix $\mathbf{X}=\left[X_{ij}\right]\in {\mathbb{R}}^{m\times n}$, the matrix completion problem can be formulated as follows:

$$min\mathrm{rank}\left(\mathbf{Y}\right)\mathrm{subject} \mathrm{to}{Y}_{ij}=X_{ij},(i,j)\in \mathsf{\Omega },$$

(15)

where $\mathbf{Y}=\left[Y_{ij}\right]\in {\mathbb{R}}^{m\times n}$ is the decision variable, and $\mathsf{\Omega }$ is the set of locations of the observed entries, with each $(i,j)\in \mathsf{\Omega }$ generated by the Bernoulli distribution, i.e., independently with probability p.

Problem (15) seeks to find the simplest explanation fitting the observed data. It is ill-posed in general. The missing entries of $\mathbf{Y}$ can be faithfully recovered with high probability under certain constraints of the matrix rank, missing rate, and sampling scheme [229,233,234].

Rank minimization-based methods [229,235,236] and matrix factorization-based methods [237] are two major categories of low-rank matrix completion methods. Matrix factorization-based methods factorize $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$ of rank-r ($r<min(m,n)$) into the products of two smaller matrices of size ${\mathbf{L}}^T\in {\mathbb{R}}^{m\times r}$ and $\mathbf{F}\in {\mathbb{R}}^{r\times n}$. The missing entries are recovered by finding such pairwise matrices [237]. The low-rank matrix completion can also be approached by the accelerated proximal gradient [238], augmented Lagrange multiplier method [239], spectral methods [240], and singular value thresholding [241].

Many weighted low-rank matrix approximation methods with missing data are presented based on $L_1$-norm and a Laplacian noise model [239,242,243]. They are computationally expensive, and it is difficult to obtain a good solution due to the non-convexity and non-smoothness of the $L_1$-norm-based cost function. In [242], convex programming and weighted median methods are derived via the alternating minimization approach of the $L_1$-norm optimization problem. Convex LP is used in [243]. In [239], a robust PCA that is based on $L_1$-norm and nuclear norm for a non-fixed rank problem is approached by the augmented Lagrange method. The procedure performs SVD at each iteration.

The matrix completion using a non-convex surrogate for the rank function, motivated by optimizing an upper bound of the rank, can be performed with closed-form solutions, such that it converges within dozens of iterations with proven convergence [244]. By exploiting the column-wise correlation, an adaptive correlation learning technique was developed.

4.1. Nuclear Norm Minimization

By nuclear norm minimization, a matrix with missing values can be exactly recovered under some general conditions [229,245,246,247]. When the observed values are noiseless, it is possible to perfectly recover a low-rank matrix [229]. For noisy measurements, recovery is constrained by an error bound that is proportional to the noise level, with high probability [246].

The nuclear norm minimization problem is formulated as in [230,231,241,246]:

$$\underset{\mathbf{Y}}{min}{\parallel \mathbf{Y}\parallel }_{*}=\sum _{k=1}^{min(m,n)}{\sigma }_k\left(\mathbf{Y}\right)\mathrm{subject} \mathrm{to}{Y}_{ij}=X_{ij},(i,j)\in \mathsf{\Omega },$$

(16)

where ${\parallel ·\parallel }_{*}$ is the nuclear norm, and ${\sigma }_k\left(\mathbf{Y}\right)$ is the kth largest singular value of $\mathbf{Y}$.

The problem can be transformed into a quadratically constrained minimization problem:

$$\underset{\mathbf{Y}}{min}{\parallel \mathbf{Y}\parallel }_{*}=\sum _{k=1}^{min(m,n)}{\sigma }_k\left(\mathbf{Y}\right)\mathrm{subject} \mathrm{to}\sum _{(i,j)\in \mathsf{\Omega }}{\left(Y_{ij}-X_{ij}\right)}^2\le \epsilon ,$$

(17)

or a regularized unconstrained problem:

$$\underset{\mathbf{Y}}{min}{\parallel \mathbf{Y}\parallel }_{*}+\lambda \sum _{(i,j)\in \mathsf{\Omega }}{\left(Y_{ij}-X_{ij}\right)}^2.$$

(18)

The nuclear norm problem (16) has to be solved iteratively and it involves SVD at each iteration, leading to high computational costs. Alternating minimization strategies are popular for matrix completion [236,248,249]. The global convergence of the gradient search method for low-rank matrix approximation is proven in [250] by optimizing the Grassmann manifold and Fubini–Study distance on this space. Some nuclear norm-based methods include singular value thresholding [241], robust PCA [17,251], and nuclear norm regularized LS [235].

Problems (16) and (17) can be formulated as semidefinite programs and then solved to global optima by standard semidefinite program solvers when the dimensions are smaller than 500. First-order algorithms, including singular value thresholding, have been proposed in (16) [241]. The proximal gradient method was implemented in (18) [235]. It has linear convergence for (18) under certain conditions [252], but the per-iteration costs for SVD and the matrix memory are high for large matrices. The alternating minimization approach can be easily parallelized, but it requires higher per-iteration computations compared to stochastic gradient descent. There are also parallelizable variants of stochastic gradient descent [253,254] and block coordinate descent [255,256].

Singular value thresholding is a gradient descent method that applies the Uzawa method in [241]

$$\underset{\mathbf{Y}}{min}{\parallel \mathbf{Y}\parallel }_{*}+\alpha {\parallel \mathbf{Y}\parallel }_F^2\mathrm{subject} \mathrm{to}{P}_{\mathsf{\Omega }}\left(\mathbf{Y}\right)=P_{\mathsf{\Omega }}\left(\mathbf{X}\right),$$

(19)

where ${\parallel \mathbf{Y}\parallel }_F={\left({\sum }_{i,j}{Y}_{ij}^2\right)}^{\frac{1}{2}}$ is the Frobenius norm (or $L_2$-norm equivalently), $P_{\mathsf{\Omega }}(·)$ is a function extracting a submatrix from a matrix, with a set of locations $\mathsf{\Omega }$, and $\alpha$ is a regularization parameter.

The nuclear norm regularized LS problem is formulated as [235,257]

$$\underset{\mathbf{Y}}{min}\frac{1}{2}\parallel P_{\mathsf{\Omega }}\left(\mathbf{Y}\right)-P_{\mathsf{\Omega }}{\left(\mathbf{X}\right)\parallel }_F^2+\mu {\parallel \mathbf{Y}\parallel }_{*},$$

(20)

where $\mu$ is a regularization parameter. This problem is solved by accelerated proximal gradient optimization [235,257]. The primal error is smaller than $\epsilon$ after $O(1/\sqrt{\epsilon })$ iterations [235,257].

In robust PCA, the nuclear norm is used for the recovery of the subspace structure from the data that are corrupted by noises or occlusions [17]. The matrix bifactorization method [258] can efficiently approximate the nuclear norm minimization problem so as to mitigate the computation costs of SVD. The method can solve a large variety of low-rank matrix recovery and completion problems, and two linearized proximal alternating optimization algorithms were developed for solving these problems [258].

The nuclear norm, however, is not an ideal approximation for the rank function. In practice, the incoherence property of the nuclear norm heuristic is difficult to meet [229]. For nuclear norm minimization, all singular values are simultaneously minimized, thus the rank cannot be suitably approximated. The truncated nuclear norm is superior to the nuclear norm since it can better approximate the rank of a matrix [259]. The truncated nuclear norm minimization method outperforms its nuclear norm counterpart in terms of convergence speed.

The truncated nuclear norm ${\parallel \mathbf{Y}\parallel }_r$ is defined as the nuclear norm ${\parallel \mathbf{Y}\parallel }_{*}$ subtracted by the r-largest singular values, i.e., the sum of the $min(m,n)-r$ minimum singular values,

$${\parallel \mathbf{Y}\parallel }_r=\sum _{i=r+1}^{min(m,n)}{\sigma }_i\left(\mathbf{Y}\right).$$

(21)

${\parallel \mathbf{Y}\parallel }_r$ is non-convex. Thus, the truncated nuclear norm minimization [259] was formulated by replacing ${\parallel \mathbf{Y}\parallel }_{*}$ in (16) with ${\parallel \mathbf{Y}\parallel }_r$. In [259], an iterative two-step scheme was implemented, where the convex subproblem in the second step was solved by ADMM with excellent convergence accuracy. The method is not robust to r, and it requires many iterations to converge. The convergence is accelerated by using an adaptive penalty parameter for ADMM.

Low-rank matrix recovery can be implemented by spectral regularization, which takes the form of regularization on the singular values of the matrix. The singular values are, in most cases, iteratively computed by applying SVD on a dense matrix. A generalized unitarily invariant gauge function for low-rank matrix recovery does not act on the singular values but generalizes some spectral functions, including the rank function, Schatten p-norm, and log-sum of singular values [260].

4.2. Matrix Factorization-Based Methods

The matrix completion problem considers a matrix $\mathbf{X}\in {\mathbb{R}}^{m\times n}$, with known elements $\mathsf{\Omega }=\left\{\right(i,j\left)\right\}$. The model is given by (15). Matrix completion can also be solved based on matrix factorization [261]. When recovering a rank-k matrix $\mathbf{Y}=\mathbf{U}{\mathbf{V}}^T$ that minimizes the distance between $\mathbf{Y}$ and $\mathbf{X}$ on the known entries of $\mathbf{X}$, we have

$$\underset{\mathbf{Y}}{min}\sum _{(i,j)\in \mathsf{\Omega }}{(Y_{ij}-X_{ij})}^2\mathrm{subject} \mathrm{to}\mathrm{rank}\left(\mathbf{Y}\right)=k.$$

(22)

This matrix factorization model has long been used in PCA.

A maximum-margin factorization method [237] solves the problem

$$\underset{\mathbf{U},\mathbf{V}}{min}\sum _{(i,j)\in \mathsf{\Omega }}{(X_{ij}-{\left(\mathbf{U}{\mathbf{V}}^T\right)}_{ij})}^2+{\beta (\parallel \mathbf{U}\parallel }_F^2+{\parallel \mathbf{V}\parallel }_F^2),$$

(23)

where $\beta$ is a regularization parameter.

The problem (23) and its extensions were investigated in [237,262,263]. In [237], the biconvex method was always stuck at a suboptimal minimum for small rank r, and the computational complexity became very high, as r, m, and n became large.

The maximum-margin matrix factorization method and nuclear norm minimization method have been combined, and the algorithm outperformed both methods for large matrix factorization and completion [264]. It is a stylized variant of the block coordinate descent. A scalable divide-and-conquer framework for noisy matrix factorization and completion achieves near-linear to superlinear speed-ups [265]. The task is randomly divided into subproblems, being solved in parallel using a nuclear norm-based matrix factorization algorithm, and the solutions are combined by using techniques from the randomized matrix approximation. In [266], two low-rank factorization methods are given for the $L_1$-based low-rank matrix approximation. By using the alternating rectified gradient method, proper projection and coefficient matrices are found at low computational and storage costs. An updated direction is first found, and then a step size is selected for updating a matrix. The weighted median algorithm is performed on the matrix at once, while in [242], it is applied column-wise.

The low-rank matrix estimation can be implemented through matrix factorization, which optimizes two low-rank factors via iterative methods, such as gradient descent and the alternating minimization approach. Despite non-convexity, these methods achieve linear convergence when initialized properly. However, for ill-conditioned matrices, the convergence of gradient descent depends linearly on the condition number $\kappa$ of the low-rank matrix, while the per-iteration cost of the alternating minimization approach is often prohibitive for large matrices. Scaled gradient descent (ScaledGD) and its alternating variants have been proposed for low-rank matrix completion [267]. ScaledGD is an adaptively preconditioned or diagonally-scaled gradient descent with a minimal computational overhead, unveiling the implicit regularization properties of ScaledGD. In [268], ScaledGD is confirmed to achieve linear convergence, independent of $\kappa$, when initialized using the standard spectral method, for solving the low-rank matrix sensing, robust PCA, and matrix completion, while maintaining the low per-iteration cost of gradient descent, which is lower than that of the projected gradient descent. ScaledGD achieves $\epsilon$-accuracy in $O(log\left(\frac{1}{\epsilon }\right))$ iterations when initialized by the spectra method.

Power factorization [269] is an efficient alternating minimization algorithm used for recovering general low-rank matrices, and its performance guarantee under a rank RIP assumption is given in [270]. When initialized with the leading right singular vector of the proxy matrix, power factorization stably recovers a rank-r matrix under rank-$2r$ RIP. Sparse power factorization [271] modifies the updates in the power factorization for compressed sensing of sparse rank-one matrices to exploit their sparsity priors. For the recovery of sparse vectors under RIP, the hard thresholding pursuit provides guarantees on both the estimation error and convergence rate [271]. Sparse power factorization converges linearly under the RIP assumption. In the rank-one case, subspace-concatenated sparse power factorization and sparse power factorization have similar near-optimal performance guarantees [271]. For rank-r matrices with a conditioning number of at most $\kappa$, the subspace-concatenated sparse power factorization succeeds with $m=O\left({\kappa }^{2}rn\right)$ measurements, substantially improving on the results of $m=O\left({\kappa }^4{r}^{3}n\right)$ for power factorization [270].

The Bayesian matrix factorization can produce low-rank representations of matrices, predict missing values, and provide confidence intervals. A distributed approach is realized by a hierarchical decomposition of the joint posterior distribution, which couples the subset inferences [272]. The Bayesian deep matrix factorization network [273] is a robust and fast low-rank matrix factorization model used for multi-image denoising. It uses a deep neural network to model the low-rank components and the model is optimized via stochastic gradient variational Bayes. A hierarchical kernelized sparse Bayesian matrix factorization model [274] integrates side information, and infers the parameters and latent variables, including the reduced rank through variational Bayesian inference. The model simultaneously achieves a low rank through sparse Bayesian learning and column-wise sparsity through an enforced constraint on latent factor matrices.

A low-rank positive semidefinite matrix can be factorized into a product of two matrices. By using a Courant penalty that penalizes the differences between certain components, the semidefinite program is formulated as a biconvex optimization problem [275]. This allows using multi-convex optimization techniques for defining simple surrogates, which can be easily minimized by using a block coordinate descent algorithm. The algorithm is as accurate as other semidefinite program algorithms but is much faster.

When factorizing a large square matrix into a number of matrices of much lower ranks, the low-rank constraint cannot be applied if the approximated matrix is intrinsically high-rank or close to full rank. In [276], a large square matrix is approximated with a product of sparse full-rank matrices, using only $N{(logN)}^2$ nonzero numbers for an $N\times N$ full matrix.

Matrix Completion with Side Information

Given all of the side information with a matrix $\mathbf{B}$, and $\mathbf{V}=\mathbf{B}\mathbf{S}$, with $\mathbf{B}={({\overrightarrow{b}}_1,{\overrightarrow{b}}_2,\dots ,{\overrightarrow{b}}_n)}^T\in {\mathbb{R}}^{n\times p}$, ${\overrightarrow{b}}_j={(B_{j1},\dots ,B_{jp})}^T\in {\mathbb{R}}^p$ ($p\ge k$), we have

$$\underset{\mathbf{S},\mathbf{U}}{min}\sum _{(i,j)\in \mathsf{\Omega }}{(Y_{ij}-X_{ij})}^2\mathrm{subject} \mathrm{to}\mathbf{Y}=\mathbf{U}{\mathbf{S}}^T{\mathbf{B}}^T,{\parallel \mathbf{S}\parallel }_F=1,$$

(24)

where the matrix of feature exposures is denoted as $\mathbf{U}\in {\mathbb{R}}^{m\times k}$ and $\mathbf{S}\in {\mathbb{R}}^{p\times k}$. In Netflix, column j corresponds to movie j; thus, ${\overrightarrow{b}}_j$ contains information on movie j.

Given perfectly predictive side information, the theoretical bound of sample complexity $O(logn)$ to retrieve the full matrix can be achieved [277]. In [278], the side information is corrupted with noise, while in [279,280], a nonlinear combination of factors in the side information is explored.

Prior knowledge of the column and row spaces of a matrix can be incorporated by minimizing a weighted nuclear norm [281]. Theoretically, reliable prior knowledge reduces the sample complexity of matrix completion by a logarithmic factor. Similar results for matrix recovery from generic linear measurements are presented in [281]. Without the incoherence assumption, a two-phase sampling algorithm does not need knowledge about the underlying structure of a matrix [233]. In the case where the observed entries are non-uniformly distributed, exact recovery guarantees for the weighted nuclear norm minimization method are provided in [233].

For matrix completion with and without side information, fastImpute [282] is a non-convex gradient descent method for the exact sparse problem. Factorization can be implemented on the matrix of the features in the side information. The method converges to a global minimum that faithfully recovers the underlying matrix and it scales well to matrices of sizes beyond ${10}^5\times {10}^5$. When a high number of entries is missing, fastImpute outperforms other methods [283] in terms of error and convergence times.

4.3. Theoretical Guarantees on the Exact Matrix Completion

Given the compact SVD of $\mathbf{Y}=\mathbf{U}\Sigma {\mathbf{V}}^T$, $\mathbf{U}\in {\mathbb{R}}^{m\times r}$, $\Sigma \in {\mathbb{R}}^{r\times r}$, $\mathbf{V}\in {\mathbb{R}}^{n\times r}$, $\mathbf{Y}$ is $\mu$-incoherent if [229]

$$\sum _{k=1}^r{U}_{ik}^2\le \frac{\mu r}{m},\sum _{k=1}^r{V}_{jk}^2\le \frac{\mu r}{n},\forall i=1,\dots ,m;j=1,\dots ,n,$$

(25)

where the coherence $\mu \in [1,\frac{max(m,n)}{r}]$ measures how spiky a matrix is. This is the standard incoherence condition for matrix completion, and it prevents information from being concentrated in a few rows or columns.

The joint or strong incoherence condition with parameter ${\mu }_1$ is defined as

$$\underset{i,j}{max}\left|{\left(\mathbf{U}{\mathbf{V}}^T\right)}_{ij}\right|\le \sqrt{\frac{{\mu }_{1}r}{mn}}.$$

(26)

Joint incoherence requires the left and right singular vectors to not be aligned. It is pointed out in [284] that the standard and joint incoherence conditions are, respectively, related to the (statistical) information and computational aspects of the matrix decomposition problem.

The algorithm and theoretical guarantees for the exact low-rank matrix completion are first given in [229]. In the case of an incoherent low-rank matrix and uniformly random sampling, nuclear norm minimization is applicable [229]. When there are more than $C{n}^{1.25}rlogn$, with a constant $C>0$, observed entries selected uniformly at random from matrix $\mathbf{X}$, with high probability, matrix $\mathbf{Y}\in {\mathbb{R}}^{n\times n}$ of rank r can be perfectly recovered by nuclear norm minimization. Similar results hold for arbitrary rectangular matrices [229].

Provable completion results for incoherent matrices and uniformly random sampling are refined via nuclear norm minimization [245,247,284], SVD followed by local descent [240], and the alternating minimization approach [285]. The case with additive noise and sparse errors was considered in [17,247,286,287,288].

Most of the existing sufficient conditions [229,247] demand a uniformly random selection of the subset of observed elements and an incoherent or non-spiky low-rank matrix (i.e., with diffuse row and column spaces). Under these conditions, the matrix is provably recoverable by convex optimization [229], alternating minimization [285], and iterative thresholding [241] methods. For the stable recovery of $\mathbf{Y}\in {\mathbb{R}}^{n\times n}$ of rank r, as few as $O(nrlogn)$ measurements suffice for a certain class of sensing systems [231,289].

If a matrix $\mathbf{Y}\in {\mathbb{R}}^{n\times n}$ of rank r satisfies certain incoherence properties, it can, with high probability, exactly reconstruct the matrix from $nr{log}^{2}n\ll n^2$ randomly sampled entries by using efficient polynomial-time algorithms for solving (16) [229,240,245,247,270,290]. This result was generalized to noisy matrix completion by solving (17) [246]. The theoretical guarantee for a variant of (18) is provided in [288]. References [17,89,287] prove the performance guarantees of some algorithms by assuming the standard and joint incoherence conditions.

If the row space can be coherent but the column space is incoherent with parameter ${\mu }_0$, the adaptive sampling algorithm proposed in [291] for matrix completion requires $O({\mu }_0{r}^{3/2}nlog(2r/\delta ))$ observed elements with a success probability of $1-\delta$ [291]. The sample complexity is improved to $O\left({\mu }_{0}rn{log}^2(r^2/\delta )\right)$ in [292].

A theoretical guarantee for the factorization-based low-rank matrix completion is given based on a regularized objective [293]. The exact recovery guarantee and linear convergence for many first-order methods, such as gradient descent and the alternating minimization approach without resampling, are proven [293].

For the exact matrix completion, the joint incoherence condition has proven to be unnecessary and can be eliminated [284]. With $\mathsf{\Omega }\left(nr{log}^{2}n\right)$ (that is, bounded below $nr{log}^{2}n$ asymptotically) uniformly sampled entries, a matrix satisfying standard incoherence but not joint incoherence (for example, a positive semidefinite matrix) can be recovered [284]. In the case of the recovery of a semidefinite matrix, the sample complexity is reduced to $O\left(nr{log}^{2}n\right)$, and the highest allowable rank is improved to $\mathsf{\Theta }(n/{log}^{2}n)$ [284]. The analysis is based on $L_{\infty ,2}$ matrix norm (i.e., the maximum of the row and column norms of a matrix). The results apply to the nuclear norm minimization approach to matrix completion.

The compressive adaptive sense and search (CASS) algorithm [294] is a simple adaptive sensing and group testing algorithm used for sparse signal recovery. To recover a k-sparse signal of dimension n, standard compressed sensing based on random Gaussian non-adaptive design matrices requires the SNR to grow, such as $logn$. Similar to standard compressed sensing, CASS requires only $klogn$ measurements, but CASS is near-optimal as it succeeds at the lowest possible signal-to-noise-ratio (SNR), which scales similar to $logk$, and is a factor $logn$ lower. CASS is substantially less computationally intensive than standard compressed sensing.

One can successfully retrieve a matrix from $O\left(\mu rn{log}^{2}n\right)$ uniform samples [233,246,284]. Low-rank matrix completion is investigated based on the leave-one-out analysis [295]. Projected gradient descent for a rank-constrained formulation is also known as the singular value projection or IHT. The projection can be efficiently computed by rank-r SVD. Projected gradient descent without regularization or sample splitting converges linearly in the infinity norm [295]. The nuclear norm minimization recovers $\mathbf{Y}\in {\mathbb{R}}^{n\times n}$ of rank r with a high probability with $O(\mu rlog(\mu r)nlogn)$ observed entries [295]. This result is better than some earlier results: $O\left(\mu rn{log}^{2}n\right)$ [284], $O({\kappa }^2\mu rnmax\{logn,\mu r{\kappa }^4\})$ [240], $O({\kappa }^2\mu rnmax\{logn,\mu r^6{\kappa }^4\})$ [293]. It is independent of the condition number $\kappa$, and matches the information-theoretic lower bound $C\mu rnlogn$ [247].

Compressed sensing suffers from a basis mismatch when imposing a discrete dictionary on the Fourier representation. This issue can be solved by enhanced matrix completion [296], which is based on the structured matrix completion without knowledge of the model order. The method arranges the data into a low-rank enhanced form that exhibits a multi-fold Hankel structure, and then implements recovery through nuclear norm minimization. Under mild incoherence conditions, the method ensures perfect recovery provided that the number of samples exceeds $O\left(r{log}^{4}n\right)$, and the performance is stable in the presence of bounded noise. Exact recovery is still possible when many samples are corrupted with noise of arbitrary magnitude, as long as the sample complexity exceeds $O\left(r^2{log}^{3}n\right)$.

For nuclear norm minimization, the observed entries are typically assumed to be sampled uniformly at random [229,247,261]. In practice, nuclear norm minimization works quite well for non-uniform data. Under distribution-free and some very mild assumptions, exact recovery is possible from $O\left(n^{3/2}\right)$ entries of a matrix $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$, $m\le n$, in a nuclear norm regime [297]. This bound is tight.

Based on a deterministic analysis, a sufficient and necessary condition for the recovery of a low-rank matrix is that the pattern of revealed entries remains dense under any relabeling of rows and columns [298]. When the missingness of the matrix entries is dependent on the unobserved values themselves, a procedure is devised to simultaneously complete the matrix and assess the covariate effect [299]. Under the assumptions of a low-rank matrix and sparse covariate effects, the statistical guarantee of the procedure and convergence is established by allowing the matrix dimensions and the number of covariates to grow ultra-high.

Suppose that a rank-r matrix $\mathbf{Y}\in {\mathbb{R}}^{n\times n}$ has up to s nonzero rows and up to s nonzero columns, and its measurements are generated as the inner products with i.i.d. Gaussian matrices. By solving a combinatorial optimization problem, faithful recovery can be guaranteed with $O(max\{rs,slog(en/s)\left\}\right)$ measurements [300], which is significantly better than the exact recovery using combinations of the nuclear norm and $L_{1,2}$-norm, $O(minrn,sn)$ measurements. When some columns are completely and arbitrarily corrupted, a matrix completion algorithm proposed in [301] combines a trimming procedure with a convex program that minimizes the sum of the nuclear norm and $L_{1,2}$-norm. As the portion of observed entries is diminishing, matrix completion is possible, even for a growing number of corrupted columns. From an information-theoretic viewpoint, the guarantees are nearly optimal with respect to the rank of the underlying matrix, the portion of corrupted columns, and the portion of sampled entries on the authentic columns [301]. When the observed samples simultaneously contain both erasures and errors, with a constant fraction of values arbitrarily corrupted, a unified performance guarantee on the exact recovery of a low-rank matrix by minimizing the nuclear norm plus $L_1$-norm is given in [287].

4.4. Discrete Matrix Completion

In recommender systems, a rating of thumbs-up or thumbs-down is denoted as a single bit for each occurrence. In the Netflix problem, movies are rated as integers from 1 to 5. Discrete matrix completion is often used in surveys.

For noisy one-bit matrix completion in a general non-uniform sampling distribution, an $L_{\infty }$-norm-constrained maximum-likelihood estimate is derived in [302]. The optimal rate of convergence for the Frobenius-norm loss is defined by the minimax upper and lower bounds together. One-bit matrix completion in the uniform sampling distribution is analyzed in [303]. The minimax optimal rate of convergence is achieved by the nuclear norm-constrained maximum-likelihood approach, and an approximately low-rank matrix $\mathbf{X}$ is recovered from a set of noisy sign (one-bit) measurements. Consider the recovery of a low-rank matrix $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$ given a subset of noisy discrete measurements. Based on the low-rank factorization of $\mathbf{Y}$, a globally convergent constrained maximum-likelihood algorithm is derived under constraints on the $L_{\infty }$-norm of $\mathbf{X}$ and exact rank [304]. The likelihood can be from any strictly log-concave distribution, including from the exponential family, which is suited for bounded discrete random variables.

In collaborative filtering, the sampling distribution is non-uniform. When certain rows or columns are sampled with probabilities that are too high, the nuclear norm minimization method might fail, and the weighted nuclear norm that takes the sampling distribution into account is viable [305]. Rigorous recovery guarantees for learning with various weighted nuclear norms are given in [306]. Approximate low-rank matrix completion using the weighted nuclear norm in the general sampling scheme is theoretically guaranteed in [288].

For discrete matrix completion, upper bounds on the error norm are derived in [302,303,307]. The results in [303] have been extended to multilevel observations for matrix completion in the case of categorical data [308].

5. Low-Rank Representation

Low-rank representation methods are robust in handling noise/corrupted data. For observed data corrupted with sparse errors, low-rank representation methods [227,228,262] jointly learn the lowest-rank representation of all data. For the recovery of a sparse and low-rank matrix from its minimal incoherent linear measurements, the exact recovery of the low-rank matrix under rank RIP is guaranteed for the minimum nuclear norm solution [231].

The low-rank matrix approximation/recovery problem [17] aims to decompose a data matrix into a sum of a low-rank matrix and a sparse matrix,

$$\mathbf{Y}=\mathbf{X}+\mathbf{E}\in {\mathbb{R}}^{m\times n},$$

(27)

where $\mathbf{X}$ is a low-rank matrix, $\mathbf{E}$ is a sparse corruption matrix with Gaussian noise. For PCA, $\mathbf{X}$ denotes a matrix of m data points in a low-dimensional subspace ${\mathbb{R}}^n$, corrupted by a sparse corruption matrix $\mathbf{E}$ of errors and Gaussian noise. Equation (27) is also known as the matrix decomposition model.

The goal is to recover these components from (27). One needs to impose conditions on the sparse and low-rank components to guarantee their identifiability. We have the following program:

$$\underset{\mathbf{X},\mathbf{E}}{min}\mathrm{rank}\left(\mathbf{X}\right)+\lambda {\parallel \mathbf{E}\parallel }_0\mathrm{subject} \mathrm{to}\mathbf{Y}=\mathbf{X}+\mathbf{E},$$

(28)

where $\lambda$ is a penalty parameter.

When the rank of $\mathbf{X}$ is not too large and $\mathbf{E}$ is sufficiently sparse, problem (28) is equivalent to the following convex version [17]:

$$\underset{\mathbf{X},\mathbf{E}}{min}{\parallel \mathbf{X}\parallel }_{*}+\lambda {\parallel \mathbf{E}\parallel }_1\mathrm{subject} \mathrm{to}\mathbf{Y}=\mathbf{X}+\mathbf{E},$$

(29)

where ${\parallel \mathbf{X}\parallel }_{*}$ is the nuclear norm, and ${\parallel \mathbf{E}\parallel }_1$ is the $L_1$-norm.

Low-rank representation is related to robust PCA [229]. Robust PCA does not have a polynomial-time algorithm under broad conditions [17]. Problem (29) is convex and can be solved in polynomial time [17]. The principal component pursuit is a method used for solving robust PCA [17]. The augmented Lagrange multipliers method [239] is used for its computational efficiency. The low-rank representation can better capture the global structure of the data while ignoring the local manifold structure. It can achieve a block-diagonal solution for independent subspaces and sufficient sampling [227]. A low-rank representation is a subspace clustering method [309].

A low-rank representation can be formulated in the same form as (28) or (29), but replacing $\mathbf{Y}=\mathbf{X}+\mathbf{E}$ with $\mathbf{Y}=\mathbf{A}\mathbf{X}+\mathbf{E}$, where $\mathbf{Y}$ is the sample set, $\mathbf{A}$ is the dictionary, and $\mathbf{E}$ is the error. In the noise-free case, the solutions to the two formulations are the same [228].

The low-rank representation simultaneously recovers the row space of a data matrix and detects outliers under mild conditions [228]. A double low-rank representation [310] learns the row and column spaces embedded in a matrix at once. A two-step approach presented in [311] achieves near-optimal sample complexity for a special measurement scheme with nested structures. When the observed matrix is noisy and contains outliers, the Huber function is used to downweight the outliers, and the method is fast and monotonically convergent [312]. A robust PCA method, referred to as the nuclear norm-based PCA [313], remedies the bias of the computed mean and low-dimensional representation of a sample. The algorithm has a closed-form solution at each iteration.

The matrix decomposition model (27) was solved using convex relaxation (29). Convex approximations lead to biased estimates. Non-convex regularizers, such as weighted nuclear norm minimization, weighted Schatten p-norm minimization, and weighted matrix gamma norm for the low-rank part, can be used. In [314], the ADMM technique is implemented on a formulation composed of the weighted minimax-concave penalty and weighted matrix gamma norm, coupled with the dynamic weight update.

For the matrix decomposition problem that separates a rank-$\mathsf{\Omega }\left(n^{\frac{1}{2}}\right)$ positive semidefinite matrix from a sparse matrix, the joint incoherence condition is unavoidable, as without it, any algorithm needs to solve the planted clique problem, which is, in general, intractable in polynomial time [284]. In [315], the decomposition of a positive semidefinite matrix (e.g., a covariance matrix) into a matrix with a given rank and a sparse matrix is implemented by using a deep neural network. The gradient descent algorithm has polynomial convergence in terms of the number of the network size.

Since nonlinear geometric structures in data are not considered in low-rank representations, locality and similarity information in data may be lost during the learning process. Manifold learning methods preserve local geometric structures; examples are ISOMAP [316], locally linear embedding [317], locality-preserving projection [318], neighborhood-preserving embedding [319], Laplacian Eigenmap [320], and nonnegative sparse hyper-Laplacian regularized low-rank representation [321].

Problem (29) can be formulated as [322]

$$\underset{\mathbf{X},\mathbf{E}}{min}{\parallel \mathbf{X}\parallel }_{*}+{\lambda \parallel \mathbf{E}\parallel }_1\mathrm{subject} \mathrm{to}{\parallel \mathbf{Y}-\mathbf{X}-\mathbf{E}\parallel }_F\le \epsilon ,$$

(30)

where perturbation $\epsilon$ is small. A regularized formulation is given by [322]

$$\underset{\mathbf{X},\mathbf{E}}{min}{\parallel \mathbf{X}\parallel }_{*}+{\lambda \parallel \mathbf{E}\parallel }_1+\frac{1}{2\mu }{\parallel \mathbf{Y}-\mathbf{X}-\mathbf{E}\parallel }_F^2,$$

(31)

where $\mu$ is a regularization parameter.

Rank-sparsity incoherence is applied to matrix decomposition. The incoherence parameters of $\mathbf{E}$ and $\mathbf{X}$ are sufficient to ensure identifiability and recovery using convex programs [286]. Following the analysis of [286], a weaker condition is given in [322]. $\mathbf{E}$ is allowed to have up to $\mathsf{\Omega }\left(mn\right)$ nonzero entries when $\mathbf{X}$ is low-rank and has non-sparse singular vectors [322]. In terms of PCA, the analysis allows for a constant fraction of the matrix entries to be corrupted by noise of arbitrary magnitude. Purely deterministic structural conditions on the sparsity pattern of $\mathbf{E}$ [322] lead to a cost of roughly a factor of $\mathrm{rank}\left(\mathbf{X}\right)$ in the allowed support size of $\mathbf{E}$, compared to the probabilistic analysis of robust PCA [17].

GoDec [18] is an efficient, robust, low-rank matrix decomposition algorithm used for solving (27). It alternatingly sets $\mathbf{X}$ as a low-rank approximation of $\mathbf{Y}-\mathbf{E}$ and sets $\mathbf{E}$ as a sparse approximation of $\mathbf{Y}-\mathbf{X}$. It can be substantially accelerated by bilateral random projections. As the error function ${\parallel \mathbf{Y}-\mathbf{X}-\mathbf{E}\parallel }_F^2$ converges to a local minimum, $\mathbf{X}$ and $\mathbf{E}$ converge linearly to their local optima [19]. GoDec+ [19] is superior to GoDec in terms of robustness and convergence speed. It maximizes a correntropy criterion by using a greedy bilateral paradigm for half-quadratic optimization. GoDec+ is robust to a variety of corruptions, including Gaussian, Laplacian, salt and pepper, and occlusion.

6. Nonnegative Matrix Factorization

The NMF problem aims to factorize a nonnegative matrix $\mathbf{X}$ into lower-rank nonnegative matrix factors,

$$\mathbf{Y}=\mathbf{A}\mathbf{X},$$

(32)

where $\mathbf{Y}=\left[{\overrightarrow{y}}_1,\dots ,{\overrightarrow{y}}_n\right]\in {\mathbb{R}}^{m\times n}$, each nonnegative data point ${\overrightarrow{x}}_i\in {\mathbb{R}}^m$, $\mathbf{A}={[{\overrightarrow{a}}_1,\dots ,{\overrightarrow{a}}_m]}^T\in {\mathbb{R}}^{m\times k}$ is a basis matrix known as the dictionary, ${\overrightarrow{a}}_i={(A_{i1},\dots ,A_{ik})}^T$ is a basic vector, the matrix of sources, denoted as $\mathbf{X}=[{\overrightarrow{x}}_1,\dots ,{\overrightarrow{x}}_n]\in {\mathbb{R}}^{k\times n}$ is the code or coefficient matrix of $\mathbf{Y}$ using the dictionary $\mathbf{A}$; all elements in these matrices are nonnegative, and the rows in $\mathbf{X}$ may be statistically dependent.

Problem (32) can be written as

$${\overrightarrow{y}}_i\approx \mathbf{A}{\overrightarrow{x}}_i.$$

(33)

That is, ${\overrightarrow{y}}_i$ is represented by a linear combination of each row of $\mathbf{A}$, multiplied by the corresponding elements of ${\overrightarrow{x}}_i$.

Usually, $k\ll m,n$; thus, NMF results in a compressed representation of a matrix. When k stands for the number of clusters, the jth row of $\mathbf{A}$, ${\overrightarrow{a}}_j^T$, is a representation of cluster j, and a feature vector ${\overrightarrow{y}}_i$ is assigned to the cluster with the largest weight.

NMF is plagued by scaling and permutation indeterminacy [323]. A unique decomposition can be obtained by imposing some constraints on the factors. The non-uniqueness problem is illustrated in two dimensions in Figure 2. The basis vectors ${\overrightarrow{h}}_1$ and ${\overrightarrow{h}}_2$ can be, respectively, put anywhere in each of the two open spaces between the coordinate axes and the data, and each data point can be exactly represented by a nonnegative linear combination of the two vectors. For some well-posed NMF problems, there are optimal and sparse solutions under the separability assumption [324].

Although both NMF and dictionary learning techniques learn sparse representations, NMF learns low-rank representations, while dictionary learning usually learns full-rank representations.

NMF is NP-hard, in general [325], and highly ill-posed. NMF is tractable under the separability assumption, under which all columns of $\mathbf{Y}$ are on the convex cone generated by some of these columns. When there exists a rank-r NMF $\mathbf{Y}=\mathbf{A}\mathbf{X}$, where columns of $\mathbf{A}$ equal some columns of $\mathbf{Y}$, i.e., the basis vectors are selected as some of the data points, this is known as separability [326]. This convex optimization problem can be solved using gradient descent and successive projection methods [327].

NMF problems are usually non-convex. NMF finds localized, parts-based representations of nonnegative multivariate data [328]. NMF is usually performed with an alternating gradient descent technique, belonging to a class of multiplicative iterative algorithms [21]. While the method has low complexity, it converges slowly and is susceptible to getting stuck in the local minima of the cost function. Alternating nonnegative LS is another popular algorithm [23]. The cost functions can be the squared Euclidean distance, Kullback–Leibler divergence, and their unified functions, such as the $\alpha$-divergence, $\beta$-divergence [329], or a broader family known as the Bregman divergence [330]. The Frobenius norm and Kullback–Leibler divergence are two examples of the Bregman divergence. The projected gradient method [331,332] is efficient for dealing with large-scale NMF problems under nonnegativity and sparsity constraints.

Algorithms for NMF can be extended to BSS by adding regularization terms for characterizing sparseness, smoothness, or effective expressions of patterns of the estimated components [333,334]. The smooth component analysis [335] imposes smoothness constraints on vectors of the factor matrix or on vectors of the mixing matrix. Boolean matrix factorization, also known as Boolean factor analysis, is a method that factorizes datasets into binary alphabets [336].

6.1. Multiplicative Update Algorithm

The NMF problem can be formulated as

$$\underset{\mathbf{A},\mathbf{X}}{min}{\parallel \mathbf{Y}-\mathbf{A}\mathbf{X}\parallel }_F,\mathrm{subject} \mathrm{to}{A}_{ij}\ge 0,X_{ij}\ge 0\forall i,j.$$

(34)

The multiplicative update rule is formulated as [21]

$$\mathbf{X}\leftarrow \mathbf{X}\otimes \frac{{\mathbf{A}}^T\mathbf{Y}}{{\mathbf{A}}^T\mathbf{A}\mathbf{X}},\mathbf{A}\leftarrow \mathbf{A}\otimes \frac{\mathbf{Y}{\mathbf{X}}^T}{\mathbf{A}\mathbf{X}{\mathbf{X}}^T},$$

(35)

where ⊗ and / are element-wise multiplication and element-wise division, respectively. Initially, $\mathbf{A}$ and $\mathbf{X}$ can be set as positive random values.

These iterates (35) guarantee monotonic convergence to a local maximum of [21]

$$F=\sum _{i=1}^m\sum _{j=1}^n(Y_{ij}ln{\left(\mathbf{A}\mathbf{X}\right)}_{ij}-{\left(\mathbf{A}\mathbf{X}\right)}_{ij}).$$

(36)

After learning $\mathbf{A}$, new data ${\mathbf{Y}}^{\prime }$ are mapped onto the k-dimensional space by fixing $\mathbf{A}$, randomizing $\mathbf{X}$, and iterating until convergence. In this procedure, $\mathbf{X}$ can be alternatively obtained by fixing $\mathbf{A}$, followed by solving ${\mathbf{Y}}^{\prime }=\mathbf{A}{\mathbf{X}}^{\prime }$ for ${\mathbf{X}}^{\prime }$ by pseudoinversion. Pseudoinversion may yield negative entries in ${\mathbf{X}}^{\prime }$. One can enforce nonnegativity by forcing negative values to zero or by using nonnegative LS. Information loss may arise from setting negative values to zero.

The multiplicative update algorithm may fail to converge to a stationary point [24,337]. The multiplicative update algorithm is guaranteed to converge to a stationary point when minimizing the Euclidean distance given by (34) [338].

6.2. Alternating Nonnegative Least Squares

Alternating nonnegative LS is a block coordinate descent method for bound-constrained optimization [23]:

$$\underset{\mathbf{A},\mathbf{X}}{min}F(\mathbf{A},\mathbf{X})=\frac{1}{2}{\parallel \mathbf{Y}-\mathbf{A}\mathbf{X}\parallel }_F^2,\mathrm{subject} \mathrm{to}{A}_{ij}\ge 0,X_{ij}\ge 0\forall i,j.$$

(37)

The algorithm is implemented as two alternating steps, each solving a convex optimization problem,

$${\mathbf{A}}_{k+1}=arg\underset{\mathbf{A}\ge 0}{min}F(\mathbf{A},{\mathbf{X}}_k),{\mathbf{X}}_{k+1}=arg\underset{\mathbf{X}\ge 0}{min}F({\mathbf{A}}_{k+1},\mathbf{X}).$$

(38)

The algorithm converges fast, but it does not warrant convergence [24]. A modified strategy presented in [337] ensures convergence to a stationary point of NMF [339].

Gradient-based algorithms converge slowly for large-scale problems. The projected gradient method [340], projected alternating LS method [24], and projected Newton method [341] speed up the convergence of NMF. The projected gradient algorithm [340] solves nonnegative LS problems for $\mathbf{A}$ and $\mathbf{X}$ in an alternating manner. By applying the projected gradient procedure, the limit point is guaranteed to be a stationary point for optimization [340]. In the multiplicative update algorithm, a fixed step size is used for gradient descent, whereas in the projected gradient algorithm, a flexible step size is used.

Hierarchical alternating LS [342] is a fast-converging method. It solves a set of column-wise nonlinear LS problems for each column. $\mathbf{A}$ and $\mathbf{X}$ are updated column-wisely.

The deficiencies of gradient descent methods, such as the multiplicative update algorithm and alternating nonnegative LS, are overcome by an improved algorithmic framework for the LS NMF problem [343]. This framework allows second-order optimization techniques to be used, such as the Newton technique, BFGS, and the conjugate gradient method [13], and allows incorporating regularization and box constraints. In a projected quasi-Newton method [341], a regularized Hessian is inverted with the Q-less QR decomposition. The Levenberg–Marquardt iterates are used for updating $\mathbf{A}$ and the fixed-point regularized LS method is used for updating $\mathbf{X}$. The quasi-Newton fixed-point algorithm slightly outperforms the conjugate gradient method with gradient projection.

6.3. Other NMF Methods

6.3.1. Sparse NMF

$\mathbf{A}$ and $\mathbf{X}$ are usually sparse matrices. Some modifications to NMF for unique decomposition impose a sparseness constraint on the mixing matrix, source matrix, or both [344]. Some circumstances for unique decomposition were investigated in [345].

Sparse NMF can be implemented by constraining the $L_1$-norm or $L_0$-norm on the factor matrices [346]. Doubly sparse NMF arises when the sparsity constraint is applied to both factor matrices, which can be useful for document–word co-clustering. For NMF with minimum volume constraint [347], which is $L_0$-norm oriented, quadratic programming is an efficient method used for small-scale problems, whereas the multiplicative update algorithm incorporating the natural gradient can be used for large-scale problems.

From a geometric perspective, NMF attempts to find an appropriate simplicial cone to contain $\mathbf{Y}$ [323]. In the large-cone NMF method [348], a large-cone penalty is applied to NMF to simultaneously reduce the reconstruction error and improve the generalization ability. The large-cone NMF produces bases that consist of large simplicial cones. The low-overlapping properties of these bases lead to sparse bases and make the algorithms robust.

6.3.2. Projective NMF

Projective NMF approximates a matrix through nonnegative subspace projection [349,350]

$$\mathbf{Y}=\mathbf{P}\mathbf{Y},$$

(39)

where $\mathbf{P}$ is a positive low-rank matrix. The Frobenius norm or a modified Kullback–Leibler divergence can be used as the cost function. Projective NMF is approached by multiplicative update rules, with proven convergence [350]. In [349], projective NMF is implemented by using a nonnegative multiplicative version of Oja’s rule.

Compared to NMF, projective NMF replaces $\mathbf{X}$ with ${\mathbf{A}}^T\mathbf{Y}$. Thus, projective NMF is close to nonnegative PCA. In projective NMF, $\mathbf{A}$ is approximately orthogonal. This leads to sparsity in the approximation, reduces the computational complexity of learning, exhibits close equivalence to clustering, and allows for easy extension to kernel methods. For highly nonlinearly distributed data, kernelization is a desirable solution. A nonlinear nonnegative component analysis method, referred to as the projected gradient kernel NMF [351], can use arbitrary positive definite kernels. The method outperforms kernel PCA and kernel ICA in terms of classification accuracy [351].

6.3.3. Graph-Regularized NMF

NMF does not include the neighborhood structure information and, thus, may fail when handling datasets with nonlinear structures. Topographic NMF [22] incorporates and organizes neighborhood connections between NMF basis functions on a topographic map. Basis functions overlap along neighboring structures due to nonnegativity. In topographic NMF, inputs are represented by multiple activity peaks, while the conventional self-organizing map model represents the input using a single activity peak on a topographic map.

Graph-regularized NMF [352] regularizes NMF formulations using a graph Laplacian matrix, which is a Euclidean distance-based similarity matrix obtained by using a fixed small subset of the neighborhood information. This problem is solved by the multiplicative update rule, which converges to one of multiple fixed points under given conditions, based on the Lyapunov indirect method [353]. The graph-regularized NMF is rather sensitive to the number of nearest neighbors and the regularization parameter.

The manifold regularization-based matrix factorization model [354] has globally optimal and closed-form solutions. It outperforms the graph-regularized NMF. The model is solved by a direct algorithm for a small data matrix and an alternating iterative algorithm with inexact inner iterations for a large data matrix. Comparatively fewer attempts have been made for the online graph-regularized NMF [355] due to the high computational requirements, making online updates of geometric weights impractical when incorporating geometric structures.

Neighborhood structure-assisted NMF [356] incorporates a neighborhood structural similarity matrix based on a minimum spanning tree. The neighborhood parameter is not used and its result is much less sensitive to the regularization parameter. Graph-regularized NMF and symmetric NMF are closely related to neighborhood structure-assisted NMF.

6.3.4. Weighted NMF

NMF attaches the same importance to all attributes of a data point. The features usually have different importance. The methods include feature-weighted NMF [357] and entropy-weighted NMF [358].

6.3.5. Bayesian NMF

NMF is not statistically regular, and the prior distribution used in variational Bayesian NMF has zero or divergence points. An analysis of the Kullback–Leibler divergence between the variational posterior and the true posterior is given in [359]. A lower bound for the approximation error of Bayesian NMF is derived, and it is dependent on the hyperparameters and the true nonnegative rank [359]. A family of NMF algorithms, including those under sparsity constraints, are derived using a statistical framework based on the generalized dual Kullback–Leibler divergence, which includes members of the exponential family of models [360].

6.3.6. Supervised or Semi-Supervised NMF

Combining label information can improve the discriminating power of the matrix decomposition. Supervised NMF methods, such as discriminant NMF [361] and max–min distance NMF [362], utilize class label information. Discriminant NMF [361] introduces Fisher’s discriminative information to NMF to enhance the classification accuracy. Max–min distance NMF [362] minimizes the maximum distance of the within-class pairs in the new NMF space and maximizes the minimum distance of the between-class pairs in an alternating way. For NMF, the learned basis is not necessarily parts-based [323]. Manifold-regularized discriminative NMF applies manifold regularization and margin maximization on NMF. It can produce parts-based bases using a Newton-based method.

Semi-supervised NMF [363] performs a joint factorization of the data and label matrices; it involves a common factor matrix $\mathbf{X}$. Constrained NMF is a semi-supervised method that incorporates the label information as additional constraints [364]. For the constrained clustering problem, domain knowledge is present in the form of must-link and cannot-link. NMF-based [365] and symmetric NMF-based [366] constrained clustering algorithms are semi-supervised NMF algorithms. They enforce the similarity between two points in a must-link constraint towards 1 and the similarity between two points in a cannot-link constraint towards 0.

6.3.7. NMF for Mixed-Sign Data

Semi-NMF, convex NMF, and cluster NMF algorithms are used to deal with mixed-sign data [367]. Semi-NMF is defined by $\mathbf{Y}=\mathbf{A}\mathbf{X}$, with the constraint that the elements of $\mathbf{X}$ are nonnegative [367]. In convex NMF, the basis vectors in $\mathbf{A}$ are constrained to be convex combinations of the data points [367]. Convex NMF can be used for a kernel extension of NMF. It is applicable to a nonnegative or mixed-sign data matrix. The generated factor matrices tend to be very sparse. Cluster NMF [367] is a particular case of convex NMF. It adopts an idea similar to projective NMF and is based on the Frobenius norm.

6.3.8. Deep NMF

By repeatedly decomposing the matrix, a hierarchical deep neural network structure for NMF [368] can provide more interpretable representations of the data.

6.3.9. NMF for BSS

Without the assumption of independence, NMF-based BSS [333] successfully estimates the original sources from the mixtures. The convergence area for NMF-based BSS is obtained using the invariant set method [369].

6.3.10. Online NMF

Online NMF [370] and incremental orthogonal projective NMF [371] are incremental NMF algorithms for data streams. In [372], online NMF in the presence of outliers is solved using projected gradient descent and ADMM, with proven convergence.

6.3.11. Coordinate Descent for NMF

Greedy coordinate descent for NMF [373] is an element-wise update algorithm. The most influential variables are selected for minimization. When a constraint, e.g., graph regularized constraint, affects all elements of one column at once, the method is invalid. It is not applicable to orthogonal NMF since orthogonality requires interactions between rows. Scalar block coordinate descent for Bregman divergence NMF [374] is a column-wise update algorithm. An element-wise algorithm is derived by using the Taylor series expansion of Bregman divergence, with complexity the same as that of a column-wise update algorithm. The scalar block coordinate descent algorithm for the Bregman divergence orthogonal NMF is a column-wise update algorithm, which incorporates the column-wise orthogonal constraint [375].

6.3.12. Robust NMF

The robust NMF using the $L_{2,1}$-norm loss function is applicable for data with Laplacian noise [376]. Truncated Cauchy NMF using truncated Cauchy loss robustly learns the subspace on noisy datasets contaminated by outliers [377].

6.4. NMF for Clustering

Many clustering methods [3] can be described as matrix factorization problems. In [21,23], nonnegative factors of matrices are interpreted as data clustering. Under certain assumptions, NMF is equivalent to clustering [26,344]. Every column of $\mathbf{A}$ represents a cluster center, while $\mathbf{X}$ denotes cluster membership.

For n feature vectors ${\overrightarrow{x}}_i\in {\mathbb{R}}^m$, $i=1,2,\dots ,n$, gathered in $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$, the C-means problem can be formulated as (34) subject to $\mathbf{X}\in {\{0,1\}}^{k\times n}$, ${\mathbf{X}}^T{\overrightarrow{1}}_k={\overrightarrow{1}}_n$. The columns of $\mathbf{A}$ can be treated as the k cluster centroids. If the jth sample belongs to the ith cluster, $X_{ij}=1$, otherwise $X_{ij}=0$.

Orthogonal NMF applies the orthogonality constraint on either $\mathbf{A}$ or $\mathbf{X}$ [26]. C-means clustering corresponds to NMF with orthogonality on $\mathbf{X}$ [8]. Orthogonal NMF with constraints on $\mathbf{A}$ (or $\mathbf{X}$) is equivalent to clustering the rows (or columns) of $\mathbf{Y}$ [378]. Cluster NMF is also close to C-means clustering. Orthogonal NMF imposes an orthogonal constraint onto NMF,

$$\underset{\mathbf{A},\mathbf{X}}{min}{\parallel \mathbf{Y}-\mathbf{A}\mathbf{X}\parallel }_F^2\mathrm{subject} \mathrm{to}\mathbf{A}\ge 0,\mathbf{X}\ge 0,{\mathbf{A}}^T\mathbf{A}=\mathbf{I},$$

(40)

where $\mathbf{I}$ is the identity matrix.

Orthogonal NMF can be solved using the multiplicative update algorithm [26,379]. In [379], zero values are replaced by small positive values to deal with the zero-lock problem.

Bregman divergence orthogonal NMF [375] is equivalent to Bregman hard clustering [380]. Let $\mathbf{Y}$ be an instance feature matrix and $\mathbf{A}$ be an indicator matrix in orthogonal NMF. Then the ith row of $\mathbf{A}$ is considered a membership vector of the i instance to all k groups (features). The ith instance is assigned to the jth cluster, $j={argmax}_j{A}_{ij}$. The Bregman divergence orthogonal NMF is formulated as

$$\underset{\mathbf{A},\mathbf{X}}{min}{\mathbf{D}}_{\varphi }(\mathbf{Y}\parallel \mathbf{A}\mathbf{X})\mathrm{subject} \mathrm{to}\mathbf{A}\ge 0,\mathbf{X}\ge 0,{\mathbf{A}}^T\mathbf{A}=\mathbf{I}.$$

(41)

Fast hierarchical alternating LS for orthogonal NMF [375] is derived from hierarchical alternating LS [342]. The scalar block coordinate descent for orthogonal NMF extends hierarchical alternating LS for orthogonal NMF but has slower convergence.

NMF and symmetric NMF can effectively cluster linearly separable data and nonlinearly separable data, respectively. Symmetric NMF is shown to be highly related to spectral clustering [8]. In [381], multiplicative update rules are used in NMF-based (or symmetric NMF-based) constrained clustering for clustering linearly (or nonlinearly) separable data. A directional clustering approach [382] is inspired by ideas from constrained low-rank matrix factorization and sparse approximation.

Symmetric NMF is defined by [383]

$$\mathbf{Y}=\mathbf{P}{\mathbf{P}}^T,$$

(42)

where $\mathbf{Y}$ is completely positive and $\mathbf{P}$ is nonnegative.

Symmetric NMF generates a symmetric and nonnegative low-rank approximation to the graph matrix. It uses the Euclidean distance-based similarity metric, considering the full graph, and can be treated as a graph clustering algorithm. The use of a dense pairwise similarity measure is computationally expensive. Symmetric NMF makes its outcome less interpretable than neighborhood structure-assisted NMF or graph-regularized NMF. Parallel multiplicative update algorithms [384] have demonstrated convergence under mild conditions, and are applied to probabilistic clustering.

Weighted symmetric NMF, also known as symmetric nonnegative tri-factorization, is defined by

$$\mathbf{Y}=\mathbf{P}\mathbf{Q}{\mathbf{P}}^T,$$

(43)

where $\mathbf{Q}$ is a symmetric nonnegative matrix.

Some alternating iterative algorithms for symmetric NMF with theoretically proven convergence include the progressive hierarchical alternating least squares method for symmetric NMF [385], symmetric NMF based on non-symmetric transformation [386], and semi-supervised structured symmetric NMF-based clustering [387] with simultaneous sparseness and smoothness constraints.

6.5. Concept Factorization

Concept factorization [388] extends NMF for data clustering. The superiority of concept factorization over NMF is demonstrated for document clustering [388].

In concept factorization, each cluster center (concept) ${\overrightarrow{x}}_c$, $c=1,\dots ,k$, is modeled as a nonnegative linear combination of data points ${\overrightarrow{y}}_j$, $j=1,\dots ,n$, and each data point ${\overrightarrow{y}}_j$ is modeled as a nonnegative linear combination of the cluster centers (concepts) ${\overrightarrow{x}}_c$,

$${\overrightarrow{x}}_c=\sum _{j=1}^n{w}_{jc}{\overrightarrow{y}}_j,{\overrightarrow{y}}_j=\sum _{c=1}^k{v}_{jc}{\overrightarrow{x}}_c,$$

(44)

where $w_{jc}\ge 0$ represents the degree of representativeness of ${\overrightarrow{y}}_j$ in concept c, and $v_{jc}\ge 0$ represents its degree belonging to concept c. Clustering is accomplished by computing the two sets of nonnegative linear coefficients that minimize the reconstruction error of the data points.

Concept factorization attempts to find the approximation

$$\mathbf{Y}=\mathbf{Y}\mathbf{W}{\mathbf{V}}^T\mathrm{subject} \mathrm{to}\mathbf{W}\ge 0,\mathbf{V}\ge 0,$$

(45)

where $\mathbf{W}=\left[w_{jc}\right]\in {\mathbb{R}}^{n\times k}$ and $\mathbf{V}=\left[v_{jc}\right]\in {\mathbb{R}}^{n\times k}$.

The model is formulated as

$$\underset{\mathbf{W},\mathbf{V}}{min}\frac{1}{2}{∥\mathbf{Y}-\mathbf{Y}\mathbf{W}{\mathbf{V}}^T∥}_F^2.$$

(46)

The multiplicative update rule is given as [388]

$$w_{jc}⟵w_{jc}\frac{{\left(\mathbf{K}\mathbf{V}\right)}_{jc}}{{\left(\mathbf{V}{\mathbf{W}}^T\mathbf{K}\mathbf{W}\right)}_{jc}},v_{jc}⟵v_{jc}\frac{{\left(\mathbf{K}\mathbf{W}\right)}_{jc}}{{\left(\mathbf{K}\mathbf{W}{\mathbf{V}}^T\mathbf{V}\right)}_{jc}},$$

(47)

where $\mathbf{K}={\mathbf{Y}}^T\mathbf{Y}$. The inner product representation of $\mathbf{K}$ allows for the easy kernelization of concept factorization.

Pairwise constrained concept factorization [389] is a semi-supervised method that incorporates pairwise constraints. Data points with pairwise must-link and cannot-link constraints have the same class labels and different class labels, respectively. The locally consistent concept factorization algorithm [390] constructs a nearest-neighbor graph to characterize the local manifold structure of the data space. Label information can be encoded into the graph directly.

7. Symmetric Positive Semi-Definite Matrix Approximation

Rank-revealing QR and truncated SVD are two standard approaches for low-rank approximation. PCA is a truncated SVD applied on recentered data. Truncated SVD attempts to minimize

$$\underset{\mathbf{A},\mathbf{X}}{min}{\parallel \mathbf{Y}-\mathbf{A}\mathbf{X}\parallel }_F^2,\mathrm{subject} \mathrm{to}{\mathbf{A}}^T\mathbf{A}={\mathbf{I}}_r,$$

(48)

where r is a pre-specified rank of approximation.

A standard eigenvalue decomposition has a computational complexity of $O\left(n^3\right)$. The symmetric positive semidefinite matrix approximation is usually used to speed up the large-scale eigenvalue computation and kernel learning. A typical low-rank decomposition is given by

$$\mathbf{K}=\mathbf{C}\mathbf{U}{\mathbf{C}}^T,$$

(49)

where $\mathbf{C}\in {\mathbb{R}}^{n\times c}$ is a sketch of $\mathbf{K}\in {\mathbb{R}}^{n\times n}$ (e.g., c randomly sampled columns of $\mathbf{K}$) and $\mathbf{U}\in {\mathbb{R}}^{c\times c}$ is obtained by

$$\underset{\mathbf{U}}{min}{\parallel \mathbf{K}-\mathbf{C}\mathbf{U}{\mathbf{C}}^T\parallel }_F^2.$$

(50)

It needs $O\left(n{c}^2\right)$ computations to approximate the eigenvalue decomposition or matrix inversion for a rank-k ($k\le c$) matrix.

The Nyström method [391,392] is a popular alternative method that originates from solving integral equations. It is a sampling-based method that can efficiently approximate large kernel matrices and their eigen systems. However, the method is inaccurate. The analysis of the Nyström method is based on the Monte Carlo algorithm. In machine learning, the Nyström method is heavily utilized to estimate large kernel matrices [391,393,394].

The Nyström method selects $c\le n$ columns from $\mathbf{K}$, and then forms a low-rank approximation of the full matrix by using the correlations between the sampled and remaining columns. The method is highly scalable. It needs to decompose a $c\times c$ matrix, formed by intersecting the selected columns and the corresponding rows. The accuracy of the approximation is decided by the number of sampled columns.

Let $\mathbf{P}\in {\mathbb{R}}^{n\times c}$ be a sketching matrix via uniform sampling. The Nyström method computes [391]

$$\mathbf{C}=\mathbf{K}\mathbf{P},\mathbf{U}={\left({\mathbf{P}}^T\mathbf{C}\right)}^{†}.$$

(51)

where † is a pseudoinverse operator. $\mathbf{C}$ can be formed by sampling $c=O(k/\epsilon )$ columns of $\mathbf{K}$, such that

$$\underset{\mathbf{U}}{min}\parallel \mathbf{K}-\mathbf{C}\mathbf{U}{\mathbf{C}}^T{\parallel }_F^2\le (1+\epsilon ){\parallel \mathbf{K}-{\mathbf{K}}_k\parallel }_F^2,$$

(52)

where ${\mathbf{K}}_k$ is the estimated $\mathbf{K}$ of rank k.

Randomized SVD for symmetric matrices [395], as a standard sketch-based method, also extends the Monte Carlo algorithm. Randomized SVD constructs a low-dimensional subspace of the input matrix, and then performs standard decomposition techniques, such as QR or SVD, on the subspace matrix. At least one pass over $\mathbf{K}$ is required. A much higher accuracy is achieved by solving $\mathbf{U}$ from (50), but at a computational complexity of $O\left(n^{2}c\right)$. Unlike the Nyström method, randomized SVD does not require c to grow with n.

The Nyström scheme proposed in [396] is as accurate as the Nyström method and has a computational complexity that is as low as the SVD of a small matrix. It first samples many columns from $\mathbf{K}$, then an approximate SVD that is based on the randomized low-rank matrix approximation is performed on the inner submatrix. The fast symmetric positive semidefinite matrix approximation model [397] approximates $\mathbf{U}$ with a computational complexity linear in n. It is as efficient as the Nyström method and as accurate as randomized SVD.

8. CX Decomposition and CUR Decomposition

In general, for large sparse matrices, the SVD or QR decomposition does not preserve sparsity. The basis vectors generated by SVD have no physical interpretation as well. It is important to find a low-rank matrix decomposition that preserves the sparse property of the original matrix. Matrix column selection and CUR decomposition are techniques used to represent a matrix by using a few columns and/or rows of the matrix. Matrix column selection is a column-based method, while CUR decomposition is a column-row-based one [398].

The matrix column selection method is also known as CX decomposition [398,399]. In CX decomposition, a data matrix $\mathbf{Y}$ composed of sample vectors is factorized as

$$\mathbf{Y}=\mathbf{C}\mathbf{X},$$

(53)

where $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$ is a large sparse matrix, $\mathbf{C}\in {\mathbb{R}}^{m\times c}$ is a submatrix of $\mathbf{Y}$ and is also sparse, and $\mathbf{X}\in {\mathbb{R}}^{c\times n}$ is a coefficient matrix, which is not sparse in general.

For term-document data and binary image data, $\mathbf{Y}$ has sparse and nonnegative columns. Due to prototype-preserving properties of the CX decomposition, the columns of $\mathbf{C}$ are also sparse and nonnegative. The CX decomposition randomly samples a few columns of $\mathbf{Y}$ by a non-uniform probability derived from SVD. The deterministic CX decomposition [400] selects columns in a deterministic manner and it approximates SVD very well. Each selected column is associated with an eigenvector of PCA.

The CUR decomposition offers a representation of data in relation to other data, facilitating the interpretation of results. It has been used in exploratory data analyses related to natural language processing [401] and subspace clustering [402]. CUR is also used as a fast approximation to SVD [394,398,403,404], and is used to accelerate algorithms for robust PCA [405,406]. The Nyström method is the CUR decomposition where the same columns and rows are selected to approximate symmetric positive semidefinite matrices.

The CUR matrix decomposition [392], sometimes referred to as the pseudoskeleton decomposition [407,408], seeks to find a subset of c columns of $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$ to form $\mathbf{C}\in {\mathbb{R}}^{m\times c}$, a subset of r rows of $\mathbf{Y}$ to form $\mathbf{R}\in {\mathbb{R}}^{r\times n}$, and computes $\mathbf{U}\in {\mathbb{R}}^{c\times r}$, such that

$$\underset{\mathbf{C},\mathbf{U},\mathbf{R}}{min}{\parallel \mathbf{Y}-\mathbf{C}\mathbf{U}\mathbf{R}\parallel }_F^2.$$

(54)

The computational complexity for the optimal $\mathbf{U}={\mathbf{C}}^{†}\mathbf{Y}{\mathbf{R}}^{†}$ is $O(mnmin\{c,r\left\}\right)$.

The CUR decomposition gives factorizations in terms of actual columns and row submatrices, and it can be cheap to form through random sampling. For matrices, uniform random sampling is known to provide good CUR approximations under incoherence assumptions [408]. Additionally, several standard randomized sampling procedures give exact CUR decompositions with high probability [409].

The CUR decomposition usually chooses many columns and rows. A standard column selection procedure is first performed. A two-stage randomized CUR algorithm, referred to as the subspace sampling algorithm [398], selects columns and rows based on their statistical leverage scores. The computational complexity is higher than that of the truncated SVD of $\mathbf{Y}$. The subspace sampling algorithm can be accelerated by a fast approximation to statistical leverage scores [410].

The Nyström method approximates a symmetric positive semidefinite matrix using a few columns. In contrast, the CUR decomposition approximates an arbitrary matrix by a certain number of columns and rows. Hence, the CUR decomposition extends the standard sketch-based method from symmetric matrices to general matrices, with the same computational problem. The upper bounds of some CUR algorithms are even much better than the lower error bounds of the Nyström method and ensemble-based Nyström method. More accurately, low-complexity CUR and Nyström algorithms have been derived from a more general error bound for adaptive column/row sampling [411]. In [397], a fast symmetric positive semidefinite matrix approximation was implemented on the CUR decompositions of general matrices. The computational complexity decreases to $O(cr\frac{1}{\epsilon }min\{m,n\}min\{c,r\})$, but with nearly the same approximation quality.

9. Conclusions

Low-rank matrix factorization is a fundamental method in signal processing, statistics, data analysis, and machine learning. It is a popular tool used for dimension reduction, feature extraction, data compression, data mining, and information retrieval when processing high-dimensional real-world data. In this paper, we provide a comprehensive, state-of-the-art review of recent matrix factorization methods, including compressed sensing, sparse coding, dictionary learning, matrix completion, low-rank approximation as an extension of sparse coding, NMF, symmetric positive semidefinite matrix approximation, CX decomposition, and CUR decomposition. All of these matrix factorization methods have primarily been developed within the past two decades, and their iterative implementations are usually viewed as unsupervised learning approaches in the machine learning field.

In this paper, the reviewed algorithms are generally derived by optimizing objective functions using traditional iterative optimization methods. The objective functions are usually defined as $L_2$-norms, subject to some constraints, such as sparse or low-rank constraints, which can be solved using a classic LP method or a Lagrange multiplier method. The derived methods can be easily implemented, and the performance is guaranteed theoretically.

Traditional matrix factorization methods, such as SVD, PCA, and ICA, are not described in this paper. There are many literature surveys on these topics [1]. This paper is unique since it brings together so many popular and related topics under the thread of matrix factorization, while previous surveys typically focused on specific topics or their applications.

Since the prevalence of compressed sensing and dictionary learning, low-rank approximation has become a primary approach for feature extraction in signal processing and machine learning. In this paper, we collected recent methods on matrix factorization and matrix decomposition, as well as their theoretical advances. Their applications to various fields are not reviewed in this paper. The idea of incorporating deep learning into low-rank matrix factorization is not dealt with in this paper. The Bayesian approach and probabilistic analyses, such as the Monte Carlo method, are also not reviewed in this paper. Interested readers can retrieve the literature or conduct their own surveys on these themes and their applications in specific fields.

This manuscript focuses on matrix factorization and matrix decomposition based on the low-rank assumption, with importance attached to algorithms and theoretical results. More general extensions to the tensor case are not treated in this manuscript. A tensor is a multidimensional array that serves as a higher-order generalization of vectors and matrices. The order of a tensor is referred to as the number of modes it possesses. High-dimensional data, especially colored images, videos, and multisensor networks, are conveniently represented by tensors. Due to space limitations, tensor factorization, tensor compressed sensing, tensor completion, and nonnegative tensor factorization are not introduced in this paper. For more on tensor factorization, refer to [4].

9.1. Optimization by Metaheuristics or Neurodynamics

In recent years, some algorithms have solved formulated optimization problems by using nature-inspired metaheuristics [412,413] or neurodynamics (i.e., recurrent neural networks), and their combinations. These methods are not reviewed in this article. Here, we describe a few representative methods published in recent years. In [414], sparse signal reconstruction was achieved through collaborative neurodynamic optimization; this approach involves using a population of recurrent neural networks operating concurrently for a scattered search of individual solutions, combined with particle swarm optimization for repeated repositioning. In [415], Boolean matrix factorization is solved through a collaborative neurodynamic approach, which uses a population of Boltzmann machines for a scattered search of factorization solutions and particle swarm optimization for re-initializing the Boltzmann machines upon local convergence. Some other examples of neurodynamics-based methods for sparse signal reconstruction include a smoothing neurodynamic neural network modeled [416] for $L_p$-norm $2\ge p\ge 1$, a projected neurodynamic neural network for $L_0$-norm [417], and a Lagrange programming neural network with $L_p$-norm [418]. Similarly, a discrete-time projection neural network for sparse NMF is presented in [419].

9.2. A Few Topics for Future Research

All of the matrix factorization or decomposition techniques discussed in this paper are linear dimensionality reduction techniques used for data analysis. A future research topic will revolve around sparse recovery from nonlinear measurements, which are prevalent in practical physical systems. This presents a challenging task that often involves discrete optimization. When the entries of a matrix are generated from nonlinear transformations of lower dimensional latent subspaces, the matrix always has a high or even a full rank. Nonlinear matrix completion can be used to recover missing entries of such data matrices. In [420], the rank of a matrix in the feature space, defined by a kernel-trick-based nonlinear mapping of the data space, is approximated using the Schatten p-norm, and is minimized.

Existing solutions to low-rank matrix completion assume uniformly random observation patterns. An open issue is how to identify matrix patterns with unique or a finite number of completions. In [421], three families of matrix patterns are presented for low-rank matrix completion (in terms of Plücker coordinates). In [422], a deterministic sampling method for matrix completion using an asymmetric Ramanujan graph and its sufficient conditions for the matrix completion are derived. For the matrix completion under deterministic sampling, compared to uniform sampling, two weaker but necessary conditions, namely isomeric conditions and relative well-conditionedness, guarantee that any arbitrary matrix can be recovered [423]. Isomeric dictionary pursuit is a method based on the Schatten p-norm for this purpose [423]. We expect more research on matrix completion with deterministic sampling.

Most of the above methods are implemented on real matrices. Extensions of these methods to the complex domains are not straightforward. In addition to amplitude, phase information has to be considered. Extending the matrices to the complex domains can lead to interesting results; this deserves further investigation. In [424], phase-only compressive sensing estimates the signal direction from only the phases of complex measurements; this is achieved by normalizing the signal by its $L_2$-norm. This is a natural extension of one-bit compressive sensing [193,196,197]. In [425], based on the RIP analysis, a low-complexity signal can be perfectly recovered with high probability from phase-only complex random observations for scenarios involving a complex Gaussian random sensing matrix and a large number of measurements, relative to the complexity level of the signal space, by using any optimal algorithm in the compressive sensing literature. This recovery requires approximately twice the number of measurements needed for its compressive sensing counterpart.

In addition to sparsity, other prior information of $\overrightarrow{x}$ can be exploited. Under some unconventional constraints, lower bounds on the probability of exact recovery of $\overrightarrow{x}$ using OMP in k iterations, and a lower bound on the number of measurements M, guaranteeing that the probability of the exact recovery of $\overrightarrow{x}$ using OMP in k iterations is greater than a given probability, are derived in [426]. The paper also addresses the recovery of a signal or matrix with unconventional constraints, along with its theoretical analysis.

The matrix factorization techniques surveyed in this paper mainly emerged from the domains of machine learning, signal processing, and statistics in the last two decades. They can be solved by using classical mathematical methods. There are more emerging topics in the field, and it is anticipated that further research will yield more results in these areas in the future.