Scalable Online Convolutional Sparse Coding

TL;DR

This paper introduces an online convolutional sparse coding algorithm that significantly improves scalability and efficiency for large datasets by reformulating the problem in the frequency domain and using ADMM, enabling faster convergence and handling more data.

Contribution

It presents a novel online learning approach for CSC that reduces computational costs and memory usage, allowing application to much larger datasets than previous batch methods.

Findings

Faster convergence compared to existing methods.

Better reconstruction performance.

Can process at least ten times more images than prior algorithms.

Abstract

Convolutional sparse coding (CSC) improves sparse coding by learning a shift-invariant dictionary from the data. However, existing CSC algorithms operate in the batch mode and are expensive, in terms of both space and time, on large datasets. In this paper, we alleviate these problems by using online learning. The key is a reformulation of the CSC objective so that convolution can be handled easily in the frequency domain and much smaller history matrices are needed. We use the alternating direction method of multipliers (ADMM) to solve the resulting optimization problem and the ADMM subproblems have efficient closed-form solutions. Theoretical analysis shows that the learned dictionary converges to a stationary point of the optimization problem. Extensive experiments show that convergence of the proposed method is much faster and its reconstruction performance is also better. Moreover, while existing CSC algorithms can only run on a small number of images, the proposed method can handle at least ten times more images.

Tables

Table 1. TABLE I: Comparing the proposed online CSC algorithm with existing batch CSC algorithms.

	batch/online	convolution operation	space	time for one data pass
DeconvNet [14]	batch	spatial	$O\left(NKP\right)$	$O\left(N{K}^2{P}^{2}M\right)$
FCSC [5]	batch	frequency	$O\left(NKP\right)$	$O\left(N{K}^{3}P+NKP\log P\right)$
FFCSC [8]	batch	frequency	$O\left(NKP\right)$	$O\left(N{K}^{2}P+NKP\log P\right)$
CBPDN [16]	batch	frequency	$O\left(NKP\right)$	$O\left(N^{2}KP+NKP\log P\right)$
CONSENSUS [17]	batch	frequency	$O\left(NKP\right)$	$O\left(NK{P}^2+NKP\log P\right)$
OCSC	online	frequency	$O\left(K^{2}P\right)$	$O\left(N{K}^{2}P+NKP\log P\right)$

Table 2. TABLE II: Summary of data sets used.

	dimension ($P$)	#training images	#testing images
Fruit	100$\times$100	10	4
City	100$\times$100	10	4
House	100$\times$100	100	4
Flower	500$\times$500	2040	6149
Dog	224$\times$224	12000	8580
CIFAR-10	32$\times$32	50000	10000

Table 3. TABLE III: Testing PSNR obtained on the small data sets. The highest and comparable PSNR (according to the pairwise t-test with 95% confidence) are in bold.

	Fruit	City	House
DeconvNet	27.41$\pm$0.13	28.11$\pm$0.25	25.39$\pm$0.27
FCSC	27.90$\pm$0.18	28.20$\pm$0.31	25.68$\pm$0.50
FFCSC	28.13$\pm$0.15	28.58$\pm$0.16	28.48$\pm$0.04
CBPDN	28.01$\pm$0.04	28.67$\pm$0.37	29.40$\pm$0.43
CONSENSUS	27.62$\pm$0.14	28.19$\pm$0.20	25.41$\pm$0.37
OCSC	28.61$\pm$0.06	28.86$\pm$0.13	29.68$\pm$0.06

Equations

$$\underset{D\in \mathcal{D},\{z_i\}}{\min }\frac{1}{N}\sum _{i=1}^N\frac{1}{2}\Vert x_i-D{z}_i{\Vert }_2^2+\beta \Vert z_i{\Vert }_1,$$

$$z_t=\arg \underset{z}{\min }\frac{1}{2}\Vert x_t-D_{t-1}z{\Vert }_2^2+\beta \Vert z{\Vert }_1,$$

$$D_t$$

$$=$$

$$A_t^{\text{(osc)}}$$

$$B_t^{\text{(osc)}}$$

$$A_t^{\text{(osc)}}$$

$$B_t^{\text{(osc)}}$$

$$\underset{x,y}{\min }f(x)+g(y):Ax+By=c,$$

$$f(x)+g(y)+{\nu }^{\top }(Ax+By-c)+\frac{\rho }{2}\Vert Ax+By-c{\Vert }^2,$$

$$x_{\tau }$$

$$y_{\tau }$$

$$u_{\tau }$$

$$\underset{D\in \mathcal{D},\{Z_i\}}{\min }\frac{1}{N}\sum _{i=1}^N\frac{1}{2}\Vert x_i-\sum _{k=1}^{K}D(:,k)\ast Z_i(:,k){\Vert }_2^2+\beta \Vert Z_i{\Vert }_1,$$

$$\underset{\tilde{D},\{{\tilde{Z}}_i\}}{\min }$$

$$\underset{{\tilde{Z}}_i,U_i}{\min }$$

$$U_i(:,k)={\mathcal{F}}^{-1}({\tilde{Z}}_i(:,k)),k=1,\dots ,K,$$

$$\underset{\tilde{D},V}{\min }$$

