Dictionary Learning with the K-SVDAlgorithm for Recovery of Highly Textured Images

1. Introduction

Motivated by studying information processed in the primary visual cortex, Field and Olshausen [1] developed a first approach to sparse dictionary learning. Technically, sparse dictionary learning is a representation learning method that aims at finding a sparse representation of input data. This is done in terms of a linear combination of basic elements called atoms. The aim of dictionary learning is to find the basic elements themselves, relying on given training data, as well as finding a useful linear combination of them. Typically, the elements of the dictionary do not make up an orthogonal basis. It is often of interest to have an overcomplete dictionary with significantly more elements than the dimensionality of the input data encourages. The representation of input data in terms of the overcomplete dictionary may be more sparse than if we had just enough elements to span the space.

In this paper, we deal with textured images as input data. For possible dictionaries in image processing, one may either employ a predefined dictionary using for example wavelets for construction, or learn it at hand of the given imagery. As indicated, we pursue the latter approach, which has lead to state-of-the-art results in many applications, for example denoising [2,3], image deblurring [4,5], or image segmentation [6]. We refer to [7,8] for more information and recent applications.

Let us turn to the learning process. We describe the dictionary in terms of a matrix $D\in {\mathbb{R}}^{n\times K}$, where $K\gg n$ so that the dictionary will be overcomplete. Thereby the atoms are the columns, and one may refer to D in terms of the atoms as ${\left\{d_i\right\}}_{i=1}^K$ with $d_i\in {\mathbb{R}}^n$. We assume that a training database ${\left\{y_i\right\}}_{i=1}^N$ with $y_i\in {\mathbb{R}}^n$ is given. Let us put together the $y_i$ into a training matrix $Y\in {\mathbb{R}}^{n\times N}$, and define a sparse representation matrix $X\in {\mathbb{R}}^{K\times N}$. The core idea is to find the dictionary D and the sparse representation matrix X such that $Y\approx DX$. This may be considered as a specific matrix decompositionof the training matrix. Finding the approximate decomposition in terms of a dictionary matrix and a sparse representation matrix is a minimization problem with two stages: dictionary updating and sparse coding.

Various methods deal with these two stages in different ways, see e.g. [7]. Among them, the K-SVDalgorithm [9] is one of the most popular methods for dictionary learning. It is a generalization of the K-means algorithm [10], combined with a step that performs the singular value decomposition (SVD) to update the atoms of the dictionary one by one. Besides, it enforces encoding of each input data by a linear combination of a limited number of non-zero elements. In the Section 2, we discuss the K-SVDmethod in some more detail.

Related work. Let us focus here on some important methods that are technically related to the K-SVDmethod and may be seen as extensions of the work [1]. The method of optimal directions (MOD) [11] is one of the first methods introduced to tackle the sparse dictionary learning problem. Its key idea is to solve the arising minimization problem in the sparse coding stage by enforcing a limited number of contributions of each atom in the training set. Concerning the latter point, let us note that the K-SVDmethod employs a similar approach. The algorithm in [12] is a process for learning an overcomplete basis by viewing it as a probabilistic model of the observed data. This method can be viewed as a generalization of the technique of independent component analysis, which is technically related to the SVD. The method in [13] proposes to learn a sparse overcomplete dictionary as unions of orthonormal bases. The dictionary update of one chosen atom is performed using SVD. Let us note that K-SVDfollows a similar approach.

There have been general experimental evaluations of the overall usefulness of the K-SVDmethod in dictionary learning [9], see also theoretical discussions in [14,15].

Our contribution. In this paper, we give a dedicated experimental analysis of the K-SVDalgorithm in its ability to recover highly textured images. Such images impose the challenges that their content is often not smooth and there are many fine structures that need to be recovered. Our investigation is motivated by ongoing research in corresponding applications.

2. Sparse Dictionary Learning and K-SVDAlgorithm

As indicated, in sparse dictionary learning, the atoms ${\left\{d_i\right\}}_{i=1}^K$ yield an overcomplete spanning set and provide an improvement in the sparsity and flexibility of the input data representation.

The sparse representation x of an input y may be exact or approximate in terms of the dictionary D. The general form of the model for exact representation in sparse coding amounts to solve

$$\underset{x}{min}{\parallel x\parallel }_0\mathrm{w}.\mathrm{r}.\mathrm{t}.y=Dx$$

(2.1)

while for the approximate representation we refer to

$$\underset{x}{min}{\parallel x\parallel }_0\mathrm{w}.\mathrm{r}.\mathrm{t}.{\parallel y-Dx\parallel }_2\le ϵ$$

