<figure><img src="image_1.png" /></figure>

HAL
open science

# Sparse signal subspace decomposition based on adaptive over-complete dictionary

Hong Sun, Cheng-Wei Sang, Didier Le Ruyet

► To cite this version:

Hong Sun, Cheng-Wei Sang, Didier Le Ruyet. Sparse signal subspace decomposition based on adaptive over-complete dictionary. EURASIP Journal on Image and Video Processing, 2017, 2017 (1), 10.1186/s13640-017-0200-7. http://hal-02470260

HAL Id: hal-02470260

https://cnam.hal.science/hal-02470260

Submitted on 26 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Sparse signal subspace decomposition based on adaptive over-complete dictionary

Hong Sun<sup>1,2*</sup> , Cheng-wei Sang<sup>1</sup> and Didier Le Ruyet<sup>3</sup>

**Abstract**

This paper proposes a subspace decomposition method based on an over-complete dictionary in sparse representation, called "sparse signal subspace decomposition" (or 3SD) method. This method makes use of a novel criterion based on the occurrence frequency of atoms of the dictionary over the data set. This criterion, well adapted to subspace decomposition over a dependent basis set, adequately reflects the intrinsic characteristic of regularity of the signal. The 3SD method combines variance, sparsity, and component frequency criteria into a unified framework. It takes benefits from using an over-complete dictionary which preserves details and from subspace decomposition which rejects strong noise. The 3SD method is very simple with a linear retrieval operation. It does not require any prior knowledge on distributions or parameters. When applied to image denoising, it demonstrates high performances both at preserving fine details and suppressing strong noise.

**Keywords:** Subspace decomposition, Sparse representation, Frequency of components, PCA, K-SVD, Image denoising

## 1 Introduction

Signal subspace methods (SSMs) are efficient techniques to reduce the dimensionality of data and to filter out noise [1]. The fundamental idea under SSM is to project the data on a basis made of two subspaces, one mostly containing the signal and the other the noise. The two subspaces are separated by a thresholding criterion associated with some measures of information.

The two most popular methods of signal subspace decomposition are wavelet shrinkage [2] and principal component analysis (PCA) [3]. Both techniques have proved to be quite efficient. However, wavelet decomposition depending on signal statistics is not equally adapted to different data and requires some knowledge on prior distributions or parameters of signals to efficiently choose the thresholds for shrinkage. A significant advantage of the PCA is its adaptability to data. The separation criterion is based on energy which may be seen as a limitation in some cases as illustrated in the next section.

In recent years, sparse coding has attracted significant interest in the field of signal denoising [4]. A sparse representation is a decomposition of a signal on a very small set of components of an over-complete basis (called dictionary) which is adapted to the processed data. A difficult aspect for signal subspace decomposition based on such a sparse representation is to define the most appropriate criterion to identify the principal components (called atoms) from the learned dictionary to build the principal subspace. The non-orthogonal property of the dictionary does not allow to use the energy criterion for this purpose, as done with PCA.

To solve this problem, we introduce a new criterion to measure the importance of atoms and propose a SSM under the criterion of the occurrence frequency of atoms. We thus make benefit both from the richness of over-complete dictionaries which preserves details of information and from signal subspace decomposition which rejects strong noise.

The remainder of this paper is organized as follows: Section 2 presents two related works to signal decomposition. Section 3 introduces the proposed sparse signal subspace decomposition based on the adaptive over-complete dictionary. Some experimental results and analysis are presented in Section 4. Finally, we draw the conclusion in Section 5.

<figure><img src="image_4.png" /></figure>

<footnote>*Correspondence: hongsun@whu.edu.cn</footnote>

<footnote n="1">¹School of Electronic Information, Wuhan University, Luojia Hill, Wuhan 430072, China</footnote>

<footnote n="2">²Signal and Image Processing Department, Telecom ParisTech, 46 rue Barrault, 75013 Paris, France</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
## 2 Review of PCA and sparse coding methods

We start with a brief description of two well-established approaches to signal decomposition that are relevant and related to the approach proposed in the next section.

### 2.1 PCA-based subspace decomposition

The basic tool of SSM is principal component analysis (PCA). PCA makes use of an orthonormal basis to capture on a small set of vectors (the signal subspace) as much energy as possible from the observed data. The other basis vectors are expected to contain noise only, and the signal projection on these vectors is rejected.

Consider a data set $\{x_m \in \mathbb{R}^{N \times 1}\}_{m=1}^M$ grouped in a matrix form $X$ of size $N \times M$: $X = \{x_m\}_{m=1}^M$. The PCA is based on singular value decomposition (SVD) with singular values $\sigma_i$ in descending order obtained from:

$$X = UA = U\Sigma V^T \quad (1)$$

