RobOMP: Robust variants of Orthogonal Matching Pursuit for sparse representations

MSE–based OMP

Orthogonal Matching Pursuit (Tropp & Gilbert, 2007) attempts to find the locally optimal solution by iteratively estimating the most correlated atom in D to the current residual. In particular, OMP initializes the residual r₀ = y, the set containing the indices of the atoms that are part of the decomposition (an active set) Λ₀ = ∅, and the iteration k = 1. In the kth iteration, the algorithm finds the predictor most correlated to the current residual: (2)${\displaystyle {\lambda }_k=\underset{i\in \Omega }{argmax}\left|〈{\mathbf{r}}_{k-1},{\mathbf{d}}_i〉\right|}$ where 〈⋅, ⋅〉 denotes the inner product operator, d_i represents the ith column of D, and Ω = {1, 2, …, n}. The resulting atom is added to the active set via Λ, i.e.: (3)${\displaystyle {\Lambda }_k={\Lambda }_{k-1}\cup {\lambda }_k.}$ The next step is the major refinement of the original Matching Pursuit algorithm (Mallat & Zhang, 1993)—instead of updating the sparse decomposition one component at the time, OMP updates all the coefficients corresponding to the active set at once according to a MSE criterion (4)${\displaystyle {\mathbf{x}}_k=\underset{\mathbf{x}\in {\mathrm{IR}}^n,\text{supp}\left(\mathbf{x}\right)\subset {\Lambda }_k}{argmin}\left|\right|\mathbf{y}-\mathbf{Dx}|{|}_2}$ where supp (x) is the support set of vector x. Equation (4) can be readily solved via OLS or Linear Regression where the predictors are the columns of D indexed by Λ_k and the response is the measurement vector y. Stopping criterion for OMP typically include a set number of iterations or compliance with a set minimum error of the residue. In the end, the estimated sparse code, x, is set as the last x_k obtained.

In practice, the true sparsity pattern, K₀, is unknown and the total number of OMP iterations, K, is treated as a hyperparameter. For a detailed analysis regarding convergence and recovery error bounds of OMP, see Donoho, Elad & Temlyakov (2006). A potential drawback of OMP is the extra computational complexity added by the OLS solver. Specifically, each incremental update of the active set affects the time complexity of the algorithm in a polynomial fashion: $\mathcal{O}\left(k^{2}n+k^3\right)$ where k is the current iteration.

Generalized Orthogonal Matching Pursuit (Wang, Kwon & Shim, 2012) refines OMP by selecting N₀ atoms per iteration. If the indices of the active set columns in the kth iteration are denoted as J_k[1], J_k[2], …, J_k[N₀], then J_k[j] can be defined recursively: (5)${\displaystyle J_k\left[j\right]=\underset{i\in \Omega \setminus \left\{J_k\left[1\right],\dots ,J_k\left[j-1\right]\right\}}{argmax}\left|〈{\mathbf{r}}_{k-1},{\mathbf{d}}_i〉\right|,1\le j\le N_0}$ The index set ${\left\{J_k\left[j\right]\right\}}_{j=1}^{N_0}$ is then added to Λ_k and, likewise OMP, GOMP exploits an OLS solver to update the current active set. Both OMP and GOMP obtain locally optimal solutions under the assumption of Gaussianity (or Normality) of the errors. Yet, if such restriction is violated (e.g., by the presence of outliers), the estimated sparse code, x, will most likely be biased.

CMP

The main drawback of MSE–based cost functions is its weighing nature in terms of influence and importance assigned to the available samples. In particular, MSE considers every sample as equally important and assigns a constant weight equal to one to all the inputs. Wang, Tang & Li (2017) proposed exploiting Correntropy (Liu, Pokharel & Príncipe, 2007) instead of MSE as an alternative cost function in the greedy sparse coding framework. Basically, the novel loss function utilizes the Correntropy Induced Metric (CIM) to weigh samples according to a Gaussian kernel g_σ(t) = exp(−t²∕2σ²), where σ, the kernel bandwidth, modulates the norm the CIM will mimic, e.g., for small σ, the CIM behaves similar to the ℓ₀-pseudonorm (aggressive non–linear weighing), if σ increases, CIM will mimic the ℓ₁–norm (moderate linear weighing), and, lastly, for large σ, the resulting cost function defaults to MSE, i.e., constant weighing for all inputs. The main conclusion here is that the CIM, unlike MSE, is robust to outliers for a principled choice of σ. This outcome easily generalizes for non–Gaussian environments with long–tailed distributions on the errors.

