On learning with shift-invariant structures

TL;DR

This paper introduces new algorithms for learning shift-invariant structures and calculating shifts between signals using circulant matrices, with applications to dictionary learning and signal alignment.

Contribution

It presents novel methods for shift-invariant dictionary learning and shift retrieval, including algorithms for convolutional and circulant matrices, and demonstrates their effectiveness on real data.

Findings

Algorithms successfully recover shift-invariant components.

Effective shift retrieval from synthetic and real signals.

Learned wavelet-like dictionaries from data.

Abstract

We describe new results and algorithms for two different, but related, problems which deal with circulant matrices: learning shift-invariant components from training data and calculating the shift (or alignment) between two given signals. In the first instance, we deal with the shift-invariant dictionary learning problem while the latter bears the name of (compressive) shift retrieval. We formulate these problems using circulant and convolutional matrices (including unions of such matrices), define optimization problems that describe our goals and propose efficient ways to solve them. Based on these findings, we also show how to learn a wavelet-like dictionary from training data. We connect our work with various previous results from the literature and we show the effectiveness of our proposed algorithms using synthetic, ECG signals and images.

Equations

$$\mathbf{C}=\text{circ}(\mathbf{c})\stackrel{\mathrm{def}}{=}$$

$$=$$

$$\mathbf{C}={\mathbf{F}}^H\mathbf{\Sigma F},\mathbf{\Sigma }=\text{diag}(\sigma )\in {\mathbb{C}}^{n\times n},$$

$$\mathbf{y}=$$

$$=$$

$${\mathbf{C}}_k^{(p)}=$$

$${\mathbf{W}}_k^{(p)}=\left[\begin{array}{cc}{\mathbf{C}}_k^{\left(\frac{p}{2^{k-1}}\right)}& 0_{\frac{p}{2^{k-1}}\times \left(p-\frac{p}{2^{k-1}}\right)}\\ 0_{\left(p-\frac{p}{2^{k-1}}\right)\times \frac{p}{2^{k-1}}}& {\mathbf{I}}_{p-\frac{p}{2^{k-1}}}\end{array}\right]\in {\mathbb{R}}^{p\times p}.$$

$$\mathbf{W}={\mathbf{W}}_1^{(p)}\cdots {\mathbf{W}}_{m-1}^{(p)}{\mathbf{W}}_m^{(p)}.$$

$$\underset{\mathbf{D},\mathbf{X}}{\mathrm{minimize}}\Vert \mathbf{Y}-\mathbf{DX}{\Vert }_F^2\text{subject to }\mathbf{X}\text{ is s–sparse}.$$

$$\underset{\mathbf{c},\mathbf{X}}{\mathrm{minimize}}$$

$$\Vert \text{vec}(\mathbf{X}){\Vert }_0\le sN,\mathbf{C}=\text{circ}(\mathbf{c}),\Vert \mathbf{c}{\Vert }_2=1.$$

$$\Vert \mathbf{Y}-\mathbf{CX}{\Vert }_F^2=\Vert \mathbf{FY}-\mathbf{\Sigma FX}{\Vert }_F^2=\Vert \tilde{\mathbf{Y}}-\mathbf{\Sigma }\tilde{\mathbf{X}}{\Vert }_F^2,$$

$${\sigma }_1=\frac{{\tilde{\mathbf{x}}}_1^H{\tilde{\mathbf{y}}}_1}{\Vert {\tilde{\mathbf{x}}}_1{\Vert }_2^2},{\sigma }_k=\frac{{\tilde{\mathbf{x}}}_k^H{\tilde{\mathbf{y}}}_k}{\Vert {\tilde{\mathbf{x}}}_k{\Vert }_2^2},{\sigma }_{n-k+2}={\sigma }_k^{\ast },k=2,\dots ,n,$$

$$\underset{\mathbf{c}}{\mathrm{minimum}}\Vert \mathbf{Y}-\mathbf{CX}{\Vert }_F^2=\sum _{k=1}^n\left(\Vert {\mathbf{y}}_k{\Vert }_2^2-\frac{|{\tilde{\mathbf{x}}}_k^H{\tilde{\mathbf{y}}}_k{|}^2}{\Vert {\tilde{\mathbf{x}}}_k{\Vert }_2^2}\right).$$

$$\Vert$$

$$=$$

$$\mathbf{D}=\left[\begin{array}{cccc}{\mathbf{C}}^{(1)}& {\mathbf{C}}^{(2)}& \dots & {\mathbf{C}}^{(L)}\end{array}\right]\in {\mathbb{R}}^{n\times nL},$$

$$\Vert \mathbf{Y}$$

$$=$$

$$\underset{{\sigma }_k^{(1)},\dots ,{\sigma }_k^{(L)}}{\mathrm{minimize}}{\Vert {\tilde{\mathbf{y}}}_k^T-\sum _{\ell =1}^L{\sigma }_k^{(\ell )}({\tilde{\mathbf{x}}}_k^{(\ell )}{)}^T\Vert }_F^2.$$

$${\mathbf{D}}_{\text{conv}}=\left[\begin{array}{cccc}{\mathbf{C}}_{\text{conv}}^{(1)}& {\mathbf{C}}_{\text{conv}}^{(2)}& \dots & {\mathbf{C}}_{\text{conv}}^{(L)}\end{array}\right]\in {\mathbb{R}}^{p\times pL},$$

$$\Vert \mathbf{Y}-$$

$$=$$

$$\mathbf{c}=({\mathbf{B}}^H\mathbf{B})\backslash {\mathbf{B}}^H\tilde{\mathbf{y}},$$

$$\Vert \mathbf{Y}-{\mathbf{C}}_{\text{conv}}{\mathbf{X}}_{\text{conv}}{\Vert }_F^2=\Vert \tilde{\mathbf{y}}-\mathbf{A}{\mathbf{F}}_{:,1:n}\mathbf{c}{\Vert }_F^2,$$

$$\mathbf{c}=(\mathbf{A}{\mathbf{F}}_{:,1:n})\backslash \tilde{\mathbf{y}}.$$

$$\mathbf{W}=\text{diag}(\left[\begin{array}{ccc}\Vert {\tilde{\mathbf{x}}}_1{\Vert }_2^2& \dots & \Vert {\tilde{\mathbf{x}}}_p{\Vert }_2^2\end{array}\right]),\mathbf{v}={\mathbf{F}}_{:,1:n}^H{\mathbf{A}}^H\tilde{\mathbf{y}}.$$

$$\mathbf{G}=\left[\begin{array}{cccc}{\mathbf{D}}^{(1)}\mathbf{C}& {\mathbf{D}}^{(2)}\mathbf{C}& \dots & {\mathbf{D}}^{(m)}\mathbf{C}\end{array}\right]\in {\mathbb{C}}^{m\times m^2},$$

$$\Vert \mathbf{Y}$$

$$=$$

Full text

I Introduction

Circulant matrices [1] are highly structured matrices where each column is a circular shift of the previous one. Because of their structure and their connection to the fast Fourier transform [2] (circulant matrices are diagonalized by the Fourier matrix [1, Section 3]), these matrices have seen many applications in the past: computing the shift between two signals (the shift retrieval problem exemplified in the circular convolution and cross-correlation theorems) for the GPS locking problem [3], time delay estimation [4], compressed shift retrieval from Fourier components [5], matching or alignment problems for image processing [6], designing numerically efficient linear transformations [7] and overcomplete dictionaries [8], matrix decompositions [9], convolutional dictionary learning [10, 11, 12] and sparse coding [13], learning shift invariant structured from data [14],[15] for medical imaging [16], EEG [17] and audio [18] signal analysis. A recent review of the methods, solutions and applications related to circulant and convolutional representations is given in [19].

In this paper, we propose several numerically efficient algorithms to extract shift-invariant components or alignments from data using several structured dictionaries related to circulant matrices.

