Multivariate Convolutional Sparse Coding with Low Rank Tensor

TL;DR

This paper proposes a novel multivariate convolutional sparse coding model using tensor algebra and low-rank constraints, improving efficiency and performance in high-dimensional signal encoding.

Contribution

It introduces a tensor-based convolutional sparse coding model with low-rank and sparsity constraints, along with an efficient optimization algorithm and theoretical guarantees.

Findings

Enhanced encoding efficiency in high-dimensional settings

Improved performance over existing methods

Theoretical links to Kruskal tensor regression

Abstract

This paper introduces a new multivariate convolutional sparse coding based on tensor algebra with a general model enforcing both element-wise sparsity and low-rankness of the activations tensors. By using the CP decomposition, this model achieves a significantly more efficient encoding of the multivariate signal-particularly in the high order/ dimension setting-resulting in better performance. We prove that our model is closely related to the Kruskal tensor regression problem, offering interesting theoretical guarantees to our setting. Furthermore, we provide an efficient optimization algorithm based on alternating optimization to solve this model. Finally, we evaluate our algorithm with a large range of experiments, highlighting its advantages and limitations.

Equations

$$\mathcal{X}=\sum _{r=1}^R{\mathbf{x}}_r^{(1)}\circ \cdots \circ {\mathbf{x}}_r^{(p)},$$

$$[[{\mathbf{X}}_1,\cdots ,{\mathbf{X}}_p]]\triangleq \sum _{r=1}^R{\mathbf{x}}_r^{(1)}\circ \cdots \circ {\mathbf{x}}_r^{(p)}.$$

$$\mathcal{Y}=\sum _{k=1}^K{\mathcal{D}}_k{\star }_{1,\cdots ,p}{\mathcal{Z}}_k+\mathcal{E},$$

$$\text{$\mathcal{C}$irc}(\boldsymbol{\mathcal{D}})(\ell_{1},k_{1},\cdots,\ell_{p},k_{p})=\left\{\begin{array}[]{ll}&0\quad\mbox{if $\exists 1\leq i\leq p$ s.t. $\ell_{i}<k_{i}$ or $\ell_{i}\geq k_{i}+w_{i}$ }\\ &\boldsymbol{\mathcal{D}}(\ell_{1}-k_{1},\cdots,\ell_{p}-k_{p})\quad\mbox{otherwise.}\end{array}\right.$$

$$\mathcal{Y}=\mathcal{A}(\mathcal{Z})+\mathcal{E},$$

$$\begin{array}[]{lclc}\boldsymbol{\mathcal{A}}:&\mathbb{R}^{K}\times\mathbb{Z}&\xrightarrow{\hskip 19.91684pt}&\mathbb{Y}\\ &\boldsymbol{\mathcal{Z}}&\xmapsto{\hskip 19.91684pt}&\Big{(}\sum_{k=1}^{K}\big{\langle}\tilde{\boldsymbol{\mathcal{D}}}_{k;i_{1},\cdots,i_{p}},\boldsymbol{\mathcal{Z}}_{k}\big{\rangle}_{F}\Big{)}^{n_{1},\ldots,n_{p}}_{i_{1}=1,\ldots,i_{p}=1},\end{array}$$

$$\displaystyle\Big{(}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}\boldsymbol{z}^{(1)}_{r}\circ\cdots\circ\boldsymbol{z}^{(p)}_{r}\Big{)}_{i_{1},\cdots,i_{p}}$$

$$\displaystyle=\Big{(}\tilde{\boldsymbol{\mathcal{D}}}_{k;i_{1},\cdots,i_{p}}\bigtimes^{p}_{m=1}\boldsymbol{z}^{(m)}_{r}\Big{)},$$

$$\operatorname*{arg\,min}_{\big{(}\boldsymbol{\mathcal{D}}_{1},\cdots,\boldsymbol{\mathcal{D}}_{K},\boldsymbol{Z}_{1,1},\cdots,\boldsymbol{Z}_{K,p}\big{)}\in\mathcal{S}}\Bigg{\lVert}\boldsymbol{\mathcal{Y}}-\sum_{k=1}^{K}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}[\![\boldsymbol{Z}_{k,1},\cdots,\boldsymbol{Z}_{k,p}]\!]\Bigg{\rVert}_{F}^{2}+\sum_{k,\ell}\alpha_{\ell}\lVert\boldsymbol{Z}_{k,\ell}\rVert_{1},$$