(2.2)

Thereby ${\parallel ·\parallel }_0$ counts the non-zero entries of a vector, ${\parallel ·\parallel }_2$ stands for the Euclidean norm; the known parameter $ϵ$ is a tolerance for the sparse representation error.

Turning now concretely to the optimization of the learning process based on the training database ${\left\{y_i\right\}}_{i=1}^N$, the aim is to estimate $D\in {\mathbb{R}}^{n\times K}$ with atoms ${\left\{d_i\right\}}_{i=1}^K$ and $X\in {\mathbb{R}}^{K\times N}$ with columns ${\left\{x_i\right\}}_{i=1}^N$, such that

$$\underset{D,x_i}{min}\sum _{i=1}^N{\parallel x_i\parallel }_0\mathrm{w}.\mathrm{r}.\mathrm{t}.{\parallel y_i-D{x}_i\parallel }_2^2\le ϵ$$

(2.3)

Solving (3) for each $y_i\in {\mathbb{R}}^n$ gives a sparse representation $x_i^{opt}\in {\mathbb{R}}^K$ over the unknown dictionary D, aiming to find the proper sparse representations and the dictionary jointly. Thus, (3) encompasses both sparse coding and dictionary updating.

One may also consider the alternative form of (3) as

$$\underset{D,x_i}{min}\sum _{i=1}^N{\parallel y_i-D{x}_i\parallel }_2^2\mathrm{w}.\mathrm{r}.\mathrm{t}.{\parallel x_i\parallel }_0\le T$$

(2.4)

where the role of the penalty and the constraint are reversed, and the parameter T controls the number of non-zero elements in the vector.

Let us now turn to the main algorithm of the paper. As the K-SVDmethod is a generalization of the K-means method, we begin by a brief explanation of the latter before discussion of K-SVD.

2.1. The K-means Algorithm

There is a relation between sparse representation and clustering [11,16]. The K-means method is a possible approach to clustering which aims to partition N inputs into K clusters, and each input belongs to the cluster with the nearest mean. The cluster centres serve as a prototype of the cluster.

The K-means method amounts to an iterative process incorporating two steps per iteration. For convenience, let us also relate to corresponding steps in dictionary learning. First, given ${\left\{d_i\right\}}_{i=1}^K$, it assigns the training samples (${\left\{y_i\right\}}_{i=1}^N$) to their nearest neighbour in terms of the $d_i$ (sparse coding). This means, the $d_i$ have the role of the closest cluster centroid. Let us note that in standard K-means this is done regarding the squared Euclidean distance, while in K-SVD[9] the non-squared Euclidean distance is employed mimicking the sparse coding problem setup. Then, the dictionary ${\left\{d_i\right\}}_{i=1}^K$ is updated to better fit the samples by use of the assignment given in the first step (dictionary update).

More precisely, in the sparse coding stage the method makes K partitions (subsets) of the training samples ${\left\{y_i\right\}}_{i=1}^N$, related to K columns of the dictionary, each holding training samples most similar to the corresponding column. At the dictionary updating stage, each dictionary column is replaced by the centroid of the corresponding subset. Let us note that K is a parameter that needs to be specified.

2.2. The K-SVDAlgorithm

The essence of the iterative K-SVDalgorithm, as a generalization of K-means , may be described as follows. First the dictionary is fixed to find the best possible sparse matrix representation X under the constraint in (4). Then the atoms of the dictionary D are updated iteratively to better fit to the training data. This is done specifically by finding a rank-1 approximation of a residual matrix via making use of the SVD while preserving the sparsity.

One may view the problem posed in (4) as a nested minimization problem, namely an inner minimization of the number of non-zeros in the representation vectors $x_i$ for a given fixed D (sparse coding stage) and an outer minimization over D (dictionary update stage). At the jth step, obtaining a sparse representation via a pursuit algorithm [7] such as e.g. orthogonal matching pursuit, we use the dictionary $D_{j-1}$ from the $(j-1)$th step and solve N instances of (4) for each $y_i$. That gives us the matrix $X_j$ containing columnwise all sparse representations (sparse representation matrix). Summarizing thus the squared Euclidean norms of vectors in (4) using the Frobenius norm, we solve for $D_j$ in terms of the following least squares problem, the solution of which has an exact representation:

$$D_j=\underset{D}{arg\; min}{\parallel Y-D{X}_j\parallel }_F^2=Y{X}_j^T{\left(X_j{X}_j^T\right)}^{-1}$$

(2.5)