Previously, several dictionary learning techniques that accommodate for shift invariance have been proposed: extending the well-known K-SVD algorithm to deal with shift-invariant structures [17, 20, 21], proposing a shift-invariant iterative least squares dictionary learning algorithm [22], extending the dictionary while solving an eigenvalue problem [23], fast online learning approach [24], research that combines shift and 2D rotation invariance [25] and proposing new algorithms that optimize directly the dictionary learning objective functions with circulant matrices [14, 15, 26]. The convolutional sparse representation model [27, 28, 29] where the dictionary is a concatenation of circulant matrices has been extensively studied in the past. Furthermore, recent work [13] uses tools developed in the sparse representations literature to provide theoretical insights into convolutional sparse coding where the dictionary is a concatenation of banded circulant matrices and its connection to convolutional neural networks [30]. Detailed literature reviews of these learning and convolutional sparse representations problems and proposed solutions have been recently described in [19, Section II].

Structured dictionaries have received a lot of attention mainly because of two reasons: the structure means that the dictionaries will be easier to store and use (lower memory footprint and lower computational complexity to perform, for example, matrix-vector multiplication or solving linear systems) and they act as regularizers modeling some property of the data that is interesting. In the case of shift-invariant dictionaries these two advantages are transparent: manipulation of circulant matrices is done via the fast Fourier transform while storing them takes only linear space (instead of quadratic) and they are able to model patterns from the data that are repetitive, as we expect real-world data (especially image data and time-series) to exhibit such patterns.

We start by outlining in Section II circulant matrices and their properties, particularly their factorization with the Fourier matrix, and other structured matrices that we use in this paper. Then, we propose algorithms to learn shift-invariant (circulant, convolutional and unions of these) and wavelet-like components from training data (Section III). Finally, in Section IV, we show experimental results with various data sources (synthetic, ECG, images) that highlight the learning capabilities of the proposed methods.

Notation: bold lowercase $\mathbf{x}\in\mathbb{R}^{n}$ is used to denote a column vector, bold uppercase $\mathbf{X}\in\mathbb{R}^{n\times m}$ is used to denote a matrix, non-bold lowercase Greek letters like $\alpha\in\mathbb{R}$ are used to denote scalar values, calligraphic letters $\mathcal{K}$ are used to denote sets and $|\mathcal{K}|$ is the cardinality of $\mathcal{K}$ (abusing this notation, $|\alpha|$ is the magnitude of a scalar). Then $\|\mathbf{x}\|_{2}$ is the $\ell_{2}$ norm, $\|\mathbf{x}\|_{0}$ is the $\ell_{0}$ pseudo-norm, $\|\mathbf{X}\|_{F}^{2}=\text{tr}(\mathbf{X}^{H}\mathbf{X})$ is the Frobenius norm, $\mathrm{tr}(\mathbf{X})$ denotes the trace, $\text{vec}(\mathbf{X})\in\mathbb{R}^{nm}$ vectorizes the matrix $\mathbf{X}\in\mathbb{R}^{n\times m}$ columnwise, $\text{diag}(\mathbf{x})$ denotes the diagonal matrix with the vector $\mathbf{x}$ on its diagonal, $\mathbf{X}^{H}$ is the complex conjugate transpose, $\mathbf{X}^{T}$ is the matrix transpose, $\mathbf{X}^{*}$ is the complex conjugate, $\mathbf{X}^{-1}$ denotes the inverse of a square matrix, $x_{kj}$ is the $(k,j)^{\text{th}}$ entry of $\mathbf{X}$. Tilde variables like $\mathbf{\tilde{X}}$ represents the columnwise Fourier transform of $\mathbf{X}$, $\mathbf{X}\otimes\mathbf{Y}$ denotes the Kronecker product and $\mathbf{X}\odot\mathbf{Y}$ is the Khatri-Rao product [31].

II The proposed structured dictionaries

II-A Circulant dictionaries

We consider in this paper circulant matrices $\mathbf{C}$. These square matrices are completely defined by their first column $\mathbf{c}\in\mathbb{R}^{n}$: every column is a circular shift of the first one. With a down shift direction the right circulant matrices are:

[TABLE]

The matrix $\mathbf{P}\in\mathbb{R}^{n\times n}$ denotes the orthonormal circulant matrix that circularly shifts a target vector $\mathbf{c}$ by one position, i.e., $\mathbf{P}=\text{circ}(\mathbf{e}_{2})$ where $\mathbf{e}_{2}$ is the second vector of the standard basis of $\mathbb{R}^{n}$. Notice that $\mathbf{P}^{q-1}=\text{circ}(\mathbf{e}_{q})$ is also orthonormal circulant and denotes a cyclic shift by $q-1$. The main property of circulant matrices (1) is their eigenvalue factorization which reads:

[TABLE]

where $\mathbf{F}\in\mathbb{C}^{n\times n}$ is the unitary Fourier matrix ($\mathbf{F}^{H}\mathbf{F}=\mathbf{F}\mathbf{F}^{H}=\mathbf{I}$) and the diagonal $\mathbf{\sigma}=\sqrt{n}\mathbf{Fc},\ \mathbf{\sigma}\in\mathbb{C}^{n}$.

Note that the multiplication with $\mathbf{F}$ on the right is equivalent to the application of the Fast Fourier Transform, i.e., $\mathbf{Fc}=\text{FFT}(\mathbf{c})$, while the multiplication with $\mathbf{F}^{H}$ is equivalent to the inverse Fourier transform, i.e., $\mathbf{F}^{H}\mathbf{c}=\text{IFFT}(\mathbf{c})$. Both transforms are applied in \operatorname{O}\bigl{(}n\log n\bigr{)} time and memory.

II-B Convolutional dictionaries

Convolutional dictionaries can be reduced to circulant dictionaries by observing that given $\mathbf{c}\in\mathbb{R}^{n}$ and $\mathbf{x}\in\mathbb{R}^{m}$ with $m\geq n$ the result of their convolution is a vector $\mathbf{y}$ of size $p=n+m-1$ as

[TABLE]

where $\text{toep}(\mathbf{c})$ is a Toeplitz matrix of size $p\times m$. Padding with zeros such that all variables are of size $p$ leads again to a matrix-vector multiplication by a circulant matrix. An alternative, but ultimately equivalent, way to write the convolution in terms of a circulant matrix is to notice that a Toeplitz matrix can be embedded into an extended circulant matrix of twice the size (see [15, Section 4]).

For our purposes, there is a fundamental difference between $\mathbf{C}$ and $\mathbf{C}_{\text{conv}}$: in the case of $\mathbf{C}$ it is exactly equivalent if we choose to operate with $\mathbf{c}$ or $\mathbf{\sigma}$ while for $\mathbf{C}_{\text{conv}}$ we necessarily have to work with $\mathbf{c}_{\text{conv}}$ in order to impose its sparse structure. As we will see, this means that in the convolutional case we cannot exploit some simplifications that occur in the Fourier domain (when working directly with $\mathbf{\sigma}$).

II-C Wavelet-like dictionaries

Multiples, powers, products, and sums of circulant matrices are themselves circulant. Therefore, extending this class of structures to include other dictionaries is not straightforward.

In order to represent a richer class of dictionaries, still based on circulants, consider now the following structured matrix

[TABLE]

where $\mathbf{G}_{k}=\text{circ}(\mathbf{g}_{k})$ and $\mathbf{H}_{k}=\text{circ}(\mathbf{h}_{k})$ are both $p\times p$ circulant matrices and $\mathbf{S}\in\mathbb{R}^{p\times\frac{p}{2}}$ is a selection matrix that keeps only every other column, i.e., $\mathbf{G}_{k}\mathbf{S}=\begin{bmatrix}\mathbf{g}_{k}&\mathbf{P}^{2}\mathbf{g}_{k}&\mathbf{P}^{4}\mathbf{g}_{k}&\dots&\mathbf{P}^{n-2}\mathbf{g}_{k}\end{bmatrix}\in\mathbb{R}^{p\times\frac{p}{2}}$ (downsampling the columns by a factor of 2). In general we assume that the filters $\mathbf{g}_{k}$ and $\mathbf{h}_{k}$ have compact support with length denoted $n\leq p$. As such, these transformations are related to the convolutional dictionaries previously described. We call $\mathbf{g}_{k}$ and $\mathbf{h}_{k}$ filters because (3) is equivalent to a filtering operation of a signal $\mathbf{x}$ where the filter coefficients are stored in the circulant matrix.