$$\{\begin{array}{rl}& \forall k,{\mathcal{D}}_k\in \mathbb{D},\Vert {\mathcal{D}}_k{\Vert }_F\le 1\\ & \forall k,\ell ,{\mathbf{Z}}_{k,\ell }\in {\mathbb{R}}^{n_{\ell }\times R}\text{ and }\Vert {\left({\mathbf{Z}}_{k,\ell }\right)}_{i,\cdot }{\Vert }_F=1\end{array}$$

$$\operatorname*{arg\,min}_{\boldsymbol{Z}_{1,\ell},\cdots,\boldsymbol{Z}_{K,\ell}}\Bigg{\lVert}\boldsymbol{\mathcal{Y}}-\sum_{k=1}^{K}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}[\![\boldsymbol{Z}_{k,1},\cdots,\boldsymbol{Z}_{k,p}]\!]\Bigg{\rVert}_{F}^{2}+\alpha_{\ell}\sum_{k=1}^{K}\lVert\boldsymbol{Z}_{k,\ell}\rVert_{1}+\beta_{\ell}\sum_{k=1}^{K}\lVert\boldsymbol{Z}_{k,\ell}\rVert^{2}_{F}.$$

$$\Big{\lVert}\boldsymbol{\mathcal{Y}}-\sum_{k=1}^{K}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}[\![\boldsymbol{Z}_{k,1},\cdots,\boldsymbol{Z}_{k,p}]\!]\Big{\rVert}_{F}^{2}=\sum^{C}_{c}\Big{\lVert}\tilde{\boldsymbol{\mathcal{Y}}}_{:,c}-\sum^{S}_{s}\tilde{\boldsymbol{\mathcal{D}}}_{s,:,c}\star\boldsymbol{z}_{s}^{(\ell)}\Big{\rVert}^{2}_{2},$$

$$\displaystyle\Big{\lVert}\boldsymbol{\mathcal{Y}}-\sum_{k=1}^{K}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}\boldsymbol{\mathcal{Z}}_{k}\Big{\rVert}_{F}^{2}=\sum_{i_{1}=1,\cdots,i_{p}=1}^{n_{1},\cdots,n_{p}}\Big{(}\boldsymbol{\mathcal{Y}}_{i_{1},\cdots,i_{p}}-\sum_{k=1}^{K}\sum_{r=1}^{R}\sum_{j_{1}=1}^{w_{1}}\bar{z}_{k,r}^{(1)}(i_{1}-j_{1})$$

$$\displaystyle\hskip 150.00023pt\sum_{j_{2}=1,\cdots,j_{p}=1}^{w_{2},\cdots,w_{p}}\boldsymbol{\mathcal{D}}_{k,j_{1},\cdots,j_{p}}\bar{z}_{k,r}^{(2)}(i_{2}-j_{2})\cdots\bar{z}_{k,r}^{(p)}(i_{p}-j_{p})\Big{)}^{2}$$

$$=$$

$$\operatorname*{arg\,min}_{\boldsymbol{z}_{s}^{(\ell)}}\sum^{C}_{c}\Big{\lVert}\tilde{\boldsymbol{\mathcal{Y}}}_{:,c}-\sum^{S}_{s}\tilde{\boldsymbol{\mathcal{D}}}_{s,:,c}\star\boldsymbol{z}_{s}^{(\ell)}\Big{\rVert}^{2}_{2}+\alpha_{\ell}\sum^{S}_{s}\lVert\boldsymbol{z}_{s}^{(\ell)}\rVert_{1}+\beta_{\ell}\sum^{S}_{s}\lVert\boldsymbol{z}_{s}^{(\ell)}\rVert^{2}_{2}.$$

$$\operatorname*{arg\,min}_{{\forall k,\boldsymbol{\mathcal{D}}_{k}\in\mathbb{D},\quad\lVert\boldsymbol{\mathcal{D}}_{k}\rVert_{F}\leq 1}}\Big{\lVert}\boldsymbol{\mathcal{Y}}-\sum^{K}_{k=1}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}[\![\boldsymbol{Z}_{k,1},\cdots,\boldsymbol{Z}_{k,p}]\!]\Big{\rVert}^{2}_{F}.$$

$$\Big{\lVert}\boldsymbol{\mathcal{Y}}-\sum_{k=1}^{K}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}[\![\boldsymbol{Z}_{k,1},\cdots,\boldsymbol{Z}_{k,p}]\!]\Big{\rVert}_{F}^{2}=\sum^{C}_{c}\Big{\lVert}\tilde{\boldsymbol{\mathcal{Y}}}_{:,c}-\sum^{S}_{s}\tilde{\boldsymbol{\mathcal{D}}}_{s,:,c}\star\boldsymbol{z}_{s}^{(\ell)}\Big{\rVert}^{2}_{2},$$

$$\displaystyle\Big{\lVert}\boldsymbol{\mathcal{Y}}-\sum_{k=1}^{K}\boldsymbol{\mathcal{D}}_{k}\star_{1,\cdots,p}\boldsymbol{\mathcal{Z}}_{k}\Big{\rVert}_{F}^{2}$$

$$=$$

Full text

1 Introduction