$$\{\begin{array}{ll}\mathcal{F}(V(:,k))=\tilde{D}(:,k),& k=1,\dots ,K,\\ \Vert \mathcal{H}((V(:,k))){\Vert }_2^2\le 1,& k=1,\dots ,K,\end{array}$$

$$\underset{D\in \mathcal{D},\{Z_i\}}{\min }\frac{1}{t}\sum _{i=1}^t\frac{1}{2}\Vert x_i-\sum _{k=1}^{K}D(:,k)\ast Z_i(:,k){\Vert }_2^2+\beta \Vert Z_i{\Vert }_1.$$

$$\underset{\tilde{D},V}{\min }$$

$$\{\begin{array}{ll}\mathcal{F}(V(:,k))=\tilde{D}(:,k),& k=1,\dots ,K,\\ \Vert \mathcal{H}(V(:,k)){\Vert }_2^2\le 1,& k=1,\dots ,K.\end{array}$$

$$\frac{1}{2tP}\sum _{i=1}^t{\Vert {\tilde{x}}_i-\dot{D}{\dot{z}}_i\Vert }_2^2,$$

$$\dot{D}=[\text{Diag}(\tilde{D}(:,1)),\dots ,\text{Diag}(\tilde{D}(:,K))],$$

$$A_t^{\text{(naive)}}$$

$$B_t^{\text{(naive)}}$$

$$\frac{1}{2P}\sum _{p=1}^P\tilde{D}(p,:)\left(\frac{1}{t}\sum _{i=1}^t{\tilde{Z}}_i^{†}(:,p){\tilde{Z}}_i(p,:)\right){\tilde{D}}^{†}(:,p)-\frac{1}{P}\sum _{p=1}^P\tilde{D}(p,:)\left(\frac{1}{t}\sum _{i=1}^t{\tilde{x}}_i(p){\tilde{Z}}_i^{†}(:,p)\right),$$

$$A_t^p$$

$$b_t^p$$

$$A_t^p$$

$$b_t^p$$

$$({\tilde{D}}_t,V_t)$$

$$\text{s.t.}\mathcal{F}(V(:,k))=\tilde{D}(:,k),k=1,\dots ,K,\Vert \mathcal{H}(V(:,k)){\Vert }_2^2\le 1,k=1,\dots ,K.$$

$${\tilde{D}}_{t,\tau }=\arg \underset{\tilde{D}}{\min }\frac{1}{2P}\sum _{p=1}^P\tilde{D}(p,:)A_t^p{\tilde{D}}^{†}(:,p)-2\tilde{D}(p,:)b_t^p+\frac{\rho }{2}\sum _{k=1}^K\Vert \tilde{D}(:,k)-{\tilde{V}}_{t,\tau -1}(:,k)+{\Theta }_{t,\tau -1}(:,k){\Vert }_2^2,$$

$${\tilde{D}}_{t,\tau }(p,:)$$

$${C}^{p}_{t}=\left\{\begin{array}[]{ll}\frac{1}{\rho P}I-\frac{1}{\rho P}\frac{\tilde{Z}_{t}^{\dagger}(:,p)\tilde{Z}_{t}(p,:)}{\rho P+\tilde{Z}_{t}(p,:)\tilde{Z}_{t}^{\dagger}(:,p)}&t=1\\ \frac{t}{t-1}\left[{C}^{p}_{t-1}-\frac{{C}^{p}_{t-1}\tilde{Z}_{t}^{\dagger}(:,p)\tilde{Z}_{t}(p,:){C}^{p}_{t-1}}{(t-1)+\tilde{Z}_{t}(p,:){C}^{p}_{t-1}\tilde{Z}_{t}^{\dagger}(:,p)}\right]&t>1\end{array}\right..$$

Full text

I Introduction

In recent years, sparse coding has been widely used in signal processing [1, 2] and computer vision [3, 4]. In sparse coding, each data sample is represented as a weighted combination of a few atoms from an over-complete dictionary learned from the data. Despite its popularity, sparse coding cannot capture shifted local patterns that are common in image samples. Often, it has to first extract overlapping image patches, which is analogous to manually convolving the dictionary with the samples. As each sample element (e.g., an image pixel) is contained in multiple overlapping patches, the separately learned representations may not be consistent. Moreover, the resultant representation is highly redundant [5].

Convolutional sparse coding (CSC) addresses this problem by learning a shift-invariant dictionary composed of many filters. Local patterns at translated positions of the samples are easily extracted by convolution, and eliminates the need for generating overlapping patches. Each sample is approximated by the sum of a set of filters convolved with the corresponding codes. The learned representations are consistent as they are obtained together. CSC has been used successfully in various image processing applications such as super-resolution image reconstruction [6], high dynamic range imaging [7], image denoising and inpainting [8]. It is also popular in biomedical applications, e.g., cell identification [9], calcium image analysis [10], tissue histology classification [11] and segmentation of curvilinear structures [12]. CSC has also been used in audio processing applications such as piano music transcription [13].

A number of approaches have been proposed to solve the optimization problem in CSC. In the pioneering deconvolutional network (DeconvNet) [14], simple gradient descent is used. As convolution is slow in the spatial domain, fast convolutional sparse coding (FCSC) [5] formulates CSC in the frequency domain, and the alternating direction method of multipliers (ADMM) [15] is used to solve the resultant optimization problem. Its most expensive operation is the inversion of a convolution-related linear operator. To alleviate this problem, convolutional basis pursuit denoising (CBPDN) [16] exploits a special structure of the dictionary, while the global consensus ADMM (CONSENSUS) [17] utilizes the matrix inverse lemma to simplify computations. Fast and flexible convolutional sparse coding (FFCSC) [8] further introduces mask matrices so as to handle incomplete samples that are common in image/video inpainting and demosaicking applications. Note that all these algorithms operate in the batch mode (i.e., all the samples/codes have to be accessed in each iteration). Hence they can become expensive, in terms of both space and time, on large data sets.

In general, online learning has been commonly used to improve the scalability of machine learning algorithms [18, 19]. While batch learning algorithms train the model after arrival of the whole data set, online learning algorithms observe the samples sequentially and update the model incrementally. Moreover, data samples need not be stored after being processed. This can significantly reduce the algorithm’s time and space complexities. In the context of sparse coding, an efficient online algorithm is proposed in [20]. In each iteration, information necessary for dictionary update is summarized in fixed-sized history matrices. The space complexity of the algorithm is thus independent of sample size. Recently, this has also been extended for large-scale matrix factorization [21].

However, though CSC is similar to sparse coding, the online sparse coding algorithm in [20] cannot be directly used. This is because convolution in CSC needs to be performed in the frequency domain for efficiency. Moreover, the sizes of history matrices depend on dimensionality of the sparse codes, which becomes much larger in CSC than in sparse coding. Storing the resultant history matrices can be computationally infeasible.

In this paper, we propose a scalable online CSC algorithm for large data sets. The algorithm, which will be called Online Convolutional Sparse Coding (OCSC), is inspired by the online sparse coding algorithm of [20]. It avoids the above-mentioned problems by reformulating the CSC objective so that convolution can be handled easily in the frequency domain and much smaller history matrices are needed. We use ADMM to solve the resultant optimization problem. It will be shown that the ADMM subproblems have efficient closed-form solutions. Consequently, to process a given number of samples, OCSC has the same time complexity as state-of-the-art batch CSC methods but requires much less space. Empirically, as OCSC updates the dictionary after coding each sample, it converges much faster than batch CSC methods. Theoretical analysis shows that the learned dictionary converges to a stationary point of the optimization problem. Extensive experiments show that convergence of the proposed method is much faster and its reconstruction performance is also better. Moreover, while existing CSC algorithms can only run on a small number of images, the proposed method can at least handle ten times more images.

The rest of the paper is organized as follows. Section II briefly reviews online sparse coding, the ADMM, and batch CSC methods. Section III describes the proposed online convolutional sparse coding algorithm. Experimental results are presented in Section IV, and the last section gives some concluding remarks.

Notations: For vector $a\in\mathbb{R}^{m}$, its $i$th element is denoted $a(i)$, its $\ell_{2}$ norm is $\|a\|_{2}=\sqrt{\sum_{i=1}^{m}a^{2}(i)}$, its $\ell_{1}$ norm is $\|a\|_{1}=\sum_{i=1}^{m}|a(i)|$, and $\text{Diag}(a)$ reshapes $a$ to a diagonal matrix with elements $a[i]$’s. Given another vector $b\in\mathbb{R}^{n}$, the convolution $a*b$ is a vector $c\in\mathbb{R}^{m+n-1}$, with $c(k)=\sum_{j=\max(1,k+1-n)}^{\min(k,m)}a(j)b(k-j+1)$. For matrix $A\in\mathbb{R}^{m\times n}$ with elements $A(i,j)$’s, $\text{vec}(A)\in\mathbb{R}^{mn}$ stacks the columns of $A$ to a vector. Given another matrix $B\in\mathbb{R}^{m\times n}$, the Hadamard product is $A\odot B=[A(i,j)B(i,j)]$. The identity matrix is denoted $I$, and the conjugate transpose is denoted $(\cdot)^{\dagger}$.

The Fourier transform that maps from the spatial domain to the frequency domain is denoted ${\cal F}$, and ${\cal F}^{-1}$ is the inverse Fourier transform. For a variable $u$ in the spatial domain, its corresponding variable in the frequency domain is denoted $\tilde{u}$.

II Related Works

II-A Online Sparse Coding

Given $N$ samples $\{x_{i},\dots,x_{N}\}$, where each $x_{i}\in\mathbb{R}^{P}$, sparse coding learns an over-complete dictionary $D\in\mathbb{R}^{P\times L}$ of $L$ atoms and sparse codes $\{z_{i}\}$ [1]. It can be formulated as the following optimization problem:

[TABLE]

where \mathcal{D}=\{D:\left\|D(:,l)\right\|_{2}\leq 1\text{ for l=1,\dots,L}\} and $\beta\geq 0$. Many efficient algorithms have been developed for solving (1). Examples include K-SVD [1] and active set method [2]. However, they require storing all the samples, which can become infeasible when $N$ is large.

To solve this problem, an online learning algorithm for sparse coding that processes samples one at a time is proposed in [20]. After observing the $t$th sample $x_{t}$, the sparse code $z_{t}$ is obtained as

[TABLE]

where $D_{t-1}$ is the dictionary obtained at the $(t-1)$th iteration. After obtaining $z_{t}$, $D_{t}$ is updated as

[TABLE]

where

[TABLE]

Each column $D_{t}(:,l)$ in (4) can be obtained by coordinate descent. $A_{t}^{\text{(osc)}}$ and $B_{t}^{\text{(osc)}}$ can also be updated incrementally as

[TABLE]

Using $A_{t}^{\text{(osc)}}$ and $B_{t}^{\text{(osc)}}$, one does not need to store all the samples and codes to update $D_{t}$. The whole algorithm is shown in Algorithm 1.

The following assumptions are made in [20].

Assumption 1**.**

1. (A)

Samples $\{x_{i}\}$ are generated i.i.d. from some distribution with $\|x_{i}\|_{2}$ bounded. 2. (B)

The code $z_{t}$ is unique w.r.t. data $x_{t}$. 3. (C)

The objective in (4) is strictly convex with lower-bounded Hessians.

Theorem 1** ([20]).**

With Assumption 1, the distance between $D_{t}$ and the set of stationary points of the dictionary learning problem converges almost surely to 0 when $t\rightarrow\infty$.

II-B Alternating Direction Method of Multipliers (ADMM)

ADMM [15] has been popularly used for solving optimization problems of the form

[TABLE]

where $f,g$ are convex functions, and $A,B$ (resp. $c$) are constant matrices (resp. vector). It first constructs the augmented Lagrangian of problem (8)

[TABLE]

where $\nu$ is the dual variable, and $\rho>0$ is a penalty parameter.

At the $\tau$th iteration, the values of $x$ and $y$ (denoted as $x_{\tau}$ and $y_{\tau}$) are updated by minimizing (9) w.r.t. $x$ and $y$ in an alternating manner. Define the scaled dual variable $u_{\tau}={\nu_{\tau}}/\rho$. The ADMM updates can be written as

[TABLE]

The above procedure converges to the optimal solution at a rate of $O(1/T)$ [22], where $T$ is the number of iterations.

II-C Convolutional Sparse Coding

Convolutional sparse coding (CSC) learns a dictionary $D\in\mathbb{R}^{M\times K}$ composed of $K$ filters, each of length $M$, that can capture the same local pattern at different translated positions of the samples. This is achieved by replacing the multiplication between dictionary and code by convolution. While each $x_{i}$ in sparse coding is represented by a single code $z_{i}\in\mathbb{R}^{K}$, each $x_{i}$ in CSC is represented by $K$ codes stored together in the matrix $Z_{i}\in\mathbb{R}^{P\times K}$.

The dictionary and codes are obtained by solving the following optimization problem:

[TABLE]

where $*$ denotes convolution in the spatial domain.

Convolution can be accelerated in the frequency domain via the convolution theorem [23]: $\mathcal{F}({D}(:,k)*{Z}(:,k))=\mathcal{F}({D}(:,k))\odot\mathcal{F}({Z}(:,k))$, where ${D}(:,k)$ is first zero-padded to $P$-dimensional. Hence, recent CSC methods [5, 8, 16, 17] choose to operate in the frequency domain. Let $\tilde{x}_{i}\equiv\mathcal{F}(x_{i})$, $\tilde{D}(:,k)\equiv\mathcal{F}(D(:,k))$ and $\tilde{Z}_{i}(:,k)\equiv\mathcal{F}(Z_{i}(:,k))$. (12) is reformulated as

[TABLE]

where the factor $\frac{1}{P}$ in the objective comes from the Parseval’s theorem [24], and $\mathcal{H}$ is the linear operation that removes the extra $P-M$ dimensions in $\mathcal{F}^{-1}(\tilde{D}(:,k))$.

Problem (II-C) can be solved by block coordinate descent [5, 8, 16, 17], which updates $\{\tilde{Z}_{i}\}$ and $\tilde{D}$ alternately.

II-C1 Updating $\{\tilde{Z}_{i}\}$

Given $\tilde{D}$, the $\{\tilde{Z}_{i}\}$ can be obtained one by one for each $i$th sample as

[TABLE]

where $U_{i}$ is introduced to decouple the loss and the $\ell_{1}$-regularizer in (II-C). This can then be solved by ADMM [5, 8, 16, 17].

II-C2 Updating $\tilde{D}$

Given $\{\tilde{Z}_{i}\}$, $\tilde{D}$ can be obtained as

[TABLE]

where $V$ is introduced to decouple the loss and constraint in (II-C). This can again be solved by using ADMM [5, 8, 16, 17].

After obtaining $\{\tilde{Z}_{i}\}$ and $\tilde{D}$, the sparse codes can be recovered as $Z_{i}(:,k)=\mathcal{F}^{-1}(\tilde{Z}_{i}(:,k))$ for $i=1,\dots,N$, and the dictionary filters as $D(:,k)=\mathcal{H}(\mathcal{F}^{-1}(\tilde{D}(:,k)))$.

The above algorithms all need $O(NPK)$ in space. They differ mainly in how to compute the linear system involved with $\tilde{D}$ in the ADMM subproblems. FCSC [5] directly solves the subproblem, which takes $O(NK^{3}P)$ time. CBPDN [16] exploits a special structure in the dictionary and reduces the time complexity to $O(N^{2}KP)$, which is efficient for small $N$. The CONSENSUS algorithm [17] utilizes the matrix inverse lemma to reduce the time complexity to $O(NKP^{2})$. The state-of-the-art is FFCSC [8], which incorporates various linear algebra techniques (such as Cholesky factorization [25] and cached factorization [25]) to reduce the time complexity to $O(NK^{2}P)$.

III Online Convolutional Sparse Coding

Existing CSC algorithms operate in the batch mode, and need to store all the samples and codes which cost $O(NPK)$ space. This becomes infeasible when the data set is large. In this section, we will scale up CSC by using online learning.111After the initial arXiv posting of our paper [26], we became aware of some very recent independent works that also consider CSC in the online setting [27, 28]. These will be discussed in Section III-E.

After observing the $t$th sample $x_{t}$, online CSC considers the following optimization problem which is analogous to (12):

[TABLE]

To solve problem (15), some naive approaches are first considered in Section III-A. The proposed online convolutional sparse coding algorithm is then presented in Section III-B. It takes the same time complexity for one data pass as state-of-the-art batch CSC algorithms, but has a much lower space complexity (Section III-C). The convergence properties of the proposed algorithm is discussed in Section III-D.

III-A Naive Approaches

As in batch CSC, problem (15) can be solved by alternating minimization w.r.t. the codes and dictionary (as in Section II-C). Given the dictionary, the codes are updated as in (13). Given the codes, the dictionary is updated by solving the following optimization subproblem analogous to (14):

[TABLE]

However, solving (16) as in Section II-C2 requires keeping all the samples and codes, and is computationally expensive on large data sets.

Alternatively, the objective in (16) can be rewritten as

[TABLE]

where

[TABLE]

and $\dot{z}_{i}=\text{vec}(\tilde{Z}_{i})$. This is of the form in (3). Hence, we may attempt to reuse the online sparse coding in Algorithm 1, and thus avoid storing all the samples and codes. However, recall that online learning the dictionary is possible because one can summarize the observed samples into the history matrices $A_{t}^{\text{(osc)}},B_{t}^{\text{(osc)}}$ in (5), (6). For (17), the history matrices become

[TABLE]

In typical CSC applications, the number of image pixels $P$ is at least in the tens of thousands, and $N$ may only be in the thousands. Hence, the $O(K^{2}P^{2})$ space required for storing ${A}_{t}^{\text{(naive)}}$ and ${B}_{t}^{\text{(naive)}}$ is even higher than the $O(NKP)$ space required for batch methods.

III-B Proposed Algorithm

Note that $\dot{D}$ in (18) is composed of a number of diagonal matrices. By utilizing this special structure, the following Proposition rewrites the objective in (16) so that much smaller history matrices can be used.

Proposition 2**.**

The objective in (16) is equivalent to the following apart from a constant:

[TABLE]

The proof is in Appendix A-A. Obviously, each $\tilde{D}(p,:)$ in (20) can then be independently optimized. This avoids directly handling the much larger $P\times PK$ matrix $\dot{D}$ in (18). Let

[TABLE]

in (20). The total space required for $\{{A}^{p}_{t},{b}^{p}_{t}\}$ is $O(K^{2}P)$, which is much smaller than the $O(K^{2}P^{2})$ space for storing ${A}_{t}^{\text{(naive)}}$ in (19). Moreover, as in (7), ${A}^{p}_{t}$ and ${b}^{p}_{t}$ can be updated incrementally as

[TABLE]

The dictionary and codes can then be efficiently updated in an alternating manner as follows.

III-B1 Updating the Dictionary

With the codes fixed, using Proposition 2, the dictionary can be updated by solving the following optimization problem:

[TABLE]

This can solved using ADMM. At the $\tau$th ADMM iteration, let $\Theta_{t,\tau}$ be the ADMM dual variable. The following shows the update equations for $\tilde{D}$ and $V$.

Updating $\tilde{D}_{t,\tau}$: From (20), $\tilde{D}_{t,\tau}$ can be updated by solving the following subproblem:

[TABLE]

where $\tilde{V}_{t,\tau-1}(:,k)\equiv\mathcal{F}({V}_{t,\tau-1}(:,k))$. Note that $\|\tilde{D}\|_{F}^{2}=\sum_{p=1}^{P}\|\tilde{D}(p,:)\|_{2}^{2}=\sum_{k=1}^{K}\|\tilde{D}(:,k)\|_{2}^{2}$. Hence, $\sum_{k=1}^{K}\|\tilde{D}(:,k)-\tilde{V}_{t,\tau-1}(:,k)+\Theta_{t,\tau-1}(:,k)\|_{2}^{2}=\sum_{p=1}^{P}\|\tilde{D}(p,:)-\tilde{V}_{t,\tau-1}(p,:)+\Theta_{t,\tau-1}(p,:)\|_{2}^{2}$. The rows of $\tilde{D}_{t,\tau}$ can then be obtained separately as

[TABLE]

where ${C}^{p}_{t}=({A}^{p}_{t}+\rho PI)^{-1}$. With $\{C^{p}_{t}\}$, we do not need to store $\{A^{p}_{t}\}$.

Computing the matrix inverse ${C}^{p}_{t}\in\mathbb{C}^{K\times K}$ takes $O(K^{3})$ time. This can be simplified by noting from (21) that ${A}^{p}_{t}+\rho PI$ is the sum of rank-1 matrices and a (scaled) identity matrix. Using the Sherman-Morrison formula222Given an invertible square matrix $A$ and vectors $u,v$, $\left(A+uv^{\top}\right)^{-1}=A^{-1}-\left(1+v^{\top}A^{-1}u\right)^{-1}\left(A^{-1}uv^{\top}A^{-1}\right)$. [25], we have

[TABLE]

This takes $O(K^{2})$, instead of $O(K^{3})$, time.

Updating $V_{t,\tau}$: From (10), each column $V_{t,\tau}(:,k)$ can be updated as

[TABLE]

It has the following closed-form solution.

Proposition 3** ([29]).**

$V_{t,\tau}(:,k)=\alpha/\max(\|\alpha\|_{2},1)$, where $\alpha=\mathcal{H}(\mathcal{F}^{-1}(\tilde{D}_{t,\tau}(:,k)+\Theta_{t,\tau-1}(:,k)))$.

Finally, the dual variables $\Theta_{t,\tau}$ is updated as in (11). The whole dictionary update procedure (DictOCSC) is shown in Algorithm 2. As (23) is convex, convergence to the globally optimal solution is guaranteed [22].

In batch CSC methods, its dictionary update in (14) is also based on ADMM (Section II-C2). However, our dictionary update step first reformulates the objective as in (20). This enables each ADMM subproblem to be solved with a much lower space complexity ($O(K^{2}P)$ vs $O(NKP)$ for the state-of-the-art [8]) but still with the same iteration time complexity (i.e., $O(NK^{2}P+NKP\log P)$).

III-B2 Updating the Code

Given the dictionary, as the codes for different samples are independent, they can be updated one by one as in batch CSC methods (Section II-C1).

The whole algorithm, which will be called Online Convolutional Sparse Coding (OCSC), is shown in Algorithm 3.

III-C Complexity Analysis

In Algorithm 3, the space requirement is dominated by $\{C^{p}_{t}\}$, which takes $O(K^{2}P)$ space. For one data pass which precesses $N$ samples, updating $\{C^{p}_{t}\}$ and $\{b^{p}_{t}\}$ takes $O(NK^{2}P)$ time, dictionary update takes $O(NK^{2}P)$ time, code update takes $O(NKP)$ time, and FFT/inverse FFT takes $O(NKP\log P)$ time. Hence, one data pass takes a total of $O(NK^{2}P+NKP\log P)$ time.

A comparison with the existing batch CSC algorithms is shown in Table I. As can be seen, the proposed algorithm takes the same time complexity for one data pass as the state-of-the-art FFCSC algorithm, but has a much lower space complexity ($O(K^{2}P)$ instead of $O(NKP)$).

III-D Convergence

In this section, we show that Algorithm 3 outputs a stationary point of the CSC problem (15) when $t\rightarrow\infty$. This is achieved by connecting Algorithm 3 to a direct application of Algorithm 1 on (15).

The convolution operation in the spatial domain can be written as matrix multiplication [5, 14]. Specifically,

[TABLE]

where $\mathcal{T}(x)$ is a linear operator which maps a vector to its associated Toeplitz matrix. Specifically,

[TABLE]

The number of columns in $\mathcal{T}(D(:,k))$ is equal to the dimension of $Z(:,k)$ (i.e., $P$).

Let $\bar{D}_{k}=\mathcal{T}(D(:,k))\in\mathbb{R}^{P\times P}$, $\mathcal{T}^{-1}$ be the inverse operator of $\mathcal{T}$ which maps $\bar{D}_{k}$ back to $D(:,k)$, $\bar{D}=\left[\bar{D}_{1},\dots,\bar{D}_{K}\right]$ and $\bar{z}_{i}=\text{vec}(Z_{i})$. Problem (15) can be rewritten as

[TABLE]

where

[TABLE]

Thus, the objective (28) is of the same form as that in (3). However, a direct use of Algorithm 1 is not feasible. First, the feasible region $\mathcal{S}$ in (29) is more complex, and coordinate descent cannot be used as there is no simple projection to $\mathcal{S}$. Second, as $\bar{z}_{i}$ is of length $KP$, the corresponding history matrices (analogous to those in (5)) require $O(K^{2}P^{2})$ space.

Though a direct application of Algorithm 1 is not practical, Theorem 1 still holds. Indeed, Theorem 1 can be further extended by relaxing its feasible region $\mathcal{C}$ on $D$. As discussed in [20], $\mathcal{C}$ can be, for example, $\{D:\|D(:,l)\|_{2}\leq 1\;\text{and}\;D(i,j)\geq 0\}$. It is mentioned in [20] that $\mathcal{C}$ has to be a union of independent constraints on each column of $D$. However, this only serves to facilitate the use of coordinate descent (step 6 in Algorithm 1), but is not required in the proof. In general, Theorem 1 holds when $\mathcal{C}$ is bounded, convex, and $\mathcal{C}\subseteq\mathcal{D}$.

The following Lemma shows that $\mathcal{S}$ in (29) satisfies the conditions. The proof is in Appendix A-B. Thus, Theorem 1 also holds for Algorithm 3, and a stationary point of problem (15) can be obtained.

Lemma 4**.**

$\mathcal{S}$* in (29) is bounded, convex, and a subset of $\mathcal{D}$.*

III-E Discussion with [27, 28]

Here, we discuss the very recent works of [27, 28] which also consider online learning of the dictionary in CSC. Extra experiments are performed in Section B, which shows our method is much faster than them.

In the online convolutional dictionary learning (OCDL) algorithm [27], convolution is performed in the spatial domain. They started with the observation that convolution is commutative. Hence, for the summation $\sum_{k=1}^{K}D(:,k)*Z_{i}(:,k)$ in (15), we have

[TABLE]

where $\bar{Z}_{i}=\left[\mathcal{T}(Z_{i}(:,1)),\dots,\mathcal{T}(Z_{i}(:,K))\right]$ with $\mathcal{T}(Z_{i}(:,k))\in\mathbb{R}^{P\times M}$ and $\bar{d}=\text{vec}(D)$. (15) can then be rewritten as

[TABLE]

where $\lambda=\beta/P$ (as each $Z_{i}(:,k)$ is repeated $P$ times in the Toeplitz matrix $\bar{Z}_{i}$). Using (31), the history matrices are constructed as

[TABLE]

Recall that $M$ is the length of filter $D(:,k)$ when CSC is solved in the spatial domain. The space complexity of [27] is dominated by $\bar{Z}_{i}$ (which takes $O(KPM)$ space) or ${A}_{t}^{\text{(ocdl1)}}$ (which takes $O(K^{2}M^{2})$ space), depending on the relative sizes of $P$ and $KM$. Though this is comparable to our $O(K^{2}P)$ space requirement, its time complexity is much larger. For one data pass, convolution in the spatial domain takes $O(NKPM)$ time and updating the history matrices above takes $O(NK^{2}PM^{2})$ time. The dictionary update in total takes $O(NK^{2}PM^{2}+NKPM)$ time. In contrast, the proposed algorithm takes $O(NK^{2}P+NKP\log P)$ time. In the experiments, $M=11\times 11$, and $P$ ranges from $32\times 32$ to $500\times 500$. Thus, the algorithm in [27] is much more expensive.

The algorithm in [28], also called online convolutional dictionary learning, considers the frequency domain and solves problem (16) as in the proposed method. First, they rewrite (16) as

[TABLE]

where $\dot{Z}_{i}=[\text{Diag}(\tilde{Z}_{i}(:,1)),\dots,\text{Diag}(\tilde{Z}_{i}(:,K))]$, and $\dot{d}=\text{vec}(\tilde{D})$. The history matrices are then constructed as

[TABLE]

They proposed to store ${A}_{t}^{\text{(ocdl2)}}$ and the history matrices in sparse format. These take $O(K^{2}P)$ space, and is the same as our $\{C^{p}_{t}\}$ in (25) which are stored in dense format. However, though their $O(\cdot)$ are the same, storing sparse matrix requires 2-3 times more space than storing dense matrix with the same number of nonzero entries, as each nonzero entry in a sparse matrix needs to be kept in the compressed sparse row (CSR)333For an introduction to the CSR format, interested readers are referred to http://www.netlib.org/utk/people/JackDongarra/etemplates/node373.html. format [30]. Moreover, though [28] and the proposed method take $O(NK^{2}P)$ time to update the history matrices, using sparse matrices as in [28] is empirically slower [30]. Preliminary experiments show that with $K=100$ filters, the proposed algorithm is $5$ times faster on an $P=100\times 100$ image and $10$ times faster on an $P=200\times 200$ image.

Moreover, we use ADMM for optimization, while [27, 28] use the FISTA algorithm [31]. Empirically, ADMM has been shown to be faster than FISTA on solving the CSC problem [16]. We also empirically verify this point in Section B-B.

IV Experiments

In this section, experiments are performed on six image data sets (Table II). The Fruit, City and House data sets444http://www.cs.ubc.ca/labs/imager/tr/2015/FastFlexibleCSC/. are standard benchmarks in CSC [5, 8, 14]. However, they are relatively small. Hence, we include three larger data sets: (i) Flower,555http://www.robots.ox.ac.uk/~vgg/data/flowers/102/. which contains images from 102 flower categories [32]; (ii) Dog,666http://vision.stanford.edu/aditya86/ImageNetDogs/. which contains images of 120 breeds of dogs [33]; and (iii) CIFAR-10,777https://www.cs.toronto.edu/~kriz/cifar.html. which contains images from 10 general object classes [34]. The training images are used to learn the dictionary, which is then used to reconstruct the test set images. We use the default training and test splits provided. The images are preprocessed as in [8, 14]. We convert each image to grayscale, and perform local contrast normalization. Edge-tapering is used to blur the edges of the samples with a random Gaussian filter.

The proposed OCSC is compared with the following (batch) CSC methods:

Deconvolutional networks (DeconvNet) [14]; 2. 2.

Fast convolutional sparse coding (FCSC) [5]; 3. 3.

Fast and flexible convolutional sparse coding (FFCSC) [8]; 4. 4.

Convolutional basis pursuit denoising (CBPDN) [16]; 5. 5.

The global consensus ADMM (CONSENSUS) algorithm [17].

All the codes are in Matlab and obtained from the respective authors (except FCSC).888DeconvNet is from http://www.matthewzeiler.com/, FFCSC from http://www.cs.ubc.ca/labs/imager/tr/2015/FastFlexibleCSC/, CBPDN from http://brendt.wohlberg.net/software/SPORCO/, and CONSENSUS is from http://zoi.utia.cas.cz/convsparsecoding. The code of FCSC is not available, so we use the code in [17] (http://zoi.utia.cas.cz/convsparsecoding). We do not compare with [27, 28], as their codes are not publicly available.

We follow the hyperparameter setting in [8], and set $\beta=1$, $M=11$ and $K=100$. The batch methods are stopped when the relative changes in $\{Z_{i}\}$ (i.e., $\frac{1}{N}\sum_{i=1}^{N}\|Z_{i}-Z_{i}^{\text{old}}\|_{2}/\|Z_{i}\|_{2}$) and $D$ (i.e., $\|D-D^{\text{old}}\|_{2}/\|D\|_{2}$) are both smaller than ${10}^{-3}$. As for the proposed online method, it is stopped when the relative changes in $Z_{t}$ (i.e., $\|Z_{t}-Z_{t}^{\text{old}}\|_{2}/\|Z_{t}\|_{2}$) and $D$ (i.e., $\|D-D^{\text{old}}\|_{2}/\|D\|_{2}$) are both smaller than ${10}^{-3}$. Experiments are run on a PC with Intel i7 4GHz CPU with 32GB memory.

To compare the empirical convergence of various methods, we monitor the test set objective as in [20, 21]. This is obtained as

[TABLE]

where $\Omega$ is the test set, and $D_{t}$ is the dictionary learned after the $t$th data pass. As for image reconstruction quality, we use the peak signal-to-noise ratio on the test set [8]:

[TABLE]

where ${x}^{\text{rec}}_{i}=\sum_{k=1}^{K}D_{T}(:,k)*Z_{i}(:,k)$ is the reconstructed image for $x_{i}$ by using the final dictionary $D_{T}$. To reduce statistical variability, we repeat the experiment five times by using different dictionary initializations and orders to present images to the algorithm. The PSNR is then averaged over these five repetitions.

IV-A Small Data Sets

In this section, we perform experiments on the small data sets (Fruit, City and House). Figure 1 shows the test set objective vs CPU time. Though FFCSC has the same time complexity as OCSC, OCSC is still much faster. This is because batch CSC methods update the dictionary only after coding all the samples, while OCSC can refine the dictionary after coding each sample. A similar behavior is also observed between online and batch SC methods [20].

Table III compares the image reconstruction performance. As OCSC converges to a lower objective than its batch counterparts (Figure 1), it also outperforms the others in terms of image reconstruction quality. Among all methods, OCSC performs the best, which is then followed by CBPDN, FCSC and FFCSC. CONSENSUS and DeconvNet perform the worst.

Figure 2 shows the difference between the reconstructed images and the corresponding ground truths of some sample test set images from the Fruit data set. As can be seen, the images reconstructed by OCSC are most similar to the ground truths, which is then followed by CBPDN, FFCSC and FCSC. Images reconstructed by CONSENSUS and DeconvNet are the least similar to the ground truths. A similar observation is also observed on City and House, which are not reported here because of the lack of space.

Figure 3 shows the dictionaries learned on Fruit (results on City and House are similar). As in [1], the filters are sorted in ascending order of the variance.999For filter $D(:,k)$, its variance is defined as $\frac{1}{M-1}\sum_{m=1}^{M}(D(m,k)-\frac{1}{M}\sum_{m=1}^{M}D(m,k))^{2}$. As can be seen, dictionaries learned by FCSC, FFCSC, CBPDN and OCSC contain Gabor-style filters, while the dictionaries learned by DeconvNet and CONSENSUS are vague. We speculate that this is due to slow convergence of these two methods, as observed in Figure 1.

IV-B Large Data Sets

As discussed in Section III-C, the space complexity of OCSC is independent of $N$. Hence, it can handle much larger training image sets than batch CSC methods. In this section, we illustrate this on the Flower, Dog and CIFAR-10 data sets. The Flower images are large ($500\times 500$), and the batch CSC methods quickly run out of memory.101010For $K=75$ or 100, FFCSC can only work with a single image, while the other batch methods can only handle 2-5 images. To allow comparison with the batch methods, we use $K$ as 50.

Figure 4(a) shows the testing PSNR with varying number of training images on the Flower data set. As can be seen, performance is improved with more training images. However, the batch CSC methods run out of memory quickly. Specifically, CBPDN can only handle 25 training images, while the other batch methods can only handle 20. Figures 4(b) and 4(c) show results on Dog and CIFAR-10, respectively. Among the batch methods, CBPDN can handle 40 training images on Dog and 2000 training images on CIFAR-10, while the others can only handle 20 training images on Dog and 1500 training images on CIFAR-10. On the other hand, OCSC can always handle the full training data set.

Figure 5 shows the CPU time with different numbers of training images. Given the same amount of time, OCSC can train on more images than the batch methods. Combined with Figure 4, OCSC obtains higher testing PSNR than batch methods given the same amount of time.

V Conclusion

In this paper, we proposed a scalable convolutional sparse coding methods in the online setting. By reformulating the CSC objective, the sizes of the history matrices required in the online setup can be significantly reduced. Moreover, the resultant optimization problem can be efficiently solved using ADMM, with closed-form solutions for the ADMM subproblems. We also provide theoretical guarantee on the convergence of the proposed algorithm. While the batch CSC methods need large space and heavy computation cost, extensive experiments show that the proposed method can efficiently learn from large data sets with much less space.

As for future work, we will consider introducing nonconvex regularizers to CSC. Traditionally, the convex $\ell_{1}$-regularizer is used to encourage sparsity among the codes. Though this leads to easier optimization, the resultant code may not be as sparse and accurate as when a nonconvex regularizer (such as the log-sum-penalty) is used [35]. Moreover, on large data sets, many filters are often needed to capture the presence of more local patterns. Though the proposed online algorithm has reduced the space consumption significantly compared to batch CSC methods, its space complexity still depends on the number of filters. A promising direction is to further reduce the storage by approximating the filters as linear combinations of a few base filters.

Appendix A Proofs

A-A Proposition 2

Proof.

(17) can be expanded as

[TABLE]

where $\dot{D}=[\text{Diag}(\tilde{D}(:,1)),\dots,\text{Diag}(\tilde{D}(:,K))]$, and $\dot{z}_{i}=\text{vec}(\tilde{Z}_{i})$. Thus,

[TABLE]

For $\tilde{D}(j,k)$, only when it is multiplied by elements of the same row index $\tilde{D}(j,l)$, $l=1,\dots,K$, will result in non-zero values. In other words, multiplications among different rows can be avoided. We then directly operate on each row of $\dot{D}$, which is exactly $\tilde{D}(p,:)$ if zeros are dropped. Consequently, we have the form in Proposition 2. ∎

A-B Lemma 4

Proof.

Obviously, $\mathcal{S}$ is a convex constraint set. To show $\mathcal{S}\subset\mathcal{D}$, recall that $\bar{D}=[T(D(:,1)),\dots,T(D(:,K))]$, where each $T(D(:,k))$ is the Toeplitz matrix associated with $D_{t}(:,k)$ as defined in (26). Because of the constraint $\|D(:,k)\|_{2}^{2}\leq 1$, each column $\bar{D}(:,l)$ of $\bar{D}$, satisfies $\|\bar{D}(:,l)\|_{2}^{2}\leq 1$. Thus $\|\bar{D}\|_{F}^{2}\leq KP$. Hence, $\mathcal{S}\subset\mathcal{D}$ and $\mathcal{S}$ is bounded. ∎

Appendix B Experimental Comparison with [27, 28]

In this section, we provide extra experimental results comparing OCSC with concurrent works OCDL1 [27] and OCDL2 [28] discussed in Section III-E. We use implementation by ourself as no pubic codes are available.

B-A Large Data Sets

We experiment on large data sets to test these online methods’ scalability. The evaluation metric used on empirical convergence is test set objective, on image reconstruction quality is PSNR. CIFAR-10 is used, which has a large number of small images ($P=32\times 32$). This is because a single image of Flower ($P=500\times 500$) makes both OCDL1 and OCDL2 run out of memory, even Dog ($P=224\times 224$) makes OCDL2 not work. Hence for complete comparison, we use CIFAR-10 only.

We still use the same training and test splits as in Section IV. Hyper-parameters are the same as Section IV. As a recap, we set $\beta=1$, $M=11$ and $K=100$. The methods are stopped when the relative changes in $Z_{t}$ (i.e., $\|Z_{t}-Z_{t}^{\text{old}}\|_{2}/\|Z_{t}\|_{2}$) and $D$ (i.e., $\|D-D^{\text{old}}\|_{2}/\|D\|_{2}$) are both smaller than ${10}^{-3}$.

Figure 6(a) shows the test set objective vs CPU time. As shown, OCSC converges the fastest, then OCDL2 follows, OCDL1 is the slowest. OCDL1 lags behinds due to its slow convolution performed in spatial domain and the costly updating of history matrices. As for OCDL2, although it shares the efficiency of performing convolution in frequency domain, it is very slow to deal with sparse matrices. They finally reach to comparable testing objective value as they are solving the same online CSC problem. Figure 6(b) shows the testing PSNR vs CPU time. For all three methods, the testing PSNR is improved given more training images, and they can reach similar final testing PSNR. However, the CPU time needed to reach it varies a lot. In other words, given same amount of time, OCSC always gets a higher testing PSNR.

B-B Synthetic Data

In this section, we perform experiments on synthetic data to compare the convergence speed of ADMM (used by our OCSC) and FISTA (used by both OCDL1 [27] and OCDL2 [28]) on CSC problem. We evaluate the convergence speed by training set objective vs time. We choose frequency domain to compare the efficiency of ADMM and FISTA. As shown in Figure 6(a) and Figure 6(b), although OCDL1 and OCDL2 both use FISTA, OCDL1 which operates in spatial domain is much slower than OCDL2 which operates in frequency domain. Specifically, we compare proposed ADMM on (23) with FISTA on (III-E) [28]. Note that the two problems are equivalent to (16), thus we compute the objective based on (16).

We set $K=100$, and the size of synthetic data $x\in\mathbb{R}^{P}$ where $P=10000$. The entries of $x$, $D$ and $Z$ are sampled i.i.d. from the standard normal distribution $\mathcal{N}(0,1)$, and each $D(:,:,k),k=1\dots K$ is projected to the $\ell_{2}$-norm unit ball.

The comparison of objective vs iterations and CPU time are shown in Figure 7(a) and Figure 7(b). In terms of the number of iterations, ADMM converges faster than FISTA. Then, in terms of time, FISTA is much slower than ADMM. This relates to the design of OCDL2 mentioned in [28]. Empirically, for our data, OCDL2’s history matrices need hundreds of GBs in dense format. Hence to fit into memory, OCDL2 has to operate on sparse matrices, which is much slower on matrix multiplications. This also verifies our discussion in Section III-E, ADMM converges faster than FISTA on CSC problem.