Now we define a new transformation that operates only on the first $\frac{p}{2^{k-1}}$ coordinates and keeps the other unchanged:

[TABLE]

Finally, we define a wavelet-like synthesis transformation that is a cascade of the fundamental stages (5) as

[TABLE]

We call this transformation wavelet-like because $\mathbf{Wx}$ applies convolutions to parts of the signal $\mathbf{x}$ at different scales (for example, see [32] for a description of discrete wavelet transformations from the perspective of matrix linear algebra).

III The proposed dictionary learning algorithms

Dictionary learning [33] provides heuristics that approximate solutions to the following problem: given a dataset $\mathbf{Y}\in\mathbb{R}^{n\times N}$, a sparsity level $s$ and the size of the dictionary $S$ we want to create a dictionary $\mathbf{D}\in\mathbb{R}^{n\times S}$ and the sparse representations $\mathbf{X}\in\mathbb{R}^{S\times N}$ then solve

[TABLE]

A classic approach to this problem, which we also use, is the iterative alternating optimization algorithm: keep the dictionary fixed and update the representations and vice-versa in a loop until convergence.

In this paper, we constrain the dictionary $\mathbf{D}$ to the structures previously discussed: circulant and convolutional (including unions in both cases) and wavelet-like. Our goal is to propose numerically efficient dictionary update rules for these structures. While for the general dictionary learning problem there are several online algorithms that have been proposed [34, 35, 36] which are computationally efficient, in this paper, we consider only batch dictionary learning and we focus on the computational complexity of the dictionary update step.

III-A Circulant dictionary learning

Given a dataset $\mathbf{Y}\in\mathbb{R}^{n\times N}$ and a sparsity level $s\geq 1$ for the representations $\mathbf{X}\in\mathbb{R}^{n\times N}$, the work in [15] introduces an efficient way of learning a circulant dictionary $\mathbf{C}\in\mathbb{R}^{n\times n}$ by approximately solving the optimization problem:

[TABLE]

For fixed $\mathbf{X}$, to update $\mathbf{c}$ we develop the objective function to

[TABLE]

and in order to minimize it we set

[TABLE]

where $\mathbf{\tilde{y}}_{k}^{T}$ and $\mathbf{\tilde{x}}_{k}^{T}$ are the rows of $\mathbf{\tilde{Y}}=\mathbf{FY}$ and $\mathbf{\tilde{X}}=\mathbf{FX}$.

Remark 1. Given $\mathbf{Y}\in\mathbb{R}^{n\times N}$ and $\mathbf{X}\in\mathbb{R}^{n\times N}$ the best circulant dictionary $\mathbf{C}$ in terms of the Frobenius norm achieves

[TABLE]

Proof. Expand the objective of (8) using the optimal (10) as

[TABLE]

In the development of (12) we use the definition of the Frobenius norm, the invariance of the Frobenius norm under unitary transformations, i.e., $\|\mathbf{F}\mathbf{X}\|_{F}=\|\mathbf{F}^{H}\mathbf{X}\|_{F}=\|\mathbf{X}\|_{F}$, and the cyclic permutation of the trace, i.e., $\text{tr}(\mathbf{ABC})=\text{tr}(\mathbf{BCA})$.

By the Cauchy-Schwarz inequality we have that $|\mathbf{\tilde{x}}_{k}^{H}\mathbf{\tilde{y}}_{k}|^{2}\leq\|\mathbf{\tilde{x}}_{k}\|_{2}^{2}\|\mathbf{\tilde{y}}_{k}\|_{2}^{2}$ which holds with equality if and only if the rows $\mathbf{\tilde{x}}_{k}$ and $\mathbf{\tilde{y}}_{k}$ are multiples of each other. In this case, the objective function (12) develops to $\|\mathbf{Y}-\mathbf{CX}\|_{F}^{2}=\sum_{k=1}^{n}\|\mathbf{y}_{k}\|_{2}^{2}-\sum_{k=1}^{n}\|\mathbf{\tilde{y}}_{k}\|_{2}^{2}$. This is the only case where the circulant dictionary can reach zero representation error.

A necessary condition that the optimal circulant dictionary $\mathbf{C}$ obeys is $\|\mathbf{\sigma}\|_{2}^{2}=\sum_{k=1}^{n}\frac{|\mathbf{\tilde{x}}_{k}^{H}\mathbf{\tilde{y}}_{k}|^{2}}{\|\mathbf{\tilde{x}}_{k}\|_{2}^{4}}=1$, i.e., the $\ell_{2}$ norm of the optimal solution of (8) is one.$\hfill\blacksquare$

In the context of dictionary learning, to obey the unit $\ell_{2}$ norm constraint on $\mathbf{c}$ we should normalize the optimal solution $\mathbf{\sigma}\leftarrow\|\mathbf{\sigma}\|_{2}^{-1}\mathbf{\sigma}$. This is avoided because we can always group this normalizing factor with the representations $\mathbf{X}$ instead of the circulant dictionary, i.e., $\|\mathbf{\sigma}\|_{2}^{-1}\mathbf{CX}=\mathbf{C}(\|\mathbf{\sigma}|\|_{2}^{-1}\mathbf{X})$. This grouping is correct because the Fourier transform preserves $\ell_{2}$ norms and we have that $\|\mathbf{\sigma}\|_{2}=\|\mathbf{c}\|_{2}$, i.e., all the columns of $\mathbf{C}$ are $\ell_{2}$ normalized after normalizing $\mathbf{\sigma}$.

The algorithm called C–DLA, first introduced in [15], has low computational complexity that is dominated by the $O(nN\log n)$ computation of $\mathbf{\tilde{Y}}$ (once) and that of $\mathbf{\tilde{X}}$ (at each iteration). The calculations in (10) take approximately $2nN$ operations: there are $\frac{n}{2}$ components in $\mathbf{\sigma}$ to be computed while $\|\mathbf{\tilde{x}}_{k}\|_{2}^{2}$ and $\mathbf{\tilde{x}}^{H}\mathbf{\tilde{y}}_{k}$ take $2N$ operations each.

We can limit which and how many of all the possible $n$ shifts of $\mathbf{c}$ are allowed. We achieve this by ensuring that rows of $\mathbf{X}$ corresponding to unwanted shifts are zero.

Remark 2 (Approximating linear operators by circulant matrices). Given $\mathbf{Y}\in\mathbb{R}^{n\times n}$, let us consider the special case of (8) when $N=n$ and we fix $\mathbf{X}=\mathbf{I}$. Now we calculate the closest, in Frobenius norm, circulant matrix to a given linear transformation $\mathbf{Y}$. Because the Frobenius norm is elementwise the optimal solution is directly $c_{k}=\frac{1}{n}\sum_{(i-j)\!\!\!\mod\!n=(k-1)}y_{ij},\ k=1,\dots,n$. Unfortunately, in general, circulant matrices do not approximate all linear transformations with high accuracy. The result is intuitive since matrices have $n^{2}$ degrees of freedom while circulants have only $n$. Furthermore, if we add other constraints, such as orthogonality for example, the situation is even worse: [37] shows that the set of orthonormal circulants is finite and constructs it. Therefore, researchers proposed approximating a linear transformation by a product of $O(n)$ circulant, or Toeplitz, and diagonal matrices [9, 38].$\hfill\blacksquare$

III-B Union of circulant dictionary learning