In recent years, dictionary learning and convolutional sparse coding techniques (CSC) have been successfully applied in a wide range of topics, including image classification (Mairal et al., 2009; Huang and Aviyente, 2007), image restoration (Aharon et al., 2006), and signal processing (Mairal et al., 2010). The main idea behind these representations is to conjointly learn a dictionary containing the patterns observed in the signal, and sparse activations that encode the temporal or spatial locations where these patterns occur. Previous works have mainly focused on the study of univariate signals or images (Garcia-Cardona and Wohlberg, 2017), solving the dictionary learning problem using sparsity constraint on atoms activations.

In many applications ranging from videos to neuroimaging, data are multivariate and therefore better encoded by a tensor structure (Zhou et al., 2013; Cichocki et al., 2015). This has led in the recent years to an increased interest in adapting statistical learning methods to the tensor framework, in order to efficiently capture multilinear relationships (Su et al., 2012). In particular, the low rank setting has been the subject of many previous works in the past few years, as an efficient tool to exploit the structural information of the data, particularly in regression problems (Zhou et al., 2013; Rabusseau and Kadri, 2016; Li et al., 2017; He et al., 2018). However, traditional low rank approaches for regression problems are difficult in this setting, as estimating or finding the rank of a tensor with order at least 3 is significantly more complex than in the traditional matrix case (Håstad, 1990), and requires new algorithms (Li et al., 2017). More recently, (He et al., 2018) have proposed a new approach –inspired by the SVD decomposition – that combines both low rank and sparsity constraints in the multivariate regression setting.

In this article, we propose to extend the classical CSC model to multivariate data by introducing a new tensor CSC model that combines both low-rank and sparse activation constraints. This model that we call Kruskal Convolutional Sparse Coding (K-CSC) (or Low-Rank Multivariate Convolutional Sparse Coding) offers the advantages of 1) taking into account the underlying structure of the data and 2) using fewer activations to decompose the data, resulting in an improved summary (dictionary) and a better reconstruction of the original multivariate signal.

Our main contributions are as follows. First, we show that under mild assumptions, this new problem can be rewritten as a Kruskal tensor regression, from which we discuss interesting properties. Then, we provide an efficient algorithm, based on an alternating procedure, that solves the K-CSC problem and scales with the number of activation parameters. Finally, we evaluate our algorithm experimentally using a large range of simulated and real tensor data, illustrating the advantages and drawbacks of our model.

2 Low-Rank Multivariate Convolutional Sparse Coding

2.1 Preliminaries on tensors

We first introduce the notation used in this paper and briefly recall some element of tensor algebra (see (Kruskal, 1977; Kolda, 2006; Kolda and Bader, 2009; Sidiropoulos et al., 2017) for a more in-depth introduction).

Across the paper, we use calligraphy font for tensors ($\boldsymbol{\mathcal{X}}$) bold uppercase letters for matrices ($\boldsymbol{X}$) bold lowercase letters for vectors ($\boldsymbol{x}$) and lowercase letter for scalars ($x$). Let $\|\boldsymbol{\mathcal{X}}\|_{F}$ and $\|\boldsymbol{\mathcal{X}}\|_{1}$ respectively denotes the Frobenius and the $\ell_{1}$ norms, and let \big{<}\cdot\big{>}_{F} be the scalar product associated with $\|\cdot\|_{F}$. The symbol $\circ$ refers to the the outer product, $\otimes$ to the Kronecker product, $\odot$ to the Khatri-Rao product, and $\times_{m}$ denotes the $m$-product. The symbol $\star_{1,\cdots,p}$ refers to the multidimensional discrete convolutional operator between scalar valued functions where the subscript indices are the dimension involved (see Figure 1). When the signal is unidimensional, $\star_{1}$ reduces to $\star,$ the 1-D discrete convolutional operator.

Proposition 1**.**

(CP Decomposition.)

For any $\boldsymbol{\mathcal{X}}\in\mathbb{X}\triangleq\mathbb{R}^{n_{1}}\times\ldots\times\mathbb{R}^{n_{p}}$, $\exists R>0,$ and, $\boldsymbol{x}^{(i)}_{r}\in\mathbb{R}^{n_{i}}$, $1\leq i\leq p$, $1\leq r\leq R$, such that

[TABLE]

where $\forall i\geq 2,\lVert\boldsymbol{x}^{(i)}_{r}\rVert_{F}=1$. The smallest $R$ for which such decomposition exists is called the Canonical Polyadic rank of $\boldsymbol{\mathcal{X}}$ ($\text{CP-rank}(\boldsymbol{\mathcal{X}})$ or rank$(\boldsymbol{\mathcal{X}})$ for short), and in this case (1) is referred to as the CP decomposition of $\boldsymbol{\mathcal{X}}$.

Definition 1**.**

(Kruskal operator.) With the notation of Proposition 1,the Kruskal operator $[\![\quad\cdot\quad]\!]$ is defined as

[TABLE]

where $\boldsymbol{X}_{i}=\left[\boldsymbol{x}^{(i)}_{1},\ldots,\boldsymbol{x}^{(i)}_{R}\right]\in\mathbb{R}^{n_{i}\times R},$ $1\leq i\leq p$.