where $U$ and $V$ are unitary matrices of sizes $N \times N$ and $M \times M$, respectively ($U^T U = I_N$, $V^T V = I_M$), and $\Sigma = \begin{bmatrix} \text{diag}[\sigma_1, \dots, \sigma_r], \mathbf{0} \end{bmatrix}$ of size $N \times M$ with $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r > 0$, $\{\sigma_i\}_{i=1}^r$ are positive real known as the singular values of $X$ with rank $r$ ($r \le N$).

Equation (1) can be rewritten in a vector form as:

$$\begin{aligned} & [\mathbf{x}_1 \quad \mathbf{x}_2 \quad \dots \quad \mathbf{x}_m \quad \dots \quad \mathbf{x}_M] \\ & = [\mathbf{u}_1 \quad \mathbf{u}_2 \quad \dots \quad \mathbf{u}_n \quad \dots \quad \mathbf{u}_N] \cdot [\boldsymbol{\alpha}_1 \quad \boldsymbol{\alpha}_2 \quad \dots \quad \boldsymbol{\alpha}_m \quad \dots \quad \boldsymbol{\alpha}_M] \end{aligned} \quad (2)$$

where $U = \{\mathbf{u}_n \in \mathbb{R}^{N \times 1}\}_{n=1}^N$ and $A = \{\boldsymbol{\alpha}_m \in \mathbb{R}^{N \times 1}\}_{m=1}^M$. Equation (2) means that the data set $\{x_m\}_{m=1}^M$ is expressed on the orthonormal basis $\{\mathbf{u}_n\}_{n=P+1}^P$ as $\{\boldsymbol{\alpha}_m\}_{m=1}^M$.

In the SVD decomposition given in Eq. (1), the standard deviation $\sigma_i$ is used as the measurement for identifying the meaningful basis vector $\mathbf{u}_i$. PCA takes the first $P(P < r)$ components $\{\mathbf{u}_n\}_{n=1}^P$ to span the signal subspace, and the remainders $\{\mathbf{u}_n\}_{n=P+1}^P$ are considered in a noise subspace orthogonal to the signal subspace. Therefore, projection on the signal subspace will hopefully filter out noise and reveal hidden structures. The reconstructed signal $\hat{\mathbf{S}}_{\text{PCA}}$ of size $N \times M$ is obtained by projecting in the signal subspace as:

$$\hat{\mathbf{S}}_{\text{PCA}} = [\mathbf{u}_1 \cdots \mathbf{u}_P \quad \mathbf{0}_{P+1} \cdots \mathbf{0}_N] \cdot [\boldsymbol{\alpha}_1 \quad \boldsymbol{\alpha}_2 \quad \dots \quad \boldsymbol{\alpha}_m \quad \dots \quad \boldsymbol{\alpha}_M] \quad (3)$$

The underlying assumption is that information in the data set is almost completely contained in a small linear subspace of the overall space of possible data vectors, whereas additive noise is typically distributed through the larger space isotropically. PCA, using the standard deviation as a criterion, implies that the components of the signal of interest in the data set have a maximum variance

and the other components are mainly due to the noise. However, in many practical cases, some components with low variances might actually be important because they carry information relative to the signal details. On the contrary, when dealing with noise with non-Gaussian statistics, it may happen that some noise components may actually have higher variances. At last, note that it is often difficult to provide a physical meaning to the orthonormal basis $\{\mathbf{u}_i\}_{i=1}^r$ of the SVD decomposition (Eq. (2)) although they have a very clear definition in the mathematical sense as orthogonal, independent, and normal. It is therefore difficult to impose known constraints on the signal features when they exist after the principal component decomposition.

### 2.2 Sparse decomposition

Recent years have shown a growing interest in research on the sparse decomposition of $M$ observations $\{x_m \in \mathbb{R}^N\}_{m=1}^M$ based on a dictionary $\mathbf{D} = \{\mathbf{d}_k\}_{k=1}^K \in \mathbb{R}^{N \times K}$. When $K > N$, the dictionary is said to be over-complete. $\mathbf{d}_k \in \mathbb{R}^N$ is a basis vector, also called an atom since it is not necessarily independent. By learning from data set $\{x_m\}_{m=1}^M$, the sparse decomposition is the solution of Eq. (4) [4]:

$$\begin{aligned} \{\mathbf{D}, \boldsymbol{\alpha}_m\} = & \operatorname*{argmin}_{\mathbf{D}, \boldsymbol{\alpha}_m} \| \boldsymbol{\alpha}_m \|_0 \\ & + \| \mathbf{D}\boldsymbol{\alpha}_m - \mathbf{x}_m \|_2^2 \le \varepsilon, \quad 1 \le m \le M \end{aligned} \quad (4)$$

