# Online Multilinear Dictionary Learning
for Sequential Compressive Sensing

Thiernithi Variddhisaï and Danilo Mandic T. Variddhisaï and D. Mandic are with
the Department of Electrical & Electronic Engineering, Imperial College
London, United Kingdom.
E-mails<EMAIL_ADDRESS>
###### Abstract

A method for online tensor dictionary learning is proposed. With the
assumption of separable dictionaries, tensor contraction is used to diminish a
$N$-way model of $\mathcal{O}\left(L^{N}\right)$ into a simple matrix equation
of $\mathcal{O}\left(NL^{2}\right)$ with a real-time capability. To avoid
numerical instability due to inversion of sparse matrix, a class of stochastic
gradient with memory is formulated via a least-square solution to guarantee
convergence and robustness. Both gradient descent with exact line search and
Newton’s method are discussed and realized. Extensions onto how to deal with
bad initialization and outliers are also explained in detail. Experiments on
two synthetic signals confirms an impressive performance of our proposed
method.

###### Index Terms:

Tensors, sparse Tucker decomposition, multilinear algebra, online dictionary
learning, stochastic gradient, least squares, quasi-Newton, forgetting factor,
sliding window.

##  1 Introduction

Dictionary learning is a type of representation leaning where the underlying
features are assumed to have sparse structure [1]. It has been an active topic
in many research fields i.e. compressive sensing, source separation,
classification and denosing [2]. Dictionary is an over-complete collection of
“atoms” which can represent a particular class of signals via their linear
combination. This dictionary can be either fixed or learned where the atoms of
former are selected from a pre-defined set of bases e.g. Fourier, Wavelet,
etc. while those of the latter are determined by the data at hand (called
traning samples), hereby the name “dictionary learning”.

###  1.1 Online 1D Dictionary Learning

Since the field of dictionary learning has been well investigated and there
exist a myriad of excellent tutorials [3, 4, 5, 6], only online methods will
be reviewed here. The very earliest one followed an LMS-inspired algorithm
with rank-$1$ stochastic gradient [7]. A more conventional gradient descent
could be employed to take into account all information propagated from the
past which was then rendered online by the block co-ordinate descent with the
most recent calculation as a warm start [8]. It is worth noting that [8] is
not a fully recursive algorithm as all past information, however already
computed, is re-calculated at every iteration to update the dictionary, not to
mention that the warm start is already explained by past calculations. This
redundancy was alleviated by a recursive least square (RLS) framework [9]
where all past calculations are stacked into a set of parameters at previous
iteration, and only the arriving data is material to the computation. Though
performing well in average cases, the RLS scheme was reported to not converge
when the coefficients are extremely sparse [10] because it is based on the
method of optimal direction which involves the unstable inversion of sparse
matrix [5]. There are more online dictionary learning routines in the
literature but most of them would simply be an extension to the abovementioned
algorithms.

###  1.2 Steps towards Tensor-Based Dictionary Learning

The increasing deployment and availability of sensing technologies have
resulted in an exponential deluge of data which is extremely highly correlated
(termed multirelational [11]), and preferable to be jointly analyzed in order
to unveil meaningful deep pattern hidden from a conventional matrix model.
Tensors have a rich century-long presence in diverse discipline from
psychometrics, chemometrics, quantum physics and, until recently, signal and
data analysis at the confluence of signal processing, machine learning and big
data [14]. For comprehensive tutorials on basic tensor decompositions and
algebra, we refer to [12, 13, 14]. The tensor-based or multilinear dictionary
learning can be considered as 1D dictionary learning with Kronecker structure
[15]. ecently, there has been many multilinear dictionary learning algorithms
proposed e.g. K-CPD [16], T-MOD [17] and HO-KSVD [18] but so far, to the best
of our knowledge, no online algorithm. It is with this spirit our intention to
generalize the existing online routines in 1D into MD cases.

##  2 Tensor-Based Learning of a Representation Dictionary

In this section, we present the general construction of online tensor-based
learning of a representation dictionary. we begin with a separability
assumption to construct a mode-wise cost function. Afterwards, a recursive
algorithm is derived from a type of descent direction called stochastic
gradient with memory a selection of past data are used to yield. Further, we
specify the choices of update scheme for our algorithm; they include steepest
descent and quasi-Newton method. Then, we give details on how to optimize and
speed up the algorithm even further by adding tuning parameters such as
forgetting factor, length of sliding window and sample-correcting weight.