2.2 Model formulation

Let $\boldsymbol{\mathcal{Y}}\in\mathbb{Y}\triangleq\mathbb{R}^{n_{1}\times\cdots\times n_{p}}$ be a multidimensional signal and $\boldsymbol{\mathcal{D}}_{1},\cdots,\boldsymbol{\mathcal{D}}_{K}$ in $\mathbb{D}\triangleq\mathbb{R}^{w_{1}\times\cdots\times w_{p}}$ a collection of $K$ multidimensional atoms such that $\forall i,1\leq w_{i}\leq n_{i}$. The Kruskal Convolutional Sparse Coding model (K-CSC) is defined as

[TABLE]

where A) $\forall 1\leq k\leq K$, $\boldsymbol{\mathcal{Z}}_{k}\in\mathbb{Z}\triangleq\mathbb{R}^{m_{1}\times\cdots\times m_{p}}$ (with $m_{i}=n_{i}-w_{i}+1$) are sparse activation tensors with CP-rank lower than $R$, with $R$ small, and B) $\boldsymbol{\mathcal{E}}\in\mathbb{Y}$ is an additive (sub)gaussian noise, whose every component are independent and centered.

Advantages of low-rank tensor. The addition of a low rank constraint offers two main advantages. First, it reduces the number of unknown activation parameters from K\big{(}\prod_{i=1}^{p}m_{i}\big{)} (unconstrained model) to K\big{(}R\sum_{i=1}^{p}m_{i}\big{)}. For instance, in the regression case of typical RGB images of size $n$-by-$n$-by-$3$, the number of parameters decrease to $R\cdot(2n+3)$ instead of $3n^{2}$, resulting in a better scalability of the problem. Second, it exploits the structural information of $\boldsymbol{\mathcal{Y}}$, and has already been proved to be effective in various contexts (e.g. (Guo et al., 2012; Liu et al., 2013)). For example, previous works have shown that the vectorization of an image removes the inherent spatial structure of it while low rank tensor regression produces more interpretable results (Zhou et al., 2013).

Relation with separable convolution. The low-rank constraint imposes that each activation $\boldsymbol{\mathcal{Z}}_{k}$ has a $\text{CP-rank}(\boldsymbol{\mathcal{Z}}_{k})\leq R$ and writes as the sum of at most $R$ separable filters (product of multiple one dimensional filters). Problem (2) is therefore a separable convolution problem, which allows to use a FFT resolution and to significantly speed up the calculus of the convolution. As an example, the complexity of filtering an $n_{1}$-by-$n_{2}$ image with a $w_{1}$-by-$w_{2}$ non-separable filter is $\mathcal{O}(n_{1}n_{2}w_{1}w_{2})$ – instead of $\mathcal{O}(n_{1}n_{2}(w_{1}+w_{2}))$ for a separable filter.

Remark 1*.*

For simplicity, we assume in this paper that, $\forall k$, $\text{CP-Rank}(\boldsymbol{\mathcal{Z}}_{k})=R$. However, it is straightforward to extend \eqref{equ:convolve_p_rank} to a model where each tensor $\boldsymbol{\mathcal{Z}}_{k}$ have a different CP-rank (i.e. $\text{CP-Rank}(\boldsymbol{\mathcal{Z}}_{k})=R_{k}$).

2.3 Kruskal CSC as Tensor Regression

CSC models are often rewritten as a regression problem by using a circulant dictionary matrix (Garcia-Cardona and Wohlberg, 2018). In this section, we adopt the same strategy. We show that the K-CSC model (2) can be reformulated as an instance of the rank-$R$ Generalized Linear tensor regression Models introduced in (Zhou et al., 2013) – also known as Kruskal tensor regression model – (Zhou et al., 2013) and then use this new formulation to discuss the properties of the model. First we introduce the notion of circulant tensor, which is a generalization of the circulant matrix.

Definition 2**.**

(Circulant tensor.) Let $w_{1}<n_{1},\ldots,w_{p}<n_{p}\in\mathbb{N_{*}}$ and $\boldsymbol{\mathcal{D}}\in\mathbb{R}^{w_{1}\times\cdots\times w_{p}}$. We define the (quasi)-circular tensor generated by $\boldsymbol{\mathcal{D}}$, \text{\mathcal{C}irc}(\boldsymbol{\mathcal{D}})\in\mathbb{R}^{(n_{1}\times m_{1})\cdots\times(n_{p}\times m_{p})}, as follows

[TABLE]

In other word, \text{\mathcal{C}irc}(\boldsymbol{\mathcal{D}}) contains the atom $\boldsymbol{\mathcal{D}}$ which is translated in every directions. The following proposition shows the equivalence between K-CSC model and Kruskal regression. The consequence of this relation are studied in the next Section.

Proposition 2**.**

(K-CSC as Kruskal regression.) Let $\boldsymbol{\mathcal{D}}_{1},\cdots,\boldsymbol{\mathcal{D}}_{K}$ in $\mathbb{D}$ be a collection of $K$ multidimensional atoms. Then, (2) is equivalent to