where $\boldsymbol{\alpha}_m = [\alpha_m(1) \ \alpha_m(2) \ \dots \ \alpha_m(K)]^T \in \mathbb{R}^{K \times 1}$ is the sparse code of the observation $\mathbf{x}_m$. The allowed error tolerance $\varepsilon$ can be chosen according to the standard deviation of the noise. An estimate of the underlying signal $\{\mathbf{s}_m\}_{m=1}^M$ embedded in the observed data set $\{x_m\}_{m=1}^M$ would be:

$$\begin{aligned} & [\hat{\mathbf{s}}_1 \quad \hat{\mathbf{s}}_2 \quad \dots \quad \hat{\mathbf{s}}_m \quad \dots \quad \hat{\mathbf{s}}_M] \\ & = [\mathbf{d}_1 \quad \mathbf{d}_2 \quad \dots \quad \mathbf{d}_K] \cdot [\boldsymbol{\alpha}_1 \quad \boldsymbol{\alpha}_2 \quad \dots \quad \boldsymbol{\alpha}_m \quad \dots \quad \boldsymbol{\alpha}_M] \quad (5) \\ & \text{or equivalently } \hat{\mathbf{S}} = \mathbf{D}\mathbf{A} \end{aligned}$$

where the matrix $\mathbf{A}$ of size $K \times M$ is composed of $M$ sparse column vectors $\boldsymbol{\alpha}_m$.

The first term on the right side of Eq. (4) is a sparsity-inducing regularization that constrains the solution with the fewest number of nonzero coefficients in each of the sparse code vectors $\boldsymbol{\alpha}_m(1 \le m \le M)$. The underlying assumption is that a meaningful signal could be represented by combining few atoms. This learned dictionary adapted to sparse signal descriptions has proved to be more effective in signal reconstruction and classification tasks than the PCA method, which is demonstrated in the next section. The second term in Eq. (4) is the residual
of the reconstruction, based on the mean-square reconstruction error estimate in the same way as in the PCA method.

On the other hand, we note that the dictionary $\mathbf{D}$, a basis in sparse decomposition, is produced by learning noisy data set $\{\mathbf{x}_m\}_{m=1}^M$, so the basis vectors $\{\mathbf{d}_k\}_{k=1}^K$ should be decomposed into a principal subspace and a residual subspace. However, it is impossible to exploit an energy-constrained subspace since $\{\mathbf{d}_k\}_{k=1}^K$ are not necessarily orthogonal or independent.

# 3 The proposed sparse subspace decomposition

In this section, we introduce a novel criterion to the subspace decomposition over a learned dictionary and a corresponding index of significance of the atoms. Then we propose a signal sparse subspace decomposition (3SD) method under this new criterion.

## 3.1 Weight vectors of learned atoms

At first, we intend to find out the weight of the atoms. In the sparse representation given in (5), coefficient matrix $\mathbf{A}$ is composed by $M$ sparse column vectors $\alpha_m$, each $\alpha_m$ representing the weight of the observation $\mathbf{x}_m$, a local parameter for the $m$-th observation. Let us consider the row vectors $\{\beta_k\}_{k=1}^K$ of coefficient matrix $\mathbf{A}$:

$$ 
\begin{aligned}
\mathbf{A} &= [\alpha_1 \ \alpha_2 \ \cdots \ \alpha_M] \\
&= \begin{bmatrix}
\alpha_1(1) & \alpha_2(1) & \cdots & \alpha_M(1) \\
\alpha_1(2) & \alpha_2(2) & \cdots & \alpha_M(2) \\
\vdots & \vdots & \ddots & \vdots \\
\alpha_1(K) & \alpha_2(K) & \cdots & \alpha_M(K)
\end{bmatrix} = \begin{bmatrix}
\beta_1 \\
\beta_2 \\
\vdots \\
\beta_K
\end{bmatrix} \\
\text{where } \beta_k &= [\alpha_1(k) \ \alpha_2(k) \ \cdots \ \alpha_M(k)] \in \mathbb{R}^{1 \times M}
\end{aligned}
\quad (6)
$$

Note that the row vector $\beta_k$ is not necessarily sparse. Then Eq. (5) can be rewritten as:

$$ 
\begin{aligned}
\hat{\mathbf{S}} &= \mathbf{D}\mathbf{A} \\
&= [\mathbf{d}_1 \cdots \mathbf{d}_k \cdots \mathbf{d}_K] \cdot [\beta_1^T \cdots \beta_k^T \cdots \beta_K^T]^T
\end{aligned}
\quad (7)
$$