Let us now consider overcomplete dictionaries (matrices that have, significantly, more columns than rows) that are unions of circulant matrices. In particular, take a dictionary which is the union of $L$ circulants:

[TABLE]

where each $\mathbf{C}^{(\ell)}=\text{circ}(\mathbf{c}^{(\ell)})=\mathbf{F}^{H}\mathbf{\Sigma}^{(\ell)}\mathbf{F},\ \mathbf{\Sigma}^{(\ell)}=\text{diag}(\mathbf{\sigma}^{(\ell)}),\ \mathbf{\sigma}^{(\ell)}=\mathbf{Fc}^{(\ell)}$, is a circulant matrix. Given training data $\mathbf{Y}\in\mathbb{R}^{n\times N}$, with this structure, the dictionary learning problem has the objective:

[TABLE]

where the tilde matrices indicate the Fourier transforms (taken columnwise) and the representations $\mathbf{X}\in\mathbb{R}^{nL\times N}$ are separated row-wise into $L$ continuous non-overlapping blocks of size $n$ denoted $\mathbf{X}^{(\ell)}\in\mathbb{R}^{n\times N}$. A way to update all circulant components using the Fourier transforms $\mathbf{\tilde{Y}}$ and $\mathbf{\tilde{X}}$ presents itself. Denote by $\mathbf{\tilde{y}}_{k}^{T}$ the $k^{\text{th}}$ row of $\mathbf{\tilde{Y}}$ and by $(\mathbf{\tilde{x}}_{k}^{(\ell)})^{T}$ the $k^{\text{th}}$ row of $\mathbf{\tilde{X}}^{(\ell)}$. The objective function of the dictionary learning problem separates into $k=1,\dots,\left\lceil\frac{n}{2}\right\rceil+1$ (given real valued training data $\mathbf{Y}$) distinct least squares problems like:

[TABLE]

Therefore, for a fixed $k$, the diagonal entries $(k,k)$ of all $\mathbf{\Sigma}^{(\ell)}$ (which are denoted $\sigma_{k}^{(\ell)}$) are updated simultaneously by solving the least squares problems (15). Given real-valued data, mirror relations $\sigma_{n-k+2}^{(\ell)}=(\sigma_{k}^{(\ell)})^{*},k=2,\dots,n,$ hold analogously to (10) for all $\ell=1,\dots,L$. Notice that this formulation is just a natural extension of the one dimensional least squares problems in (10). To compute all the components of all $\mathbf{\sigma}^{(L)}$ we solve this least squares problem $\frac{n}{2}$ times – the computational complexity is $O(nL^{2}N)$.

In comparison, the union of circulants dictionary learning method presented in [15], UC–DLA, updates each circulant block $\mathbf{C}^{(\ell)}$ sequentially and separately (this can be seen as a block coordinate descent approach). The new proposed learning method, called Union of Circulant Dictionary Learning Algorithm with Simultaneous Updates (UCirc–DLA–SU), is described in Algorithm 1.

**Remark 3 (Updating an unused circulant component). **Assuming that $\mathbf{X}^{(\ell)}=\mathbf{0}_{n\times N}$, i.e., the $\ell^{\text{th}}$ circulant matrix is never used in the representations, then we use the update $\mathbf{c}^{(\ell)}=\underset{\mathbf{z};\ \|\mathbf{z}\|_{2}^{2}=1}{\arg\max}\quad\left\|\left(\mathbf{Y}-\sum_{i=1,i\neq\ell}^{L}\mathbf{C}^{(i)}\mathbf{X}^{(i)}\right)\mathbf{z}\right\|_{2}^{2}.$This is the block update method use in UC–DLA [15]. Furthermore, similarly to [39], this update could be used also when atoms of block $\ell$ have a lower contribution to the reconstruction than atoms from other blocks, i.e., $\|\mathbf{X}^{(\ell)}\|_{F}^{2}\ll\|\mathbf{X}^{(i)}\|_{F}^{2},\forall\ i\neq\ell$.$\hfill\blacksquare$

III-C Union of convolutional dictionary learning

Analogously to (13), we define the dictionary which is a union of $L$ convolutional matrices as

[TABLE]

and the objective function to minimize with respect to this union dictionary given the fixed representations $\mathbf{X}_{\text{conv}}\in\mathbb{R}^{pL\times N}$ (separated into $L$ continuous non-overlapping blocks denoted $\mathbf{X}_{\text{conv}}^{(\ell)}=\begin{bmatrix}\mathbf{X}^{(\ell)}\\ \mathbf{0}_{(n-1)\times N}\end{bmatrix}\in\mathbb{R}^{p\times N}$) is developed as

[TABLE]

where $\mathbf{B}=\begin{bmatrix}\mathbf{A}^{(1)}\mathbf{F}_{:,1:n}&\dots&\mathbf{A}^{(L)}\mathbf{F}_{:,1:n}\end{bmatrix}\in\mathbb{R}^{pN\times nL}$, $\mathbf{A}^{(\ell)}=\begin{bmatrix}\mathbf{\tilde{x}}_{1}^{(\ell)}\otimes\mathbf{e}_{1}&\dots&\mathbf{\tilde{x}}_{p}^{(\ell)}\otimes\mathbf{e}_{p}\end{bmatrix}\in\mathbb{C}^{pN\times p},$ with the rows of $\mathbf{\tilde{X}}_{\text{conv}}^{(\ell)}$, $\mathbf{c}=\begin{bmatrix}\mathbf{c}^{(1)}&\mathbf{c}^{(2)}&\dots&\mathbf{c}^{(L)}\end{bmatrix}^{T}\in\mathbb{R}^{nL}$ and $\mathbf{F}_{:,1:n}\in\mathbb{C}^{p\times n}$ is the $p\times p$ Fourier matrix restricted to its first $n$ columns. The solution here is given by the least squares

[TABLE]

where $\mathbf{B}^{H}\mathbf{B}\in\mathbb{R}^{nL\times nL}$ is a positive definite block symmetric matrix, where each $n\times n$ block is a Toeplitz matrix like $\mathbf{T}_{\ell_{1}\ell_{2}}=\mathbf{F}_{:,1:n}^{H}(\mathbf{A}^{(\ell_{1})})^{H}\mathbf{A}^{(\ell_{2})}\mathbf{F}_{:,1:n}$ for the $(\ell_{1},\ell_{2})^{\text{th}}$ block – the diagonal blocks are symmetric positive definite Toeplitz. Therefore, $\mathbf{B}^{H}\mathbf{B}$ is determined by $nL+(2n-1)\frac{L(L-1)}{2}$ parameters – the first term covers the parameters of the $L$ symmetric Toeplitz diagonal blocks and the second terms covers the parameters of all $\frac{L(L-1)}{2}$ non-diagonal Toeplitz blocks. The computational burden is highest in order to calculate the diagonals $\mathbf{W}^{(\ell_{1},\ell_{2})}=(\mathbf{A}^{(\ell_{1})})^{H}\mathbf{A}^{(\ell_{2})}$ with entries $w_{k}^{(\ell_{1},\ell_{2})}=(\mathbf{\tilde{x}}_{k}^{(\ell_{1})})^{H}\mathbf{\tilde{x}}_{k}^{(\ell_{2})}$, i.e., the inner product of the corresponding rows from $\mathbf{\tilde{X}}_{\text{conv}}^{(\ell_{1})}$ and $\mathbf{\tilde{X}}_{\text{conv}}^{(\ell_{2})}$ respectively. These calculations take $O(pL^{2}N)$ operations. The inverse Fourier transforms of $\mathbf{W}^{(\ell_{1},\ell_{2})}$ to recover the entries of $\mathbf{T}_{\ell_{1}\ell_{2}}$ take only $O(p\log_{2}p)$ operations.