[TABLE]

where $\boldsymbol{\mathcal{Z}}$ is in $\mathbb{R}^{K}\times\mathbb{Z}$ such that $\boldsymbol{\mathcal{Z}}_{k}\triangleq\boldsymbol{\mathcal{Z}}(k,\cdot,\cdots,\cdot)=\sum_{r=1}^{R}z_{k,r}^{(1)}\circ\cdots\circ z_{k,r}^{(p)}$ and $\boldsymbol{\mathcal{A}}$ is the linear map defined by

[TABLE]

where \forall k,\tilde{\boldsymbol{\mathcal{D}}}_{k;i_{1},\cdots,i_{p}}=\text{\mathcal{C}irc}(\boldsymbol{\mathcal{D}}_{k})(i_{1},\cdot,i_{2},\cdot,\cdots,i_{p},\cdot).

Proof.

The proposition is a consequence of the following equality

[TABLE]

where $\tilde{\boldsymbol{\mathcal{D}}}_{k;i_{1},\cdots,i_{p}}\bigtimes^{p}_{m=1}\boldsymbol{z}^{(m)}_{r}\triangleq\tilde{\boldsymbol{\mathcal{D}}}_{k;i_{1},\cdots,i_{p}}\times_{1}\boldsymbol{z}^{(1)}_{r}\times_{2}\boldsymbol{z}^{(2)}_{r}\cdots\times_{p}\boldsymbol{z}^{(p)}_{r}.$ ∎

3 Model estimation

In order to solve (2), we minimize the associated negative log-likelihood with an additional sparsity constraint – a common tool in CSC – on the activation tensors. This leads to the following minimization problem:

[TABLE]

where $\forall\ell$ $\alpha_{\ell}>0$ and \big{(}\boldsymbol{\mathcal{D}}_{1},\cdots,\boldsymbol{\mathcal{D}}_{K},\boldsymbol{Z}_{1,1},\cdots,\boldsymbol{Z}_{K,p}\big{)}\in\mathcal{S} i.f.f.

[TABLE]

Notice that the low rank constraint is embedded in $\mathcal{S}$ by the use of the Kruskal operator $[\![\quad\cdot\quad]\!]$ and the fact that $\forall k,\ell,\quad\boldsymbol{Z}_{k,\ell}\in\mathbb{R}^{n_{\ell}\times R}$. This use of the Kruskal operator can be seen as a generalization of the Burer-Monteiro heuristic for matrix (Burer and Monteiro, 2003).

Multiple CP Decompositions.

It should be noted that the CP decomposition is known to be unique when it satisfies the Kruskal condition (Kruskal, 1989), but only up to permutation of the normalized factor matrices. Therefore, the $\boldsymbol{Z}_{k,\ell}$ that solve (4) may not be unique, but they are isolated equivalent minimizers, and thus this problem does not negatively impact the optimization.

Regularizations.

In (4), the sparsity constraint enforces the sparsity of each element of the CP-decomposition for every activation tensors independently. Hence, the sparsity of each mode can be controlled. In addition, it is possible to add a specific ridge penalization (i.e. $\sum_{\ell=1}^{p}\beta_{\ell}\sum_{k=1}^{K}\lVert\boldsymbol{Z}_{k,\ell}\rVert^{2}_{F}$, $(\beta_{1},\cdots,\beta_{p})\succeq 0$) to (4) to limit numerical instabilities, decrease risks of overfitting and ensure the unity of the solution (Zhou et al., 2013).

Sparse Low Rank Regression.

A consequence of Proposition 2 is that (4) can be seen as an instance of a low-rank sparse regression problem. It has been shown that while both properties are desirable, a balance between the two constraints has to be found, as the two regularizations may have adversarial influence on the minimizer (Richard et al., 2012). This is achieved in our setting by using a Ivanov regularization for the rank (CP-rank $\leq R$) and a Tykhonov regularization for the sparsity ($\sum_{k,\ell}\alpha_{\ell}\lVert\boldsymbol{Z}_{k,\ell}\rVert_{1}$): the solution should be as sparse as possible while having a CP-rank lower than $R$.

Solving the optimization problem with AK-CSC

The non-convex problem (4)(due to the rank constraint) is convex with respect to each $Z$ block $([\boldsymbol{Z}_{k,i},\cdots,\boldsymbol{Z}_{k,i}])_{i},$ and $([\boldsymbol{\mathcal{D}}_{1},\cdots,\boldsymbol{\mathcal{D}}_{K}])$. Hence, we use a block coordinate strategy to minimize it. Our algorithm, called Alternated K-CSC (AK-CSC), splits the main non-convex problem into several convex subproblems; 1) by freezing $\boldsymbol{\mathcal{D}}$ and all except one $Z$ block at a time ($\boldsymbol{\mathcal{Z}}$-step) 2) by freezing only the activation tensor ($\boldsymbol{\mathcal{D}}$-step). Algorithm 1 presents this process.