###  2.1 Separable Dictionaries

Let $\mathcal{X}_{t}\in\mathbb{R}^{J_{1}\times J_{2}\times\cdots\times J_{N}}$
be an observation of a tensor sequence at time instant $t$. Its sparse
representation model can be written in the following Tucker form

$\mathcal{X}_{t}=\mathcal{S}_{t}\times_{1}\bm{\Psi}^{(1)}\times_{2}\bm{\Psi}^{(2)}\times_{3}\cdots\times_{N}\bm{\Psi}^{(N)}+\mathcal{W}_{t},$
(1)

where $\bm{\Psi}^{(n)}\in\textbf{\em C}^{(n)}\subset\mathbb{R}^{J_{n}\times
L_{n}},\,n\in\textbf{\em N}$, $\textbf{\em N}=\\{1,2,...,N\\}$, is the
mode-$n$ over-complete dictionary (i.e. $J_{n}<L_{n},\ \forall n\in\textbf{\em
N}$) and $\mathcal{S}_{t}\in\mathbb{R}^{L_{1}\times L_{2}\times\cdots\times
L_{N}}$ is $K$-sparse coefficient tensors associated with $\mathcal{X}_{t}$,
with $\mathcal{W}_{t}$ representing white Gausssian noise. In congruence with
the 1-D case, the number of the non-zero entries in $\mathcal{S}_{t}$, $K$,
should be fewer than the total dimension of the observation, that is,
$K<{\prod\limits_{n=1}^{N}}J_{n}$.

Assuming that the $K$-sparse coefficient tensors $\mathcal{S}_{t}$ is fixed
and $\left\\{\bm{\Psi}\right\\}=\left\\{\bm{\Psi}^{(n)},\,\ \forall
n\in\textbf{\em N}\right\\}$ is a set containing dictionaries of all modes,
let the objective function at time instance t,
$\mathcal{J}_{t}(\left\\{\bm{\Psi}\right\\})$, be defined as follows

$\mathcal{J}_{t}(\left\\{\bm{\Psi}\right\\})=\mathcal{J}_{t}={\sum\limits_{\tau=1}^{t}}\left\|\mathcal{E}_{\tau}\right\|_{F}^{2},$
(2)

where
$\mathcal{E}_{\tau}=\mathcal{X}_{\tau}-\mathcal{S}_{\tau}\times_{1}\bm{\Psi}^{(1)}\times_{2}\bm{\Psi}^{(2)}\cdots\times_{N}\bm{\Psi}^{(N)}$
The aim of our tensor OLD (T-OLD) is to find $\left\\{\bm{\Psi}\right\\}$
which minimizes $\mathcal{J}_{t}$; this can be expressed as

$\bm{\Psi}_{t}^{(1)},\bm{\Psi}_{t}^{(2)},\ldots,\bm{\Psi}_{t}^{(N)}=\underset{\left\\{\bm{\Psi}\right\\}}{\textrm{arg}\,\min}\
\mathcal{J}_{t}.$ (3)

Solving eq. 3 is a non-convex (in fact multilinear) optimization problem and
known to be NP-complete. One popular method is to assume that all mode-$n$
dictionaries are separable (i.e. each multivariate atom is a rank-1 tensor
made up of the outer product of atoms from each mode-n dictionaries) [19].
With this assumption held, we can solve for each mode-$n$ dictionary by fixing
the other modes. Without loss of generality, we consider the case of mode-$n$
dictionary to derive our algorithm. Let’s modify eq. 2 to yield

$\mathcal{J}_{t}(\bm{\Psi}^{(n)})={\sum\limits_{\tau=1}^{t}}\left\|\mathcal{X}_{\tau}-\tilde{\mathcal{S}}_{\tau}^{(n)}\times_{n}\bm{\Psi}^{(n)}\right\|_{F}^{2},$
(4)