Correntropy Matching Pursuit (CMP) exploits the CIM robustness to update the active set in the sparse coding solver. The algorithm begins in a similar fashion as OMP, i.e., r₀ = y, Λ₀ = ∅, and k = 1. Then, instead of the MSE–based update of Eq. (4), CMP proceeds to minimize the following CIM–based expression: (6)${\displaystyle {\mathbf{x}}_k=\underset{\mathbf{x}\in {\mathrm{IR}}^n,\text{supp}\left(\mathbf{x}\right)\subset {\Lambda }_k}{argmin}{L}_{\sigma }\left(\mathbf{y}-\mathbf{Dx}\right)}$ where $L_{\sigma }\left(\mathbf{e}\right)=\frac{1}{m}{\sum }_{i=1}^m{\sigma }^2\left(1-g_{\sigma }\left(\mathbf{e}\left[i\right]\right)\right)$ is the simplified version (without constant terms independent of e) of the CIM loss function and e[i] is the ith entry of the vector e. The Half–Quadratic (HQ) technique is utilized to efficiently optimize the invex CIM cost function (Geman & Yang, 1995; Nikolova & Ng, 2005). The result is a local minimum of Eq. (6) alongside a weight vector w that indicates the importance of the components of the observation vector y: (7)${\displaystyle {\mathbf{w}}^{\left(t+1\right)}\left[i\right]=g_{\sigma }\left(\mathbf{y}\left[i\right]-\left({\mathbf{D}\mathbf{x}}^{\left(t\right)}\right)\left[i\right]\right),i=1,2,\dots ,m}$ where t is the iteration in the HQ subroutine. In short, the HQ optimization performs block coordinate descent to separately optimize the sparse code, x, and the weight vector, w, in order to find local optima. The hyperparameter σ is iteratively updated without manual selection according to the following heuristic: (8)${\displaystyle {\sigma }^{\left(t+1\right)}={\left(\frac{1}{2m}{∥\mathbf{y}-{\mathbf{D}\mathbf{x}}^{\left(t+1\right)}∥}_2^2\right)}^{1∕2}.}$ In Wang, Tang & Li (2017), the authors thoroughly illustrate the advantage of CMP over many MSE–based variants of OMP when dealing with non-Gaussian error distributions and outliers in computer vision applications. And even though they mention the improved performance of the algorithm when σ is iteratively updated versus manual selection scenarios, they fail to explain the particular heuristic behind Eq. (8) or its statistical validity. In addition, the HQ optimization technique is succinctly reduced to a weighted Least Squares problem than can be solved explicitly. Therefore, more principled techniques that exploit weighted Least Squares and robust estimators for linear regression can readily provide the needed statistical validity, while at the same time, generalize the concepts of CMP under the umbrella of M–estimators.

M–Estimators