Activations update, $\boldsymbol{\mathcal{Z}}$-step. In order to solve (4) with respect to $\boldsymbol{\mathcal{Z}}$, we proceed as follows: we assume that the dictionary $\boldsymbol{\mathcal{D}}$ is fixed and we iteratively solve the problem where all mode except the $\ell$-th one of each activation tensor are constant, for $\ell$ varying between $1$ and $p$. In other words, for each value of $\ell$, we solve the problem

[TABLE]

where the $\lVert\cdot\rVert^{2}_{F}$ is added to improve the minimization process, as previously discussed. Without any loss of generality, we set $\ell=1$ in the rest of this section, as the other values of $\ell$ can be treated similarly.

Proposition 3**.**

The first term of minimization problem (5) can be rewritten as

[TABLE]

with $C=\prod_{i=2}^{p}n_{i}$, $S=KR$ and \boldsymbol{z}_{s}^{(\ell)}=\big{[}z_{1,1}^{(\ell)},\cdots,z_{1,R}^{(\ell)},z_{2,1}^{(\ell)},\cdots,z_{K,R}^{(\ell)}\big{]}.

Proof.

In the following, for all $k$, we denote by $\bar{\boldsymbol{\mathcal{Z}}_{k}}=\sum_{r=1}^{R}\bar{z}_{k,r}^{(1)}\circ\cdots\circ\bar{z}_{k,r}^{(p)}\in\mathbb{Y}$ the tensor where we add [math] on each dimension to reach the one of $\boldsymbol{\mathcal{Z}}$.

[TABLE]

[TABLE]

where \tilde{\boldsymbol{\mathcal{D}}}_{k;r,j_{1},i_{2}\cdots,i_{p}}=\Big{(}\boldsymbol{\mathcal{D}}_{k;j_{1},:,\cdots,:}\star_{2,\cdots,p}z_{k,r}^{(2)}\circ\cdots\circ z_{k,r}^{(p)}\Big{)}_{i_{2},\cdots,i_{p}}. ∎

From the previous proposition, it is clear that in (5), each subproblem is a CSC with multichannel dictionary filters and single-channel activation maps (Wohlberg, 2016a), i.e. we need to solve

[TABLE]

Therefore, this Z block step can be solved using standard multi-channel CSC algorithms (see (Garcia-Cardona and Wohlberg, 2018) for a complete review).

Dictionary update, $\boldsymbol{\mathcal{D}}$-step. Given the $K$ activation tensors $(\boldsymbol{\mathcal{Z}}_{k})_{k}$, the dictionary update aims at improving how the model reconstructs $\boldsymbol{\mathcal{Y}}$ by solving

[TABLE]

This step presents no significant difference with existing methods. This problem is smooth and convex and can be solved using classical algorithms (Mairal et al., 2010; Yellin et al., 2017; Chalasani et al., 2013).

4 Related work

Related models. With specific choices on the parameters or on the dimension values, the K-CSC model reduces to well-known CSC problems. In the following, we enumerate some of them which also use the (multidimensional) convolutional product.

• •

Univariate CSC: For vector-valued atoms and signals ($p=1$), our model reduces to the $1$-D Convolutional Dictionary Learning model (CDL). The sparse coding step – i.e. the minimization on $(\boldsymbol{z}_{k})_{k}$ – is commonly referred as Convolutional Basis Pursuit DeNoising where the two leading approaches are based on the the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) (Beck and Teboulle, 2009) and on ADMM (Boyd et al., 2011). The dictionary update step – i.e. the minimization on $(\boldsymbol{d}_{k})_{k}$ – is referred as Method of Optimal Directions (with a constraint on the filter normalization) where, again, the most efficient solutions are based on FISTA and ADMM (Garcia-Cardona and Wohlberg, 2018).

• •

Multivariate CSC: If $p>1$ and $R=+\infty$ (i.e. no low-rank constraint), our model reduces to the multivariate CDL model also referred as multi-channel CDL. To the best of our knowledge, the only reference to multivariate CSC is (Wohlberg, 2016a), where the author proposes two models for $3$-channel images. More recently, (Garcia-Cardona and Wohlberg, 2018) propose scalable algorithms for hyper-spectral images (tensor of order $3$ with a high dimension on the third mode).

• •

Multivariate CSC with rank-1 constraint: If $p=2$, $R=1$ and $w_{2}=1$ (i.e. vector-valued atoms), our model reduces to the one recently presented in (La Tour et al., 2018). In this work, the authors strongly rely on the rank one constraint to solve the CSC problem and only consider spatio-temporal data.