The inverse problem in (18) can be solved exactly in $O(n^{3}L^{3})$ via block Cholesky factorization [40, Chapter 4.2]. When $nL$ is large, an alternative approach is to use some iterative procedure like the Conjugate Gradient approach (already used in convolutional problems [41]) where the computational burden falls on computing matrix-vector products with the matrix $\mathbf{B}^{H}\mathbf{B}$ which take only $O(pL^{2})$ operations. The vector $\mathbf{v}=\mathbf{B}^{H}\mathbf{\tilde{y}}$ has the structure $\mathbf{v}=\begin{bmatrix}\mathbf{v}^{(1)}&\dots&\mathbf{v}^{(L)}\end{bmatrix}^{T}$ where $\mathbf{v}^{(\ell)}=\mathbf{F}_{:,1:n}^{H}\mathbf{z}^{(\ell)}\in\mathbb{R}^{n}$ and $\mathbf{z}^{(\ell)}\in\mathbb{C}^{p}$ with entries $z_{j}^{(\ell)}=(\mathbf{\tilde{x}}_{j}^{(\ell)})^{H}\mathbf{\tilde{y}}_{j}$, i.e., the $j^{\text{th}}$ entry of the $\ell^{\text{th}}$ component is the inner product of the corresponding rows from $\mathbf{\tilde{X}}_{\text{conv}}^{(\ell)}$ and $\mathbf{\tilde{Y}}$ respectively. The cost of computing $\mathbf{v}\in\mathbb{R}^{nL}$ is $O(pLN)$ since the inverse Fourier transforms are only $O(p\log_{2}p)$. Also, we need to compute once the Fourier transform of the dataset ($\mathbf{\tilde{Y}}$) and at each iteration all the $L$ Fourier transforms of the sparse representations $\mathbf{\tilde{X}}_{\text{conv}}^{(\ell)}$ which take $O(Np\log_{2}p)$ and $O(LNp\log_{2}p)$ overall operations respectively.

The new proposed learning method, called Union of Convolutional Dictionary Learning Algorithm with Simultaneous Updates (UConv–DLA–SU), is described in Algorithm 2.

The theoretical and algorithmic importance of dictionaries that are unions of convolutional matrices where the first columns have different (and potentially overlapping) supports has been highlighted in [13] (see in particular the clear presentation of the convolutional sparse model in [13, Figure 1]): group together the first columns of each $\mathbf{C}_{\text{conv}}^{(\ell)}$ into the local dictionary $\mathbf{D}^{(1)}$ of size $p\times L$, do the same with the second columns into the local dictionary $\mathbf{D}^{(2)}$ and so on until $\mathbf{D}^{(p)}$. These dictionaries are called local because they are localized to the reconstruction of only a few (depending on the size of the support $n$) grouped entries from a signal, as compared to a global signal model.

Notice that (17) makes use of the matrix $\mathbf{F}_{:,1:n}$ due to the sparsity pattern in $\mathbf{c}_{\text{conv}}$, i.e., only the first $n$ entries are non-zero. This can be generalized so that if the support of $\mathbf{c}_{\text{conv}}$ is denoted by $\mathcal{S}(\mathbf{c}_{\text{conv}})$ then the least squares problem (17) can be solved only on this support my making use of $\mathbf{F}_{:,\mathcal{S}(\mathbf{c}_{\text{conv}})}\in\mathbb{C}^{p\times|\mathcal{S}(\mathbf{c}_{\text{conv}})|}$, i.e., the Fourier matrix of size $p\times p$ restricted to the columns indexed in $\mathcal{S}(\mathbf{c}_{\text{conv}})$. Alternatively, if we do not want to (or cannot) decide the support a priori, we can add an $\ell_{1}$ penalty to (17) to promote sparse solutions.

Without the sparsity constraints in $\mathbf{c}_{\text{conv}}$, UConv–DLA–SU essentially reduces to Ucirc–DLA–SU but with a major numerical disadvantage: the decoupling that allows for (15) is no longer valid and therefore the least squares problem (18) is significantly larger and harder to solve. This extra computational effort seems unavoidable if we want to impose the sparsity in $\mathbf{c}_{\text{conv}}$. This motivates us to discuss the possibility of updating each convolutional dictionary sequentially, just like [15] does for circulant dictionaries.

Remark 4 (The special case of $L=1$ – the computational simplifications of constructing a single convolutional dictionary). Following a similar development as [15, Remark 2], consider again the objective function of (8) and develop

[TABLE]

where we have defined $\mathbf{A}=\begin{bmatrix}\mathbf{\tilde{x}}_{1}\otimes\mathbf{e}_{1}&\dots&\mathbf{\tilde{x}}_{p}\otimes\mathbf{e}_{p}\end{bmatrix}\in\mathbb{C}^{pN\times p}$, with $\{\mathbf{e}_{k}\}_{k=1}^{p}$ the standard basis for $\mathbb{R}^{p}$, i.e., $\mathbf{A}$ is composed of columns from $(\mathbf{\tilde{X}}_{\text{conv}}^{T}\otimes\mathbf{I})$ corresponding to the non-zero diagonal entries of $\mathbf{\Sigma}_{\text{conv}}$ (see also the Khatri-Rao product [31]).

By the construction in (3) only the first $n$ entries of $\mathbf{c}_{\text{conv}}$ are non-zero and therefore the optimal minimizer of (19) is

[TABLE]

The large matrix $\mathbf{AF}_{:,1:n}$ is never explicitly constructed in the computation of $\mathbf{c}=(\mathbf{F}_{:,1:n}^{H}\mathbf{A}^{H}\mathbf{A}\mathbf{F}_{:,1:n})^{-1}\mathbf{F}_{:,1:n}^{H}\mathbf{A}^{H}\mathbf{\tilde{y}}=(\mathbf{F}_{:,1:n}^{H}\mathbf{W}\mathbf{F}_{:,1:n})^{-1}\mathbf{v}$, where we have defined

[TABLE]

Observe that $\mathbf{T}=\mathbf{F}_{:,1:n}^{H}\mathbf{W}\mathbf{F}_{:,1:n}$ is a real-valued symmetric positive definite Toeplitz matrix – it is the upper left (the leading principal) submatrix of the $p\times p$ circulant matrix $\mathbf{F}^{H}\mathbf{W}\mathbf{F}$. The matrix $\mathbf{T}$ is never explicitly computed, but its first column is contained in the $n$ entries of the vector $\mathbf{F}^{H}\text{diag}(\mathbf{W})$. Also, notice that $\mathbf{A}^{H}\mathbf{\tilde{y}}$ is a vector whose $j^{\text{th}}$ entry is the inner product of the corresponding rows from $\mathbf{\tilde{X}}$ and $\mathbf{\tilde{Y}}$ respectively, i.e., $\mathbf{\tilde{x}}_{j}^{H}\mathbf{\tilde{y}}_{j}$.

The computation of $\mathbf{W}$ and $\mathbf{v}$ take approximately $O(pN)$ operations each – because $\mathbf{X}$ is real-valued, the Fourier transforms exhibit symmetries and only half the $\ell_{2}$ norms in $\mathbf{W}$ and of the entries in $\mathbf{A}^{H}\mathbf{\tilde{y}}$ need to be computed. These calculations dominate the computational complexity since they depend on the size of the dataset $N\gg p$ – together with the computations of the Fourier transforms $\mathbf{\tilde{X}}_{\text{conv}}$ (at each iteration) and $\mathbf{\tilde{Y}}$ (only once) which take $O(LNp\log_{2}p)$ and $O(Np\log_{2}p)$ operations respectively. The least squares problem (20) is solved via the Levinson-Durbin algorithm [40, Section 4.7.3], whose complexity is $4n^{2}$ instead of the regular $O(n^{3})$ computational complexity for unstructured inverse problems.