The K-SVDmethod handles the atoms in D sequentially. Keeping all the columns fixed except the $i_0$th one, $\left(d_{i_0}\right)$, it updates $d_{i_0}$ along with the coefficients that multiply it in the sparse representation matrix. To this end, we may rewrite the penalty term in (5), omitting the iteration number j, as

$${\parallel Y-DX\parallel }_F^2={∥Y-\sum _{i=1}^K{d}_i{x}_i^r∥}_F^2={∥E_{i_0}-d_{i_0}{x}_{i_0}^r∥}_F^2$$

(2.6)

to isolate the dependency on $d_{i_0}$. Here, $E_{i_0}=Y-{\sum }_{i\ne i_0}{d}_i{x}_i^r$ is a precomputed residual matrix and $x_i^r$ stands for the row r of X, so that $x_i^r$ indicates the contribution of $d_i$ in the training set. Which means it includes sparse coefficients of elements of the training set that currently use $d_i$.

The optimal values of $d_{i_0}$ and $x_{i_0}^r$ minimizing (6) are given by the rank one approximation of $E_{i_0}$, which one can obtain using the SVD. That typically would yield a dense vector $x_{i_0}^r$, increasing the number of non-zeros in X. Keeping the cardinalities of all the representations fixed during the minimization, one may restrict $E_{i_0}$ to those columns where the entries in the row $x_{i_0}^r$ are non-zero. In this manner, only the existing non-zero coefficients in $x_{i_0}^r$ may vary, preserving in this way the sparsity.

Considering the SVD of the restricted $E_{i_0}$ as $E_{i_0}^R=U\Sigma V$, the updated dictionary column is the first column of U, and the updated non-zero coefficient of $x_{i_0}^r$ is the first column of V multiplied by $\Sigma (1,1)$. The K-SVDalgorithm performs an SVD for each of the K different residuals related to the K columns of the dictionary. This method is a direct generalization of the K-means process. When the model orders $T=1$ in (4), K-SVDexactly acts like K-means [9].

3. Data Set of Textured Images and Preprocessing

In this section, we want to explain which data set we consider, which transformation process is applied to the data set, and how we learn the dictionary.

Our experiments are based on a subset of the Brodatz texture database [18] obtained via the website [17], see Figure 1. The complete data set contains 112 highly textured gray scale images. For our studywe select images with a dot-like texture. The images in the Brodatz texture database have a resolution of $640\times 640$ pixels.

3.1. Preprocessing

We begin by describing the basic steps to transform the resource images $I_{\mathrm{orig}}$ into a useful format.In doing this, we mostly follow the basic steps of the Matlab code [19]. However, we also incorporate a slicing technique, as explained below.

Our starting data set contains 12 images $I_{\mathrm{orig}}\in {\mathbb{R}}^{640\times 640}$, which we want to extend by slicing the images. In doing this, we aim for two goals: First, we increase the number of elements in the data set. Second, working with smaller images will reduce computation time. By this motivation, we divide the images into $4\times 4$ blocks equal in size. Each block, with a resolution of $160\times 160$ pixels, will be saved as a separate image, yielding in total a new larger data set. This enlarged data set now contains 192 images and is used to learn the dictionary D.

We now pursue to explain the preprocessing step $\mathcal{P}:I\to W$. In the first step of the preprocessing, we are randomly choosing 20 out of these 192 images to assemble the foundation for training our dictionary D. At this point, we actually need the enlarged data set. The images I will be completely divided into small image patches $B_i\in {\mathbb{R}}^{b\times b}$, $i\in \{1,2,\dots ,n_b\}$ of a resolution $b\times b$. As our focus will now be on numerical computing rather than on interpretation as image patches, we will call $B_i$ batch, b the batch size, and $n_b$ is the number of batches in the image I. This dividing is not overlapping, and a single batch is visualized in Figure 2.

Next, in the third step, we will concatenate these batches $B_i$ to create the vector $w_i\in {\mathbb{R}}^{b^2}$.

$$B_i=\left[v_1,\dots ,v_b\right]\in {\mathbb{R}}^{b\times b}\to {\left[v_1^T,\dots ,v_b^T\right]}^T=w_i\in {\mathbb{R}}^{b^2}$$

(3.1)

Concatenating all $n_b$ batches $B_i$ of I and pinning them together will construct the matrix $W=[w_1,\dots w_{n_b}]\in {\mathbb{R}}^{b^2\times n_b}$.

Finally, we have transformed the images $I\in {\mathbb{R}}^{160\times 160}$ into the matrix $W\in {\mathbb{R}}^{b^2\times n_b}$. We can reverse this process ${\mathcal{P}}^{-1}:W\to I$ to construct an image I from a matrix W. This is used to visualize the recovered images $w_i=D{x}^{\mathrm{opt}}$, $\forall i\in [1,n_b]$, where $x^{\mathrm{opt}}$ is the solution from Equation 9.