Equation (7) means that the row vector $\beta_k$ is the weight of the atom $\mathbf{d}_k$, which is a global parameter over the data set $\mathbf{X}$. Denoting $\|\beta_k\|_0$ the $l^0$ zero pseudo-norm of $\beta_k$. $\|\beta_k\|_0$ is the number of occurrences of atom $\mathbf{d}_k$ over the data set $\{\mathbf{x}_m\}_{m=1}^M$. We call it the frequency of the atom $\mathbf{d}_k$ denoted by $f_k$:

$$ f_k \triangleq \text{Frequency}(\mathbf{d}_k|\mathbf{X}) = \|\beta_k\|_0 \quad (8) $$

In the sparse decomposition, basis vectors $\{\mathbf{d}_k\}_{k=1}^K$ are prototypes of signal segments. That allows us to take them

as a signal patterns. Thereupon, some important features of this signal pattern could be considered as a criterion to identify significant atoms. It is demonstrated [5] that $f_k$ is a good description of the signal texture. Intuitively, a signal pattern must occur in meaningful signals with higher frequency even with a lower energy. On the contrary, a noise pattern would hardly be reproduced in observed data even with a higher energy.

It is reasonable to take this frequency $f_k$ as a relevance criterion to decompose the over-complete dictionary into a principal signal subspace and a remained noise subspace. Here, we use the word “subspace,” but in fact, these two subspaces are not necessarily independent.

## 3.2 Subspace decomposition based on over-complete dictionary

Taking vectors $\{\beta_k\}_{k=1}^K$, we calculate their $l^0$-norms $\{\|\beta_k\|_0\}_{k=1}^K$ and rank them in descending order as follows. The index $k$ of vectors $\{\beta_k\}_{k=1}^K$ are belonging to the set $\mathbf{C} = \{1, 2, \cdots, k, \cdots, K\}$. A one-to-one index mapping function $\pi$ is defined as:

$$ 
\begin{aligned}
\pi(\mathbf{C} \rightarrow \mathbf{C}) : k = \pi(\tilde{k}), \quad k, \tilde{k} \in \mathbf{C} \\
\text{s.t.} \quad \|\beta_{\pi(1)}\|_0 \ge \|\beta_{\pi(2)}\|_0 \ge \cdots \ge \|\beta_{\pi(\tilde{k})}\|_0 \ge \cdots \ge \|\beta_{\pi(K)}\|_0
\end{aligned}
\quad (9)
$$

By the permutation $\pi$ of the row index $k$ of matrix $\mathbf{A} = [\beta_1^T \cdots \beta_k^T \cdots \beta_K^T]^T$, the reordered coefficient matrix $\tilde{\mathbf{A}}$ becomes

$$ \tilde{\mathbf{A}} = [\beta_{\pi(1)}^T \ \beta_{\pi(2)}^T \ \cdots \ \beta_{\pi(k)}^T \ \cdots \ \beta_{\pi(K)}^T]^T \quad (10) $$

With corresponding reordered dictionary $\tilde{\mathbf{D}} = \{\mathbf{d}_{\pi(k)}\}_{k=1}^K$, Eq. (7) can be written as:

$$ 
\begin{aligned}
\hat{\mathbf{S}} &= \tilde{\mathbf{D}}\tilde{\mathbf{A}} \\
&= [\mathbf{d}_{\pi(1)} \cdots \mathbf{d}_{\pi(k)} \cdots \mathbf{d}_{\pi(K)}] \cdot [\beta_{\pi(1)}^T \cdots \beta_{\pi(k)}^T \cdots \beta_{\pi(K)}^T]^T
\end{aligned}
\quad (11)
$$

Then, the span of the first $P$ atoms can be taken as a principal subspace $\mathbf{D}_P^{(S)}$, and the remaining atoms span a noise subspace $\mathbf{D}_{K-P}^{(N)}$ as:

$$ 
\begin{aligned}
\mathbf{D}_P^{(S)} &= \operatorname{span}\{\mathbf{d}_{\pi(1)}, \mathbf{d}_{\pi(2)}, \cdots, \mathbf{d}_{\pi(P)}\} \\
\mathbf{D}_{K-P}^{(N)} &= \operatorname{span}\{\mathbf{d}_{\pi(P+1)}, \mathbf{d}_{\pi(P+2)}, \cdots, \mathbf{d}_{\pi(K)}\}
\end{aligned}
\quad (12)
$$