$\mathcal{J}_{t}(\bm{\Psi}^{(n)})=tr\left(\frac{1}{2}\bm{\Psi}^{(n)}\left({\sum\limits_{\tau=1}^{t}}\tilde{\mathcal{S}}_{\tau}^{(n)}\times_{\textbf{\em
N}/\\{n\\}}\tilde{\mathcal{S}}_{\tau}^{(n)}\right){\bm{\Psi}^{(n)}}^{T}-\left({\sum\limits_{\tau=1}^{t}}\mathcal{X}_{\tau}\times_{\textbf{\em
N}/\\{n\\}}\tilde{\mathcal{S}}_{\tau}^{(n)}\right){\bm{\Psi}^{(n)}}^{T}\right)$
(5)

where
$\tilde{\mathcal{S}}_{\tau}^{(n)}=\mathcal{S}_{\tau}\times_{1}\bm{\Psi}^{(1)}\times_{2}\bm{\Psi}^{(2)}\cdots\times_{n-1}\bm{\Psi}^{(n-1)}\times_{n+1}\bm{\Psi}^{(n+1)}\cdots\times_{N}\bm{\Psi}^{(N)}$
By using properties of tensor contraction and Frobenius inner product, it is
not difficult to reformulate the right-hand side of eq. 4 as eq. 5 shown at
the bottom of the page where the notation $\times_{\textbf{\em N}/\\{n\\}}$
means that the contraction is performed on “all-except-$n^{th}$” modes.

###  2.2 Fast Recursion via Dual Update

In order to make the online learning possible, a recursive formula of eq. 5 is
required. So let’s define the arriving data $\mathbf{A}_{t}^{(n)}$ and
$\mathbf{B}_{t}^{(n)}$ as

$\mathbf{A}_{t}^{(n)}=\tilde{\mathcal{S}}_{t}^{(n)}\times_{\textbf{\em
N}/\\{n\\}}\tilde{\mathcal{S}}_{t}^{(n)},$ (6a)
$\mathbf{B}_{t}^{(n)}=\mathcal{X}_{t}\times_{\textbf{\em
N}/\\{n\\}}\tilde{\mathcal{S}}_{t}^{(n)}.$ (6b)

With eqs. 5 and 6, we arrive at the following recursion

$\mathcal{J}_{t}(\bm{\Psi}^{(n)})=tr\left(\frac{1}{2}\bm{\Psi}^{(n)}\mathbf{R}_{t}^{(n)}{\bm{\Psi}^{(n)}}^{T}-\mathbf{P}_{t}^{(n)}{\bm{\Psi}^{(n)}}^{T}\right),$
(7)

where

$\mathbf{R}_{t}^{(n)}=\mathbf{R}_{t-1}^{(n)}+\mathbf{A}_{t}^{(n)},$ (8a)
$\mathbf{P}_{t}^{(n)}=\mathbf{P}_{t-1}^{(n)}+\mathbf{B}_{t}^{(n)}.$ (8b)

Now that the cost function $\mathcal{J}_{t}(\bm{\Psi}^{(n)})$ is quadratic
with respect to $\bm{\Psi}^{(n)}$, standard descent direction schemes will
suffice to solve for $\bm{\Psi}^{(n)}$. If we define a priori gradient matrix
$\mathbf{G}_{t}^{(n)}$ to be
$\partial\mathcal{J}_{t}(\bm{\Psi}^{(n)})/\partial\bm{\Psi}^{(n)}|_{\bm{\Psi}^{(n)}=\bm{\Psi}_{t-1}^{(n)}}$,
we will have

$\mathbf{G}_{t}^{(n)}=\bm{\Psi}_{t-1}^{(n)}\mathbf{R}_{t}^{(n)}-\mathbf{P}_{t}^{(n)}.$
(9)

Consequently, we can express a mode-$n$ dictionary update equation as

$\bm{\Psi}_{t}^{(n)}=\bm{\Psi}_{t-1}^{(n)}+\alpha_{t}\mathbf{D}_{t}^{(n)},$
(10)

where the step size $\alpha_{t}>0$ and the descent direction
$\mathbf{D}_{t}^{(n)}$ are respectively any positive scalar and matrix which
make
$\mathcal{J}_{t}(\bm{\Psi}_{t}^{(n)})<\mathcal{J}_{t}(\bm{\Psi}_{t-1}^{(n)})$.
It is straightforward to show that $\mathbf{D}_{t}^{(n)}$ satisfies