If the errors are not normally distributed, the estimator of Eq. (12) will be suboptimal. Hence, a different function is utilized to model the statistical properties of the errors. Following the same premises of independence and equivalence of the optimum under the log–transform, the following estimator arises: (13)${\displaystyle {\hat{\beta }}_{M-Est}=\underset{\beta }{argmin}\sum _{i=1}^m\rho \left(\right.\frac{\mathbf{e}\left[i\right]}{s}\left)\right.=\underset{\beta }{argmin}\sum _{i=1}^m\rho \left(\right.\frac{\left(\mathbf{y}\left[i\right]-\left(\Phi \beta \right)\left[i\right]\right)}{s}\left)\right.}$ where ρ(e) is a continuous, symmetric function (also known as the objective function) with a unique minimum at e = 0 (Andersen, 2008). Clearly, ρ(e) reduces to half the sum of squared errors for the Gaussian case. s is an estimate of the scale of the errors in order to guarantee scale–invariance of the solution. The usual standard deviation cannot be used for s due to its non–robustness; thus, a suitable alternative is usually the “re–scaled MAD”: (14)${\displaystyle s=1.4826\times MAD}$ where the MAD (median absolute deviation) is highly resistant to outliers with a breakdown point (BDP) of 50%, as it is based on the median of the errors ($\tilde{\mathbf{e}}$) rather than their mean (Andersen, 2008): (15)${\displaystyle MAD=\text{median}|\mathbf{e}\left[i\right]-\tilde{\mathbf{e}}|.}$ The re–scaling factor of 1.4826 guarantees that, for large sample sizes and $\mathbf{e}\left[i\right]\sim \mathcal{N}\left(0,{\sigma }^2\right)$, s reduces to the population standard deviation (Hogg, 1979). M–Estimation then, likewise OLS, finds the optimal coefficients vector via partial differentiation of Eq. (13) with respect to each of the k parameters in question, resulting in a system of k equations: (16)${\displaystyle \sum _{i=1}^m{\Phi }_{ij}\psi \left(\right.\frac{\mathbf{y}\left[i\right]-{\varphi }_i^T\beta }{s}\left)\right.=\sum _{i=1}^m{\Phi }_{ij}\psi \left(\right.\frac{\mathbf{e}\left[i\right]}{s}\left)\right.=0,j=1,2,\dots ,k}$ where ϕ_i represents the ith row of the matrix Φ while Φ_ij accesses the jth component of the ith row of $\Phi .\psi \left(\right.\frac{\mathbf{e}\left[i\right]}{s}\left)\right.=\frac{\partial \rho }{\partial \frac{\mathbf{e}\left[i\right]}{s}}$ is known as the score function while the weight function is derived from it as: (17)${\displaystyle \mathbf{w}\left[i\right]=\mathbf{w}\left(\right.\frac{\mathbf{e}\left[i\right]}{s}\left)\right.=\frac{\psi \left(\right.\frac{\mathbf{e}\left[i\right]}{s}\left)\right.}{\frac{\mathbf{e}\left[i\right]}{s}}.}$ Substituting Eqs. (17) into (16) results in: (18)${\displaystyle \sum _{i=1}^m{\Phi }_{ij}\mathbf{w}\left[i\right]\frac{\mathbf{e}\left[i\right]}{s}=\sum _{i=1}^m{\Phi }_{ij}\mathbf{w}\left[i\right]\left(\mathbf{y}\left[i\right]-{\varphi }_i^T\beta \right)\frac{1}{s}=0j=1,2,\dots ,k\sum _{i=1}^m{\Phi }_{ij}\mathbf{w}\left[i\right]\left(\mathbf{y}\left[i\right]-{\varphi }_i^T\beta \right)\frac{1}{s}=0j=1,2,\dots ,k\sum _{i=1}^m{\Phi }_{ij}\mathbf{w}\left[i\right]{\varphi }_i\beta =\sum _{i=1}^m{\Phi }_{ij}\mathbf{w}\left[i\right]\mathbf{y}\left[i\right]j=1,2,\dots ,k}$ which can be succinctly reduced in matrix form as: (19)${\displaystyle {\Phi }^T\mathbf{W}\Phi \beta ={\Phi }^T\mathbf{W}\mathbf{y}}$ by defining the weight matrix, W, as a square diagonal matrix with non–zero elements equal to the entries in w, i.e.: W = diag({w[i]:i = 1, 2, …, m}). Lastly, if Φ^TWΦ is well-conditioned, the closed form solution of the robust M–Estimator is equal to: (20)${\displaystyle {\hat{\beta }}_{M-Est}={\left({\Phi }^T\mathbf{W}\Phi \right)}^{-1}{\Phi }^T\mathbf{Wy}.}$ Equation (20) resembles its OLS counterpart (Eq. (12)), except for the addition of the matrix W that weights the entries of the observation vector and mitigates the effect of outliers according to a linear fit. A wide variety of objective functions (and in turn, weight functions) have been proposed in the literature (for a thorough review, see Zhang (1997)). For the present study, we will focus on five different variants that are detailed in Table 1. Every M–Estimator weighs its entries according to a symmetric, decaying function that assigns large weights to errors in the vicinity of zero and small coefficients to gross contributions. Consequently, the estimators downplay the effect of outliers and samples, in general, that deviate from the main mode of the residuals.

Table 1:

Comparison between OLS estimator and M–Estimators.

Objective ρ(e) and weight w(e) functions of OLS solution and five different M–Estimators. For M–Estimators, error entries are standardized, i.e. divided by the scale estimator, s. Each robust variant comes with a hyperparameter c. Exemplary plots in the last column utilize the optimal hyperparameters detailed in Table 2.