3.2. Revisiting Sparse Coding

In accordance with Section 2, the K-SVDmethod is based on the sparse coding problem that may be formulated in general in terms of Equations (1) and (2). After preprocessing, let us now reformulate them as one problem mainly used for learning as

$$\forall i\in [1,M]\underset{D,x_i}{min}{\parallel w_i-D{x}_i\parallel }_2^2\mathrm{w}.\mathrm{r}.\mathrm{t}.{\parallel x_i\parallel }_0\le T_{\mathrm{train}}$$

(3.2)

and one for testing a computed reconstruction:

$$\forall i\in [1,n_b]x^{\mathrm{opt}}=\underset{x_i}{min}{\parallel w_i-D{x}_i\parallel }_2^2\mathrm{w}.\mathrm{r}.\mathrm{t}.{\parallel x_i\parallel }_0\le T_{\mathrm{test}}$$

(3.3)

3.3. Solving the Minimization Problem

Solving Equation 8 will follow as shown in [9] (Figure 2), which we will briefly recall now. In the sparse coding stage, the orthogonal matching pursuit algorithm is employed to solve Equation 8 only for $x_i$. In the dictionary update stage, we can use the solution of the previous stage to enhance the dictionary. This iterative process of using D to find a better x and then use x to find a better D is done, until the process converges, or it is stopped.

We want to give a short overview on the update process of the dictionary D now. However, first we need to discuss the initialization of D. We use $k\in \mathbb{N}$ random batches $B_i$ from our 20 images used for training. But, before the concatenating step, we normalize the matrix $\parallel B_i\parallel =1$. Therefore, the dimension of the dictionary will be $D\in {\mathbb{R}}^{b^2\times k}$. The main part of the updating stage is, to control if the atoms in the dictionary are used or not. If we do not use an atom, we simply remove it from the dictionary and put in another normalized concatenated batch.

4. Experimental Results and Discussion

Let us start with the dictionaries, and continue with the recovery of selected images. After discussing these, we proceed with a parameter study and the impact of the texture to the recovery process.

4.1. Learning the Dictionaries

We learn the dictionaries $D_b\in {\mathbb{R}}^{b^2\times k}$ for different batch sizes $b=\{4,8,16\}$. The batch size of $b=4$ can be interpreted as the minimum batch size, since batches of size $2\times 2$ are too small and end up in learning almost every pixel by itself. Doubling the batch size to $b=8$, and again to get $b=16$, appears to be a natural choice. With $D_b$, we will denote the learned dictionaries corresponding to batch size. The number of atoms in the dictionary will be set to $k=128$. Additionally, we fix the number $T_{\mathrm{train}}=5$. The corresponding learning results can be seen in Figure 3. These are the dictionaries we will use to recover images. Inspecting these dictionaries, we notice, that with an increasing batch size b, the atoms more and more look like each other.

4.2. Generalization Experiment

After learning the dictionary $D_b$, we can solve Equation 9 to recover an image, via

$$I_{\mathrm{rec}}={\mathcal{P}}^{-1}\left(W_{\mathrm{rec}}\right)\mathrm{with}{W}_{\mathrm{rec}}=[w_1^{\mathrm{rec}},\dots ,w_{n_b}^{\mathrm{rec}}]\mathrm{and}{w}_{n_b}^{\mathrm{rec}}=D_b{x}^{\mathrm{opt}}$$

(4.1)

For this experiment, we have chosen a $T_{\mathrm{rec}}=10$. The results can be seen in Figure 4. The 16 images on the left (cf. Figure 4) present the results of two images (D30 and D73) from the subset of the Brodatz data set also used for training. On the right (cf. Figure 4) we illustrate the results using two images (D42 and D44) that are not part of our training set. The left and the right figures in Figure 4 are organized the same. In the first and third row, we show the image in the original resolution of $640\times 640$. In the rows two and four, we are presenting a zoomed part of the image above. The first column is always presenting the original image. The columns two, three, and four reveal the recoveries using $D_4$, $D_8$ and $D_{16}$, in this order.

Visual inspection of Figure 4 reveals that using $D_4$ (2nd column) will recover the images almost perfectly. The recovery using $D_8$ (3rd column) works pretty well, if only focusing on the full resolution image. Examining zoomed parts, we notice that block-like artefacts emerge. In recovery with $D_{16}$, these artefacts are even more dominant, even in the full resolution image.