$tr\left(\mathbf{D}_{t}^{(n)}{\mathbf{G}_{t}^{(n)}}^{T}\right)<0.$ (11)

Given $\mathbf{D}_{t}^{(n)}$, one legitimate way to find $\alpha_{t}$ is the
exact line search, which results in

$\alpha_{t}=\left.-{\left\langle\mathbf{D}_{t}^{(n)},\mathbf{G}_{t}^{(n)}\right\rangle}_{F}\middle/{\left\langle\mathbf{D}_{t}^{(n)},\mathbf{H}_{t}^{(n)}\right\rangle}_{F}\right.,$
(12)

where

$\mathbf{H}_{t}^{(n)}=\mathbf{D}_{t}^{(n)}\mathbf{R}_{t}^{(n)}$ (13)

and ${\langle\cdot,\cdot\rangle}_{F}$ denotes the Frobenius inner product. To
accelerate the algorithm, stochastic gradient matrix
$\hat{\mathbf{H}}_{t}^{(n)}$ and a posteriori gradient matrix
$\hat{\mathbf{G}}_{t}^{(n)}$ are defined respectively as

$\hat{\mathbf{H}}_{t}^{(n)}=\bm{\Psi}_{t-1}^{(n)}\mathbf{A}_{t}^{(n)}-\mathbf{B}_{t}^{(n)}$
(14)

and

$\hat{\mathbf{G}}_{t}^{(n)}=\bm{\Psi}_{t}^{(n)}\mathbf{R}_{t}^{(n)}-\mathbf{P}_{t}^{(n)}.$
(15)

With some manipulation of eqs. 9, 10, 12, 13, 14 and 15, a dual recursion
between $\mathbf{G}_{t}^{(n)}$ and $\hat{\mathbf{G}}_{t}^{(n)}$ are given by

$\mathbf{G}_{t}^{(n)}=\hat{\mathbf{G}}_{t-1}^{(n)}+\hat{\mathbf{H}}_{t}^{(n)},$
(16a)
$\hat{\mathbf{G}}_{t}^{(n)}=\mathbf{G}_{t}^{(n)}+\alpha_{t}\mathbf{H}_{t}^{(n)}.$
(16b)

###  2.3 Choice of Descent Direction

It is logical to begin with the simplest form of unconstrained optimization
techniques - the steepest descent, therefore, $\mathbf{D}_{t}^{(n)}$ is given
by

$[\mathbf{D}_{t}^{(n)}]_{SD}=-\mathbf{G}_{t}^{(n)},$ (17)

Clearly, eq. 17 satisfies eq. 11. Another possibility is to choose a second-
order or Newton’s direction. From the construction of the cost function in eq.
7 and gradient matrix in eq. 9, it is not difficult to see that the Newton-
type alternative of $\mathbf{D}_{t}^{(n)}$ will be

$[\mathbf{D}_{t}^{(n)}]_{QN}=-\mathbf{G}_{t}^{(n)}\mathbf{C}_{t-1}^{(n)},$
(18)

$\mathbf{C}_{t}^{(n)}=\left(\mathbf{I}-{\mathbf{F}_{t}^{(n)}}^{T}\mathbf{H}_{t}^{(n)}\right)\mathbf{C}_{t-1}^{(n)}\left(\mathbf{I}-{\mathbf{H}_{t}^{(n)}}^{T}\mathbf{F}_{t}^{(n)}\right)+{\mathbf{D}_{t}^{(n)}}^{T}\mathbf{F}_{t}^{(n)}$
(19)

where $\mathbf{C}_{t-1}^{(n)}$ is the approximate of
${\mathbf{R}_{t-1}^{(n)}}^{-1}$. A variety of quasi-Newton methods can be
utilized to update $\mathbf{C}_{t}^{(n)}$ For example, the BFGS method will
render the iterative correction of $\mathbf{C}_{t}^{(n)}$ given in eq. 19
where

$\mathbf{F}_{t}^{(n)}={\left(\mathbf{D}_{t}^{(n)}{\mathbf{H}_{t}^{(n)}}^{T}\right)}^{-1}\mathbf{D}_{t}^{(n)}.$
(20)