An estimate $\hat{\mathbf{S}}_P$ of the underlying signal $\mathbf{S}$ embedded in the observed data set $\mathbf{X}$ can be obtained on the principal subspace $\mathbf{D}_P^{(S)}$ simply by linear combination:
$$ 
\begin{aligned}
\hat{\mathbf{S}}_P &= \mathbf{D}_p^{(S)} \cdot \mathbf{A}_p^{(S)} \\
&= [\mathbf{d}_{\pi(1)} \cdots \mathbf{d}_{\pi(k)} \cdots \mathbf{d}_{\pi(P)}] \cdot [\beta_{\pi(1)}^T \cdots \beta_{\pi(k)}^T \cdots \beta_{\pi(P)}^T]^T
\end{aligned}
\quad (13) $$

### 3.3 Threshold of atom's frequency

Determining the number $P$ of atoms spanning the signal subspace $\mathbf{D}_p^{(S)}$ is always a hard topic especially for wide-band signals. Here, $P$ is the threshold of atom's frequency $f_k$ to distinguish a signal subspace and a noise subspace. One of the advantages of 3SD is that this threshold $P$ can be easily chosen without any prior parameter.

For a noiseless signal even with some weak details, such as the image example in Fig. 1a, the atoms' frequencies $f_{\pi(k)}^{\text{image}}$s shown in Fig. 1d (in black line) are almost always high except the zero value. For a signal with strong noise, such as the example in Fig. 1b, the atoms' frequencies $f_{\pi(k)}^{\text{noise}}$s shown in Fig. 1d (in red line) are almost always equal to 1 without zero and very few with a value 2 or 3. It is easy to set a threshold $P$ of $f_k$ (dotted line in the Fig. 1d) to separate the signal's atoms from the noise's atoms. By contrast, using the values of atom's energies $\|\beta_k\|_2$s for the two images shown in Fig. 1c, it is rather a puzzle to identify principal bases.

For a noisy signal, such as an image example in Fig. 2a, its adaptive over-complete dictionary (Fig. 2b) consists of atoms of principal signal patterns, strong noise patterns, and noisy signal patterns. Principal signal atoms should have higher frequencies, strong noise atoms lower frequencies and noisy signal atoms moderate frequencies. Intuitively, the red line (Fig. 2c) should be a suitable threshold $P$ of the frequencies $f_k$s. In practical implementation, the value of $P$ could be simply decided relying on the histogram of $f_k$. As shown in Fig. 2d, one can set

the value of $f_k$ associated with the maximum point of its histogram to $P$ as follows:

$$ P = \arg \max_k \text{Hist}(\|\beta_k\|_0) \quad (14) $$

In fact, the performances in signal analyses by 3SD method are not sensitive to the threshold $P$, owed to the dependence of the atoms. To illustrate this point, we take three images, Barbara, Lena, and Boat. Their histograms of $f_k$ are shown in Fig. 3a with the maximum points in dotted lines, 121, 97, and 92, respectively. Figure 3b reports the peak signal-to-noise ratio (PSNR) of the retrieved images $\hat{\mathbf{S}}_P$ on the signal principal subspace $\mathbf{D}_p^{(S)}$ with respect to $P$. We can see that the PSNRs of the results remain the same in a large range around the maximum points (in dotted lines). Consequently, taking the value of $f_k$ associated to the maximum point of its histogram as the threshold $P$ is a reasonable solution.

## 4 Results and discussion

### 4.1 Signal decomposition methods

Taking a part of the noisy Barbara image (Fig. 4a), we show an example of the sparse signal subspace decomposition (3SD) and the corresponding retrieved image (Fig. 4b). For comparison, the traditional sparse decomposition and the PCA-based subspace decomposition are shown in Fig. 4c, d.

We use the PSNR to assess the noise removal performance:

$$ 
\begin{aligned}
\text{PSNR} &= 20 \cdot \log_{10} [\text{MAX}\{\mathbf{S}(i,j)\}] - 10 \cdot \log_{10} [\text{MSE}] \\
\text{MSE} &= \frac{1}{IJ} \sum_{i=0}^{I-1} \sum_{j=0}^{J-1} [\mathbf{S}(i,j) - \hat{\mathbf{S}}(i,j)]^2
\end{aligned}
\quad (15) $$

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 2 The threshold P of the frequencies f_k's. a Noisy image. b Over-complete dictionary D. c Frequency of d_k. d Histogram of d_k's frequency</figcaption></figure>

and the structural similarity index metric (SSIM) between the denoised image and the pure one to evaluate the preserving detail performance:

$$SSIM(\hat{S}, \hat{S}) = \frac{(2u_S u_{\hat{S}} + c_1)(2\sigma_{S\hat{S}} + c_2)}{(u_{\hat{S}}^2 + u_{\hat{S}}^2 + c_1)(\sigma_{S}^2 + \sigma_{\hat{S}}^2 + c_2)} \quad (16)$$

where $u_x$ is the average of x, $\sigma_x^2$ is the variance of x, $\sigma_{xy}$ is the covariance of x and y, and $c_1$ and $c_2$ are small variables to stabilize the division with a weak denominator.