We will use the mean squared error (MSE) $\mathrm{MSE}={\parallel I-I_{\mathrm{rec}}\parallel }_2/N$ with N the number of pixels in the image, to quantify the quality of the recovered image. The MSE values for the images in Figure 4 are listed in Table 1. There we see, that if we switch from a batch size $b=4$ to $b=8$, we will increase the error by almost a power of ten. Then again, switching from $b=8$ to $b=16$ only roughly doubles the error. Generalizing from images that are part of our training set, to images that are not, seems to have no significant impact on recovery quality.

The amount of details/texture in an image appears to be impactful on MSE and the quality. The bottom image in Figure 4 (D73) and the top image in Figure 4 (D42) are highly textured images and have higher MSE values than the other two, less detailed, images in Figure 4.

4.3. Binary Image Experiment

We ended the last section with the hypothesis, that highly textured images produce larger MSE values. For investigation of this, we notice that the Brodatz images contain some visually not apparent noise that could influence the observed MSE values, as the method will attempt to recover it.

To study this issue, we created a binary (black and white) version of our training dataset. Then we learned a dictionary in the same way as before. No parameters were changed or adjusted. From the learning process we obtained the three dictionaries $D_4^{\mathrm{bw}}$, $D_8^{\mathrm{bw}}$ and $D_{16}^{\mathrm{bw}}$.

An example of original image and recoveries can be seen in Figure 5, keeping the positioning as in Figure 4. We notice that the recovery using $D_4^{\mathrm{bw}}$ leads to almost perfect recoveries. However, on the top right, there are some problems with the recovery. In the original image, we see there a hook-like structure in the zoomed image. This feature is missing in the recovery. Considering the recovery using $D_8^{\mathrm{bw}}$, we notice that the full-size image still looks very similar to the original image, on one hand. On the other hand, in the zoomed image, a blurring occurs that gets even more dominant in recovery with $D_{16}^{\mathrm{bw}}$. The recovery from $D_{16}^{\mathrm{bw}}$ can no longer be visually declared as a black-and-white image. This gray blurring comes from the normalization process when learning the dictionary. Therefore, using only binary images will not lead to a binary dictionary.

Computing the MSE for the recovered images, we get $\mathrm{MSE}=1.6891·{10}^{-5}$ for $b=4$, $\mathrm{MSE}=1.6699·{10}^{-4}$ for $b=8$, and $\mathrm{MSE}=4.8149·{10}^{-4}$ for $b=16$. The MSE belonging to $D_{16}^{\mathrm{bw}}$ is still smaller than the MSE of the gray valued recovery with $D_4$.

In this figure, we plot the difference between the original image in black-and-white to the recovery using $D_8^{\mathrm{bw}}$. On the right plot, we present a zoomed in image. The zoomed part was indicated left plot with a black square. Red dots indicate negative and blue dots a positive difference between the original and the recovery.

4.4. Studying the Parameter $T_{\mathrm{rec}}$

In the previous experiments, we always fixed $T_{\mathrm{rec}}=10$. This parameter controls the number of used atoms in the recovery. To study its influence, we compute the MSE of the recovery obtained with different values of $T_{\mathrm{rec}}\in [3,20]$ in Equation 9. The results are visualized in Figure 7.

In this figure, we gave each batch a plot. On the left we presented the MSE from a recovery using $D_4$ with a batch size of $b=4$. In the middle, we plot the results for $b=8$ and on the right plot the MSE for $b=16$ is presented. The blue line indicates the D30 image (cf. Figure 4 top), the red line belongs to the image D73 (cf. Figure 4 bottom), and the black denotes the black-and-white version of D30 (cf. Figure 5).

As a main observation, the batch size seems to have a bigger impact on increasing the recovery quality, than increasing the number of used atoms $T_{\mathrm{rec}}$. However, reducing the batch size will lead to an increase in computation time. In the end, we have to balance low MSE versus computation time.

In the left plot ($b=4$) we see a considerable drop in the MSE for $T_{\mathrm{rec}}\ge 16$. The same drop could be noticed if we use $T_{\mathrm{rec}}\ge 64$ or $T_{\mathrm{rec}}\ge 256$ for $8\times 8$ and $16\times 16$ batches, respectively. Thus, using values $T_{\mathrm{rec}}\ge b^2$ appears to be beneficial in terms of the MSE, but they are not a reasonable option for the K-SVDprocedure itself.

5. Conclusion and Future Work

Concerning recovery of highly textured images, we found that the biggest part in reconstruction error is by the recovery of edges. Additionally, we have shown that increasing the number of used atoms has a limited effect on increasing the quality of the recovered image.