###  2.4 Constraint Space of Dictionary

Last but not least, eq. 10 does not guarantee
$\bm{\Psi}_{t}^{(n)}\subset\textbf{\em C}^{(n)},\ \forall t$, the projection
operator $\bm{\mathcal{P}}_{\textbf{\em C}^{(n)}}$ is used on
$\bm{\Psi}_{t}^{(n)}$ such that

$\bm{\Psi}_{t}^{(n)}\leftarrow\bm{\mathcal{P}}_{\textbf{\em
C}^{(n)}}\left(\bm{\Psi}_{t}^{(n)}\right)\subset\textbf{\em C}^{(n)}.$ (21)

For simplicity, this paper consider the case where the convex set $\textbf{\em
C}^{(n)}$ is a linear map which conserves a space spanned by the dictionary
atoms i.e. the column space. We can express such a map with a diagonal matrix
$\bm{\Pi}_{t}^{(n)}\in\mathbb{R}^{L_{n}\times L_{n}}$ and re-write eq. 21 as

$\bm{\Psi}_{t}^{(n)}\leftarrow\bm{\Psi}_{t}^{(n)}\bm{\Pi}_{t}^{(n)}.$ (22)

If, for example, $\textbf{\em C}^{(n)}$ is the space where the dictionary has
a unit norm, each diagonal element in $\bm{\Pi}_{t}^{(n)}$ will equal the
inverse of the 2-norm of the corresponding column of $\bm{\Psi}_{t}^{(n)}$.
Since the mapping in eq. 22 will require $\mathbf{R}_{t}^{(n)}$,
$\mathbf{P}_{t}^{(n)}$ and $\mathbf{G}_{t}^{(n)}$ to be adjusted accordingly,
complicating the algorithm and contributing to large complexity in each
iteration, and since it does not alter the dictionaries’ column space, it will
not affect the desired algebraic structure of the dictionaries, can be skipped
completely from the algorithm, and can be done only to yield the dictionaries
at time instant required.

###  2.5 Additional Settings

The recursion in eq. 8 takes all past values of eq. 6 with equal weight. Since
the sparse coefficients $\mathcal{S}_{t}$ at initial iterations would be far
off its actual values due to the randomly initialized dictionaries,
$\mathbf{R}_{t}^{(n)}$ and, thus, $\mathbf{G}_{t}^{(n)}$ will be biased.
Moreover in real-world setting, it is possible for the data at some time
instants to be overwhelmed by outliers which could propagate through the whole
sequence of algorithm.

A. Decimating Bad Dictionaries

Two techniques can be utilized to alleviate the coefficient mismatch resulting
from bad initial dictionaries, the forgetting factor $\lambda_{t}$ \-
forgetting bad past data - and sliding window of length $T$ \- retaining a few
$T$ good recent samples. Combining two techniques, the cost function in eq. 2
is given by

$\mathcal{J}_{t}={\sum\limits_{\tau=t-T+1}^{t-1}}{\prod\limits_{\tau_{0}=1}^{t-\tau}}\lambda_{t-\tau_{0}}\left\|\mathcal{E}_{\tau}\right\|_{F}^{2}+\left\|\mathcal{E}_{t}\right\|_{F}^{2}.$
(23)

Note that eq. 23 is for $t>T$. For $t\leq T$, the starting sum index $\tau$
will equal $1$, instead of $t-T+1$.

As the iterations grow, the dictionaries would get closer to their real
values. It is hence sensible to begin with a small value of $\lambda_{t}$ \-
to quickly forget the not so good dictionaries - and then gradually increase
$\lambda_{t}$ towards $1$ \- to remember more the more accurate dictionaries.
There are several adaptive schemes for $\lambda_{t}$ e.g. linear, quadratic,
cubic, exponential and hyperbola, etc. Here, we resort to a quartic function
which sets $\lambda_{t}$ to increase from $\lambda_{0}$ to $1$ in $\tau$
iterations

$\lambda_{t}=\begin{cases}1-(1-\lambda_{0})(1-t/\tau)^{4}&\textnormal{if}\
i<\tau\\\ 1&\textnormal{if}\ i\geq\tau\end{cases}$ (24)