Differences with recent tensor based approaches. Some previous works have proposed different approaches to CSC in a tensor setting, albeit without a low rank setting – a key component of our approach. In (Jiang et al., 2018), the authors introduce a new formulation based on a t-linear combination (related to the t-product). In (Bibi and Ghanem, 2017) they propose a generic CSC model based on tensor factorization strategy called t-SVD which also use the t-product. Notice that, while this new formulation reduces back to 2-D CSC when the tensor order is set to particular sizes, the notion of low-rank on the activations is not considered. Other tensor-based CSC models enforce the low-rank constraint on the dictionary instead of the activations using a Tucker or a CP decomposition (Zhang et al., 2017; Tan et al., 2015) or tensor factorization techniques (Huang and Anandkumar, 2015).

5 Experiments

In this section we evaluate our tensor-based dictionary learning framework K-CSC on both synthetic and real-world datasets. We compare our algorithm AK-CSC to state-of-the-art dictionary learning algorithm based on ADMM. Results are for tensor-structured output regression problems on which we report $\ell_{2}$-distance. All experiments are conducted on a linux laptop with $4$-core $2.5$GHz Intel CPUs using Tensorly (Kossaifi et al., 2019), Sporco (Wohlberg, 2017) and standard python libraries.

$\boldsymbol{\mathcal{Z}}$-step solver. As shown in the previous section, the main problem can be rewriten as a regression problem. Hence, to solve the $\boldsymbol{\mathcal{Z}}$-step, it is possible to use standard tensor regression solvers (Zhou et al., 2013; Li et al., 2017; He et al., 2018). However, it necessitates the construction of an enormous circulant tensor which is not tractable in practice due to memory limitation. In section 3, we show that the $\boldsymbol{\mathcal{Z}}$-step necessitates to solve $p$ multi-channel CSC problem with a large amount of channels ($C=\prod_{i=2}^{p}n_{i}$). While this problem has received only limited attention for many more than three channels, in our experiment, we used the algorithm proposed in (Garcia-Cardona and Wohlberg, 2018) which turns out to be scalable regarding the value of $C$.

Initialization. Each activation subproblem is regularized with a $\ell_{1}$-norm, which induces sparsity on a particular mode. As an example, consider the update of the $Z$ block $([\boldsymbol{Z}_{k,i},\cdots,\boldsymbol{Z}_{k,i}]),$ during the $\boldsymbol{\mathcal{Z}}$-step (see equation (5)). From Tibshirani et al. (2015), we know that there exists a value $\alpha^{\max}_{i}$ above which the subproblem solution is always zeros. As $\alpha^{\max}_{i}$ depends on the "atoms" $\tilde{\boldsymbol{\mathcal{D}}}_{s,:,c}$ and on the multichannel signal $\tilde{\boldsymbol{\mathcal{Y}}}_{:,c}$, its value changes after a complete loop. In particular, its value might change a lot between the initialization and the first loop. This is problematic since we cannot use a regularization $\alpha_{i}$ above this initial $\alpha^{\max}_{i}$. The standard strategy to initialize univariate CSC methods is to generate Gaussian white noise atoms. However, as these atoms generally poorly correlate with the signals, the initial value of $\alpha^{\max}_{i}$ are low compared to the following ones. To fix this problem, we use the strategy propose in (La Tour et al., 2018) and initialize the dictionary with random parts of the signal. Furthermore, in all the experiments, a zero-padding is added to the initial tensor $\boldsymbol{\mathcal{Y}}$.

Synthetic data.

To highlight the behavior of AK-CSC and the influence of the hyperparameters, we consider two different scenarios; 1) The rank is unknown 2) The rank is known. In both cases, the true dictionary is given and only the CSC task is evaluated. We generate $K=10$ random tensor atoms of size $\mathbb{R}^{2\times 4\times 8}$ where entries follow a Gaussian distribution with mean [math] and standard deviation in $[1,10]$. Each atom is associated to a sparse activation tensor of CP-rank $R^{*}=4$. The third order tensor $\boldsymbol{\mathcal{Y}}$ generated by model (2) is in $\mathbb{R}^{16\times 32\times 64}$. First, we illustrate the convergence of AK-CSC by plotting the loss (4) as a function of the rank parameter $R$ (Figure 2 A)). Convergence is reached after a reasonable number of full loops. It appears that choosing a rank $R$ greater than the true one permits a faster convergence. However, as depicted in Figure 2 B), an over estimation of the rank does not increase nor decrease the reconstruction error. Also note that, as expected, the error after convergence drastically decreases when we reach the true rank $R_{*}$. The research of the best set of hyperparameters for $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3})$ is done in $[0,100]^{3}$. Figure 2 C) illustrates how the three hyperparameters related to the sparsity influence the $\ell_{2}$-distance, i.e. the reconstruction of $\boldsymbol{\mathcal{Y}}$. We see that our method is robust to small modification of the hyperparameters. This is an important property which facilitates the search of the best set of hyperparameters.

Color Animated Pictures.