Let us look at the proposed sparse signal subspace decomposition on the top of Fig. 4b. The 128 atoms $d_k$ of the learned over-complete dictionary D are shown in descending order of their energies measured by $\|\beta_k\|_2$. The 32 principal signal atoms are chosen from the dictionary D under the frequency criterion. They are shown in descending order of their frequencies measured by $\|\beta_k\|_0$ composing a signal subspace $D_{32}^{(S)}$. We can see that some of the principal atoms are not among the first 32 atoms with the largest energy in the over-complete

<figure><img src="image_7.png" /><figcaption>Fig. 3 The insensitivity of the threshold P. a Histograms of f_k. b PSNR of Ŝ_P with respect to P</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Signal decompositions. **a** Image sample. **b** Sparse Subspace decomposition. **c** Sparse decomposition. **d** Subspace decomposition</figcaption></figure>

dictionary **D**. The retrieved images are shown at the bottom of Fig. 4b. The image **S** on **D** is apparently denoised. The image **ŷ** on the signal subspace $D_{32}^{(S)}$ improves

obviously by preserving fine details with a high SSIM = 0.86 and at suppressing strong noise with a high PSNR = 36.41. On the other hand, the residual image on noise
subspace $D_{96}^{(N)}$ contains some very noisy information. This is because the atoms of the over-complete dictionary are not independent.

For the same example, the classical sparse decomposition is shown in Fig. 4c, using the K-SVD algorithm [6] in which the allowed error tolerance $\epsilon$ (in Eq. (4)) is set to a larger value to filter out noise. The retrieved image **S** has a high PSNR = 29.62, but it has obviously lost the weak information with SSIM = 0.82. This is because signal distortion and residual noise cannot be minimized simultaneously at dictionary learning by Eq. (4).

In another comparison, the PCA-based subspace decomposition is shown in Fig. 4d. The 64 components are orthonormal and the 32 principal components are of the largest variance. The retrieved image by projecting on the signal subspace is rather noisy with PSNR = 29.62. This is because it cannot suppress strong noise and preserve weak details of information only using the variance criterion.

## 4.2 Application to image denoising

The application of 3SD to image denoising is presented here. A major difficulty of denoising is to separate the underlying signal from the noise. The proposed 3SD method could win this challenge. In the 3SD method, the important components are selected from the over-complete dictionary relying on their occurrence number over the noisy image set. Evidently, the occurrence numbers would be large for the signal, even for weak details, such as edges or textures. On the other hand, the occurrence numbers would be low for different kinds of white Gaussian or non-Gaussian noises, even strong at intensity.

The 3SD algorithm for image denoising is presented as follows:

**Input:** Noisy image **X**

**Output:** Denoised image $\hat{\mathbf{S}}$

- Sparse representation {$\mathbf{D}, \mathbf{A}$}: using K-SVD algorithm [6] by (4)

- Identify principal atoms from **D** based on **A**:

- Compute the frequencies of atoms $\{\|\beta_k\|_0\}_{k=1}^K$ according to (6) and (8)

- Get the permutation $\pi$ sorting the index $k$ of $\{\|\beta_k\|_0\}_{k=1}^K$ by (9)

- Compute the threshold $P$ by (14)

- Obtain the signal principal atoms $\{\mathbf{d}_{\pi(k)}\}_{k=1}^P$ by (12)

- Reconstruct image $\hat{\mathbf{S}}_P$ by (13)

In this application, we intend to preserve faint signal details under a situation of strong noise.

In the experiments, dictionaries used $\mathbf{D}_s$ of size $64 \times 256$ ($K = 256$ atoms), designed to handle image patches $\mathbf{x}_m$ of size $N = 64 = 8 \times 8$ pixels.

## 4.3 Image denoising

A noisy Lena image $\mathbf{X} = \mathbf{S}+\mathbf{V}$ with an additive zero-mean white Gaussian noise $\mathbf{V}$ is used. The standard deviation of noise is $\sigma = 35$. A comparison is made between the 3SD method and the K-SVD method [6] which is one of the best denoising methods reported in the recent literatures.

From the results shown in Fig. 5, the 3SD method outperforms the K-SVD method by about 1 dB in PSNR and by about 1% in SSIM (depending on how much details in the images and how faint the details). In terms of subjective visual quality, we can see that the corner of the mouth and the nasolabial fold with weak intensities are much better recovered by the 3SD method.

## 4.4 SAR image despeckling

In the second experiment, a simulated SAR image with speckle noise is used. Speckle is often modeled as multiplicative noise as $x(i,j) = s(i,j)v(i,j)$ where $x, s$, and $v$ correspond to the contaminated intensity, the original intensity, and the noise level, respectively.