To select a proper size of sliding window ($T$), its effect on
$\mathbf{R_{t}^{(n)}}$ and $\mathbf{P_{t}^{(n)}}$ in the cost function eq. 7
must be clearly understood. In the rank-$1$ stochastic gradient method [7],
only the most recent samples are used to estimate $\mathbf{R_{t}^{(n)}}$ and
$\mathbf{P_{t}^{(n)}}$, thereby the term ”rank-$1$”. In [8],
$\mathbf{R_{t}^{(n)}}$ is approximated by all past information, making it a
better estimate. As explicated above, the data at early iterations would
probably contain and even with forgetting factor will take some time for those
erroneous to vanish completely. Since the maximal rank the matrix
$\mathbf{R_{t}^{(n)}}$ can attain is $L_{n}$ and if the rank of
$\mathbf{A_{t}^{(n)}}$ is at least $1$ for all $t$, then only $L_{n}$ samples
or $T\geq L_{n}$ are required to construct a full rank $\mathbf{R_{t}^{(n)}}$.

B. Correcting The Arriving Sample

Even though the dictionaries get closer to their real values after some
iterations, the samples may not be consistent with the dictionaries in the
presence of outliers which could make $\left\|\mathcal{E}_{t}\right\|_{F}$,
the error of the arriving data, undesirably large. The influence of the
arriving data can be curbed by tuning it with some weight $\mu_{t}$, which
should be disproportionate to its consistency with the existing dictionaries.
Extending eq. 23 yields

$\mathcal{J}_{t}={\sum\limits_{\tau=t-T+1}^{t-1}}{\prod\limits_{\tau_{0}=1}^{t-\tau}}\lambda_{t-\tau_{0}}\mu_{\tau}\left\|\mathcal{E}_{\tau}\right\|_{F}^{2}+\mu_{t}\left\|\mathcal{E}_{t}\right\|_{F}^{2}.$
(25)

Such disproportionate weight can be given by

$\mu_{t}=\left(1+\frac{\left\|\hat{\mathbf{H}}_{t}^{(n)}\right\|_{F}}{\sqrt{tr\left(\mathbf{A}_{t}^{(n)}\right)}}\right)^{-1}$
(26)

From eq. 25, in order to keep the form of the mode-$n$ cost function in eq. 7
unchanged, the recursion in eq. 8 will become

$\mathbf{R}_{t}^{(n)}=\lambda_{t-1}\mathbf{R}_{t-1}^{(n)}+\mathbf{S}_{t}^{(n)}$
(27a)
$\mathbf{P}_{t}^{(n)}=\lambda_{t-1}\mathbf{P}_{t-1}^{(n)}+\mathbf{Q}_{t}^{(n)}$
(27b)

where