Also, observe that this least squares solution when applied to minimizing (9) leads to the same optimal solution from (10): $\mathbf{c}=(\mathbf{F}^{H}\mathbf{A}^{H}\mathbf{A}\mathbf{F})^{-1}\mathbf{F}^{H}\mathbf{A}^{H}\mathbf{\tilde{y}}=\mathbf{F}^{H}\mathbf{W}^{-1}\mathbf{F}\mathbf{F}^{H}\mathbf{A}^{H}\mathbf{\tilde{y}}=\mathbf{F}^{H}\mathbf{W}^{-1}\mathbf{A}^{H}\mathbf{\tilde{y}}=\mathbf{F}^{H}\mathbf{\sigma}$, with $\mathbf{W}=\text{diag}(\begin{bmatrix}\|\mathbf{\tilde{x}}_{1}\|_{2}^{2}&\dots&\|\mathbf{\tilde{x}}_{n}\|_{2}^{2}\end{bmatrix})$. The approach in (10) is preferred to (20) since in the Fourier domain the overall problem is decoupled into a series of smaller size independent subproblems that can be efficiently solved in parallel.

Therefore, each single convolutional dictionary can be updated efficiently and an algorithm in the style of UC–DLA [15] can be proposed whenever running time is of concern or the dataset is large. In this case, because the convolutional components would not be updated simultaneously, we would expect worse results on average.$\hfill\blacksquare$

There are several places where structures related to unions of circulant and convolutional dictionaries appear. We briefly discuss next Gabor frames and then, in the following section, wavelet-like dictionaries.

Remark 5 (A special case of time-frequency synthesis dictionary). Union of circulant matrices also show up when studying discrete time-frequency analysis/synthesis matrices [42]. Consider the Gabor synthesis matrix given by

[TABLE]

where $\mathbf{C}=\text{circ}(\mathbf{g})$ for $\mathbf{g}\in\mathbb{C}^{m}$ which is called the Gabor window function, i.e., the matrix $\mathbf{G}$ contains circular shifts and modulations of $\mathbf{g}$. The matrices $\mathbf{D}^{(\ell)},\ell=1,\dots,m,$ are diagonal with entries $d^{(\ell)}_{kk}=\omega^{(\ell-1)(k-1)}$ and $\omega=e^{2\pi i/m}$.

In the context of compressed sensing with structured random matrices, Gabor measurement matrices have been used for sparse signal recovery [43]: the Gabor function of length $m$ is chosen with random entries independently and uniformly distributed on the torus $\{z\in\mathbb{C}\ |\ |z|=1\}$ [43, Theorem 2.3]. Here, our goal is to learn the Gabor function $\mathbf{g}$ from a given dataset $\mathbf{Y}$ such that the dictionary $\mathbf{G}$ allows good sparse representations. Now the objective function develops to:

[TABLE]

where we have $\mathbf{A}\in\mathbb{C}^{mN\times m}$ and we have denoted $\mathbf{y}=\text{vec}(\mathbf{Y})$, $\mathbf{D}=\begin{bmatrix}\mathbf{D}^{(1)}&\dots&\mathbf{D}^{(m)}\end{bmatrix}$, $(\mathbf{\tilde{x}}^{(\ell)}_{k})^{T}$ is the $k^{\text{th}}$ row of $\mathbf{\tilde{X}}^{(\ell)}$ and $\mathbf{\tilde{d}}^{(\ell)}_{k}$ is the $k^{\text{th}}$ column of $\mathbf{D}^{(\ell)}\mathbf{F}^{H}$. Similarly to the convolutional dictionary learning case, Gabor atoms typically have compact support and we can add the sparse structure to $\mathbf{g}$ (the support of size $n\leq m$) and find the minimizer by solving the least squares problem which this time is unstructured (there is no Toeplitz structure).$\hfill\blacksquare$

III-D Wavelet-like dictionary learning

Starting from (6), in the spirit of dictionary learning, our goal is to learn a transformation from given data such that it has sparse representations. The strategy we use is to update each $\mathbf{W}_{k}^{(p)}$ (actually, the $\mathbf{C}_{k}^{(p)}$ component) while keeping all the other transformations fixed. Therefore, for the $k^{\text{th}}$ component we want to minimize

[TABLE]

where $\mathbf{F}$ is the Fourier matrix of size $\frac{p}{2^{k-1}}$, we denoted $\mathbf{W}_{\text{A}}=\mathbf{W}_{1}^{(p)}\cdots\mathbf{W}_{k-1}^{(p)}$, $\mathbf{W}_{\text{B}}=\mathbf{W}_{k+1}^{(p)}\cdots\mathbf{W}_{m}^{(p)}$ and $\mathbf{\bar{X}}=\mathbf{W}_{\text{B}}\mathbf{X}$, $\mathbf{\bar{Y}}=\mathbf{Y}-\mathbf{W}_{\text{A},2}\mathbf{\bar{X}}_{2}$, $\mathbf{\bar{y}}=\text{vec}(\mathbf{\bar{Y}})$, $\mathbf{W}_{\text{A},1}$ are the first $\frac{p}{2^{k-1}}$ columns of $\mathbf{W}_{\text{A}}$, $\mathbf{\bar{X}}_{1}$ are the first $\frac{p}{2^{k-1}}$ rows of $\mathbf{\bar{X}}$. We have also denote $\mathbf{\Sigma}_{\mathbf{g}_{k}}=\text{diag}(\mathbf{Fg}_{k})$ and similarly for $\mathbf{\Sigma}_{\mathbf{g}_{k}}$. The matrix $\mathbf{A}\in\mathbb{C}^{pN\times\frac{2p}{2^{k-1}}}$ is made up of a subset of the columns from $((\mathbf{I}_{2}\otimes\mathbf{FS})\mathbf{\bar{X}}_{1})^{T}\otimes\mathbf{W}_{\text{A},1}\mathbf{F}^{H}$ corresponding to the non-zero entries from $\text{vec}(\begin{bmatrix}\mathbf{\Sigma}_{\mathbf{g}_{k}}&\mathbf{\Sigma}_{\mathbf{h}_{k}}\end{bmatrix})$. Minimizing the quantity in (24) leads to a least squares problem where both $\mathbf{g}_{k}$ and $\mathbf{h}_{k}$ have a fixed non-zero support of known size $n$. It obeys $n\leq\frac{p}{2^{m-1}}$ such that the circulants for $\mathbf{C}_{m}^{(p)}$ can be constructed.