Type	ρ(e)	w(e)	w(e)
OLS	$\frac{1}{2}{e}^2$	1
Cauchy	$\frac{c^2}{2}log\left(1+{\left(\frac{e}{c}\right)}^2\right)$	$\frac{1}{1+{\left(e∕c\right)}^2}$
Fair	$c^2\left(\right.\frac{\left|e\right|}{c}-log\left(\right.1+\frac{\left|e\right|}{c}\left)\right.\left)\right.$	$\frac{1}{1+\left|e\right|∕c}$
Huber  ${\displaystyle \left\{\right.\begin{array}{c}\text{if}\left|e\right|\le c\hfill \\ \text{if}\left|e\right|\ge c\hfill \\ \hfill \end{array}}$	${\displaystyle \left\{\right.\begin{array}{c}\frac{e^2}{2}\hfill \\ c\left(\right.\left|e\right|-\frac{c}{2}\left)\right.\hfill \\ \hfill \end{array}}$	${\displaystyle \left\{\right.\begin{array}{c}1\hfill \\ \frac{c}{\left|e\right|}\hfill \\ \hfill \end{array}}$
Tukey ${\displaystyle \left\{\right.\begin{array}{c}\text{if}\left|e\right|\le c\hfill \\ \text{if}\left|e\right|>c\hfill \\ \hfill \end{array}}$	${\displaystyle \left\{\right.\begin{array}{c}\frac{c^2}{6}\left(\right.1-\left(\right.1-\left(\right.\frac{e}{c}{\left)\right.}^2{\left)\right.}^3\left)\right.\hfill \\ \frac{c^2}{6}\hfill \\ \hfill \end{array}}$	${\displaystyle \left\{\right.\begin{array}{c}\left(\right.1-{\left(\frac{e}{c}\right)}^2{\left)\right.}^2\hfill \\ 0\hfill \\ \hfill \end{array}}$
Welsch	$\frac{c^2}{2}\left(1-exp\left(-{\left(\frac{e}{c}\right)}^2\right)\right)$	$exp\left(-{\left(\frac{e}{c}\right)}^2\right)$

DOI: 10.7717/peerjcs.192/table-1

Solving the M-Estimation problem is not as straightforward as the OLS counterpart. In particular, Eq. (20) assumes the optimal W is readily available, which, in turn, depends on the residuals, which, again, depends on the coefficient vector. In short, the optimization problem for M–Estimators can be posed as finding both ${\hat{\beta }}_{M-Est}$ and ${\hat{\mathbf{w}}}_{M-Est}$ that minimize Eq. (13). Likewise Half–Quadratic, the common approach is to perform block coordinate descent on the cost function with respect to each variable individually until local optima are found. In the robust regression literature, this optimization procedure is the well known Iteratively Reweighted Least Squares or IRLS (Andersen, 2008). The procedure is detailed in Algorithm 1. In particular, the routine runs for either a fixed number of iterations or until the estimates change by less than a selected threshold between iterations. The main hyperparameter is the choice of the robust M–Estimator alongside its corresponding parameter c. However, it is conventional to select the value that provides a 95% asymptotic efficiency on the standard Normal distribution (Zhang, 1997). Throughout this work, we exploit such optimal values to avoid parameter tuning by the user (see Table 2). In this way, the influence of outliers and non-Gaussian errors are expected to be diminished in the OMP update stage of the coefficients corresponding to the active set.

 
__________________________________________________________________________________________________ 
Algorithm 1 IRLS–based M–Estimation 
__________________________________________________________________________________________________ 
  1:    function IRLS(y ∈ IRm, Φ ∈ IRm × k, wc(u))   ⊳ Weight function w(u) with 
     hyperparameter c 
  2:       t ← 0 
  3:       β(0) = βOLS ← (ΦTΦ)−1ΦTy                          ⊳ OLS initialization 
  4:       e(0) ← y − Φβ(0) 
  5:       MAD ← median|e(0)[i] − ˜ e(0)| 
  6:       s(0) ← 1.4826 × MAD 
  7:       w(0)[i] ← wc(e(0)[i] s(0)  )       i = 1,2,...,m ⊳ Initial weight vector 
  8:       W(0) ← diag(w(0)) 
  9:       t ← 1 
 10:       while NO CONVERGENCE do 
 11:            β(t) ← (ΦTW(t−1)Φ)−1ΦTW(t−1)y       ⊳ Block coordinate descent 
 12:            e(t) ← y − Φβ(t) 
 13:            MAD ← median|e(t)[i] − ˜ e(t)| 
 14:            s(t) ← 1.4826 × MAD 
 15:            w(t)[i] ← wc(e(t)[i] s(t)  )i = 1,2,...,m ⊳ Block coordinatedescent 
 16:            W(t) ← diag( w(t)) 
 17:            t ← t + 1 
 18:       return ˆ βM–Est ← β(t)     ˆ wM–Est ← w(t)                         ⊳ Final estimates    

M–Estimators–based OMP

Here, we combine the ideas behind greedy approximations to the sparse coding problem and robust M–Estimators; the result is RobOMP or Robust Orthogonal Matching Pursuit. We propose five variants based on five different M–Estimators (Table 1). We refer to each RobOMP alternative according to its underlaying M–Estimator; for instance, Fair–Estimator–based–OMP is simply referred to as Fair. As with OMP, the only hyperparameter is the stopping criterion: either K as the maximum number of iterations (i.e., sparseness of the solution), or ϵ, defined as a threshold on the error norm.

For completeness, Algorithm 2 details the RobOMP routine for the case of set maximum number of iterations (the case involving ϵ is straightforward). Three major differences are noted with respect to OMP:

1. The robust M–Estimator–based update stage of the active set is performed via IRLS,

2. The updated residuals are computed considering the weight vector ${\hat{\mathbf{w}}}_k$ from IRLS, and

3. The weight vector constitutes an additional output of RobOMP.

The last two differences are key for convergence and interpretability, respectively. The former guarantees shrinkage of the weighted error in its first and second moments, while the latter provides an intuitive, bounded, m–dimensional vector capable of discriminating between samples from the main mode and potential outliers at the tails of the density.

 
_______________________________________________________________________________________________________ 
Algorithm 2 RobOMP 
_______________________________________________________________________________________________________ 
  1:  function RobOMP(y ∈ IRm,D ∈ IRm×n,wc(u),K) 
  2:       k ← 1                                                                  ⊳ Initializations 
  3:       r0 ← y 
  4:       Λ0 ←∅ 
  5:       while k < K do 
  6:            λk = argmaxi∈Ω |〈rk−1,di〉|        Ω = {1,2,⋅⋅⋅ ,n} 
  7:            Λk = Λk−1 ∪{λk} 
  8:            Φ = [dΛk[1],dΛk[2],⋅⋅⋅ ,dΛk[k]] 
  9:            {ˆβM–Est, ˆ wk}← IRLS(y,Φ,wc(u))                    ⊳ Robust linear fit 
  10:            xk[Λk[i]] ← ˆ βM–Est[i]                        i = 1,2,...,k       ⊳ Update 
     active set 
  11:            rk[i] ← ˆ wk[i] × (y[i] − (Dxk)[i])                  i = 1,2,...,m        ⊳ 
    Update residual 
  12:            k ← k + 1 
  13:       return xRobOMP ← xK,w ← ˆ wK                           ⊳ Final Estimates    

Sparse coding with synthetic data

The dictionary or observation matrix, D ∈ IR^100×500, is generated with independent entries drawn from a zero–mean Gaussian random variable with variance equal to one. The ideal sparse code, x₀ ∈ IR⁵⁰⁰, is generated by randomly selecting ten entries and assigning them independent samples from a zero–mean, unit–variance Gaussian distribution. The rest of the components are set equal to zero, i.e., K₀ = 10. The resulting observation vector y ∈ IR¹⁰⁰ is computed as the linear combination of the dictionary with weights from the ideal sparse code plus a noise component n ∈ IR¹⁰⁰: (21)${\displaystyle \mathbf{y}={\mathbf{Dx}}_0+\mathbf{n}.}$

The first set of experiments considers different noise distributions. In particular, five noise cases are analyzed: Gaussian ($\mathcal{N}\left(0,2\right)$), Laplacian with variance equal to 2, Student’s t–distribution with 2 degrees of freedom, Chi–squared noise with 1 degree of freedom, and Exponential with parameter λ = 1. Then, OMP, GOMP, CMP, and the 5 variants of RobOMP estimate the sparse code with parameter K = 10. For the active set update stage of CMP and RobOMP, the maximum allowed number of HQ/IRLS iterations is set to 100. For GOMP, N₀ ∈ {2, 3, 4, 5} where the best results are presented.

The performance measure is defined as the normalized ℓ₂–norm of the difference between the ground truth sparse code, x₀, and its estimate. The average results for 100 independent runs are summarized in Table 3. As expected, most of the algorithms perform similar under Gaussian noise, which highlights the adaptive nature of CMP and RobOMP. For the non–Gaussian cases, CMP and Tukey are major improvements over ordinary OMP. The rest of the RobOMP flavors consistently outperform the state of the art OMP and GOMP techniques. This confirms the optimality of MSE–based greedy sparse decompositions when the errors are Normally distributed; yet, they degrade their performance when such assumption is violated.

The second set of results deals with non–linear additive noise or instance–based degradation. Once again, D and x₀ are generated following the same procedure of the previous set of results (K₀ = 10). Yet now, noise is introduced by means of zeroing randomly selected entries in y. The number of missing samples is modulated by a rate parameter ranging from 0 to 0.5. Figure 2 summarizes the average results for K = 10 and 100 independent runs. As expected, the performance degrades when the rate of missing entries increases. However, the five variants of RobOMP are consistently superior than OMP and GOMP until the 0.4–mark. Beyond that point, some variants degrade at a faster rate. Also, CMP achieves small sparse code error norms for low missing entries rate; however, beyond the 0.25–mark, CMP seems to perform worse than OMP and even GOMP. This experiment highlights the superiority of RobOMP over MSE–based and Correntropy–based methods.

Now, the effect of the hyperparameter K is studied. Once again, 100 independent runs are averaged to estimate the performance measure. The rate of missing entries is fixed to 0.2 while K is the free variable. Figure 3 shows how the average error norm is a non–increasing function of K for the non–MSE–based variants of OMP (slight deviation in some cases beyond K = 8 might be due to estimation uncertainty and restricted sample size). On the other hand, both OMP and GOMP seem to stabilize after a certain number of iterations, resulting in redundant runs of the algorithm. These outcomes imply that RobOMP is not only a robust sparse code estimator, but also a statistically efficient one that exploits the available information in the data in a principled manner. It is also worth noting that CMP underperforms when compared to most flavors of RobOMP.

Impulsive noise is the other extreme of instance–based contamination. Namely, a rate of entries in y are affected by aggressive high–variance noise while the rest of the elements are left intact. The average performance measure of 100 independent runs is reported for K = K₀ = 10. Figure 4A details the results for varying rates of entries affected by −20 dB impulsive noise. Again, RobOMP and CMP outperform OMP and GOMP throughout the entire experiment. Tukey and Welsch seem to handle this type of noise more effectively; specifically, the error associated to the algorithms in question seem to be logarithmic or radical for OMP and GOMP, linear for Fair, Cauchy, Huber and CMP, and polynomial for Tukey and Welsch with respect to the percentage of noisy samples. On the other hand, Figure 4B reflects the result of fixing the rate of affected entries to 0.10 and modulating the variance of the impulsive noise in the range [−25,0]. RobOMP and CMP again outperform MSE–based methods (effect visually diminished due to log–transform of the performance measure for plotting purposes). For this case, CMP is only superior to the Fair version of RobOMP.

In summary, the experiments concerning sparse coding with synthetic data confirm the robustness of the proposed RobOMP algorithms. Non–Gaussian errors, missing samples and impulsive noise are handled in a principled scheme by all the RobOMP variants and, for most cases, the results outperform the Correntropy–based CMP. Tukey seems to be the more robust alternative that is able to deal with a wide spectrum of outliers in a consistent, efficient manner.

RobOMP–based classifier

We introduce a novel robust variant for sparse representation–based classifiers (SRC) fully based on RobOMP. Let ${\mathbf{A}}_i=\left[{\mathbf{a}}_1^i,{\mathbf{a}}_2^i,\dots ,{\mathbf{a}}_{n_i}^i\right]\in {\mathrm{IR}}^{m\times n_i}$ be a matrix with n_i examples from the ith class for i = 1, 2, …, N. Then, denote the set N = {1, 2, …, N} and the dictionary matrix of all training samples A = [A₁, A₂, …, A_N] ∈ IR^m×n where $n={\sum }_{i=1}^N{n}_i$ is the number of training examples from all N classes. Lastly, for each class i, the characteristic function δ_i:IRⁿ → IRⁿ extracts the coefficients associated with the ith label. The goal of the proposed classifier is to assign a class to a test sample y ∈ IR^m given the generative “labeled” dictionary A.

The classification scheme proceeds as follows: N different sparse codes are estimated via Algorithm 2 given the subdictionaries A_i for i = 1, 2, …, N. The class–dependent residuals, rⁱ(y) are computed and the test example is assigned to the class with minimal residual norm. To avoid biased solutions based on the scale of the data, the columns of A are set to have unit–ℓ₂–norm. The result is a robust sparse representation–based classifier or RSRC, which is detailed in Algorithm 3.

 
_______________________________________________________________________________________________________ 
Algorithm 3 RSRC 
______________________________________________________________________________________________________ 
Inputs: Normalized matrix of training samples A = [A1,A2,...,AN] ∈ IRm×n 
Test Example, y ∈ IRm 
M–Estimator weight function, wc(u) 
Stopping criterion for RobOMP, K 
Output: class(y) 
  1:  (xRobOMP,w) ← RobOMP(y,A,wc(u),K)  ⊳ Compute robust sparse code 
     and weight vector 
  2:  ri(y) = ||diag(w)(y − Aiδi(xRobOMP))||2,   i ∈ N      ⊳ Calculate norm of 
     class–dependent residuals 
  3:  class(y) ← argmini∈N ri(y)                                           ⊳ Predict label    

Similar algorithms can be deployed for OMP, GOMP and CMP (Wang, Tang & Li, 2017). In particular, the original SRC (Wright et al., 2009) exploits a ℓ₁–minimization approach to the sparse coding problem; however the fidelity term is still MSE, which is sensitive to outliers. In this section we opt for greedy approaches to estimate the sparse representation. Moreover for RobOMP, the major difference is the computation of the residual—we utilize the weight vector to downplay the influence of potential outlier components and, hence, reduce the norm of the errors under the proper dictionary. CMP utilizes a similar approach, but the weight matrix is further modified due to the HQ implementation (see Wang, Tang & Li (2017) for details). We compare the 7 SRC variants under two different types of noise on the Extended Yale B Database.

Occlusions and missing pixels

Two different types of noise are simulated: blocks of salt and pepper noise, i.e., occlusions, and random missing pixels. In all the following experiments, the sparsity parameter for learning the sparse code is set to K = 5 (for GOMP, N₀ ∈ {2, 3} and the best results are presented). Also, 10 different independent runs are simulated for each noise scenario.

For the occlusion blocks, a rate of affected pixels is selected beforehand in the range [0, 0.5]. Then, as in the original SRC (Wright et al., 2009), we downsampled the inputs mainly for computational reasons. In particular, we utilized factors of 1/2, 1/4, 1/8, and 1/16 resulting in feature dimensions of 8232, 2058, 514, and 128, respectively. Next, every test example is affected by blocks of salt and pepper noise (random pixels set to either 0 or 255). The location of the block is random and its size is determined by the rate parameter. Every sample is assigned a label according to SRC variants based on OMP and GOMP, CMP–based classifier (coined as CMPC by Wang, Tang & Li (2017)), and our proposed RSRC. For simplicity, we use the same terminology as before when it comes to the different classifiers. The performance metric is the average classification accuracy in the range [0,1]. Figure 5 highlights the superiority of RSRC over OMP and GOMP. Particularly, Huber, Tukey and Welsch are consistently better than CMP while Fair and Cauchy seem to plateau after the 0.3–mark.

Next, the effects of the feature dimension and the sparsity parameter are investigated. Figure 6 confirms the robustness of the proposed discriminative framework. As expected, when the feature dimension increases, the classification accuracy increases accordingly. However, the baselines set by OMP and GOMP are extremely low for some cases. On the other hand, CMP and RSRC outperform both MSE–based approaches, and even more, the novel M–Estimator–based classifiers surpass their Correntropy–based counterpart. When it comes to the sparsity parameter, K, it is remarkable how OMP and GOMP do not improve their measures after the first iteration. This is expected due to the lack of principled schemes to deal with outliers. In contrast, RSCR shows a non–decreasing relation between classification accuracy and K, which implies progressive refinement of the sparse code over iterations. To make these last two findings more evident, Table 4 illustrates the classification accuracy for a very extreme case: 0.3 rate of occlusion and feature dimension equal to 128, i.e., each input image is roughly 12 × 11 pixels in size (the downsampling operator introduces rounding errors in the final dimensionality). This scenario is very challenging and, yet, most of RSRC variants achieve stability and high classification after only four iterations. On the other hand, OMP and GOMP degrade their performance over iterations. This confirms the robust and sparse nature of the proposed framework.

Table 4:

Average classification accuracy on the Extended Yale B Database over K for a fixed rate of 0.3 pixels affected by blocks of salt and pepper noise.

Best result for each classifier is marked bold. Feature dimension = 128.

K	OMP	GOMP	CMP	Cauchy	Fair	Huber	Tukey	Welsch
1	0.38	0.36	0.59	0.53	0.55	0.53	0.51	0.51
2	0.39	0.34	0.90	0.81	0.83	0.87	0.86	0.87
3	0.39	0.30	0.97	0.88	0.88	0.98	0.98	0.98
4	0.37	0.28	0.97	0.88	0.89	0.99	0.99	0.99
5	0.36	0.28	0.97	0.88	0.88	0.99	0.99	0.99
6	0.34	0.28	0.97	0.88	0.88	0.99	0.99	0.99
7	0.34	0.28	0.97	0.88	0.88	0.98	0.99	0.99
8	0.33	0.28	0.97	0.88	0.88	0.98	0.99	0.99
9	0.32	0.28	0.97	0.88	0.88	0.98	0.99	0.99
10	0.31	0.28	0.97	0.88	0.88	0.98	0.99	0.98

DOI: 10.7717/peerjcs.192/table-4

For the missing pixels case, a rate of affected pixels is selected beforehand in the range [0, 1]. Then, every test example is affected by randomly selected missing pixels—the chosen elements are replaced by samples drawn from a uniform distribution over the range [0, y_max] where y_max is the largest possible intensity of the image in question. Figures 7 and 8 summarize similar experiments as in the occlusion case. Again, the RSRC are superior than MSE–based methods and consistently increase the performance measure as the sparsity parameter grows. The extreme case here involves a rate of 0.4 affected pixels by distorted inputs and a feature dimension of 128. Table 5 reinforces the notion that robust methods achieve higher classification accuracy even in challenging scenarios.

Lastly, it is worth noting that CMP performs better in the missing pixels case; yet, it fails to surpass the Welsch variant of RSRC which is its equivalent in terms of weight function of errors. Once again, Tukey is the algorithm with overall best results that is able to handle both kinds of noise distributions in a more principled manner.

Image denoising via robust, sparse and redundant representations

The last set of results introduces a preliminary analysis of image denoising exploiting sparse and redundant representations over overcomplete dictionaries. The approach is based on the seminal paper by Elad & Aharon (2006). Essentially, zero–mean white and homogeneous Gaussian additive noise with variance σ² is removed from a given image via sparse modeling. A global image prior that imposes sparsity over patches in every location of the image simplifies the sparse modeling framework and facilitates its implementation via parallel processing. In particular, if the unknown image Z can be devised as the spatial (and possibly overlapping) superposition of patches that can be effectively sparsely represented given a dictionary D, then, the optimal sparse code, ${\hat{\mathbf{x}}}_{ij}$, and estimated denoised image, $\hat{\mathbf{Z}}$, are equal to: (22)${\displaystyle \left\{{\hat{\mathbf{x}}}_{ij},\hat{\mathbf{Z}}\right\}=\underset{{\mathbf{x}}_{ij},\mathbf{Z}}{argmin}\lambda \left|\right|\mathbf{Z}-\mathbf{Y}|{|}_2^2+\sum _{ij}{\mu }_{ij}|\left|{\mathbf{x}}_{ij}\right|{|}_0+\sum _{ij}\left|\right|\mathbf{D}{\mathbf{x}}_{ij}-{\mathbf{R}}_{ij}\mathbf{Z}|{|}_2^2}$ where the first term is the log–likelihood component that enforces close resemblance (or proximity in an ℓ₂ sense) between the measured noisy image, Y, and its denoised (and unknown) counterpart Z. The second and third terms are image priors that enforce that every patch, z_ij = R_ijZ, of size $\sqrt{n}\times \sqrt{n}$ in every location of the constructed image Z has a sparse representation with bounded error. λ and μ_ij are regularization parameters than can easily be reformulated as constraints.

Block coordinate descent is exploited to solve Eq. (22). In particular, ${\hat{\mathbf{x}}}_{ij}$ is estimated via greedy approximations of the sparse code of each local block or patch. The authors suggest OMP with stopping criterion set by $\left|\right|\mathbf{D}{\mathbf{x}}_{ij}-{\mathbf{R}}_{ij}\mathbf{Z}|{|}_2^2\le {\left(C\sigma \right)}^2$ for all {ij} combinations (sequential sweep of the image to extract all possible $\sqrt{n}\times \sqrt{n}$ blocks). Then, the estimated denoised image has the following closed form solution: (23)${\displaystyle \hat{\mathbf{Z}}=\left(\right.\lambda I+\sum _{ij}{\mathbf{R}}_{ij}^T{\mathbf{R}}_{ij}{\left)\right.}^{-1}\left(\right.\lambda \mathbf{Y}+\sum _{ij}{\mathbf{R}}_{ij}^T\mathbf{D}{\mathbf{x}}_{ij}\left)\right.}$ where I is the identity matrix. The authors go one step further and propose learning the dictionary, D, as well; this is accomplished either from a corpus of high–quality images or the corrupted image itself. The latter alternative results in a fully generative sparse modeling scheme. For more details regarding the denoising mechanisms, refer to Elad & Aharon (2006).

For our case, we focus on the sparse coding subproblem alone and utilize an overcomplete Discrete Cosine Transform (DCT) dictionary, D ∈ IR^64×256, and overlapping blocks of size 8 × 8. The rest of the free parameters are set according to the heuristics presented in the original work: λ = 30∕σ and C = 1.15. Our major contribution is the robust estimation of the sparse codes via RobOMP in order to handle potential outliers in a principled manner. Two types of zero–mean, homogeneous, additive noise (Gaussian and Laplacian) are simulated with different variance levels on 10 independent runs. Each run comprises of separate contaminations of 4 well known images (Lena, Barbara, Boats and House) followed by the 7 different denoising frameworks, each one based on a distinct variant of OMP. As before, every algorithm is referred to as the estimator exploited in the active set update stage.

Tables 6 and 7 summarize the average performance measures (PSNR in dB) for 5 different variance levels of each noise distribution. As expected, OMP is roughly the best denoising framework for additive Gaussian noise. However, in the Laplacian case, Cauchy achieves higher PSNR levels throughout the entire experiment. This suggests the Cauchy M–Estimator is more suitable for this type of non–Gaussian environment. It is worth noting though that the averaging performed in Eq. (23) could easily blur the impact of the sparse code solvers for this particular joint optimization. Also, no attempt was made to search over the hyperparameter space of λ and C, which we suspect have different empirical optima depending on the noise distribution and sparse code estimator. These results are simply preliminary and highlight the potential of robust denoising frameworks based on sparse and redundant representations.