$\mathbf{S}_{t}^{(n)}=\left\\{\begin{array}[]{cl}\mu_{t}\mathbf{A}_{t}^{(n)}-\mathit{\Lambda}_{t-T}\mu_{t-T}\mathbf{A}_{t-T}^{(n)}&t>T\\\
\mu_{t}\mathbf{A}_{t}^{(n)}&t\leq T\end{array}\right.$ (28)

and

$\mathbf{Q}_{t}^{(n)}=\left\\{\begin{array}[]{cl}\mu_{t}\mathbf{B}_{t}^{(n)}-\mathit{\Lambda}_{t-T}\mu_{t-T}\mathbf{B}_{t-T}^{(n)}&t>T\\\
\mu_{t}\mathbf{B}_{t}^{(n)}&t\leq T\end{array}\right.$ (29)

with $\mathit{\Lambda}_{t-T}={\prod\limits_{\tau=1}^{T}}\lambda_{t-\tau}$.
Finally, $\hat{\mathbf{H}}$ in eq. 14 is redefined as

$\hat{\mathbf{H}}_{t}^{(n)}=\bm{\Psi}_{t-1}^{(n)}\mathbf{S}_{t}^{(n)}-\mathbf{Q}_{t}^{(n)}$
(30)

and $\mathbf{G}_{t}^{(n)}$ in eq. 16a will be accordingly modified as

$\mathbf{G}_{t}^{(n)}=\lambda_{t-1}\hat{\mathbf{G}}_{t-1}^{(n)}+\hat{\mathbf{H}}_{t}^{(n)}$
(31)

##  3 Numerical Experiments

In this section, two experiments are conducted to show the benefits of the
proposed method.

The performance is tested against two criteria, the mean square error and the
percentage of atoms recovered. Here, if
$|\mathbf{\psi}^{T}\mathbf{\psi}|>0.95$, then it is said the atom is
recovered. 100 trials are run to compare both T-MOD algorithm and our online
multilinear dictionary learning (shorten for OMDL) method. For fair
comparison, the forgetting factor $\lambda_{t}$ for all experiments will
follow eq. 24 with $\lambda_{0}=0.8$. The tests are run on 3D tensor model
where each mode-$n$ dictionary $\bm{\Psi}^{(n)}\in\textbf{\em
C}^{(n)}\subset\mathbb{R}^{10\times 20}$ for $n=1,2,3$ is generated by
Gaussian random variable with mean and variance equal $0$ and $1$
respectively. $\textbf{\em C}^{(n)}$ is a set of matrices with unit columns.
The sparse core $\mathcal{S}_{t}\in\mathbb{R}^{20\times 20\times 20}$ at each
$t$ has $10$ non-zero elements ($K=10$) randomly selected from Gaussian
distribution. The SNR is set to $50$.

In the first experiment, we test the ability of the proposed method against
other tensor-based dictionary learning algorithms like T-MOD and HO-KSVD. The
result is presented in figs. 1 and 2 and our method shows edges over the T-MOD
algorithm. This is due to the T-MOD involves the inversion of the matrix
derived from contracted sparse core tensor which is highly probably not
invertible; this leads the algorithm to diverge. Unlike T-MOD, our method
avoids matrix inversion by reverting to a more classical method of
optimization. Also unlike LMS-type method where the gradient is a product of
an instantaneous sample, the gradient of our method is built in a similar way
to RLS algorithm where a number of past samples are strategically utilized.

Figure 1: MSE learning curves of T-MOD and OMDL, averaged over 100
realisations, for the identification of 3D tensor signal with sparsity $=10$
and SNR $=50$ Figure 2: Percentage of atoms recovered from T-MOD and OMDL,
averaged over 100 realisations, for the identification of 3D tensor signal
with sparsity $=10$ and SNR $=50$

For the second experiment, we compare how our proposed method behaves under
different level of sparsity which is expressed as a ratio between the number
of non-zero elements and total number of elements of output tensor which is,
in this case, $10\times 10\times 10=1000$ (e.g. $k=10$ means sparsity
$=0.01$). It is clear from figs. 3 and 4 that as sparsity increases, it will
take longer time for the algorithm to converge because the less sparse the
core, the more complex the observed data, and consequently the longer time it
takes to resolve the mixture.

Figure 3: MSE learning curves of OMDL for different levels of sparsity,
averaged over 100 realisations, for the identification of 3D tensor signal
with SNR $=50$ Figure 4: Percentage of atoms recovered from OMDL for different
level of sparsity, averaged over 100 realisations, for the identification of
3D tensor signal with sparsity $=10$ and SNR $=50$

##  4 Conclusion

We have introduced an online dictionary learning for the processing of tensor
signals the underlying structure of which is sparse. This has been achieved by
extending the online 1-D dictionary learning routines to the tensor domain in
the light of separable dictionaries, and then finalising the recursive
solution via dual update of gradient matrices. Unlike its 1-D counterpart,
this method follows the least square pathway to obtain the recursion, so that
the past calculations are utilized efficiently which reduces computation
further. The method then further improved by incorporating forgetting factor,
sliding window and correcting weight to eliminate erroneous dictionaries,
accelerate convergence and dispose of outliers, respectively. Simulations on
benchmark 3-D synthetic signals supported the analysis.

## References

  * [1] R. Rubinstein, A. M. Bruckstein and M. Elad, “Dictionaries for Sparse Representation Modeling,” in Proceedings of the IEEE, vol. 98, no. 6, pp. 1045-1057, June 2010.
  * [2] J. M. Duarte-Carvajalino and G. Sapiro, “Learning to Sense Sparse Signals: Simultaneous Sensing Matrix and Sparsifying Dictionary Optimization,” in IEEE Transactions on Image Processing, vol. 18, no. 7, pp. 1395-1408, July 2009.
  * [3] Kjersti Engan, Sven Ole Aase, John Håkon Husøy, “Multi-frame compression: theory and design,” in Signal Processing, Volume 80, Issue 10, October 2000, Pages 2121-2140.
  * [4] M. Aharon, M. Elad and A. Bruckstein, “K -SVD: An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation,” in IEEE Transactions on Signal Processing, vol. 54, no. 11, pp. 4311-4322, Nov. 2006.
  * [5] Kjersti Engan, Karl Skretting, John Håkon Husøy, “Family of iterative LS-based dictionary learning algorithms, ILS-DLA, for sparse signal representation,” in Digital Signal Processing, Volume 17, Issue 1, 2007, Pages 32-49.
  * [6] M. Yaghoobi, T. Blumensath and M. E. Davies, “Dictionary Learning for Sparse Approximations With the Majorization Method,” in IEEE Transactions on Signal Processing, vol. 57, no. 6, pp. 2178-2191, June 2009.
  * [7] M. Aharon and M. Elad, “Sparse and Redundant Modeling of Image Content Using an Image-Signature-Dictionary,” in SIAM Journal on Imaging Sciences, vol. 1, no. 3, pp. 228-247, 2008.
  * [8] Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro. 2009. “Online dictionary learning for sparse coding”. In Proceedings of the 26th Annual International Conference on Machine Learning (ICML ’09). ACM, New York, NY, USA, 689-696.
  * [9] K. Skretting and K. Engan, “Recursive Least Squares Dictionary Learning Algorithm,” in IEEE Transactions on Signal Processing, vol. 58, no. 4, pp. 2121-2130, April 2010.
  * [10] Y. Naderahmadian, S. Beheshti and M. A. Tinati, “Correlation Based Online Dictionary Learning Algorithm,” in IEEE Transactions on Signal Processing, vol. 64, no. 3, pp. 592-602, Feb.1, 2016.
  * [11] G. Zhou, Q. Zhao, Y. Zhang, T. Adalı, S. Xie and A. Cichocki, “Linked Component Analysis From Matrices to High-Order Tensors: Applications to Biomedical Data,” in Proceedings of the IEEE, vol. 104, no. 2, pp. 310-331, Feb. 2016.
  * [12] A. Anandkumar, R. Ge, D. Hsu, S. M. Kakade, and M. Telgarsky. 2014. “Tensor decompositions for learning latent variable models,” J. Mach. Learn. Res. 15, 1 (January 2014), 2773-2832.
  * [13] A. Cichocki et al., “Tensor Decompositions for Signal Processing Applications: From two-way to multiway component analysis,” in IEEE Signal Processing Magazine, vol. 32, no. 2, pp. 145-163, March 2015.
  * [14] Sidiropoulos, N. D. and De Lathauwer, L. and Fu, X. and Huang, K. and Papalexakis, E. E. and Faloutsos, C., “Tensor Decomposition for Signal Processing and Machine Learning,” in arXiv e-prints, 1607.01668, Jul., 2016.
  * [15] C. F. Caiafa and A. Cichocki, “Computing Sparse Representations of Multidimensional Signals Using Kronecker Bases,” in Neural Computation, vol. 25, no. 1, pp. 186-220, Jan. 2013.
  * [16] G. Duan, H. Wang, Z. Liu, J. Deng and Y. W. Chen, “K-CPD: Learning of overcomplete dictionaries for tensor sparse coding,” in Proceedings of the 21st International Conference on Pattern Recognition (ICPR2012), Tsukuba, 2012, pp. 493-496.
  * [17] F. Roemer, G. Del Galdo and M. Haardt, “Tensor-based algorithms for learning multidimensional separable dictionaries,” 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Florence, 2014, pp. 3963-3967.
  * [18] Zhang, Zemin and Shuchin Aeron. “Denoising and Completion of 3D Data via Multidimensional Dictionary Learning.” in CoRR abs/1512.09227, 2016.
  * [19] S. Hawe, M. Seibert and M. Kleinsteuber, “Separable Dictionary Learning,” 2013 IEEE Conference on Computer Vision and Pattern Recognition, Portland, OR, 2013, pp. 438-445.