Similarly to the union of circulants cases described before, some computational benefits arise when minimizing (24), i.e., computing $(\mathbf{I}_{2}\otimes\mathbf{F}^{H})\mathbf{A}^{H}\mathbf{A}(\mathbf{I}_{2}\otimes\mathbf{F})\backslash(\mathbf{I}_{2}\otimes\mathbf{F}^{H})\mathbf{A}^{H}\mathbf{\bar{y}}$. For convenience we will denote $\mathbf{Q}=(((\mathbf{I}_{2}\!\otimes\!\mathbf{FS})\mathbf{\bar{X}}_{1})^{T}$ and $\mathbf{R}=\mathbf{W}_{\text{A},1}\mathbf{F}^{H}$, such that $\mathbf{A}=\mathbf{Q}\odot\mathbf{R}$. Notice that $\mathbf{A}^{H}\mathbf{A}=\frac{p}{2^{k-1}}\begin{bmatrix}\mathbf{D}^{(1)}&\mathbf{D}^{(2)*}\\ \mathbf{D}^{(2)}&\mathbf{D}^{(3)}\end{bmatrix}$, where the blocks are diagonal with entries $d_{ii}^{(1)}=\|\mathbf{q}_{i}\|_{2}^{2}\|\mathbf{r}_{i}\|_{2}^{2}$, $d_{ii}^{(3)}=\|\mathbf{q}_{\frac{p}{2^{k-1}}+i}\|_{2}^{2}\|\mathbf{r}_{i}\|_{2}^{2}$ and $d_{ii}^{(2)}=\mathbf{q}_{\frac{p}{2^{k-1}}+i}^{H}\mathbf{q}_{i}\|\mathbf{r}_{i}\|_{2}^{2}$ where $\mathbf{q}_{i}$ and $\mathbf{r}_{i}$ are columns of $\mathbf{Q}$ and $\mathbf{R}$, respectively, $i=1,\dots,\frac{p}{2^{k-1}}$. Therefore, $(\mathbf{I}_{2}\otimes\mathbf{F}^{H})\mathbf{A}^{H}\mathbf{A}(\mathbf{I}_{2}\otimes\mathbf{F})$ is a $2\times 2$ block matrix whose blocks are real-valued circulant matrices (and the diagonal blocks are also symmetric). Also, because $\mathbf{A}^{H}\mathbf{A}$ is symmetric positive definite it allows for a Cholesky factorization $\mathbf{L}\mathbf{L}^{T}$ where the matrix $\mathbf{L}$ has only the main diagonal and the secondary lower diagonal of size $\frac{p}{2^{k-1}}$ of non-zero values. A further computational benefit comes from when $n\ll p$ and we solve a least squared problem in $2n$ variables, as compared to $2p$, i.e., $\mathbf{A}(\mathbf{I}_{2}\otimes\mathbf{F})\begin{bmatrix}\mathbf{g}_{k}\\ \mathbf{h}_{k}\end{bmatrix}=\mathbf{A}(\mathbf{I}_{2}\otimes\mathbf{F}_{:,1:n})\begin{bmatrix}\mathbf{\bar{g}}_{k}\\ \mathbf{\bar{h}}_{k}\end{bmatrix}$, with both $\mathbf{\bar{g}}_{k},\mathbf{\bar{h}}_{k}\in\mathbb{R}^{n}$. Finally, notice that $(\mathbf{I}_{2}\otimes\mathbf{F}_{:,1:n}^{H})\mathbf{A}^{H}\mathbf{A}(\mathbf{I}_{2}\otimes\mathbf{F}_{:,1:n})$ has Toeplitz blocks, like in the case of UConv–DLA–SU.

The linear transformation (6) has two major advantages: i) the computational complexity of matrix-vector multiplications $\mathbf{Wx}$ with a fixed $\mathbf{x}\in\mathbb{R}^{p}$ is controlled by the number of stages $m$ and the length of the filters $n$, instead of the fixed matrix-vector multiplication complexity which is $O(p^{2})$ and ii) it allows for learning atoms that capture features from the data at different scales, i.e., atoms of different sparsity levels.

The new proposed procedure, called Wavelet-like Dictionary Learning Algorithm (W–DLA), is described in Algorithm 3.

Remark 6 (Extending and constraining the wavelet-like structure). Unlike wavelets that use the same filters $\mathbf{g}$ and $\mathbf{h}$ (known as the low and high pass filters, respectively) at each stage $k$ of the transformation, we learn different filters $\mathbf{g}_{k}$ and $\mathbf{h}_{k}$. Also, the support of the filters (and their size) can be decided dynamically at every stage and the downsampling can be replaced with a general column selection matrix $\mathbf{S}$. Note that if we allow full support then the least squares problem can be solved in the Fourier domain with the complex-valued unknowns $\mathbf{\tilde{g}}_{k}=\mathbf{Fg}_{k}$ and $\mathbf{\tilde{h}}_{k}=\mathbf{Fh}_{k}$. Furthermore, each stage in (5) applies only to the decomposition of the left-most (so-called low-frequency) components from the previous stage. In the spirit of optimal sub-band tree structuring (also known as wavelet packet decompositions) [44], we can propose a transformation $\mathbf{W}=\mathbf{C}_{1}^{(p)}\cdots\mathbf{C}_{m-1}^{(p)}\mathbf{C}_{m}^{(p)}$ factored into $m$ stages all of the form (4) (where again each stage has its own filters $\mathbf{g}_{k}$ and $\mathbf{h}_{k}$ which now can have support $n\leq p$) and use the same optimization procedure. Lastly, we could also propose a structured transformation similar to (4) but based on more than two filters, e.g., $\mathbf{C}_{k}^{(p)}=\begin{bmatrix}\mathbf{G}_{k}\mathbf{S}&\mathbf{H}_{k}\mathbf{S}&\mathbf{J}_{k}\mathbf{S}\end{bmatrix}\in\mathbb{R}^{p\times p}$, for a new selection (downsampling by 3) matrix $\mathbf{S}\in\mathbb{R}^{p\times\frac{p}{3}}$. Heuristics to choose $m$, $n$ (maybe at each stage $m$, i.e., having $n_{k}$), the location of the $n$ non-zero entries in each filter, the structure and sizes of $\mathbf{C}_{k}$ and $\mathbf{S}$ may be proposed to further improve the accuracy of the algorithm (or the trade-off between numerical efficiency and accuracy in terms of the representation error).

Orthonormal wavelets, i.e., in our case meaning that $\mathbf{W}$ and all $\mathbf{C}_{k}^{(p)}$ are orthonormal, are also extensively used in many application. With these orthogonality constraints the objective function develops into a simpler form than (24) as

[TABLE]

where we have now denoted $\mathbf{\bar{y}}=\text{vec}(\mathbf{F}(\mathbf{W}_{\text{A}}^{T}\mathbf{Y}-\mathbf{\bar{X}}_{2}))$ and of course the Kronecker products are never explicitly built. We minimize this quadratic objective with the additional orthogonality constraint $\mathbf{g}_{k}^{T}\mathbf{h}_{k}=0$ (in order to keep $\mathbf{C}_{k}^{(p)}$ orthogonal). Minimizing quadratic functions under quadratic equality constraints has been studied in the past and numerical procedures are available [45]. This constraint ensures orthogonality between the columns of $\mathbf{GS}$ and $\mathbf{HS}$. To ensure orthogonality among the columns of $\mathbf{GS}$ (and $\mathbf{HS}$, respectively) we can explicitly add symmetry constraints to the filter coefficients, e.g., if $\mathbf{g}_{k}$ has a non-zero support of size four we have $g_{k3}=g_{k1}$ and $g_{k4}=-g_{k2}$. Alternatively, both orthogonality constraints can be added leading to a quadratic optimization problem with two quadratic constraints [46].$\hfill\blacksquare$

IV Experimental results

We now discuss numerical results that show how well the proposed methods extract shift-invariant structures from data. Since we are dealing with the dictionary learning problem our goal is to build good representations of the datasets we consider. In our proposed algorithms, in the sparse recovery phase, we use the orthogonal matching pursuit (OMP) algorithm [47], but any other sparse approximation method could be chosen. Also, because the sparse approximation steps are numerically expensive (at least quadratic complexity in general and applied for all $N\gg n$ data points) and repetitive operations (done at each iterative step of the alternating optimization algorithms), OMP is chosen from practical considerations, as numerically efficient implementations are available (for example [48]).

Notice from the description of all the proposed algorithms that each dictionary update step necessarily decrease the objective function but, unfortunately, the overall algorithms may not be monotonically convergent to a local minimum since OMP is not guaranteed in general to always reduce the objective function. As such, in this section, we also provide some experimental results where we empirically observe the converge of the proposed methods.

For datasets with a strong DC component the learned circulant dictionary might be $\mathbf{C}\approx\frac{1}{\sqrt{n}}\mathbf{1}_{n\times n}$. Therefore, preprocessing the dataset $\mathbf{Y}$ by removing the mean component is necessary and we have $\sigma_{1}=0$ in (10) since $\mathbf{\tilde{y}}_{1}=\mathbf{0}_{N\times 1}$. This operation is assumed performed before the application of any algorithm developed in this paper.

IV-A Synthetic experiments