Figure 6 shows the despeckling results of a simulated one-look SAR scenario with a fragment of the Barbara image. A comparison is made with 3SD method and a probabilistic patch-based (PPB) filter based on nonlocal means approach [7] which can cope with non-Gaussian noise. We can see that PPB can well remove speckle noise. However, it also removes the low-intensity details. The 3SD method shows advantages at preserving fine details and at suppressing strong noise.

## 4.5 Comparison with BM3D method

With a spatial complicated image scene, we make a comparison of the 3SD-based denoising method with the BM3D algorithm [8], one of the best methods especially for image denoising reported in many recent literatures.

The effectiveness of any signal analysis method depends on the different conditions in different applications. For the image denoising application, the signals involved should be homogeneous. Therefore, a procedure of grouping is generally adopted to select homogeneous pixels. In the BM3D method, a block-matching grouping is taken before filtering. We adopt the same grouping technique and then filter each homogeneous group by the proposed 3SD method.

Firstly, we take a 256 × 256 Barbara image (Fig. 7a) with a strong additive zero-mean white Gaussian noise where $\sigma = 70$ (Fig. 7b). The denoising result by the BM3D algorithm is shown in Fig. 7c. It displays a quite high performance. The denoising result by the 3SD-based method is shown in Fig. 7d. It demonstrates a higher PSNR and a higher SSIM and a better subjective visual quality over the BM3D algorithm.

Secondly, we take a simulated one-look SAR image with this 256 × 256 Barbara image (Fig. 8a) where PSNR =
<figure><img src="image_9.png" /><figcaption>**Fig. 5** Image denoising comparing the proposed 3SD method with the K-SVD method</figcaption></figure>

−0.1042. Figure 8b shows the despeckling result by the PPB method [7], in which grouping is realized based on nonlocal similarity and filtering is implemented by averaging each homogeneous group. The despeckled image is too smooth due to average filtering. Figure 8c shows the despeckling result by the BM3D-based method [9], in which grouping is realized based on similar 2-D fragments and filtering is implemented by Wiener shrinkage

coefficients from the energy of the 3-D transform coefficients. The despeckled image is much better in PSNR and in SSIM, but it seems a little noisy still due to the used energy criterion which is not effective enough to separate noise elements from the principal elements. Figure 8d shows the despeckling result by the proposed 3SD-based method, in which grouping is realized based on nonlocal similarity [7] and filtering is implemented by the

<figure><img src="image_10.png" /><figcaption>**Fig. 6** SAR image despeckling comparing the proposed 3SD method with the PPB method</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7** Denoising for spatial complicated image scene comparing BM3D method with 3SD-based method. **a** Original: Barbara (a fragment) 256 * 256.</figcaption></figure>

**b** Noisy: $\sigma = 70$; PSNR = 11.22; SSIM = 0.142. **c** BM3D denoising: PSNR = 24.08; SSIM = 0.7026. **d** 3SD denoising: PSNR = 24.21; SSIM = 0.765

proposed sparse subspace decomposition. The despeckled image demonstrates some advantages of the 3SD method at preserving fine details and at suppressing speckle noise, attributed to the principal subspace decomposition.

## 5 Conclusions

We proposed a method of sparse signal subspace decomposition (3SD). The central idea of the proposed 3SD is to identify principal atoms from an adaptive over-complete dictionary relying on the occurrence frequency of atoms over the data set (Eq. (8)). The atom frequency is measured by zero pseudo-norms of weight vectors of atoms (Eqs. (6) and (8)). The principal subspace is spanned by the maximum frequency atoms (Eq. (12)).

The 3SD method combines the variance criterion, the sparsity criterion, and the component's frequency criterion into a uniform framework. As a result, it can identify more effectively the principal atoms with the three important signal features. On the contrary, PCA uses only variance criterion and sparse coding method uses the variance and the sparsity criterions. In those ways, it is more difficult to distinguish weak information from strong noise.

Another interesting asset of the 3SD method is that it takes benefits from using an over-complete dictionary which reserves details of information and from subspace decomposition which rejects strong noise. On the contrary, some undercomplete dictionary methods [10] and

<figure><img src="image_12.png" /><figcaption>**Fig. 8** Despeckling for spatial complicated SAR image scene. **a** Noisy: 1-look; PSNR = -0.1042; SSIM = 0.21. **b** PPB despeckling: PSNR = 7.8943; SSIM = 0.63. **c** SAR-BM3D despeckling: PSNR = 9.9312; SSIM = 0.73. **d** 3SD-based despeckling: PSNR = 10.3316; SSIM = 0.74</figcaption></figure>
some sparse shrinkage methods [11, 12] might lose weak information when suppressing noise.