We consider a RGB-animated picture composed of $20$ images of size $30\times 30\times 3$ that make up the animation of Mario running. Hence, the animated picture is a $4$-th order tensor of size $\boldsymbol{\mathcal{Y}}\in\mathbb{R}^{30\times 30\times 3\times 20}$. The objective is to learn a maximum of $K=20$ RGB-animated atoms in $\mathbb{R}^{20\times 10\times 3\times 3}$ in order to reconstruct $\boldsymbol{\mathcal{Y}}$. We compare the reconstruction of our algorithm AK-CSC when $R=1$ to the ADMM approach without low-rank constraint, which is a classical and efficient CSC solver (see comparative section (Wohlberg, 2016b)). Figure 3 A) shows that for the same level of sparsity $\boldsymbol{\alpha}$ our method always uses fewer non-zero coefficients and yet provides a better reconstruction. Indeed, the rank constraint on the activations allows to choose more accurate non-zero coefficients. For instance, Figure 3B) shows that even with $2.5$ times less non-zero coefficients, AK-CSC provides a visual better reconstruction.

Functional Magnetic Resonance Imaging dataset.

In a last example, we analyze functional Magnetic Resonance Imaging data (fMRI) from the Autism Brain Imaging Data Exchange (ABIDE) (Di Martino et al., 2008). The fMRI data is a third-order tensor in $\mathbb{R}^{31\times 37\times 31}$. We learn $K=20$ atoms of size $\mathbb{R}^{10\times 10\times 10}$ and of rank $R=1$. As in the previous example, we compare the reconstruction of our algorithm AK-CSC to the ADMM approach. The resulting reconstruction is displayed on Figure 4 and the learned atoms are in the supplementary materials. It appears that for the same level of sparsity, the reconstruction performance is more efficient with AK-CSC. Furthermore, the learned atoms are more informative than those with the classical method.

6 Conclusion

In this paper, we introduced K-CSC, a new multivariate convolutional sparse coding model using tensor algebra that enforces both element-wise sparsity and low-rankness of the activations tensors. We provided an efficient algorithm based on alternate minimization called AK-CSC to solve this model, and we showed that K-CSC can be rewritten as a Kruskal Tensor regression problem. Finally, we showed that AK-CSC achieves good performances on both synthetic and real data.

7 Detailed proof

We detail the proof of the following proposition from the main paper.

Proposition 4**.**

The first term of minimization problem ($4$) can be rewritten as

[TABLE]

with $C=\prod_{i=2}^{p}n_{i}$, $S=KR$ and \boldsymbol{z}_{s}^{(\ell)}=\big{[}z_{1,1}^{(\ell)},\cdots,z_{1,R}^{(\ell)},z_{2,1}^{(\ell)},\cdots,z_{K,R}^{(\ell)}\big{]}.

Proof.

In the following, for all $k$, we denote by $\bar{\boldsymbol{\mathcal{Z}}_{k}}=\sum_{r=1}^{R}\bar{z}_{k,r}^{(1)}\circ\cdots\circ\bar{z}_{k,r}^{(p)}\in\mathbb{Y}$ the tensor where we add [math] on each dimension to reach the one of $\boldsymbol{\mathcal{Z}}$.

[TABLE]

where \tilde{\boldsymbol{\mathcal{D}}}_{k;r,j_{1},i_{2}\cdots,i_{p}}=\Big{(}\boldsymbol{\mathcal{D}}_{k;j_{1},:,\cdots,:}\star_{2,\cdots,p}z_{k,r}^{(2)}\circ\cdots\circ z_{k,r}^{(p)}\Big{)}_{i_{2},\cdots,i_{p}}. ∎

8 Additional results

8.1 Synthetic data

We provide the full heatmap obtained in the Synthetic data section on Figure 5. This Figure shows how the three hyperparameters related to the sparsity influence the $\ell_{2}$-distance, i.e. the reconstruction of $\boldsymbol{\mathcal{Y}}$. We see that our method is robust to small modification of the hyperparameters. This is an important property which facilitates the search of the best set of hyperparameters.

8.2 Functional Magnetic Resonance Imaging dataset

The atoms learned from the main article for functional Magnetic Resonance Imaging (fMRI) are exposed Figure LABEL:fig:fMRI_atoms. We can see that our method learns interesting atoms (on the left) while standard methods failed to use the full dictionary (on the right). Indeed, more than half of the atoms learned by ADMM remain noise.

8.3 Color Animated Picture

To highlight the behavior of AK-CSC with regard to its non-convexity, we performed $20$ trials of the same minimization program. Results of such an experiments is provided on Figure LABEL:fig:converge_mario. This illustration shows that the convergence is globally identical on each trials. Hence, our algorithm seems to be robust.

8.4 Black and White Animated Picture

In this section, we focus on a black and white version of the animated pictures of Mario. Interestingly, the number of used atoms and of activations is much more smaller for our method. Indeed, it appears that only important atoms are kept as for the other one, the resulting dictionary is redundant. Quantitatively, for a similar error of reconstruction, the number of activations is much smaller. Figure 8 shows one example of reconstruction with associated important atoms.