We create a synthetic dataset generated from a fixed number of atoms and their shifts and the proceed to measure how well we recover them from the dictionary learning perspective. The experimental setup follows: generate $N=2000$ signals of length $n=20$, that are linear combinations (with sparsity $s=4$) of $L=45$ randomly generated kernel columns which are allowed to have only $q=3$ circular shifts (out of the possible $n=20$), i.e., $\mathbf{Y}\in\mathbb{R}^{n\times N}$ where each columns is $\mathbf{y}_{i}=\sum_{\ell=1}^{L}\alpha_{i\ell}\mathbf{P}^{q_{i\ell}}\mathbf{c}_{\ell}+\mathbf{n}_{i}$ for $i=1,\dots,N$ with fixed $\|\mathbf{c}_{\ell}\|_{2}=1$, $\|\mathbf{\alpha}_{i}\|_{0}=s$ where $\alpha_{i\ell}\in[-10,10]$ and $q_{i\ell}\in\{0,\dots,q-1\}$ are randomly uniformly distributed and $\mathbf{n}_{i}$ is a random Gaussian vector representing noise.

First, using the synthetic dataset, we show in Figure 1 how the UCirc–DLA–SU outperforms previously proposed methods in the task of recovering the atoms used in creating the dataset. This shows the benefit (as compared to UC–DLA) of updating all the circulant components simultaneously with each step of the algorithm. We observe that UCirc–DLA–SU achieves lower error approximately $75\%$ of the time. The typical counter-example is one where UC–DLA converges slower (in more iterations) to a slightly lower representation error, i.e., sub-optimal block calculations ultimately lead to a better final result. This observation is not surprising since both heuristic methods only approximately solve the overall original dictionary learning problem (with unknowns both $\mathbf{C}$ and $\mathbf{X}$). To show this, with the same synthetic dataset for noise level $\text{SNR}=30\text{dB}$ in Figure 2 we calculate how many times on average each atom in all circulant components (from all the $L=45$) is used in the sparse representations. On average, UCirc–DLA–SU recovers the correct supports (in effect, the indices of the shifts used) more often than UC–DLA.

Figure 3 shows the learning times for the union of circulants algorithms, with blocks and with simultaneous updates, for a fixed number of $K=100$ iterations. For this test, we created a synthetic dataset $\mathbf{Y}$ of size $n=64$ with $N=8192$ data points and we vary the sparsity level $s$ and the number of bases in the union $L$. For $s\in\{8,12\}$, in our experiments, UCirc–DLA–SU is always faster than UC–DLA [15] and the speedup is on average $20\%$ while for $s=4$ the average speedup is only $10\%$ and for large $L$ there are cases where UCirc–DLA–SU is slower than UC–DLA. In principle, UCirc–DLA–SU should always be faster but in practice, the algorithm involves memory manipulations in order to build all the matrices for all the subproblems (15). This is the reason why the running time difference is not larger. Furthermore, for small $s$ the blocks used by UC–DLA are calculated rapidly (calculations are similar to the matrix operations in Remark 5) because $\mathbf{X}$ is very sparse and the overall running time is, therefore, lower.

IV-B Experiments on ECG data

Electrocardiography (ECG) signals [49] have many repetitive sub-structures that could be recovered by shift-invariant dictionary learning. Therefore, in this section, we use the proposed UConv–DLA–SU to find in ECG signals short (compact support) features that are repeated. We use the MIT-BIH arrhythmia database111https://www.physionet.org/physiobank/database/mitdb/ from which we extract a normal sinus rhythm signal composed of equality length samples from five different patients, all sampled at 128 Hz. This signal is reshaped into a matrix of centered, non-overlapping sections of length $p=64$, leading to the training dataset $\mathbf{Y}\in\mathbb{R}^{64\times 101000}$.

Because we are searching for sub-signals with limited support, we use the UConv–DLA–SU with parameters $n=12$, $s=2$ and $L=2$ to recover the shift-invariant structure. In Figure 4 we show the original ECG signal and its reconstruction in the union of convolutional dictionaries. Of course, the reconstruction is not perfect but it is able to accurately capture the spikes in the data and remove some of the high-frequency features, i.e., the signal looks filtered (denoised). The second plot, Figure 5, shows the $L=2$ learned atoms from the data which capture the spiky nature of the training signal.

IV-C Experiments on image data

The training data $\mathbf{Y}$ that we consider are taken from popular test images from the image processing literature (pirate, peppers, boat, etc.). The test dataset $\mathbf{Y}\in\mathbb{R}^{p\times N}$ consists of $8\times 8$ non-overlapping image patches with their means removed. We consider $N=12288$ and we have $p=64$. To evaluate the learning algorithms, in this section we consider the relative representation error of the dataset $\mathbf{Y}$ in the dictionary $\mathbf{D}$ given the sparse representations $\mathbf{X}$ as

[TABLE]

We consider image data because there are well-known wavelet transforms that efficiently encode such data. We will use the filters of the Haar and Daubechies D4 wavelet transforms, i.e., with $n=2$, $m=\log_{2}p$, $\mathbf{\bar{h}}_{k}=\begin{bmatrix}\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}\end{bmatrix}$, $\mathbf{\bar{g}}_{k}=\begin{bmatrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\end{bmatrix}$ and $n=4$, $m=-1+\log_{2}p$, $\mathbf{\bar{h}}_{k}=\begin{bmatrix}\frac{1-\sqrt{3}}{4\sqrt{2}}&-\frac{3-\sqrt{3}}{4\sqrt{2}}&\frac{3+\sqrt{3}}{4\sqrt{2}}&-\frac{1+\sqrt{3}}{4\sqrt{2}}\end{bmatrix}$, $\mathbf{\bar{g}}_{k}=\begin{bmatrix}\frac{1+\sqrt{3}}{4\sqrt{2}}&\frac{3+\sqrt{3}}{4\sqrt{2}}&\frac{3-\sqrt{3}}{4\sqrt{2}}&\frac{1-\sqrt{3}}{4\sqrt{2}}\end{bmatrix}$, respectively, for all $k$. The filters are chosen such that the resulting $\mathbf{W}$ is orthonormal.

First, we show in Figure 6 the experimental convergence of the proposed W–DLA. If we also impose an orthogonality constraint on $\mathbf{W}$ then we can avoid OMP for the sparse representations and just a projection operation guarantees optimal sparse representations. The figure shows that, when available, it is convenient to initialize the W–DLA with well-known wavelet filters since the algorithm converges faster and to slightly lower representation errors. Still, the differences are not significant and wavelets do not exist for every $n,m$.

Then, in Figure 7 we show how the representation error varies with the sparsity level $s$. We also run W–DLA initialized with the well-known wavelet filter coefficients Haar and D4, respectively. W–DLA is always able to improve the representation performance, even when starting with the wavelet coefficients. In the Haar case the improvement is small due to the small number of free filter parameters to learn, i.e., only 24: $m=\log_{2}p=6$ stages each with 2 filters and each with 2 coefficients. In the D4 case, the representation error is improved significantly. Both figures show that, when available, wavelet coefficients provide an excellent initialization even slightly better than the proposed W–DLA (also confirmed in Figure 6). Since these wavelets are not available for all choices $n,m$ the purpose of this plot is to show that the proposed initialization provides very good results in general.

Finally, in Figure 8 we show the effect that parameters $n$ and $m$ have on the representation error. For reference, we show the representation error of C–DLA which has $p=64$ free parameters to learn. The performance of this dictionary is approximately matched by W–DLA with $n=4$ and $m=-1+\log_{2}p$ which has $40$ free parameters to learn: $m=5$ stages each with 2 filters of support 4 each. Notice that the representation error plateaus after $n=8$. In general, dictionaries built with W–DLA have $2nm$ degrees of freedom. We also show a version of W–DLA where we keep $m=1$ and vary only $m$ in which case, of course, the representation error decreases. Note that in this case each run of W–DLA is initialized with a random set of coefficients. To show monotonic convergence it would help to initialize the filters of size $n$ with those previously computed of support $n-1$.

V Conclusions

In this paper, we propose several algorithms that learn, under different constraints, shift-invariant structures from data. Our work is based on using circulant matrices and finding numerically efficient closed-form solutions to the dictionary update steps, by least-squares. We analyze the behavior of the algorithms on various data sources, we compare and show we outperform previously proposed algorithms from the literature.