Moreover, the 3SD method is very simple with a linear retrieval operation (Eq. (13)). It does not require any prior knowledge on distribution or parameter to determine a threshold (Eq. (14)). On the contrary, some sparse shrinkage methods, such as [11], necessitate non-linear processing with some prior distributions of signals.

The proposed 3SD could be interpreted as a PCA in sparse decomposition, so it admits straightforward extension to applications of feature extraction, inverse problems, or machine learning.

**Acknowledgements**

The idea of the sparse signal subspace decomposition here arises through a lot of deep discussions with Professor Henri Maitre at Telecom ParisTech in France; he also gave suggestion on the structure of the manuscript.

**Funding**

This work was supported by the National Natural Science Foundation of China (Grant No. 60872131).

**Authors' contributions**

HS proposed the idea, designed the algorithms and experiments, and drafted the manuscript. CWS gave suggestion on the design of the algorithm, realized the algorithms, and carried out the experiments. DLR gave suggestions on the mathematical expressions of the manuscript and experiment analysis as well as explanation and helped draft the manuscript. All authors read and approved the final manuscript.

**Competing interests**

The authors declare that they have no competing interests.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Author details**

¹School of Electronic Information, Wuhan University, Luojia Hill, Wuhan 430072, China. ²Signal and Image Processing Department, Telecom ParisTech, 46 rue Barrault, 75013 Paris, France. ³CEDRIC Laboratory, CNAM, 292 rue Saint Martin, 75003 Paris, France.

Received: 13 September 2016 Accepted: 17 July 2017

Published online: 28 July 2017

**References**

1. K Hermus, P Wambacq, HV Hamme, A review of signal subspace speech enhancement and its application to noise robust speech recognition. EURASIP J. Adv. Signal Process. 888–896 (2007). doi:10.1155/2007/45821

2. DL Donoho, IM Johnstone, G Kerkyacharian, D Picard, Wavelet shrinkage: asymptopia? J. R. Stat. Soc. Ser. B. 57, 301–369 (1995). http://citeseer ist. psu.edu/viewdoc/summary?doi=10.1.1.162.1643

3. DW Tufts, R Kumaresan, I Kirsteins, Data adaptive signal estimation by singular value decomposition of a data matrix. IEEE Proc. 70(6), 684–685 (1982). doi:10.1109/PROC.1982.12367

4. M Elad, M Aharon, Image denoising via sparse and redundant representations over learned dictionaries. IEEE Trans. Image Process. 15(12), 3736–3745 (2006). doi:10.1109/TIP.2006.881969

5. G Tartavel, Y Gousseau, G Peyré, Variational texture synthesis with sparsity and spectrum constraints. J. Math. Imaging Vis. 52(1), 124–144 (2015). doi:10.1007/s10851-014-0547-7

6. M Aharon, M Elad, A Bruckstein, K-SVD: an algorithm for designing overcomplete dictionaries for sparse representation. IEEE Trans. Signal Process. 54(11), 4311–4322 (2006). doi:10.1109/TSP.2006.881199

7. CA Deledalle, L Denis, F Tupin, Iterative weighted maximum likelihood denoising with probabilistic patch-based weights. IEEE Trans. Image Process. 18(12), 2661–2672 (2009). doi:10.1109/TIP.2009.2029593

8. K Dabov, A Foi, V Katkovnik, K Egiazarian, Image denoising by sparse 3-D transform-domain collaborative filtering. IEEE Trans. Image Process. 16, 2080–2095 (2007). doi:10.1109/TIP.2007.901238

9. S Parrilli, M Poderico, CV Angelino, L Verdoliva, A nonlocal SAR image denoising algorithm based on LLMNSE wavelet shrinkage. IEEE Trans. Geosci. Remote Sens. 50, 606–616 (2012). doi:10.1109/TGRS.2011.2161586

10. F Porikli, R Sundaresan, K Suwa, SAR depeckling by sparse reconstruction on affinity nets. EUSAR2012. 18, 796–799 (2012). https://www.research gate.net/publication/232905473_SAR_Despeckling_by_Sparse_ Reconstruction_on_Affinity_Nets_SRAN

11. A Hyvarinen, P Hoyer, E Oja, Sparse code shrinkage for image denoising. IEEE World Congr. Comput. Intell. 2, 59–864 (1998). doi:10.1109/IJCNN. 1998.685880

12. R Malutan, R Terebes, C Germain, Speckle noise removal in ultrasound images using sparse code shrinkage. 5th IEEE Conf. E-Health Bioeng. Conf. 2, 1–4 (2015). doi:10.1109/EHB.2015.7391394

Submit your manuscript to a SpringerOpen journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com