#  Online Orthogonal Dictionary Learning Based on Frank-Wolfe Method

Ye Xue, and Vincent LAU Y. Xue and V. Lau are with the Department of
Electronic and Computer Engineering, Hong Kong University of Science and
Technology, Hong Kong (E-mail: yxueaf, eeknlau@ust.hk).

###### Abstract

Dictionary learning is a widely used unsupervised learning method in signal
processing and machine learning. Most existing works on dictionary learning
adopt an offline approach, and there are two main offline ways of conducting
it. One is to alternately optimize both the dictionary and the sparse code,
while the other is to optimize the dictionary by restricting it over the
orthogonal group. The latter, called orthogonal dictionary learning, has a
lower implementation complexity, and hence, is more favorable for low-cost
devices. However, existing schemes for orthogonal dictionary learning only
work with batch data and cannot be implemented online, making them
inapplicable for real-time applications. This paper thus proposes a novel
online orthogonal dictionary scheme to dynamically learn the dictionary from
streaming data, without storing the historical data. The proposed scheme
includes a novel problem formulation and an efficient online algorithm design
with convergence analysis. In the problem formulation, we relax the orthogonal
constraint to enable an efficient online algorithm. We then propose the design
of a new Frank-Wolfe-based online algorithm with a convergence rate of
$\mathcal{O}(\ln t/t^{1/4})$. The convergence rate in terms of key system
parameters is also derived. Experiments with synthetic data and real-world
internet of things (IoT) sensor readings demonstrate the effectiveness and
efficiency of the proposed online orthogonal dictionary learning scheme.

###### Index Terms:

online learning, dictionary learning, Frank-Wolfe, sensor, convergence
analysis.

##  I Introduction

Sparse representation of data has been widely used in signal processing,
machine learning and data analysis, and shows a highly expressive and
effective representation ability [1, 2, 3, 4]. It represents the data
$\bm{y}\in\mathbb{R}^{N}$ by a linear combination
$\bm{y}\approx\bm{D}^{true}\bm{x}$, where $\bm{x}\in\mathbb{R}^{M}$ is a
sparse code (i.e., the number of non-zero entries of $\bm{x}$ is much smaller
than $M$), and $\bm{D}^{true}\in\mathbb{R}^{N\times M}$ is the dictionary that
contains the compact information of $\bm{y}$. At the initial stage in this
line of research, predefined dictionaries based on a Fourier basis and various
types of wavelets were successfully used for signal processing. However, using
a learned dictionary instead of a generic one has been shown to dramatically
improve the performance on various tasks, e.g., image denoising and
classification.

One method to learn a dictionary is to alternately optimize (AO) problems with
both the dictionary and the sparse code as the variables [1, 2, 3, 4, 5]. In
this approach, the dictionary usually has no constraints or has a bounded norm
constraint on each atom111One column of the dictionary matrix $\bm{D}$ is
called an atom. of the dictionary to prevent trivial solutions [3, 4]. Another
method is to restrict the dictionary over the orthogonal group
$\mathbb{O}(N,\mathbb{R})$ and solve the following optimization problem only
for the dictionary [6, 7, 8]:

$\displaystyle\underset{\bm{D}\in\mathbb{O}(N,\mathbb{R})}{\text{minimize}}$
$\displaystyle\frac{1}{T}\sum^{T}_{t=1}{\text{S}p}(\bm{D}^{\text{T}}\bm{y}_{t}),$
(1)

where ${\text{S}p}(\cdot)$ is a sparsity-promoting function. This formulation
is motivated by the fact that the sparse code can be obtained as
$\\{\bm{x}_{t}=({\bm{D}}^{true})^{T}\bm{D}^{true}\bm{x}_{t}\approx({\bm{D}}^{true})^{\text{T}}\bm{y}_{t}\\}^{T}_{t=1}$
if the dictionary $\bm{D}^{true}$ is orthogonal.

To enable efficient data processing, the latter method, called orthogonal
dictionary learning (ODL), is more favorable than the AO method for the
following reasons. First, ODL has a lower computational complexity. It updates
only the dictionary at each iteration, while the AO method requires solving
two sub-problems respectively for the dictionary and the sparse code. Second,
ODL has a lower sample complexity since it restricts the dictionary over a
smaller optimization space.222Though the orthogonal constraint seemingly
brings performance loss as it narrows the optimization space, ODL has a
performance competitive with that of the AO method in real applications [9].
Third, ODL allows efficient transmission of the dictionary. This is because an
$N\times N$ orthogonal matrix can be represented by $\frac{N(N-1)}{2}$
statistically independent angles via the Givens rotation,333Although matrices
belonging to orthogonal group $\mathbb{O}(N)$, but not the special orthogonal
group $\mathbb{SO}(N)$, cannot be directly factorized by the Givens rotations,
a factorization can be obtained up to a permutation with a negative sign,
e.g., by flipping two columns [10]. and the angles can be quantized
efficiently to achieve a minimum quantization loss [11]. This angle-based
transmission has been adopted in many wireless communication standards [12].

Despite the benefits of ODL, however, the existing ODL methods only work with
batch data. In other words, they require the whole data set to run the
algorithm. Hence, they are not applicable in many real-time applications,
including real-time network monitoring [13], sensor networks [14], and Twitter
analysis [15], where data arrives continuously in rapid, unpredictable, and
unbounded streams. To deal with streaming data, we propose an online ODL
approach that processes the data in a single sample or in a mini-batch. The
online ODL problem can be regarded as taking ${T}\to\infty$ in Problem (1),
and is formally formulated as

$\displaystyle\underset{\bm{D}\in\mathbb{O}(N,\mathbb{R})}{\text{minimize}}$
$\displaystyle F(\bm{D})=\mathbb{E}_{\bm{y}\sim
P}[{\text{S}p}(\bm{D}^{\text{T}}\bm{y})],$ (2)

where $\bm{y}$ is the realization of the random variable $Y$ drawn from a
distribution $P$. Since the distribution $P$ is usually unknown and the cost
of computing the expectation is prohibitive, the main challenge is to solve
Problem (2) without the accessibility of $F(\bm{D})$ or its gradient $\nabla
F(\bm{D})$. In this case, we can only rely on the sampled data to calculate
the approximation of the objective function or the gradient. These
approximations will jeopardize the performance and the convergence of the
algorithms compared to the case where the exact objective value and gradient
are available [16].

Problem (2) is an online constraint optimization problem and many general
algorithms are available for solving this type of problem, such as regularized
dual averaging (RDA) [17] and stochastic mirror descent (SMD) [18]. However,
none can be directly applied to Problem (2) due to the nonconvex constraint
set and possibly nonconvex sparsity-promoting function, e.g.,
$\text{Sp}(\cdot)=-\|\cdot\|_{p}^{p},(p\in\mathbb{N},\quad p>2)$ [8, 19, 20].
In this work, to enable an efficient online algorithm with a convergence
guarantee, we propose a novel online ODL solution with a convex relaxation of
the orthogonal constraint in Problem (2). After the relaxation, the problem
becomes a nonconvex optimization over a convex set. One could transform this
relaxed problem into an unconstrained problem with a composite nonconvex
objective $\bar{F}(\bm{D})$ (see Example 3 in [17, Section 1.1]) and solve the
transformed problem by ProxSGD [21]. However, the use of ProxSGD gives rise to
o issues: first, the proximal operator in ProxSGD creates a high per-iteration
computational complexity; and second, ProxSGD can only be guaranteed to
converge to an $\epsilon$-stationary point (a point $\bm{D}^{*}$, such that
$\mathbb{E}[\|\nabla\bar{F}(\bm{D}^{*})\|^{2}_{F}]\leq\epsilon$) when the
mini-batch size is increasing by $1/\epsilon$ [21]. To achieve small error,
the mini-batch size needs to be large, which is not suitable for most online
processors with limited memory.

In this work, we propose a Frank-Wolfe-based [22] algorithm, the Nonconvex
Stochastic Frank-Wolfe (NoncvxSFW) method, to solve the relaxed problem
directly without the problem transformation. NoncvxSFW can achieve low-
complexity per-iteration computation thanks to the linear minimization oracle
(LMO) in the Frank-Wolfe method. We also prove that the proposed algorithm
with a single sample or a fixed mini-batch size can be guaranteed to converge
to a stationary point of the relaxed online ODL problem. The main
contributions are summarized as follows.

  * •

Novel Online ODL Formulation: We propose an online ODL problem with an
$\ell_{3}$-norm-based sparsity-promoting function and a convex relaxation of
the orthogonal constraint. We prove that all the optimal solutions of the
original problem are also the optimal solutions of the relaxed problem, which
enables an efficient online algorithm with guaranteed convergence.

  * •

Online Frank-Wolfe-Based Algorithm: We develop an online algorithm, the
NoncvxSFW method, to solve the relaxed optimization problem with a single
sample or a fixed mini-batch size. The convergence is analyzed, and its rate
is shown to be $\mathcal{O}(\ln t/t^{1/4})$, where $t$ is the number of
iterations. As far as we are aware, this is the first non-asymptotic
convergence rate for online nonconvex optimization using the Frank-Wolfe-based
method. The proposed algorithm and the corresponding theoretical results can
also be generalized into general online nonconvex problems with convex
constraints.

  * •

Effective and Efficient Application on IoT Sensor Data Compression: We provide
extensive simulations with both synthetic data and a real-world IoT sensor
data set. The simulation results demonstrate the effectiveness and efficiency
of our proposed online ODL method. They also verify the correctness of our
theoretical results. For the synthetic data, the proposed scheme can achieve
superb performance in terms of the convergence rate and the recovery error,
while for the real-world sensor data, it can achieve a better root-mean-square
error (RMSE) with a higher compression ratio for sensor data compression
compared to the state-of-the-art baselines [8, 23, 24, 4, 25].

The rest of the paper is organized as follows. In Section II we illustrate the
signal model and present the problem formulation for the online ODL. In
Section III, we present the Frank-Wolfe-based online algorithm with non-
asymptotic convergence analysis. Application examples and numerical simulation
results are provided in Section IV and Section V, respectively. Finally,
Section VI summarizes the work.

##  II  Signal Model and Problem Formulation

In this section, we introduce the signal model and the proposed online ODL
problem formulation.

###  II-A Signal Model

In the online ODL, we consider that the data samples arrive in streams. At
time $t$, there is one mini-batch of samples
$\bm{Y}_{t}=[\bm{y}^{1}_{t},\ldots,\bm{y}_{t}^{M_{t}}]\in\mathbb{R}^{N\times{M_{t}}}$
arriving at the processor, where $M_{t}$ is the mini-batch size at time $t$
and $\bm{y}^{j}_{t}$ is the $j$-th sample in the $t$-th mini-batch. Each
sample is assumed to be generated by

$\bm{y}^{j}_{t}=\bm{D}^{true}\bm{x}^{j}_{t},\forall j=1,\ldots,M_{t},\forall
t,$ (3)

where $\bm{D}^{true}$ is the orthogonal dictionary and $\bm{x}^{j}_{t}$ is a
realization of the random variable $X$ drawn from some distribution that
induces sparsity. The basic goal of the online ODL processor is to dynamically
learn the dictionary in an online manner without storing all the historical
samples. That is to say, the processor updates the dictionary $\bm{D}_{t}$ at
time $t$ only according to the samples arriving at that time, i.e.,
$\bm{Y}_{t}$.

###  II-B $\ell_{3}$-Norm-Based Formulation

For the online ODL scheme, the dynamic updating of the dictionary can be done
by solving the generic Problem (2) in an online manner. In Problem (2), the
sparsity-promoting function $Sp(\cdot)$ needs to be carefully designed since
it determines the performance of the online ODL and the complexity of the
algorithm. In this work, we use $Sp(\cdot)=-\|\cdot\|^{3}_{3}$, which results
in the following optimization problem:

$\displaystyle\underset{\bm{D}\in\mathbb{O}(N,\mathbb{R})}{\text{minimize}}$
$\displaystyle F(\bm{D})=\mathbb{E}_{\bm{y}\sim
P}[-\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}].$ (4)

The choice of $-\|\cdot\|^{3}_{3}$ is inspired by the recent result that
minimizing the negative $p$-th power of the $\ell_{p}$-norm
($p\in\mathbb{N},\quad p>2$) with the unit $\ell_{2}$-norm constraint leads to
sparse (or spiky) solutions [8, 19, 26]. An illustration is given in Fig. 1.

Figure 1: Unit spheres of the of $\ell_{p}$ in $\mathbb{R}^{2}$, where
$p=1,2,3,4$. The sparsest points on the unit $\ell_{2}$-norm-sphere, e.g.,
points $(0,1),(0,-1),(1,0)$ and $(-1,0)$ in $\mathbb{R}^{2}$, have the largest
$\ell_{p}$-norm ($p\in\mathbb{N},\quad p>2$).

Compared to the widely used $\ell_{1}$-norm, the negative $\ell_{p}$-norm
formulation allows a convex relaxation of the orthogonal constraint. Hence it
provides flexibility for the algorithm design, as we will illustrate in
Section II-C. Also, the differentiability of this formulation enables a faster
convergence of algorithms. In this paper, we choose $p=3$ since the sample
complexity and the total computation complexity achieve the minimum when $p=3$
among all the choices of $p\quad(p\in\mathbb{N},p>2)$ for maximizing the
$\ell_{p}$-norm over the orthogonal group [19].

###  II-C Orthogonal Constraint Relaxation

After determining the sparsity-promoting function, to facilitate an efficient
online solver with a convergence guarantee, we propose a convex relaxation of
the orthogonal constraint in Problem (4) based on the following Lemma 1.

######  Lemma 1 ([27, 3.4] Convex Hull of the Orthogonal Group).

The convex hull of the orthogonal group is the (closed) unit spectral ball:

$conv(\mathbb{O}(N,\mathbb{R}))=\mathbb{B}_{sp}(N,\mathbb{R}),$

where $\mathbb{B}_{sp}(N,\mathbb{R}):=\\{\bm{X}\in\mathbb{R}^{N\times
N}:\|\bm{X}\|\leq 1\\}$ is the unit spectral ball.

From Lemma 1, we know that the unit spectral ball is the minimal convex set
containing the orthogonal group. Then we propose the following relaxed problem
for the online ODL:

$\boxed{\mathscr{P}:\underset{\bm{D}\in\mathbb{B}_{sp}(N,\mathbb{R})}{\text{minimize}}\quad
F(\bm{D})=\mathbb{E}_{\bm{y}\sim P}[-\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}].}$
(5)

Problem $\mathscr{P}$ is a proper relaxation of Problem (4) since all the
optimal solutions of Problem (4) belong to the set of the optimal solutions of
Problem $\mathscr{P}$ under some general statistical model for $\bm{y}$. We
will show this relationship formally in the following.

We first introduce the following definition of the sign-permutation matrix for
characterizing the optimal solutions.

######  Definition 1.

(Sign-permutation Matrix) The $N$-dimensional sign-permutation matrix
$\bm{\Xi}\in\mathbb{R}^{N\times N}$ is defined as

$\bm{\Xi}=\bm{\Sigma}\bm{\Pi},$ (6)

where $\bm{\Sigma}=diag(\pm\bm{1}_{N})$ with $\bm{1}_{N}$ the $N$-dimensional
all-one vector, and
$\bm{\Pi}=[\bm{e}_{\pi(1)},\bm{e}_{\pi(2)},\ldots,\bm{e}_{\pi(N)}]$, with
$\bm{e}_{n}$ being a standard basis vector and $[\pi(1),\pi(2),\ldots,\pi(N)]$
being any permutations of the $N$ elements.

Next, the relationship between the optimal solutions of Problem (4) and
Problem $\mathscr{P}$ is formally presented in the following Theorem 1.

######  Theorem 1 (Consistency of the Relaxation).

If $\bm{y}$ follows the distribution $P$ such that
$\bm{y}=\bm{D}^{true}\bm{x}$ with $\bm{D}^{true}\in\mathbb{O}(N,\mathbb{R})$
and the entries of $\bm{x}$ being i.i.d Bernoulli Gaussian,444The Bernoulli
Gaussian is a typical statistical model for the sparse coefficient, and it is
widely used for the analysis of dictionary learning [6, 7, 8].
$\bm{x}\overset{i.i.d}{\sim}\mathcal{BG}(\theta)$, then the optimal solution
of Problem (4) is

$\bm{D}^{opt}=\bm{D}^{true}\bm{\Xi}^{\text{T}},$ (7)

which is also an optimal solution of Problem $\mathscr{P}$. In (7), $\bm{\Xi}$
is a sign-permutation matrix.

###### Proof.

See Appendix B. ∎

######  Remark 1.

This consistency of the relaxation no longer holds when $-\|\cdot\|^{3}_{3}$
in Problem $\mathscr{P}$ is replaced by the widely used sparsity-promoting
function $\|\cdot\|_{1}$. When $\|\cdot\|_{1}$ is used, the optimal solution
over the orthogonal group is still $\bm{D}^{true}\bm{\Xi}^{\text{T}}$, but the
optimal solution over the unit spectral ball is $\bm{0}$. Hence, the
relaxation no longer maintains the optimality of the minimizers of the
original problem if $\|\cdot\|_{1}$ is used. The minimizers of the relaxed
Problem $\mathscr{P}$ form a larger set compared to the minimizers of Problem
(4). Hence, Problem $\mathscr{P}$ maintains the optimality of the minimizers
of Problem (4), but the minimizers of Problem $\mathscr{P}$ are not
necessarily the minimizers of Problem (4). However, we can show that if the
dictionary in Problem $\mathscr{P}$ is further restricted to be fullrank, then
the minimizers of Problem $\mathscr{P}$ are also the minimizers of Problem (4)
(see Appendix B). Fortunately, the full-rank condition holds, as shown from
the extensive experiments in Section V. We believe the probability of
encountering a rank-deficient instance is very low when adopting algorithms
with full-rank dictionary initialization and randomness in updating the
dictionary, e.g., the proposed Algorithm 2.

After the convex relaxation, we then focus on solving Problem $\mathscr{P}$,
which is a nonconvex optimization problem over a convex set. The convex set
has a key property: it contains the convex combination of any two points in
the set. That is to say, if we have
$\bm{A},\bm{B}\in\mathbb{B}_{sp}(N,\mathbb{R})$, then

$\eta\bm{A}+(1-\eta)\bm{B}\in\mathbb{B}_{sp}(N,\mathbb{R}),\eta\in(0,1).$ (8)

This property enables an efficient online algorithm with a convergence
guarantee, as we will illustrate in the next section.

##  III Online Nonconvex Frank-Wolfe-Based Algorithm

In this section, we first outline the proposed Frank-Wolfe-based algorithm,
NoncvxSFW, for general online convex-constraint nonconvex problems, and then
specialize it to solve Problem $\mathscr{P}$.

###  III-A NoncvxSFW for General Online Non-Convex Optimization

To solve a general online convex-constraint nonconvex problem,

$\underset{\bm{X}\in\underbrace{\mathcal{C}}_{\text{convex}}}{\text{min}}\quad
F_{gen}(\bm{X})=\mathbb{E}_{\bm{y}\sim
P}[\underbrace{f(\bm{X},\bm{y})}_{\text{nonconvex in $\bm{X}$}}],$ (9)

we require the algorithm to have the following properties:

  * •

Computational Efficiency: Efficient per-iteration computation.

  * •

Theoretical Effectiveness: Theoretical guarantee of convergence to a
stationary point.

To fulfill the above properties, we propose NoncvxSFW, as shown in Algorithm
1, which is a variant of the Stochastic Frank-Wolfe method (SFW) in [23]. The
SFW method and the corresponding analysis can only be applied to solve convex
problems. However, NoncvxSFW and the analysis we propose in this paper are
also applicable to nonconvex problems. In the following, we will elaborate on
how the proposed NoncvxSFW algorithm satisfies the required properties.

Data: $\\{\bm{Y}_{t}\\}_{t=1}^{\infty}$ with
$\bm{Y}_{t}=[\bm{y}^{1}_{t},\ldots,\bm{y}_{t}^{M_{t}}]$

Result: $\\{\bm{X}_{t}\\}_{t=1}^{\infty}$

Initialization: $\bm{G}_{0}=\bm{0}$ and random $\bm{X}_{0}\in\mathcal{C}$

for  _$t=1,2,\ldots$_ do

$\rho_{t}=4(t+1)^{-1/2}$, $\gamma_{t}=2(t+2)^{-3/4}$

1\. Gradient Approximation:
$\bm{G}_{t}=(1-\rho_{t})\bm{G}_{t-1}+\frac{\rho_{t}}{M_{t}}\sum_{j\in[M_{t}]}\nabla
f(\bm{X}_{t-1},\bm{y}^{j}_{t})$

2\. LMO:
$\bm{S}_{t}=\arg\min_{\bm{S}\in\mathcal{C}}\langle\bm{G}_{t},\bm{S}\rangle$

3\. Variable Update:
$\bm{X}_{t}=\mathcal{P}[(1-\gamma_{t})\bm{X}_{t-1}+\gamma_{t}\bm{S}_{t}]$.

Algorithm 1 NoncvxSFW for General Nonconvex Problem

####  III-A1 Computational Efficiency of General Non-Convex Optimization

Algorithm 1 comprises three main steps.

  * •

Step 1 (Gradient Approximation) approximates the true gradient $\nabla
F_{gen}(\bm{X})$ with $\bm{G}_{t}$ in a recursive way. In the calculation,
$M_{t}$ can be fixed along all $t$. Hence, compared to the methods in [28] and
[29] that require an increasing number of stochastic gradient evaluations as
the number of iterations $t$ grows, NoncvxSFW is more computationally
efficient.

  * •

Step 2 (LMO) is a procedure to handle the constraint. It can be regarded as
solving a linear approximation of the objective function over the constraint
set $\mathcal{C}$ using the approximated gradient produced by Step 1. Compared
to the Quadratic Minimization Oracle (QMO) in the proximal-based methods,
e.g., ProxSGD [21], the LMO can be more computationally efficient for many
constraint sets, such as the trace norm and the $\ell_{p}$ balls [30].

  * •

Step 3 (Variable Update) updates the variable $\bm{X}_{t}$ by a simple convex
combination of $\bm{S}_{t}\in\mathcal{C}$ and $\bm{X}_{t-1}\in\mathcal{C}$,
$\mathcal{P}[\bm{X}_{t}=(1-\gamma_{t})\bm{X}_{t-1}+\gamma_{t}\bm{S}_{t}]$,
where $\mathcal{P}[\bm{X}]$ is any operation that satisfies
$\mathcal{P}[\bm{X}]\in\mathcal{C}$ and $F_{gen}(\mathcal{P}[\bm{X}])\leq
F_{gen}(\bm{X})$. According to (8), we have
$\bm{X}_{t}=(1-\gamma_{t})\bm{X}_{t-1}+\gamma_{t}\bm{S}_{t}\in\mathcal{C}$.
This step only requires the output from Step 2 and the variable of the last
iteration and automatically ensures the feasibility of the output
$\bm{X}_{t}$.

####  III-A2 Theoretical Effectiveness of General Non-Convex Optimization

To show the convergence of the NoncvxSFW, we first introduce the following
Frank-Wolfe gap as the measure for the first-order stationarity.

######  Definition 2 (Frank-Wolfe Gap).

The Frank-Wolfe gap at the $t$-th iteration, $g^{gen}_{t}$, is defined as

$g^{gen}_{t}:=\underset{\bm{S}\in\mathcal{C}}{\max}\langle-\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}-\bm{X}_{t-1}\rangle.$ (10)

The Frank-Wolfe gap is a valid first-order stationary measure because of the
following Lemma 2.

######  Lemma 2 (Frank-Wolfe Gap is a Measure for Stationarity).

A point $\bm{X}_{t-1}$ is a stationary point for the optimization problem (9)
if and only if $g^{gen}_{t}=0$.

###### Proof.

See Appendix C. ∎

Then we assume the following conditions hold for the general problem (9).

######  Assumption 1.

  1. 1.

(Bounded Constraint Set) The constraint set $\mathcal{C}$ is bounded with
diameter $diam(\mathcal{C})$ in terms of the Frobenius norm for matrices,
i.e.,

$\|\bm{X}-\bm{Y}\|_{F}\leq
diam(\mathcal{C}),\quad\forall\bm{X},\bm{Y}\in\mathcal{C}.$ (11)

  2. 2.

(Lipschitz Smoothness) $F_{gen}(\bm{X})$ is $L$-smooth over the set
$\mathcal{C}$, i.e.,

$\displaystyle\|\nabla F_{gen}(\bm{X})-\nabla F_{gen}(\bm{Y})\|_{F}\leq
L\|\bm{X}-\bm{Y}\|_{F},$ (12)
$\displaystyle\forall\bm{X},\bm{Y}\in\mathcal{C}.$

  3. 3.

(Unbiased Mini-batch Gradient) The mini-batch gradient
$\frac{1}{M_{t}}\sum_{j\in[M_{t}]}\nabla f(\bm{X}_{t-1},\bm{y}^{j}_{t})$ is an
unbiased estimation of the true gradient $\nabla F_{gen}(\bm{X})$, i.e.,

$\mathbb{E}\Big{[}\frac{1}{M_{t}}\sum_{j\in[M_{t}]}\nabla
f(\bm{X}_{t-1},\bm{y}^{j}_{t})\Big{]}=\nabla F_{gen}(\bm{X}),\quad\forall t.$
(13)

  4. 4.

(Bounded Variance of the Mini-batch Gradient) The variance of the mini-batch
gradient $\frac{1}{M_{t}}\sum_{j\in[M_{t}]}\nabla
f(\bm{X}_{t-1},\bm{y}^{j}_{t})$ is bounded; i.e., for all $t$, we have

$\mathbb{E}\Big{[}\|\frac{1}{M_{t}}\sum_{j\in[M_{t}]}\nabla
f(\bm{X}_{t-1},\bm{y}^{j}_{t})-\nabla
F_{gen}(\bm{X})\|_{F}^{2}\Big{]}\leq\frac{V}{M_{t}}.$ (14)

The above assumptions ensure the convergence of NoncvxSFW, which is formally
stated in the following Theorem 2.

######  Theorem 2 (Convergence of NoncvxSFW for the General Nonconvex
Problem).

If the conditions in Assumption 1 hold, using Algorithm 1 to solve the general
problem (9), we have that the expected Frank-Wolfe gap converges to zero, in
the sense that

$\displaystyle\underset{1<s\leq t}{\inf}\mathbb{E}\Big{[}g^{gen}_{s}\Big{]}$
(15) $\displaystyle\leq$
$\displaystyle\frac{c_{1}(\sqrt{\max\\{C_{0},C_{1}\\}}diam(\mathcal{C})+Ldiam(\mathcal{C})^{2})\ln(t+2)}{(t+3)^{\frac{1}{4}}},$

where $C_{0}=\|\nabla F_{gen}(\bm{X}_{0})\|^{2}_{F}$,
$C_{1}=\frac{4V}{\underset{t}{\min}\\{M_{t}\\}}+2L^{2}diam(\mathcal{C})^{2}$
and $c_{1}$ is some positive constant. In other words, Algorithm 1 is
guaranteed to converge to a stationary point of the general problem (9) at a
rate of $\mathcal{O}(\ln(t)/t^{1/4})$ in expectation.

###### Proof.

See Appendix D. ∎

###  III-B NoncvxSFW for the Proposed Online ODL Problem

In this section, we will apply the proposed NoncvxSFW to solve the proposed
online ODL Problem $\mathscr{P}$. The algorithm is summarized in Algorithm 2.

Data: $\\{\bm{Y}_{t}\\}_{t=1}^{\infty}$ with
$\bm{Y}_{t}=[\bm{y}^{1}_{t},\ldots,\bm{y}_{t}^{M_{t}}]$

Result: $\\{\bm{D}_{t}\\}_{t=1}^{\infty}$

Initialization: $\bm{G}_{0}=\bm{0}$ and random
$\bm{D}_{0}\in\mathbb{O}(N,\mathbb{R})\subset\mathbb{B}_{sp}(N,\mathbb{R})$

for  _$t=1,2,\ldots$_ do

$\rho_{t}=4(t+1)^{-1/2}$, $\gamma_{t}=2(t+2)^{-3/4}$

1\. Gradient Approximation:
$\bm{G}_{t}=(1-\rho_{t})\bm{G}_{t-1}+\frac{\rho_{t}}{M_{t}}\sum_{j\in[M_{t}]}-\nabla\|\bm{D}_{t-1}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}$

2\. LMO: $\bm{U},\bm{\Sigma},\bm{V}^{\text{T}}=\text{SVD}(-\bm{G}_{t})$

$\quad\quad\quad\quad\bm{S}_{t}=\bm{U}\bm{V}^{\text{T}}$

3\. Variable Update:
$\bm{D}_{t}=Polar((1-\gamma_{t})\bm{D}_{t-1}+\gamma_{t}\bm{S}_{t})$.

Algorithm 2 NoncvxSFW for the proposed Online ODL problem

Similar to Section III-A, we will illustrate the NoncvxSFW for the proposed
Online ODL problem from the computational and theoretical aspects.

####  III-B1 Computational Efficiency of the Proposed Online ODL Problem

When adopting NoncvxSFW for Problem $\mathscr{P}$, we can obtain a
computationally efficient ODL algorithm whose complexity will not increase as
$t$ grows.

Step 1 (Gradient Approximation) remains unchanged in Algorithm 1. The sampled
gradient can be expressed as
$-\nabla\|\bm{D}^{T}\bm{y}^{j}_{t}\|^{3}_{3}=-\bm{y}_{t}^{j}(|(\bm{D}^{(t-1)})^{\text{T}}\bm{y}_{t}^{j}|\odot(\bm{D}^{(t-1)})^{\text{T}}\bm{y}_{t}^{j})^{\text{T}}$.
Hence, at each iteration, the time complexity and memory complexity of Step 1
are $\mathcal{O}(N^{2}M_{t})$ ( $M_{t}$ can be fixed along all $t$) and
$\mathcal{O}(N^{2})$, respectively.

Step 2 (LMO) is calculated based on the following Lemma 3.

######  Lemma 3 (LMO for the Unit Spectral Ball).

The minimum value of
$\langle\bm{G},\bm{S}\rangle,\forall\bm{S}\in\mathbb{B}_{sp}(N,\mathbb{R})$ is
the nuclear norm of $-\bm{G}$, i.e.,

$\min_{\bm{S}\in\mathbb{B}_{sp}(N,\mathbb{R})}\langle\bm{G},\bm{S}\rangle=\|-\bm{G}\|_{*}.$
(16)

The minimum is achieved when $\bm{S}$ belongs to the subdifferential of
$\|-\bm{G}\|_{*}$, i.e.,

$\displaystyle\bm{S}^{*}=$
$\displaystyle\arg\min_{\bm{S}\in\mathbb{B}_{sp}(N,\mathbb{R})}\langle\bm{G},\bm{S}\rangle=\bm{U}\bm{V}^{\text{T}}\in\partial\|-\bm{G}\|_{*},$
(17)

where $\bm{U},\bm{\Sigma},\bm{V}^{\text{T}}=\text{SVD}(-\bm{G})$.

###### Proof.

We can prove Lemma 3 by simply using the definition of the dual norm and the
subdifferential of the norm.555 The detailed deduction can be found in
https://stephentu.github.io/blog/convex-analysis/2014/10/01/subdifferential-
of-a-norm.html. ∎

Using the LMO to deal with the constraint is much more computationally
friendly than using the QMO in the proximal-based method. In the QMO, the
proximal operator for the matrix spectral norm requires a proximal operator
for the $\ell_{\infty}$-norm of the singular vector, which has no closed-form
solution [31]. The calculation of the LMO can be simplified via the fact that
$-\bm{G}=\bm{U}\bm{\Sigma}\bm{V}^{\text{T}}=\bm{U}\bm{V}^{\text{T}}\bm{V}\bm{\Sigma}\bm{V}^{\text{T}}=\bm{S}^{*}\bm{V}\bm{\Sigma}\bm{V}^{\text{T}}$.
This indicates that $\bm{S}^{*}$ can be calculated directly from the polar
decomposition of $-\bm{G}$, which has many efficient calculations [32]. Using
Coppersmith-Winograd matrix multiplication [4] in the calculation of polar
decomposition, the time complexity and memory complexity of Step 2 are
$\mathcal{O}(N^{2.38})$ and $\mathcal{O}(N^{2})$, respectively.

In Step 3 (Variable Update), we further adopt polar decomposition, inspired by
its projection property [33, Proposition 3.4], which has
$\mathcal{O}(N^{2.38})$ time complexity and $\mathcal{O}(N^{2})$ memory
complexity. The following lemma shows that this is a valid specification of
$\mathcal{P}$ in Algorithm 1.

######  Lemma 4.

(Validation of Variable Update Step) Let
$\bm{X}\in\mathbb{B}_{sp}(N,\mathbb{R})$. Then we have
$Polar(\bm{X})\in\mathbb{B}_{sp}(N,\mathbb{R})$ and $F(Polar(\bm{X}))\leq
F(\bm{X})$.

###### Proof.

See Appendix E. ∎

####  III-B2 Theoretical Effectiveness of the Proposed Online ODL Problem

In this subsection, we will adapt the convergence theory in Theorem 2 to the
proposed online ODL problem. We first show in the following Lemma 5 that the
conditions in Assumption 1 are satisfied by Problem $\mathscr{P}$.

######  Lemma 5 (The ODL Problem Satisfies the Convergence Condition).

If $\bm{y}$ follows the distribution $P$ such that, for all $t$ and all $j$,
$\bm{y}^{j}_{t}=\bm{D}^{true}\bm{x}^{j}_{t}$, with
$\bm{D}^{true}\in\mathbb{O}(N,\mathbb{R})$ and the entries of $\bm{x}^{j}_{t}$
being i.i.d Bernoulli Gaussian,
$\bm{x}^{j}_{t}\overset{i.i.d}{\sim}\mathcal{BG}(\theta)$, then Problem
$\mathscr{P}$ satisfies the conditions in Assumption 1. Specifically, we have:

  1. 1.

(Bounded Constraint Set)

$\displaystyle\|\bm{D}_{1}-\bm{D}_{2}\|_{F}\leq\sqrt{2N},\forall\bm{D}_{1},\bm{D}_{2}\in\mathbb{B}_{sp}(N,\mathbb{R}).$
(18)

  2. 2.

(Lipschitz Smoothness) $F(\bm{D})$ is $L$-smooth over the set
$\mathbb{B}_{sp}(N,\mathbb{R})$, i.e.,

$\displaystyle\|\nabla F(\bm{D}_{1})-\nabla F(\bm{D}_{2})\|_{F}$ (19)
$\displaystyle\leq$
$\displaystyle\sqrt{\frac{2}{\pi}}N^{3/2}(N+1)\theta\|\bm{D}_{1}-\bm{D}_{2}\|_{F},$
$\displaystyle\forall\bm{D}_{1},\bm{D}_{2}\in\mathbb{B}_{sp}(N,\mathbb{R}).$

  3. 3.

(Unbiased Mini-batch Gradient) The mini-batch gradient
$\frac{1}{M_{t}}\sum_{j\in[M_{t}]}-\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}$
is an unbiased estimation of the true gradient $\nabla F(\bm{D})$, i.e.,

$\mathbb{E}\Big{[}\frac{1}{M_{t}}\sum_{j\in[M_{t}]}-\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}\Big{]}=\nabla
F(\bm{D}),\quad\forall t.$ (20)

  4. 4.

(Bounded Variance of the Mini-batch Gradient) The variance of the mini-batch
gradient
$\frac{1}{M_{t}}\sum_{j\in[M_{t}]}-\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}$
is bounded; i.e., for all $t$, we have

$\mathbb{E}\Big{[}\|\frac{1}{M_{t}}\sum_{j\in[M_{t}]}-\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}-\nabla
F(\bm{D})\|_{F}^{2}\Big{]}\leq\frac{3\theta N^{2}}{M_{t}}.$ (21)

###### Proof.

See Appendix F. ∎

Based on Lemma 5 and Theorem 2, we have the convergence result for Algorithm
2, as shown in the following Theorem 3.

######  Theorem 3 (Convergence of NoncvxSFW for the Proposed Online ODL
Problem).

Using Algorithm 2 to solve Problem $\mathscr{P}$, we have that the expected
Frank-Wolfe gap

$\mathbb{E}\Big{[}g_{t}\Big{]}:=\mathbb{E}\Big{[}\underset{\bm{S}\in\mathbb{B}_{sp}(N,\mathbb{R})}{\max}\langle-\nabla
F(\bm{D}_{t-1}),\bm{S}-\bm{D}_{t-1}\rangle\Big{]}$

converges to zero, in the sense that

$\displaystyle\underset{1<s\leq t}{\inf}\mathbb{E}\Big{[}g_{s}\Big{]}$ (22)
$\displaystyle\leq$ $\displaystyle\frac{c_{2}\Big{(}\sqrt{C_{0}N+\frac{\theta
N^{3}}{\underset{t}{\min}\\{M_{t}\\}}+\theta^{2}N^{5}(N+1)^{2}}+\theta
N^{\frac{5}{2}}(N+1)\Big{)}\ln(t+2)}{(t+3)^{\frac{1}{4}}},$

where $C_{0}=\|\nabla F(\bm{D}_{0})\|^{2}_{F}$ and $c_{2}$ is a positive
constant.

###### Proof.

The Proof can be made by substituting the results in Lemma 5 into Theorem 2. ∎

######  Remark 2 (Impact of the Key System Parameters).

Theorem 3 suggests that a larger value of the smallest mini-batch size
$\underset{t}{\min}\\{M_{t}\\}$, a smaller value of dictionary size $N$, and a
smaller sparsity level $\theta$ (data becomes more sparse with a smaller
$\theta$) will lead to a faster convergence speed. These conclusions are
consistent with the simulation results in Section V-B. Though the above
theorem only proves the convergence to stationary points, we have observed in
experiments that the algorithm actually converges to the global optimum under
very broad conditions, as shown in Section V-B. Similar phenomena have also
appeared in many other works on offline ODL [8, 19, 26, 7, 6].

##  IV  Application Examples

In this subsection, we give two important application examples where the
online ODL method should be adopted.

####  IV-1 Example 1 (Online Data Compression on Edge Devices in the IoT
Network [34])

Consider an IoT network architecture shown in Fig. 2, where the data are
collected from smart objects, such as wearables and industrial sensor devices,
and are sent periodically to an edge device using short-range communication
protocols (e.g., WiFi and Bluetooth). The edge device is responsible for low-
level processing, filtering, and sending the data to the cloud.

Figure 2:  Illustration of edge data compression in an IoT network using the
proposed online ODL scheme.

We assume that an edge device is connected to a total number of $N$
geographically distributed IoT sensors. When sensor measurements (temperature,
humidity, concentration, etc.) are required by the cloud from the edge for the
data analytics, the edge device transmits a compressed version of the data to
save communication resources. At the $t$-th time slot, the ${M_{t}}$ samples
of the sensor measurements from all $N$ sensors,
$\bm{Y}_{t}=[\bm{y}^{1}_{t},\ldots,\bm{y}_{t}^{M_{t}}]\in\mathbb{R}^{N\times{M_{t}}}$,
are transmitted to an edge device. When the $j$-th sample from all $N$
sensors, i.e., $\bm{y}^{j}_{t},j\in[M_{t}]$, is required by the cloud, the
edge data compression is executed using the following steps:

  1. 1.

Preprocessing (Sparse Coding on the Edge): The edge device calculates the
sparse code $\hat{\bm{x}}^{j}_{t}$ based on the latest edge dictionary
$\bm{D}_{t-1}$ and the input $\bm{y}^{j}_{t}$.

  2. 2.

Preprocessing (Transmission Content Decision on the Edge): Upon obtaining the
sparse code $\hat{\bm{x}}^{j}_{t}$, the edge device calculates the error
between $\bm{y}^{j}_{t}$ and $\bm{D}_{cloud}\hat{\bm{x}}^{j}_{t}$ using a
certain error metric $l(\bm{y}^{j}_{t},\bm{D}_{cloud}\hat{\bm{x}}^{j}_{t})$,
where $\bm{D}_{cloud}$ is a local copy of the latest cloud dictionary in the
cloud. Then, the edge decides the content to transmit:

     * •

If the error metric $l(\bm{y}^{j}_{t},\bm{D}_{cloud}\hat{\bm{x}}^{j}_{t})$ is
larger than a predetermined threshold, the edge device updates its local cloud
dictionary copy as $\bm{D}_{cloud}=\bm{D}_{t-1}$, and transmits the updated
$\bm{D}_{cloud}$ to the cloud in a compressed format. It also transmits the
sparse code $\hat{\bm{x}}^{j}_{t}$ to the cloud in a compressed format;

     * •

Otherwise, the edge transmits the sparse code $\hat{\bm{x}}^{j}_{t}$ to the
cloud in a compressed format.

  3. 3.

Core Procedure (Online ODL on the Edge): The edge device runs the online ODL
method to produce $\bm{D}_{t}$ using $\bm{D}_{t-1}$ and the input
$\bm{Y}_{t}$.

  4. 4.

Postprocessing (Sensor Data Recovery on the Cloud): The cloud recovers the
required data $\bm{y}^{j}_{t}$ by
$\hat{\bm{y}}_{t}^{j}=\bm{D}_{cloud}\hat{\bm{x}}^{j}_{t}$, where
$\hat{\bm{x}}^{j}_{t}$ is the sparse code received from the edge and
$\bm{D}_{cloud}$ is the latest dictionary in the cloud.

The proposed ODL module produces the $\bm{D}_{t}$, $\bm{D}_{t-1}$ and
$\bm{D}_{cloud}$, which play critical roles in the above example of data
compression on IoT edge devices.

####  IV-2 Example 2 (Real-time Novel Document Detection [35])

Novel document detection can be used to find breaking news or emerging topics
on social media. In this application,
$\bm{Y}_{t}=[\bm{y}^{1}_{t},\ldots,\bm{y}_{t}^{M_{t}}]\in\mathbb{R}^{N\times{M_{t}}}$
denotes the mini-batch of documents arriving at time $t$, where each column of
$\bm{Y}_{t}$ represents a document at that time, as shown in Fig. 3.

Figure 3:  Illustration of real-time novel document detection using the
proposed online ODL scheme.

Each document is represented by a conventional vector space model such as TF-
IDF [36]. For the mini-batch of documents $\bm{Y}_{t}$ arriving at time $t$,
the novel document detector operates using the following steps:

  1. 1.

Preprocessing (Sparse Coding): For all $\bm{y}^{j}_{t}$ in $\bm{Y}_{t}$, the
detector calculates the sparse code $\hat{\bm{x}}^{j}_{t}$ based on the latest
dictionary $\bm{D}_{t-1}$ and the input $\bm{y}^{j}_{t}$.

  2. 2.

Preprocessing (Novel Document Detection): For all $\bm{y}^{j}_{t}$ in
$\bm{Y}_{t}$, the detector calculates the error between $\bm{y}^{j}_{t}$ and
$\bm{D}_{t-1}\hat{\bm{x}}^{j}_{t}$ with an error metric
$l(\bm{y}^{j}_{t},\bm{D}_{t-1}\hat{\bm{x}}^{j}_{t})$.

     * •

If the error $l(\bm{y}^{j}_{t},\bm{D}_{t-1}\hat{\bm{x}}^{j}_{t})$ is larger
than some predefined threshold, the detector marks the document
$\bm{y}^{j}_{t}$ as novel;

     * •

Otherwise, the detector marks the document $\bm{y}^{j}_{t}$ as non-novel.

  3. 3.

Core Procedure (Online ODL): The detector runs the online ODL method to
produce the new dictionary $\bm{D}_{t}$ using $\bm{D}_{t-1}$ and the input
$\bm{Y}_{t}$.

##  V Experiments

This section provides experiments demonstrating the effectiveness and the
efficiency of our scheme compared to the state-of-the-art prior works. All the
experiments are conducted in Python 3.7 with a 3.6 GHz Intel Core I7
processor.

###  V-A List of Baseline Methods

The baseline methods are listed as follows: 666The proposed method, Baselines
1 and 2 have no hyperparameters. For Baselines $3\sim 5$, grid search is
adopted for tuning the hyperparameters. The grids are drawn around the
hyperparameter values given in the baseline papers [24, 4, 25]. In the
experiments with real-world sensor data, we extend the grids for $\lambda$ to
$[1,20]$ for Baselines 3 and 4 when $\eta_{0}=2,8,10$. We pick the
hyperparameter with the best performance in hindsight for the online method,
Baseline 3, in both the synthetic data and the real-world data experiments.
For the offline methods, Baselines 4 and 5, we adopt the walk-forward
validation method [37] with $4593$ testing data as the holdout data to pick
the hyperparameters in the real-world data experiment.

  * •

Baseline 1 (SFW) [23]: This baseline adopts the recently proposed SFW
algorithm to solve the online ODL problem $\mathscr{P}$ in (5). We compare the
proposed NoncvxSFW to this baseline in terms of the convergence property to
show the effectiveness and efficiency of the proposed NoncvxSFW algorithm for
solving the online ODL problem $\mathscr{P}$.

  * •

Baseline 2 ($\ell_{4}$_NoncvxSFW) [8]: In this baseline, we replace the
sparsity-promoting function $-\|\cdot\|^{3}_{3}$ in the ODL problem
$\mathscr{P}$ with $-\|\cdot\|^{4}_{4}$. Then, we solve the problem by the
NoncvxSFW algorithm. This baseline is adopted to demonstrate the advantage of
the choice of the negative $\ell_{3}$-norm objective in the online ODL
formulation.

  * •

Baseline 3 (Online AODL) [24]: This baseline alternately learns the sparse
code and the dictionary by solving the following optimization problem in an
online manner:

$\displaystyle\underset{\bm{D}\in\mathbb{R}^{N\times
N},\\{\bm{x}_{t}\in\mathbb{R}^{N}\\}_{t=1}^{T}}{\text{minimize}}$
$\displaystyle\sum_{t=1}^{T}\|\bm{y}_{t}-\bm{D}\bm{x}_{t}\|_{F}^{2}+\lambda\|\bm{x}_{t}\|_{1}.$
(23)

The online AODL is introduced to show the advantage of the online ODL scheme.
The Python SPAMS toolbox is used to implement this baseline.777Python code and
documents are available at http://spams-devel.gforge.inria.fr/downloads.html.

  * •

Baseline 4 (Offline-DeepSTGDL) [4]: This baseline is a recently proposed
offline dictionary learning method for behind-the-meter load and photovoltaic
forecasting. In DeepSTGDL, a deep encoder first transforms the load
measurements into a latent space which captures the spatiotemporal patterns of
the data. Then, a spatiotemporal dictionary and a sparse code are alternately
learned to capture the significant spatiotemporal patterns for forecasting.

  * •

Baseline 5 (Offline-NCBDL) [25]: This baseline is an offline nonparametric
Bayesian approach in which the dictionary is inferred from a hierarchical
Bayesian model based on the Beta-Bernoulli process and Gibbs sampling.

###  V-B Experiments with Synthetic Data

####  V-B1 Experiment Settings

We evaluate the convergence property of different online dictionary learning
methods with synthetic data. For all the experiments, we conduct $100$
independent Monte Carlo trials. At the $l$-th trial, we generate the
measurements
$\bm{y}^{j}_{t}(l)={\bm{D}}^{true}(l)\bm{x}^{j}_{t}(l)\quad(j\in[M_{t}],1\leq
t\leq T)$, with the ground truth dictionary ${\bm{D}}^{true}(l)$ drawn
uniformly randomly from the orthogonal group $\mathbb{O}(N,\mathbb{R})$, and
with sparse signals $\bm{x}^{j}_{t}(l)\in\mathbb{R}^{N}$ drawn from an i.i.d.
Bernoulli-Gaussian distribution, i.e.,
$\bm{x}^{j}_{t}(l)\overset{i.i.d}{\sim}\mathcal{BG}(\theta)$. For a fair
comparison, all the methods share the same random initial point at each trial.
Without loss of generality, we set $M_{t}=B,$ $\forall t$ and $T=3\times
10^{3}$ for data generation.

To evaluate the convergence property, the error metric at time index $t$ is
calculated by

$\text{Error}_{t}=\frac{1}{100}\sum_{l=1}^{100}|1-\frac{\|\bm{D}^{\text{T}}_{t}(l){\bm{D}}^{true}(l)\|_{4}^{4}}{N}|.$
(24)

This metric is an averaged measure for the difference between the dictionary
learning result and the ground truth dictionary, since the true dictionary at
the $l$-th trial will be perfectly recovered if
$\frac{\|\bm{D}^{\text{T}}_{t}(l){\bm{D}}^{true}(l)\|_{4}^{4}}{N}=1$ [8, 19].
We compare the $\text{Error}_{t}$ versus the number of iterations of the
proposed method and the online baseline methods as follows.

####  V-B2 Convergence Comparison with Different System Parameters

In Fig. 4, we show the convergence properties under different mini-batch sizes
$B$ with a fixed dictionary size $N=10$ and sparsity level $\theta=0.3$. We
tune $\lambda$ for Baseline 3 at $\lambda=0.1$. The results show that a larger
mini-batch size $B$ can accelerate the convergence of all methods, and the
proposed method with $B=10$ can converge to an error at $10^{-3}$ at around
the $1000$-th time index, which is faster than the baselines.888 If
$t=\infty$, there should be no gap between the results for Baseline 2 and the
proposed method under the same mini-batch size. However, in the simulation,
$t=\infty$ is prohibitive. Therefore, when $t$ is finite, the difference of
the curves comes from the fact that the $\ell_{3}$-norm-based formulation has
a lower sample complexity than the $\ell_{4}$-norm-based formulation [19].
Hence, when we have a finite number of samples ($t$ is finite), the optimal
point of the proposed formulation will be closer to the ground truth than the
formulation in Baseline 2, and therefore shows a better performance.

Figure 4:  Convergence comparison under different mini-batch sizes $B$ with
dictionary size $N=10$ and sparsity level $\theta=0.3$.

In Fig. 5, the convergence curves with different dictionary sizes $N$ are
plotted. We fix the mini-batch size at $B=10$ and the sparsity level at
$\theta=0.3$. We tune $\lambda$ for Baseline 3 at $\lambda=0.1$. The results
show that all the methods have a faster convergence rate under a smaller
dictionary size $N$. In addition, the proposed method can achieve a smaller
error with fewer number of iterations than the baselines.

Figure 5:  Convergence comparison under different dictionary sizes $N$ with
mini-batch size $B=10$ and sparsity level $\theta=0.3$.

In Fig. 6, we compare the convergence properties with different sparsity
levels $\theta$. The mini-batch size and dictionary size are fixed at $B=10$
and $N=10$. We tune $\lambda$ for Baseline 3. Specifically, we set
$\lambda=0.1$ with $\theta=0.3$ and $\lambda=0.05$ with $\theta=0.5$. The
results show that both the proposed method and the baselines have a faster
convergence rate with a more sparse $\bm{x}$ (smaller $\theta$). Moreover, the
proposed method achieves $10^{-3}$ error at around the $2000$-th time index
with $\theta=0.5$, which is faster than the baselines.

Figure 6:  Convergence comparison under different sparsity levels $\theta$
with dictionary size $N=10$ and mini-batch size $B=10$.

###  V-C Experiments with a Real-World Sensor Data Set

####  V-C1 Experiment Data Set

We evaluate the performance of the IoT sensor data compression task with
different dictionary learning methods. The experiments are carried out on the
Airly network air quality data set [38], which records temperature, air
pressure, humidity, and the concentrations of particulate matter from 00:00,
Jan. 1, 2017 to 00:00, Dec. 25, 2017. The sensor readings are measured by a
network of $N=56$ low-cost sensors located in Krakow, Poland, and each sensor
has its own location.999Detailed location information with latitude and
longitude can be found in the data set [38]. There are $8593$ readings for
each item from each sensor sampled per hour. In this work, we use the
temperature readings from all the sensors as the input to the dictionary
learning schemes. Since there is a missing data issue in the raw data, we
replace the missing data with the mean readings over all the sensors at the
times that data are missing. Fig. 7 displays the temperature readings that the
dictionary learning schemes process.

Figure 7:  $56\times 8593$ temperature readings from the Airly network air
quality data set [38]. The readings are from a sensor network with $56$
sensors deployed in different locations in Krakow, Poland. The readings are
collected every hour from Jan 1, 2017 to Dec. 25, 2017.

####  V-C2 Performance Metrics

At each time $t$, $M_{t}$ readings,
$\bm{y}^{j}_{t}\in\mathbb{R}^{56}\quad(j\in[M_{t}],1\leq t\leq T)$, are
uploaded to the dictionary learning processors, and will be approximated by
$\tilde{\bm{y}}^{j}_{t}\quad(j\in[M_{t}],1\leq t\leq T)$ through calculations
of the learned dictionary and the sparse code. To show the compression
performance of the dictionary learning methods, two performance metrics are
calculated under different compression ratios. The compression ratio is
defined as

$\text{compression ratio}=\Big{\lfloor}\frac{56}{\eta_{0}}\Big{\rfloor},$ (25)

where $\eta_{0}$ is the number of nonzero values in the sparse code.

The first performance metric is the RMSE, which is defined by

RMSE
$\displaystyle=\sqrt{\frac{\sum^{T}_{t=0}\sum^{M_{t}}_{j=1}\|\tilde{\bm{y}}^{j}_{t}-\bm{y}^{j}_{t}\|_{2}^{2}}{\sum^{T}_{t=0}\sum^{M_{t}}_{j=1}\|\bm{y}^{j}_{t}\|_{2}^{2}}},$
(26)
$\displaystyle\text{with}\quad\tilde{\bm{y}}^{j}_{t}=\phi(\bm{D}^{est}_{t},\tilde{\bm{x}}^{j}_{t}),\|\tilde{\bm{x}}^{j}_{t}\|_{0}=\eta_{0},$

where $\phi(\cdot)$ is calculations of the dictionary and the sparse code
determined by each method, as we will specify in Section V-C3,
$\bm{D}^{est}_{t}$ is the dictionary for compression at time $t$, and
$\tilde{\bm{x}}^{j}_{t}$ is the sparse code for the $j$-th reading at time $t$
with $\eta_{0}$ nonzero values.

To provide more persuasive results, we have also calculated the HLN-corrected
Diebold-Mariano (HLNDM) [39] test results to quantitatively evaluate the
compression accuracy of different methods from a statistical point of view.
Specifically, the HLNDM statistic is calculated by

$\displaystyle
HLNDM=\sqrt{\frac{T-1-2h+h(h-1)}{T}}\frac{\bar{d}}{\sqrt{\frac{\hat{f}_{d}(0)}{T}}},$
(27)

where we set $h=4$ as the time horizon in the experiments, and
$\bar{d}=\sum_{t=1}^{T}d_{t}=\sum_{t=1}^{T}\Big{(}\text{RMSE}_{t}^{2}(1)-\text{RMSE}_{t}^{2}(2)\Big{)}$
is the average of the distance between the instantaneous RMSE produced by two
different methods (method $1$ and method $2$) at time index $t$. Specifically,
we have $d_{t}=\text{RMSE}_{t}^{2}(1)-\text{RMSE}_{t}^{2}(2)$ and

$\displaystyle\text{RMSE}_{t}(1)$
$\displaystyle=\sqrt{\frac{\sum^{M_{t}}_{j=1}\|\tilde{\bm{y}}^{j}_{t}(1)-\bm{y}^{j}_{t}\|_{2}^{2}}{\sum^{M_{t}}_{j=1}\|\bm{y}^{j}_{t}\|_{2}^{2}}},$
(28)
$\displaystyle\text{with}\quad\tilde{\bm{y}}^{j}_{t}(1)=\phi(\bm{D}^{est}_{t}(1),\tilde{\bm{x}}^{j}_{t}(1)),\|\tilde{\bm{x}}^{j}_{t}(1)\|_{0}=\eta_{0},$

where $\tilde{\bm{y}}^{j}_{t}(1)$ is the estimated readings calculated from
method $1$ for the $j$-th reading at time $t$. Furthermore, in (27), we have

$\displaystyle\hat{f_{d}}(0)=$
$\displaystyle\sum_{k=-(T-1)}^{T-1}\Big{(}\mathbb{I}(\frac{k}{h-1})\frac{1}{T}$
(29)
$\displaystyle\times\sum_{t=|k|+1}^{T}(d_{t}-\bar{d})(d_{t-|k|}-\bar{d})\Big{)},$

where the indicator function is expressed as

$\mathbb{I}(\frac{k}{h-1})=\begin{cases}\begin{array}[]{c}1\\\
0\end{array}&\begin{array}[]{c}|\frac{k}{h-1}|\leq 1,\\\
otherwise.\end{array}\end{cases}$ (30)

To interpret the HLNDM statistic, we define the null hypothesis
$H_{0}:\mathbb{E}(d_{t})=0,\forall t$, which means that the two methods have
the same compression accuracy in terms of statistics. The null hypothesis of
no difference will be rejected if the computed HLNDM statistic falls outside
the range of $-z_{\alpha/2}$ to $z_{\alpha/2}$, i.e.,

$|HLNDM|>z_{\alpha/2},$ (31)

where $z_{\alpha/2}$ is the upper (or positive) $z$-value from the standard
normal table corresponding to half of the desired $\alpha$ significance level
of the test. In the experiment, we set the proposed method as the reference
method $1$ to calculate the HLNDM values for the baselines.

####  V-C3 Experiment Procedures

In the experiments, both online and offline methods are considered. We first
introduce the experiment procedures for the offline methods (Baseline 4 and
Baseline 5).

  * •

Offline Training: The first $4000$ temperature readings are used to train the
weights and dictionaries in Baseline 4 and Baseline 5. For Baseline 4 [4], the
first three terms in the training loss [4, Eq.(4)] are considered, and a
spatiotemporal long short-term memory (ST-LSTM) is trained as the deep encoder
$f_{enc}(\cdot)$ and a $4$-layer rectified linear unit (ReLU) neural network
is trained as the node decoder $f_{n}(\cdot)$ with the Adam optimizer using
$200$ epochs. The edge decoder $f_{e}(\cdot)$ follows the expression in [4,
Eq.(10)]. The hyper-parameters are fine-tuned to be $m=6$, $K=50$, $d_{h}=20$,
$\lambda_{e}=0.25$ and $\lambda_{n}=0.25$ for better performance. In addition,
the choices of $\lambda$ are listed in Table I. For Baseline 5 [25], the upper
bound of the number of atoms for the dictionary is set to $K=80$ and the Beta
distribution parameters are set to $a_{0}=b_{0}=\frac{4000}{8}$ for better
performance. Other parameters have the same values as those in [25, Section
V].

  * •

Compression: The remaining $4593$ readings are grouped into $766$ mini-batches
with mini-batch size $M_{t}=6$.101010The last mini-batch of readings has the
mini-batch size $M_{766}=3$. Then each mini-batch of readings is provided to
the offline dictionary learning methods at one specific time index to test the
performance of the offline learned dictionary. Specifically, for Baseline 4,
the sparse code for the $j$-th test reading at time $t$,
$\bm{y}^{j}_{t}\quad(j\in[M_{t}])$, is calculated by solving the following
sparse coding problem provided in [4]:

$\tilde{\bm{x}}^{j}_{t}=\mathcal{T}_{\eta_{0}}(\arg\min\|f_{enc}(\mathcal{G}(\bm{y}^{j}_{t}))-\bm{D}^{est}_{t}\bm{x}\|_{F}^{2}+\lambda\|\bm{x}\|_{1}),$
(32)

with the sklearn Lasso toolbox.111111Implementation details can be found at
http://scikit-
learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html The
estimated reading $\tilde{\bm{y}}^{j}_{t}$ is calculated by
$\tilde{\bm{y}}^{j}_{t}=\phi(\bm{D}^{est}_{t},\tilde{\bm{x}}^{j}_{t})=f_{n}(\bm{D}^{est}_{t}\tilde{\bm{x}}^{j}_{t}),$
where $\bm{D}^{est}_{t}$ is equal to the offline learned dictionary for all
$t$, $f_{enc}(\cdot)$ and $f_{n}(\cdot)$ are the offline learned encoder and
decoder, and $\mathcal{G}(\bm{y}^{j}_{t})$ is the graph embedding for
$\bm{y}^{j}_{t}$ given by [4, Section II-B]. $\mathcal{T}_{\eta_{0}}(\bm{a})$
selects $\eta_{0}$ elements in $\bm{a}$ with the largest $\ell_{2}$-norm and
sets the other elements to zero. For Baseline 5, the sparse code for the
$j$-th test reading at time $t$, $\bm{y}^{j}_{t}\quad(j\in[M_{t}])$, is
calculated by solving the following problem provided by [25]:

$\tilde{\bm{x}}^{j}_{t}=\underset{\|\bm{x}\|_{0}\leq\eta_{0}}{\arg\min}\|\bm{y}^{j}_{t}-\bm{D}^{est}_{t}\bm{x}\|_{2}^{2},$
(33)

with the sklearn OMP toolbox121212Implementation details can be found at
https://scikit-
learn.org/stable/modules/generated/sklearn.linear_model.OrthogonalMatchingPursuit.html,
and the estimated reading $\tilde{\bm{y}}^{j}_{t}$ is calculated by
$\tilde{\bm{y}}^{j}_{t}=\phi(\bm{D}^{est}_{t},\tilde{\bm{x}}^{j}_{t})=\bm{D}^{est}_{t}\tilde{\bm{x}}^{j}_{t}$,
where $\bm{D}^{est}_{t}$ is equal to the offline learned dictionary for all
$t$.

  * •

Performance Metric Calculation: The RMSE values and the HLNDM test results
under different compression ratios are calculated for Baseline 4 and Baseline
5 according to (26) and (27).

Next, we elaborate on the experiment procedures for the online methods (the
proposed method, Baseline 1, Baseline 2, Baseline 3).

  * •

Initialization: Since there is no offline training in the online methods, the
last $4593$ readings are provided to the online methods. The first $100$
readings are utilized to initialize the dictionary $\bm{D}_{0}$ by the batch
version of each online method with $20$ iterations.

  * •

Online Learning: The remaining $4493$ readings are grouped into $749$ mini-
batches with mini-batch size $M_{t}=6$.131313The last mini-batch of readings
has the mini-batch size $M_{749}=5$. The readings
$\bm{y}^{j}_{t}\quad(j\in[M_{t}])$ at time $t$ are provided to the online
methods to produce $\bm{D}_{t}$, sequentially along $t=1,\ldots,749$.

  * •

Online Compression: Each reading in the mini-batch is compressed by the sparse
code after obtaining the dictionary $\bm{D}_{t}$. For the proposed method,
Baseline 1, and Baseline 2, the sparse code for the $j$-th reading at time
$t$, $\bm{y}^{j}_{t}\quad(j\in[M_{t}])$, is calculated using

$\tilde{\bm{x}}^{j}_{t}=\mathcal{T}_{\eta_{0}}(({\bm{D}_{t}})^{\text{T}}\bm{y}^{j}_{t}).$
(34)

For Baseline 3, the sparse code is

$\tilde{\bm{x}}^{j}_{t}=\mathcal{T}_{\eta_{0}}(\arg\min\|\bm{y}^{j}_{t}-\bm{D}_{t}\bm{x}\|_{F}^{2}+\lambda\|\bm{x}\|_{1}).$
(35)

For all online methods, the estimated reading $\tilde{\bm{y}}^{j}_{t}$ is
calculated by
$\tilde{\bm{y}}^{j}_{t}=\phi(\bm{D}^{est}_{t},\tilde{\bm{x}}^{j}_{t})=\bm{D}_{t}\tilde{\bm{x}}^{j}_{t}$.

  * •

Performance Metric Calculation: Then, the RMSE values and the HLNDM test
results under different compression ratios are calculated for the proposed
method, Baseline 1, Baseline 2, and Baseline 3 according to (26) and (27).

####  V-C4 Performance Comparison and Discussion

The RMSE, HLNDM and CPU time comparison among different dictionary learning
schemes under different compression ratios are given in Table I. The results
demonstrate that the proposed online ODL scheme achieves a lower RMSE than the
other methods for the data compression task under different compression
ratios. The proposed online ODL scheme keeps tracking the input readings and
adapts the dictionary learning within a small optimization space. Hence it is
capable of finding a good dictionary with fewer streaming data received at
each time index.

To interpret the HLNDM results, we adopt a widely used significance level of
$\alpha=0.05$ with $z_{\alpha/2}=1.96$. The HLNDM statistic results and the
RMSE results show that the proposed method has better sensor data compression
performance in most cases. Though the HLNDM between the proposed method and
Baseline 2 is $1.51$ ($<1.96$) when $\eta_{0}=17$, the probability that these
two methods have the same level of accuracy is only $6.55\%$.

We count the per-time-slot CPU time to show the computational efficiency of
the proposed method. For the offline methods, we count the per-time-slot CPU
time on the test reading compression stage; i.e., we ignore the offline
training cost. In this case, the proposed online ODL method still costs fewer
computational resources. This is because it only executes one iteration at
each time index with a simple thresholding of a matrix-vector product to
reconstruct the sensor readings, while the offline methods require solving an
optimization problem.

Fig. 8 shows the compression results of the proposed online ODL method for
$1000$ temperature readings of the $45$-th wireless sensor from 21:00, Oct.
27, 2017 to 12:00, Dec. 08, 2017 under different compression ratios. As shown
in this figure, the proposed method can compress the temperature readings with
a satisfactory accuracy. The maximum compression errors are $0.92$
${}^{\circ}C$ with $\eta_{0}=8$ and $0.23$ ${}^{\circ}C$ with $\eta_{0}=17$.

TABLE I:  The performance of different methods in compressing sensor readings
of temperature in 2017 in Krakow, Poland [38]

Compression ratio |  $28$ $(\eta_{0}=2)$ |  $7$ $(\eta_{0}=8)$ |  $5$ $(\eta_{0}=10)$
---|---|---|---
Methods * | Time (ms) | RMSE | HLNDM | Time (ms) | RMSE | HLNDM | Time (ms) | RMSE | HLNDM
Proposed | 1.018 | 4.82% | Ref | 1.021 | 2.74% | Ref | 1.032 | 2.53% | Ref
Baseline 1 (SFW) | 1.022 | 22.53% | 7.56 | 1.030 | 14.78% | 3.57 | 1.034 | 13.90% | 3.48
Baseline 2 ($\ell_{4}$_NoncvxSFW) | 1.076 | 5.22% | 6.212 | 1.083 | 3.16% | 2.65 | 1.088 | 3.04% | 2.70
Baseline 3 (Online AODL) | 8.932 | 47.45% | 4.38 | 19.58 | 25.06% | 3.29 | 39.65 | 24.71% | 2.11
Baseline 4 (Offline-DeepSTGDL) | 4.606 | 24.8% | 2.42 | 10.14 | 10.66% | 5.47 | 18.04 | 7.27% | 3.41
Baseline 5 (Offline-NCBDL) | 2.717 | 34.72% | 5.94 | 4.028 | 24.87% | 6.17 | 4.347 | 24.72% | 6.51
Compression ratio |  $3$ $(\eta_{0}=17)$ |  $2$ $(\eta_{0}=25)$ |  $1$ $(\eta_{0}=35)$
Methods * | Time (ms) | RMSE | HLNDM | Time (ms) | RMSE | HLNDM | Time (ms) | RMSE | HLNDM
Proposed | 1.036 | 1.97% | Ref | 1.040 | 1.20% | Ref | 1.104 | 0.68% | Ref
Baseline 1 (SFW) | 1.040 | 9.68% | 3.04 | 1.063 | 9.40% | 2.37 | 1.163 | 8.81% | 7.54
Baseline 2 ($\ell_{4}$_NoncvxSFW) | 1.089 | 2.40% | 1.51 | 1.096 | 1.48% | 2.91 | 1.106 | 1.08% | 2.17
Baseline 3 (Online AODL) | 105.9 | 19.54% | 4.14 | 111.2 | 19.49% | 5.34 | 137.7 | 17.90% | 6.15
Baseline 4 (Offline-DeepSTGDL) | 38.69 | 7.00% | 5.28 | 39.29 | 6.74% | 5.21 | 53.08 | 6.15% | 6.28
Baseline 5 (Offline-NCBDL) | 6.057 | 16.86% | 6.80 | 7.550 | 13.51% | 5.33 | 9.675 | 13.42% | 6.24

  * *

For $\eta_{0}=2,8,10,17,25,35$, the hyper-parameter $\lambda$ for Baseline 3
is set to $18,5,1,0.3,0.2,0.1$, and $\lambda$ for Baseline 4 is set to
$20,8,4,0.4,0.4,0.2$, for better performance.

Figure 8:  Sensor temperature reading compression results of the proposed
method for the $45$-th wireless sensor from 21:00, Oct. 27, 2017 to 12:00,
Dec. 08, 2017. Top to bottom rows: $\eta_{0}=8,17$, compression ratio $=7,3$.

###  V-D Generalization Study

To demonstrate the generalization of the proposed Algorithm 1 based on the
generic problem (9), we evaluate the performance of Algorithm 1 on the Sparse
Principal Component Analysis (SPCA) problem [40, 27] which is a powerful
framework to find a sparse principal component for seeking a reasonable trade-
off between the statistical fidelity and interpretability of data.
Specifically, for the single-unit SPCA problem, one can have the following
formulation:

$\displaystyle\underset{\bm{\bm{z}}\in\mathbb{B}^{N}}{min}\mathbb{E}_{\bm{y}\sim
P}[-\bm{\bm{z}}^{T}\bm{y}\bm{y}^{T}\bm{\bm{z}}+\lambda H_{\mu}(\bm{\bm{z}})],$
(36)

where ${\mathbb{B}^{N}}$ is the $N$ dimensional closed unit ball, which is a
convex set, $H_{\mu}(\cdot)$ is the Huber loss with parameter $\mu$ to promote
the sparsity of the principal component, and $\lambda$ is a regularization
parameter. The objective of Problem (36) is non-convex. Problem (36) can be
solved by the proposed generic Algorithm 1, with the LMO over the unit ball as
$\bm{s}_{t}=\arg\min_{\bm{s}\in\mathbb{B}^{N}}\langle\bm{g}_{t},\bm{s}\rangle=\frac{-\bm{g}_{t}}{\|\bm{g}_{t}\|_{2}}$.

We test the convergence property of the proposed Algorithm 1 for Problem (36)
in a number of experimental settings with synthetic data. In the experiments,
we follow the procedures in [40, 27] to generate random data with a covariance
matrix $\bm{\text{$\Sigma$}}=\bm{V\Upsilon V}^{T}\in\mathbb{R}^{N\times N}$
containing a sparse leading eigenvector. To do so, $\Sigma$ is synthesized by
constructing $\bm{\Upsilon}$ and $\bm{V}$. Specifically, we have
$\bm{\Upsilon}=diag(\bm{\gamma})$, where $\bm{\gamma}=[100,1,1,\ldots,1]$. To
synthesize $\bm{V}=[\bm{v}_{1},\bm{v}_{2},\ldots,\bm{v}_{N}]$, we first
construct the sparse leading eigenvector $\bm{v}_{1}$. Specifically, the
$i$-th component of $\bm{v}_{1}$ is constructed by the following expression:

$v_{1,i}=\begin{cases}\begin{array}[]{c}\frac{1}{\text{$\sqrt{q}$}}\\\
0\end{array}&\begin{array}[]{c}i=1,\ldots,q,\\\
otherwise.\end{array}\end{cases}$ (37)

Next, vectors $\bm{v}_{2},\ldots,\bm{v}_{N}$ are randomly drawn from
$\mathbb{S}^{N-1}$ until a full-rank $\bm{V}$ is obtained. Then, we apply the
Gram-Schmidt orthogonalization method to the full-rank $\bm{V}$ to obtain an
orthogonal matrix $\bm{V}$. Finally, $T$ mini-batches of data with mini-batch
size $B$ are generated by drawing independent samples from a zero-mean normal
distribution with covariance matrix $\bm{\Sigma}=\bm{V\Upsilon V}^{T}$. The
resultant data samples are denoted by $\bm{y}_{t}^{j}\quad(j\in[B],1\leq t\leq
T)$, where $t$ is the index for the mini-batch and $j$ is the index for a
sample within one mini-batch.

The target for the algorithm is to find the sparse leading principal component
$\bm{v}_{1}$ with the streaming input $\bm{y}_{t}^{j}\quad(j\in[B],1\leq t\leq
T)$ by solving Problem (36). We set $\mu=0.2$, $\lambda=1$ for the problem and
generate $T=3\times 10^{3}$ mini-batches of data from the above procedures. In
the experiments, we test the recovery error by varying the dimensions $N$, the
mini-batch size $B$, and the number of non-zero elements in the principal
component $q$. The recovery error is defined by

$\text{Error}_{t}=|1-|\bm{z}_{t}^{T}\bm{v}_{1}||$ (38)

since the maximum of $|\bm{z}_{t}^{T}\bm{v}_{1}|$ is $1$, and it is achieved
when $\bm{z}_{t}=\bm{v}_{1}$. The results with $\text{Error}_{t}$ against the
number of iterations are shown in Fig. 9, where $\text{Error}_{t}$ at each $t$
is averaged over $100$ independent Monte Carlo trials. As indicated by the
results under various experimental settings, the proposed method for generic
problem (9) generalizes well to this SPCA problem.

Figure 9:  Performance of the proposed Algorithm 1 in solving the SPCA
problem.

##  VI Conclusion

We have proposed in this paper a novel online orthogonal dictionary learning
scheme with a new relaxed problem formulation, a Frank-Wolfe-based online
algorithm, and the convergence analysis. Experiments on synthetic data and a
real-world data set show the superior performance of our proposed scheme
compared to the baselines and verify the correctness of our theoretical
results. Future focus on application will include a distributed or federated
version of the proposed scheme to further leverage the spatially distributed
data. As for the theory, the sample complexity analysis of the proposed scheme
will be an important direction since it characterizes how close the result
returned by the proposed scheme at a given iteration is to the ground truth
dictionary.

##  Appendix A Auxiliary Lemmas

######  Lemma 6.

Let $\langle\bm{A},\bm{B}\rangle$ be the Frobenius inner product of real
matrices $\bm{A}\in\mathbb{R}^{N\times M}$ and $\bm{B}\in\mathbb{R}^{N\times
M}$, then we have for any $\epsilon>0$

$\langle\bm{A},\bm{B}\rangle\leq\frac{\epsilon^{p}}{p}\|vec(\bm{A})\|_{p}^{p}+\frac{\epsilon^{-q}}{q}\|vec(\bm{B})\|_{q}^{q}.$
(39)

###### Proof.

According to the definition of the Frobenius norm, we have

$\displaystyle\langle\bm{A},\bm{B}\rangle=$ $\displaystyle
vec(\bm{A})^{T}vec(\bm{B})=\sum_{i=1}^{M}\sum_{j=1}^{N}A_{j,i}B_{j,i}$ (40)
$\displaystyle=$ $\displaystyle\sum_{i=1}^{M}\sum_{j=1}^{N}\epsilon
A_{j,i}\epsilon^{-1}B_{j,i}\leq\sum_{i=1}^{M}\sum_{j=1}^{N}(\frac{\epsilon^{p}}{p}A_{j,i}^{p}+\frac{\epsilon^{-q}}{q}B_{j,i}^{q}),$

where the last inequality is according to Young’s Inequality [41]. ∎

######  Lemma 7.

Let $f(s)$ be any decreasing function of $s$. We have

$\int_{\mu=a}^{b+1}f(\mu)d\mu\leq\sum_{s=a}^{b}f(s)\leq\int_{\mu=a-1}^{b}f(\mu)d\mu.$
(41)

##  Appendix B Proof Theorem 1

Let $\bm{W}=\bm{D}^{T}\bm{D}^{true}\in\mathbb{R}^{N\times N}$, we have

$\displaystyle\mathbb{E}_{\bm{y}\sim\mathcal{BG}(\theta)}[-\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}]=-\sum_{n=1}^{n=N}\mathbb{E}[|\bm{W}_{n,:}\bm{x}|^{3}]$
(42) $\displaystyle=$
$\displaystyle-\sum_{n=1}^{n=N}\mathbb{E}[|(\bm{W}_{n,:}\odot\bm{b}^{T})\bm{g}|^{3}],$

where we denote $\bm{b}\sim_{i.i.d}\mathbb{B}_{sp}(\theta)$ and
$\bm{g}\sim_{i.i.d}\mathcal{CN}(0,1)$. Using the rotation-invariant property
of Guassian random variables, we have
$-\mathbb{E}[|\bm{W}_{n,:}\odot\bm{b}^{T}\bm{g}|^{3}]=-\gamma_{1}\mathbb{E}[||\bm{W}_{n,:}\odot\bm{b}^{T}||_{2}^{3}]$
with $\gamma_{1}=\frac{2^{3/2}}{\sqrt{\pi}}$ calculated by the 3rd-order non-
central moment of the Gaussian distribution.

For Problem (4), we have $\bm{D}\in\mathbb{O}(N,\mathbb{R})$ and
$||\bm{W}_{n,:}||_{2}=||\bm{D}_{:,n}^{T}\bm{D}^{true}||_{2}\leq 1$. Hence, we
have
$0\leq\mathbb{E}[||\bm{W}_{n,:}\odot\bm{b}^{T}||_{2}^{3}]\leq\mathbb{E}[||\bm{W}_{n,:}\odot\bm{b}^{T}||_{2}^{2}]=\theta$.
The equality holds if and only if
$||\bm{W}_{n,:}\odot\bm{b}^{T}||_{2}\in\\{0,1\\}$ for all $\bm{b}$, which is
only satisfied at $\bm{W}_{n,:}\in\\{\pm\bm{e}^{T}_{i}:i\in[N]\\}$[7].
Therefore, we have
$\mathbb{E}_{\bm{y}\sim\mathcal{BG}(\theta)}[-\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}]=-\sum_{n=1}^{n=N}\mathbb{E}[||\bm{W}_{n,:}\odot\bm{b}^{T}||_{2}^{3}]\geq-N\gamma_{1}\theta$
and the minimum is achieved when
$\bm{W}_{n,:}\in\\{\pm\bm{e}^{T}_{i}:i\in[N]\\}$. Furthermore, if
$\bm{W}_{n_{1},:}=\pm\bm{e}^{T}_{i}$ and $\bm{W}_{n_{2},:}=\pm\bm{e}^{T}_{i}$,
then $Tr(\bm{W}_{n_{1},:}\bm{W}^{T}_{n_{2},:})=\pm 1$. However, we have
$Tr(\bm{W}_{n_{1},:}\bm{W}^{T}_{n_{2},:})=Tr((\bm{D}_{:,n_{1}}^{T}\bm{D}^{true})(\bm{D}_{:,n_{2}}^{T}\bm{D}^{true})^{T})=Tr(\bm{D}_{:,n_{1}}^{T}\bm{D}_{:,n_{2}})=0$.
This indicates that two different columns of $\bm{W}$ cannot simultaneously
equal the same standard basis vector. Hence,
$\mathbb{E}_{\bm{y}\sim\mathcal{BG}(\theta)}[-\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}]$
achieves the minimum $-N\gamma_{1}\theta$ when $\bm{D}=\bm{D}^{opt}$ with
$(\bm{D}^{opt})^{T}\bm{D}^{true}=\bm{\Xi}$. In other words, the optimal
solution of Problem (4) is $\bm{D}^{opt}=\bm{D}^{true}\bm{\Xi}^{T}$.

For Problem $\mathscr{P}$, we have $\bm{D}\in\mathbb{B}_{sp}(N,\mathbb{R})$,
but $||\bm{W}_{n,:}||_{2}=||\bm{D}_{:,n}^{T}\bm{D}^{true}||_{2}\leq 1$ still
holds. Hence, we also have
$\mathbb{E}_{\bm{y}\sim\mathcal{BG}(\theta)}[-\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}]=-\sum_{n=1}^{n=N}\mathbb{E}[||\bm{W}_{n,:}\odot\bm{b}^{T}||_{2}^{3}]\geq-N\gamma_{1}\theta$
and the minimum is achieved when
$\bm{W}_{n,:}\in\\{\pm\bm{e}^{T}_{i}:i\in[N]\\}$. Since
$\bm{D}^{opt}=\bm{D}^{true}\bm{\Xi}^{T}\in\mathbb{B}_{sp}(N,\mathbb{R})$ is
feasible for Problem $\mathscr{P}$, $\bm{D}^{opt}$ is also an optimal solution
for Problem $\mathscr{P}$. This finishes the proof of Theorem 1.

Then we will demonstrate that if the minimizers of Problem $\mathscr{P}$ are
restricted to be full rank, then these minimizers are also the minimizers of
Problem (4). The minimum of Problem $\mathscr{P}$ is achieved when
$\bm{W}_{n,:}\in\\{\pm\bm{e}^{T}_{i}:i\in[N]\\}$. If there are two columns of
$\bm{W}$ that can simultaneously equal the same standard basis vector, say
$\bm{W}_{n_{1},:}=\pm\bm{e}^{T}_{i}$ and $\bm{W}_{n_{2},:}=\pm\bm{e}^{T}_{i}$,
then $|Tr(\bm{W}_{n_{1},:}\bm{W}^{T}_{n_{2},:})|=1$. However, if we we have
$\bm{D}\in\mathbb{B}_{sp}(N,\mathbb{R})$ and $\bm{D}$ being full rank, then
$|Tr(\bm{W}_{n_{1},:}\bm{W}^{T}_{n_{2},:})|=|Tr((\bm{D}_{:,n_{1}}^{T}\bm{D}^{true})(\bm{D}_{:,n_{2}}^{T}\bm{D}^{true})^{T})|=|Tr(\bm{D}_{:,n_{1}}^{T}\bm{D}_{:,n_{2}})|<1$.

##  Appendix C Proof of Lemma 2

The definition of the stationary points for a constraint problem is

######  Definition 3 (Definition of Stationary Points).

We will say that $\bm{X}^{*}\in\mathcal{C}$ is a stationary point if

$\langle\nabla F_{gen}(\bm{X}^{*}),\bm{S}-\bm{X}^{*}\rangle\geq
0,\forall\bm{S}\in\mathcal{C}.$

Then, we first show the sufficient condition. If $g^{gen}_{t}=0$, we have
$\underset{\bm{S}\in\mathcal{C}}{\min}\langle\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}-\bm{X}_{t-1}\rangle=0.$ According to Definition
3, obviously $\bm{X}_{t-1}$ is a stationary point.

Next, we show the necessary condition. If $\langle\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}-\bm{X}_{t-1}\rangle\geq
0,\forall\bm{S}\in\mathcal{C}$, then we have $\langle-\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}-\bm{X}_{t-1}\rangle\leq
0,\forall\bm{S}\in\mathcal{C},$ which indicates

$\displaystyle\underset{\bm{S}\in\mathcal{C}}{\max}\langle-\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}-\bm{X}_{t-1}\rangle=g^{gen}_{t}\leq 0.$ (43)

Let $\bm{S}^{*}=\underset{\bm{S}\in\mathcal{C}}{\arg\max}\langle-\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}\rangle$, then we have $g^{gen}_{t}=\langle-\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}^{*}-\bm{X}_{t-1}\rangle\geq 0$. Combining the
result in (43), we have that if $\langle\nabla
F_{gen}(\bm{X}_{t-1}),\bm{S}-\bm{X}_{t-1}\rangle\geq
0,\forall\bm{S}\in\mathcal{C}$, then $g^{gen}_{t}=0$.

The sufficient and necessary conditions complete the proof.

##  Appendix D  Proof of Theorem 2

For simplicity, we let $F(\cdot)$ represent $F_{gen}(\cdot)$. To prove Theorem
2, we first prove the Iterates Contraction by the following Lemma.

######  Lemma 8 (Iterates Contraction).

Under Assumptions 1, the Frank-Wolfe gap for Problem (9) using Algorithm (1)
satisfies

$\displaystyle\gamma_{t}g^{gen}_{t}$ $\displaystyle\leq
F(\bm{X}_{t})-F(\bm{X}_{t+1})+\gamma_{t}\sqrt{\epsilon_{t}}diam(\mathcal{C})$
(44) $\displaystyle+\frac{L}{2}\gamma^{2}_{t}Ndiam(\mathcal{C})^{2},$

where $N$ is the dimension of the subspace that $\mathcal{C}$ belongs to and

$\epsilon_{t}:=\|\bm{G}_{t}-\nabla F(\bm{X}_{t-1})\|^{2}_{F}$ (45)

is the gradient estimation error.

###### Proof.

We introduce the auxiliary variable

$\bm{S}^{*}_{t}:=\arg\underset{\bm{S}\in\mathcal{C}}{\max}\langle-\nabla
F(\bm{X}_{t-1}),\bm{S}\rangle$ (46)

for the proof. At the $t$-th iteration, we have

$\displaystyle F(\bm{X}_{t})\leq$ $\displaystyle
F(\bm{X}_{t-1}+\gamma_{t}(\bm{S}_{t}-\bm{X}_{t-1}))$ (47)
$\displaystyle\overset{(a)}{\leq}$ $\displaystyle
F(\bm{X}_{t-1})+\gamma_{t}\langle\nabla
F(\bm{X}_{t-1}),\bm{S}_{t}-\bm{X}_{t-1}\rangle$
$\displaystyle+\frac{L\gamma^{2}_{t}}{2}\|\bm{S}_{t}-\bm{X}_{t-1}\|^{2}_{F}$
$\displaystyle=$ $\displaystyle
F(\bm{X}_{t-1})+\gamma_{t}\langle\bm{G}_{t},\bm{S}_{t}-\bm{X}_{t-1}\rangle$
$\displaystyle+\gamma_{t}\langle\nabla
F(\bm{X}_{t-1})-\bm{G}_{t},\bm{S}_{t}-\bm{X}_{t-1}\rangle$
$\displaystyle+\frac{L\gamma^{2}_{t}}{2}\|\bm{S}_{t}-\bm{X}_{t-1}\|^{2}_{F}$
$\displaystyle\overset{(b)}{\leq}$ $\displaystyle
F(\bm{X}_{t-1})+\gamma_{t}\langle\bm{G}_{t},\bm{S}^{*}_{t}-\bm{X}_{t-1}\rangle+\frac{L\gamma^{2}_{t}}{2}\|\bm{S}_{t}-\bm{X}_{t-1}\|^{2}_{F}$
$\displaystyle+\gamma_{t}\langle\nabla
F(\bm{X}_{t-1})-\bm{G}_{t},\bm{S}_{t}-\bm{X}_{t-1}\rangle$ $\displaystyle=$
$\displaystyle F(\bm{X}_{t-1})+\gamma_{t}\langle\nabla
F(\bm{X}_{t-1})-\bm{G}_{t},\bm{S}_{t}-\bm{S}^{*}_{t}\rangle$
$\displaystyle+\gamma_{t}\langle\nabla
F(\bm{X}_{t-1}),\bm{S}^{*}_{t}-\bm{X}_{t-1}\rangle+\frac{L\gamma^{2}_{t}}{2}\|\bm{S}_{t}-\bm{X}_{t-1}\|^{2}_{F}$
$\displaystyle\overset{(c)}{\leq}$ $\displaystyle
F(\bm{X}_{t-1})+\gamma_{t}\|\nabla
F(\bm{X}_{t-1})-\bm{G}_{t}\|_{F}\|\bm{S}_{t}-\bm{S}^{*}_{t}\|_{F}$
$\displaystyle-$
$\displaystyle\gamma_{t}g^{gen}_{t}+\frac{L\gamma^{2}_{t}}{2}\|\bm{S}_{t}-\bm{X}_{t-1}\|^{2}_{F}$
$\displaystyle=$ $\displaystyle
F(\bm{X}_{t-1})+\gamma_{t}\sqrt{\epsilon_{t}}diam(\mathcal{C})-\gamma_{t}g_{t}+\frac{L\gamma^{2}_{t}}{2}diam(\mathcal{C})^{2}$
$\displaystyle\Rightarrow$ $\displaystyle\gamma_{t}g_{t}\leq
F(\bm{X}_{t-1})-F(\bm{X}_{t})+\gamma_{t}\sqrt{\epsilon_{t}}diam(\mathcal{C})$
$\displaystyle+\frac{L}{2}\gamma^{2}_{t}diam(\mathcal{C})^{2},$

where $(a)$ is from Assumption 1, $(b)$ is from Algorithm (1), and $(c)$ is
from the Cauchy–Schwarz inequality. ∎

The next key ingredient for the proof is the diminishing gradient estimation
error as formally shown in the following Lemma.

######  Lemma 9 (Diminishing Gradient Estimation Error).

Let the gradient estimation error at the $t$-th iteration in Algorithm (1)
defined as

$\epsilon_{t}:=\|\bm{G}_{t}-\nabla F(\bm{X}_{t-1})\|^{2}_{F}.$ (48)

Using the updating rule of Algorithm (1), we have

$\displaystyle\mathbb{E}\big{[}\epsilon_{t}\big{]}\leq\frac{\max\\{C_{0},C^{t}_{1}\\}}{(t+2)^{1/2}},$
(49)

where $C_{0}=\|\nabla F(\bm{X}_{0})\|^{2}_{F}$, and
$C^{t}_{1}={\frac{16V}{M_{t}}+2L^{2}diam(\mathcal{C})^{2}}$.

###### Proof.

To prove Lemma 9 we have the following bound with $\mathbb{E}_{t}$ the
conditional expectation w.r.t. the randomness sampled at the $t$-th iteration,
conditioned on all randomness up to the $t$-th iteration.

$\displaystyle\mathbb{E}_{t}[\|\bm{G}_{t}-\nabla F(\bm{X}_{t-1})\|^{2}_{F}]$
(50) $\displaystyle=$ $\displaystyle\mathbb{E}_{t}[\|\rho_{t}(\nabla
F(\bm{X}_{t-1})-\frac{1}{|M_{t}|}\sum_{i\in M_{t}}\nabla
f(\bm{X}_{t-1},\bm{y}_{i}))$ $\displaystyle+(1-\rho_{t})(\nabla
F(\bm{X}_{t-1})-\bm{G}_{t-1})\|^{2}_{F}]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{t}[\|\rho_{t}(\nabla
F(\bm{X}_{t-1})-\frac{1}{|M_{t}|}\sum_{i\in M_{t}}\nabla
f(\bm{X}_{t-1},\bm{y}_{i}))$ $\displaystyle+(1-\rho_{t})(\nabla
F(\bm{X}_{t-1})-\nabla F(\bm{X}_{t-2}))$ $\displaystyle+(1-\rho_{t})(\nabla
F(\bm{X}_{t-2})-\bm{G}_{t-1})\|^{2}_{F}]$
$\displaystyle(\ref{graderror1})=\rho^{2}_{t}\mathbb{E}_{t}[\|\nabla
F(\bm{X}_{t-1})-\frac{1}{|M_{t}|}\sum_{i\i nM_{t}}\nabla
f(\bm{X}_{t-1},\bm{y}_{i})]\|^{2}_{F}$ (51)
$\displaystyle+(1-\rho_{t})^{2}\|\nabla F(\bm{X}_{t-1})-\nabla
F(\bm{X}_{t-2})\|^{2}_{F}$ $\displaystyle+(1-\rho_{t})^{2}\|\nabla
F(\bm{X}_{t-2})-\bm{G}_{t-1}\|^{2}_{F}$
$\displaystyle+2\rho_{t}(1-\rho_{t})\mathbb{E}_{t}[\langle\nabla
F(\bm{X}_{t-1})$ $\displaystyle-\frac{\sum_{i\in M_{t}}\nabla
f(\bm{X}_{t-1},\bm{y}_{i}),\nabla F(\bm{X}_{t-1})-\nabla
F(\bm{X}_{t-2})\rangle}{|M_{t}|}]$
$\displaystyle+2\rho_{t}(1-\rho_{t})\mathbb{E}_{t}[\langle\nabla
F(\bm{X}_{t-1})$ $\displaystyle-\frac{\sum_{i\in M_{t}}\nabla
f(\bm{X}_{t-1},\bm{y}_{i}),F(\bm{X}_{t-2})-\bm{G}_{t-1}\rangle}{|M_{t}|}]+2(1-\rho_{t})^{2}$
$\displaystyle\times$ $\displaystyle\langle\nabla F(\bm{X}_{t-1})-\nabla
F(\bm{X}_{t-2}),\nabla F(\bm{X}_{t-2})-\bm{G}_{t-1}\rangle$
$\displaystyle(\ref{graderror2})=\rho^{2}_{t}\mathbb{E}_{t}[\|\nabla
F(\bm{X}_{t-1})-\frac{1}{|M_{t}|}\sum_{i\in M_{t}}\nabla
f(\bm{X}_{t-1},\bm{y}_{i})\|^{2}_{F}]$ (52)
$\displaystyle+(1-\rho_{t})^{2}\|\nabla F(\bm{X}_{t-1})-\nabla
F(\bm{X}_{t-2})\|^{2}_{F}$ $\displaystyle+$
$\displaystyle(1-\rho_{t})^{2}\|\nabla F(\bm{X}_{t-2})-\bm{G}_{t-1}\|^{2}_{F}$
$\displaystyle+$ $\displaystyle 2(1-\rho_{t})^{2}\langle\nabla
F(\bm{X}_{t-1})-\nabla F(\bm{X}_{t-2}),\nabla
F(\bm{X}_{t-2})-\bm{G}_{t-1}\rangle$ $\displaystyle\overset{(a)}{\leq}$
$\displaystyle\rho_{t}^{2}\frac{V}{|M_{t}|}+(1-\rho_{t})^{2}\|\nabla
F(\bm{X}_{t-1})-\nabla F(\bm{X}_{t-2})\|^{2}_{F}$ $\displaystyle+$
$\displaystyle(1-\rho_{t})^{2}\|\nabla F(\bm{X}_{t-2})-\bm{G}_{t-1}\|^{2}_{F}$
$\displaystyle+$ $\displaystyle 2(1-\rho_{t})^{2}\langle\nabla
F(\bm{X}_{t-1})-\nabla F(\bm{X}_{t-2}),\nabla
F(\bm{X}_{t-2})-\bm{G}_{t-1}\rangle$ $\displaystyle\overset{(b)}{\leq}$
$\displaystyle\rho_{t}^{2}\frac{V}{|M_{t}|}+(1-\rho_{t})^{2}\|\nabla
F(\bm{X}_{t-1})-\nabla F(\bm{X}_{t-2})\|^{2}_{F}$ $\displaystyle+$
$\displaystyle(1-\rho_{t})^{2}\|\nabla F(\bm{X}_{t-2})-\bm{G}_{t-1}\|^{2}_{F}$
$\displaystyle+$
$\displaystyle(1-\rho_{t})^{2}\Big{(}\frac{2}{\rho_{t}}\|\nabla
F(\bm{X}_{t-1})-\nabla F(\bm{X}_{t-2})\|_{F}$
$\displaystyle+\frac{\rho_{t}}{2}\|\nabla
F(\bm{X}_{t-2})-\bm{G}_{t-1}\|_{F}\Big{)}$ $\displaystyle\overset{(c)}{\leq}$
$\displaystyle\rho_{t}^{2}\frac{V}{|M_{t}|}+(1-\rho_{t})^{2}(1+\frac{2}{\rho_{t}})\gamma_{t-1}^{2}L^{2}diam(\mathcal{C})^{2}$
$\displaystyle+(1-\rho_{t})^{2}(1+\frac{\rho_{t}}{2})\|\nabla
F(\bm{X}_{t-2})-\bm{G}_{t-1}\|^{2}_{F}$ $\displaystyle\overset{(d)}{\leq}$
$\displaystyle\rho_{t}^{2}\frac{V}{|M_{t}|}+\frac{2}{\rho_{t}}\gamma_{t-1}^{2}L^{2}diam(\mathcal{C})^{2}$
$\displaystyle+(1-\frac{\rho_{t}}{2})\|\nabla
F(\bm{X}_{t-2})-\bm{G}_{t-1}\|^{2}_{F},$

where $(a)$ is obtained using Assumption 1, $(b)$ is according to Lemma 6 with
$\epsilon=\frac{\rho_{t}}{2}$ and $p=q=2$ , and $(c)$ is obtained from

$\displaystyle\|\nabla F(\bm{X}_{t})-\nabla F(\bm{X}_{t-1})\|^{2}_{F}\leq
L^{2}\|\bm{X}_{t}-\bm{X}_{t-1}\|^{2}_{F}$ (53) $\displaystyle=$ $\displaystyle
L^{2}\gamma^{2}_{t}\|\bm{S}_{t}-\bm{X}_{t-1}\|^{2}_{F}\leq
L^{2}\gamma^{2}_{t}diam(\mathcal{C})^{2}$

via Assumption (1) and Algorithm (1). In $(d)$, we use the inequality
$(1-\rho_{t})^{2}(1+\frac{2}{\rho_{t}})\leq(1-\rho_{t})(1+\frac{2}{\rho_{t}})\leq\frac{2}{\rho_{t}}$
and
$(1-\rho_{t})^{2}(1+\frac{\rho_{t}}{2})\leq(1-\rho_{t})(1+\frac{\rho_{t}}{2})\leq(1-\frac{2}{\rho_{t}}).$

Therefore, we have

$\displaystyle\mathbb{E}_{t}\big{[}\epsilon_{t}\big{]}\leq$
$\displaystyle\big{(}1-\frac{\rho_{t}}{2})\|\bm{G}_{t-1}-\nabla
F(\bm{X}_{t-2})\|^{2}_{F}$ (54)
$\displaystyle+\frac{2L^{2}diam({\mathcal{C}})^{2}\gamma^{2}_{t}}{\rho_{t}}+\rho^{2}_{t}\frac{V}{M_{t}},$

and

$\displaystyle\mathbb{E}[\epsilon_{t}]=\mathbb{E}_{0,1,\ldots,t}[\epsilon_{t}]$
(55) $\displaystyle\leq$
$\displaystyle\big{(}1-\frac{4(t+1)^{-1/2}}{2})\mathbb{E}_{0,1,\ldots,t-1}[\epsilon_{t-1}]$
$\displaystyle+\big{(}2L^{2}diam({\mathcal{C}})^{2}+16\frac{V}{M_{t}}\big{)}(\eta_{0}+t)^{-1}.$

Then, using the result in [42, Lemma 17], we have finished the proof. ∎

We are now able to to prove Theorem 2. Based on Lemma 8, Lemma 9, and Jensen’s
inequality, we have

$\displaystyle\mathbb{E}[\gamma_{t}g^{gen}_{t}]$
$\displaystyle\leq\mathbb{E}[F(\bm{X}_{t-1})-F(\bm{X}_{t})]+\gamma_{t}\sqrt{\mathbb{E}[\epsilon_{t}]}diam(\mathcal{C})$
(56) $\displaystyle+\frac{L}{2}\gamma^{2}_{t}diam(\mathcal{C})^{2}.$

Define $C_{*}=\max F(\bm{X})-\min F(\bm{X})$, we have

$\displaystyle\underset{1\leq s\leq
t}{\inf}\mathbb{E}[g_{s}]\sum^{t}_{s=1}\gamma_{s}\leq$ $\displaystyle
C_{*}+\sum^{t}_{s=1}\big{(}\gamma_{s}\sqrt{\frac{\max\\{C_{0},C^{t}_{1}\\}}{(s+\eta_{0}+1)^{1/2}}}diam(\mathcal{C})$
(57) $\displaystyle+\frac{L}{2}\gamma^{2}_{s}diam(\mathcal{C})^{2}\big{)}.$

Hence, we have

$\displaystyle\underset{1\leq s\leq t}{\inf}\mathbb{E}[g_{s}]$
$\displaystyle\leq\Big{(}C_{*}+\sum^{t}_{s=1}\big{(}\gamma_{s}\sqrt{\frac{\max\\{C_{0},C^{t}_{1}\\}}{(s+\eta_{0}+1)^{1/2}}}diam(\mathcal{C})$
(58)
$\displaystyle+\frac{L}{2}\gamma^{2}_{s}diam(\mathcal{C})^{2}\big{)}\Big{)}/(\sum^{t}_{s=1}\gamma_{s})$
$\displaystyle\leq$
$\displaystyle\Big{(}C_{*}+\sum^{t}_{s=1}\big{(}2(s+2)^{-1}\sqrt{\max\\{C_{0},C^{t}_{1}\\}}diam(\mathcal{C})$
$\displaystyle+\frac{L}{2}4(s+2)^{-3/2}diam(\mathcal{C})^{2}\big{)}\Big{)}/\Big{(}\sum^{t}_{s=1}2(s+2)^{-3/4}\Big{)}$
$\displaystyle\leq$
$\displaystyle\Big{(}C_{*}+\sum^{t}_{s=1}\big{(}(s+2)^{-1}(2\sqrt{\max\\{C_{0},C^{t}_{1}\\}}diam(\mathcal{C})$
$\displaystyle+2Ldiam(\mathcal{C})^{2})\big{)}\Big{)}/\Big{(}\sum^{t}_{s=1}2(s+2)^{-3/4}\Big{)}$

Since both $(s+2)^{-1}$ and $(s+2)^{-3/4}$ are decreasing in terms of $s$,
from Lemma 7, we have

$\sum_{s=1}^{t}(s+2)^{-1}\leq\ln(2+t)-\ln(2),$ (59)

and

$\sum_{s=1}^{t}(s+2)^{-3/4}\geq
4((t+3)^{1/4}-3^{1/4})\geq\frac{2}{5}(t+3)^{1/4},\quad\forall s>1.$ (60)

Hence, we have

$\displaystyle\underset{1<s\leq t}{\inf}\mathbb{E}[g_{s}]$ $\displaystyle\leq
5\Big{(}C_{*}+10(\ln(2+t)-\ln 2)$ (61)
$\displaystyle\times(\sqrt{\max\\{C_{0},C_{1}\\}}diam(\mathcal{C})+Ldiam(\mathcal{C})^{2})\Big{)}$
$\displaystyle/\Big{(}4(t+3)^{1/4}\Big{)}$ $\displaystyle\leq$
$\displaystyle\Big{(}5C_{*}+10\ln(t+2)(\sqrt{\max\\{C_{0},C_{1}\\}}diam(\mathcal{C})$
$\displaystyle+Ldiam(\mathcal{C})^{2})\Big{)}/\Big{(}4(t+3)^{1/4}\Big{)}.$

This finishes the proof.

##  Appendix E Proof of Lemma 4

Since $Polar(\bm{D})\in\mathbb{O}(N,\mathbb{R})$ for
$\bm{D}\in\mathbb{B}_{sp}(N,\mathbb{R})$, we know that
$\mathbb{B}_{sp}(N,\mathbb{R})$ from Lemma 1.

Define $Polar(\bm{D})$ as $\bm{D}^{O}$. Let $\bm{W}=\bm{D}^{T}\bm{D}_{true}$
and $\bm{W}^{O}=(\bm{D}^{O})T\bm{D}_{true}$. It is easy to know that
$\bm{W}^{O}$ has orthonormal rows while rows of $\bm{W}^{O}$ lie in the unit
ball. Using similar calculations as in Eq.(42) and the discussion below, we
have

$-\sum_{n=1}^{N}\mathbb{E}[||\bm{W}^{O}_{n,:}\odot\bm{b}^{T}||_{2}^{3}]\leq-\sum_{n=1}^{N}\mathbb{E}[||\bm{W}_{n,:}\odot\bm{b}^{T}||_{2}^{3}].$
(62)

Hence, we have $F(Polar(\bm{D}))\leq F(\bm{D})$ from Eq.(42). This completes
the proof.

##  Appendix F Proof of Lemma 5

We first show condition (1) is satisfied.
$\forall\bm{D}_{1},\bm{D}_{2}\in\mathbb{B}_{sp}(N,\mathbb{R})$, we have

$\displaystyle\|\bm{D}_{1}-\bm{D}_{2}\|_{F}=\sqrt{trace\big{(}(\bm{D}_{1}-\bm{D}_{2})^{T}(\bm{D}_{1}-\bm{D}_{2})\big{)}}$
(63) $\displaystyle\leq$
$\displaystyle\sqrt{trace\big{(}\bm{D}^{T}_{1}\bm{D}_{1}+\bm{D}^{T}_{2}\bm{D}_{2}\big{)}}\leq\sqrt{N+N}=\sqrt{2N}.$

For condition (2), we define $\bm{d}=vec(\bm{D})$ and
$f_{v}(\bm{d})=f_{v}(vec(\bm{D}))=\nabla
F(\bm{D})=\nabla\mathbb{E}_{\bm{y}\sim\mathcal{P}}[-\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}]$.
Since we have $\|\bm{D}_{1}-\bm{D}_{2}\|_{2}=\|\bm{d}_{1}-\bm{d}_{2}\|_{2}$,
the following holds

$\displaystyle\|\nabla f_{v}(\bm{d}_{1})-\nabla
f_{v}(\bm{d}_{2})\|_{2}\leq\|\mathrm{D}_{\bm{d}}[f_{v}(\bm{d}_{1})]\|\|\bm{d}_{1}-\bm{d}_{2}\|_{2}$
(64)
$\displaystyle\leq\|\mathrm{D}_{\bm{d}}[f_{v}(\bm{d}_{1})]\|_{F}\|\bm{d}_{1}-\bm{d}_{2}\|_{2}$
$\displaystyle=\sqrt{\sum_{i=1}^{N^{2}}\sum_{j=1}^{N^{2}}\mathrm{D}_{\bm{d}}[f_{v}(\bm{d}_{1})]_{i,j}}\|\bm{d}_{1}-\bm{d}_{2}\|_{2},$

where $\mathrm{D}_{\bm{d}}[f_{v}(\bm{d})]$ is the differential of
$f_{v}(\bm{d})$ in terms of $\bm{d}$. To prove condition (2), we only need to
show that

$\sqrt{\sum_{i=1}^{N^{2}}\sum_{j=1}^{N^{2}}\big{(}\mathrm{D}_{\bm{d}}[f_{v}(\bm{d}_{1})]_{i,j}}\big{)}^{2}\leq\sqrt{\frac{2}{\pi}}N^{3/2}(N+1)\theta.$
(65)

Letting $\bm{W}=\bm{D}_{1}^{T}\bm{D}_{true}$ and using the strategy in [7,
B.1], we have the $j$-th element of $f_{v}(\bm{d}_{1})$ being

$f_{v}(\bm{d}_{1})_{j}=\mathbb{E}_{\Omega}\Big{[}\bm{D}_{true([j]\\_N,:)}\big{(}\|\bm{W}^{\Omega}_{(\lceil\frac{j}{N}\rceil,:)}\|\bm{W}^{\Omega}_{(\lceil\frac{j}{N}\rceil,:)}\big{)}^{T}\Big{]},$
(66)

where $\Omega$ is used to denote the generic support set of the Bernoulli
random viable contained in $\bm{x}$ with $\bm{y}=\bm{D}_{true}\bm{x}$. Then,
we have the $i,j$-th element in $\mathrm{D}_{\bm{d}}[f_{v}(\bm{d}_{1})]$ as

$\big{(}\mathrm{D}_{\bm{d}}[f_{v}(\bm{d}_{1})]_{i,j}\big{)}^{2}\begin{cases}\leq\Big{(}\mathbb{E}_{\Omega}\big{[}N\|\bm{W}^{\Omega}_{(\lceil\frac{j}{N}\rceil,:)}\|+\|\bm{W}^{\Omega}_{(\lceil\frac{j}{N}\rceil,:)}\|\big{]}\Big{)}^{2},\\\
\text{if}\quad N(\lceil\frac{j}{N}\rceil-1)+1\leq i\leq
N(\lceil\frac{j}{N}\rceil)\\\ =0,\text{ortherwise}.\\\ \end{cases}$ (67)

From [19, B.1], we know that

$\mathbb{E}_{\Omega}\Big{[}\|\bm{W}^{\Omega}_{(\lceil\frac{j}{N}\rceil,:)}\|\big{]}\leq\sqrt{\frac{2}{\pi}}\theta.$
(68)

Combining (68) with (67), we have (65). This finishes the proof.

Condition (3) holds obviously due to the i.i.d assumption of $\bm{x}_{t}$.

To prove that condition (4) is satisfied, we have

$\displaystyle\mathbb{E}\big{[}\|\frac{1}{M_{t}}\sum_{j\in[M_{t}]}-\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}-\nabla
F(\bm{D})\|_{F}^{2}\big{]}$ (69) $\displaystyle\leq$
$\displaystyle\mathbb{E}\big{[}\|\frac{1}{M_{t}}\sum_{j\in[M_{t}]}-\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}\|_{F}^{2}\big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{M^{2}_{t}}\sum_{j\in[M_{t}]}\mathbb{E}\big{[}trace\big{(}(\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3})^{T}(\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3})\big{)}\big{]}$
$\displaystyle\overset{(a)}{\leq}$
$\displaystyle\frac{1}{M_{t}}\mathbb{E}\big{[}\sigma^{2}\big{(}\nabla\|\bm{D}^{\text{T}}\bm{y}\|^{3}_{3}\big{)}\big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{M_{t}}\mathbb{E}\big{[}\|(\bm{x}^{T}\bm{D}^{T}_{true}\bm{D}\odot|\bm{x}^{T}\bm{D}^{T}_{true}\bm{D}|)\bm{D}_{true}\bm{x}\|^{2}\big{]}$
$\displaystyle\leq$
$\displaystyle\frac{1}{M_{t}}\mathbb{E}\big{[}\||\bm{x}^{T}\bm{D}^{T}_{true}\bm{D}|^{\odot
2}\bm{D}_{true}\bm{x}\|^{2}\big{]}$ $\displaystyle\overset{(b)}{\leq}$
$\displaystyle\frac{1}{M_{t}}\sqrt{\big{(}\mathbb{E}\big{[}\||\bm{x}^{T}\bm{D}^{T}_{true}\bm{D}|^{\odot
2}\|^{4}\big{]}\big{)}}\sqrt{\big{(}\mathbb{E}\big{[}\|\bm{D}_{true}\bm{x}\|^{4}\big{]}\big{)}},$

where (a) is due to that matrix
$\nabla\|\bm{D}^{\text{T}}\bm{y}^{j}_{t}\|^{3}_{3}$ is rank one, (b) is due to
the Cauchy-Schwarz inequality. Let $\bm{D}^{T}_{true}\bm{D}=\bm{W}$, and
$\bm{W}_{:,i}$ represent the $i$-th column vector in $\bm{W}$. We have

  1. 1.

$\displaystyle\mathbb{E}\big{[}\||\bm{x}^{T}\bm{D}^{T}_{true}\bm{D}|^{\odot
2}\|^{4}\big{]}=\mathbb{E}\big{[}(\sum_{i=1}^{N}\|\bm{x}^{T}\bm{W}_{:,i}\|^{2})^{2}\big{]}$
(70) $\displaystyle=$
$\displaystyle\mathbb{E}\big{[}\sum_{i=1}^{N}\|\bm{x}^{T}\bm{W}_{:,i}\|^{4}\big{]}+\mathbb{E}\big{[}\sum_{i^{\prime}\neq
i}\|\bm{x}^{T}\bm{W}_{:,i^{\prime}}\|^{2}\|\bm{x}^{T}\bm{W}_{:,i}\|^{2}\big{]}$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{N}\mathbb{E}\big{[}\|\bm{x}^{T}\bm{W}_{:,i}\|^{4}\big{]}$
$\displaystyle+\sum_{i^{\prime}\neq
i}\big{(}\mathbb{E}\big{[}\|\bm{x}^{T}\bm{W}_{:,i^{\prime}}\|^{4}\big{]}\big{)}^{\frac{1}{2}}\big{(}\mathbb{E}\big{[}\|\bm{x}^{T}\bm{W}_{:,i}\|^{4}\big{]}\big{)}^{\frac{1}{2}}$
$\displaystyle\leq$ $\displaystyle\sum_{i=1}^{N}3\theta+\sum_{i^{\prime}\neq
i}3\theta=N^{2}3\theta,$

where the last inequality is according to [19, Lemma B.1].

  2. 2.

Let $\bm{D}_{true(i,:)}$ represent the $i$-th row vector in $\bm{D}_{true}$.

$\displaystyle\mathbb{E}\big{[}\|\bm{D}_{true}\bm{x}\|^{4}\big{]}=\mathbb{E}\big{[}(\sum_{i=1}^{N}\|\bm{D}_{true(i,:)}\bm{x}\|^{2})^{2}\big{]}$
(71) $\displaystyle\leq$
$\displaystyle\sum_{i=1}^{N}3\theta+\sum_{i^{\prime}\neq
i}3\theta=N^{2}3\theta.$

Substituting (70) and (71) into (69) we finish the proof.

## References

  * [1] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” _IEEE Trans. Signal Process._ , vol. 54, no. 11, pp. 4311–4322, 2006.
  * [2] J. Mairal, M. Elad, and G. Sapiro, “Sparse representation for color image restoration,” _IEEE Trans. on Imag. Process._ , vol. 17, no. 1, pp. 53–69, 2007.
  * [3] M. Khodayar, J. Wang, and Z. Wang, “Energy disaggregation via deep temporal dictionary learning,” _IEEE Trans. on Neural Netw. and Learning Syst._ , vol. 31, no. 5, pp. 1696–1709, 2020.
  * [4] M. Khodayar, G. Liu, J. Wang, O. Kaynak, and M. E. Khodayar, “Spatiotemporal behind-the-meter load and PV power forecasting via deep graph dictionary learning,” _IEEE Trans. on Neural Netw. and Learning Syst._ , pp. 1–15, 2020\.
  * [5] M. Khodayar, M. E. Khodayar, and S. M. J. Jalali, “Deep learning for pattern recognition of photovoltaic energy generation,” _The Electricity J._ , vol. 34, no. 1, p. 106882, 2021.
  * [6] J. Sun, Q. Qu, and J. Wright, “Complete dictionary recovery over the sphere I: Overview and the geometric picture,” _IEEE Trans. Inf. Theory_ , vol. 63, no. 2, pp. 853–884, 2017.
  * [7] Y. Bai, Q. Jiang, and J. Sun, “Subgradient descent learns orthogonal dictionaries,” in _Proc. Int. Conf. Learn. Representations_ , 2019.
  * [8] Y. Zhai, Z. Yang, Z. Liao, J. Wright, and Y. Ma, “Complete dictionary learning via l4-norm maximization over the orthogonal group,” _J. Mach. Learn. Res._ , vol. 21, no. 165, pp. 1–68, 2020.
  * [9] C. Bao, J.-F. Cai, and H. Ji, “Fast sparsity-based orthogonal dictionary learning for image restoration,” in _IEEE Int. Conf. Comput. Vision_ , 2013, pp. 3384–3391.
  * [10] T. Frerix and J. Bruna, “Approximating orthogonal matrices with effective givens factorization,” in _Int. Conf. on Mach. Learn._ PMLR, 2019, pp. 1993–2001.
  * [11] C. Yuen, S. Sun, M. M. S. Ho, and Z. Zhang, “Beamforming matrix quantization with variable feedback rate,” _EURASIP J. Wireless Commun. Netw._ , vol. 2012, no. 1, pp. 1–9, 2012.
  * [12] M. S. Gast, _802.11 ac: a survival guide: Wi-Fi at Gigabit and beyond_. " O’Reilly Media, Inc.", 2013\.
  * [13] S. Babu and J. Widom, “Continuous queries over data streams,” _ACM Sigmod Record_ , vol. 30, no. 3, pp. 109–120, 2001.
  * [14] G. De Francisci Morales, A. Bifet, L. Khan, J. Gama, and W. Fan, “IoT big data stream mining,” in _ACM SIGKDD Int. Conf. Knowl. Discovery and Data Mining_ , 2016, pp. 2119–2120.
  * [15] A. Bifet and E. Frank, “Sentiment knowledge discovery in twitter streaming data,” in _Int. Conf. on Discovery Science_. Springer, 2010, pp. 1–15.
  * [16] R. Johnson and T. Zhang, “Accelerating stochastic gradient descent using predictive reduction,” _Adv. in Neural Inf. Process. Syst._ , vol. 26, pp. 315–323, 2013.
  * [17] L. Xiao, “Dual averaging methods for regularized stochastic learning and online optimization,” _J. Mach. Learn. Res._ , vol. 11, p. 2543–2596, Dec. 2010.
  * [18] Z. Zhou, P. Mertikopoulos, N. Bambos, S. Boyd, and P. W. Glynn, “Stochastic mirror descent in variationally coherent optimization problems,” in _Adv. in Neural Inf. Process. Syst._ , 2017, pp. 7040–7049.
  * [19] Y. Shen, Y. Xue, J. Zhang, K. Letaief, and V. Lau, “Complete dictionary learning via $\ell_{p}$-norm maximization,” in _Proc. Conf. Uncertainty in Artif. Intell. (UAI)_ , vol. 124. Virtual: PMLR, 03–06 Aug 2020, pp. 280–289.
  * [20] Y. Xue, Y. Shen, V. Lau, J. Zhang, and K. B. Letaief, “Blind data detection in massive MIMO via $\ell_{3}$-norm maximization over the Stiefel manifold,” _IEEE Trans. on Wireless Commun._ , pp. 1–1, 2020.
  * [21] S. Ghadimi, G. Lan, and H. Zhang, “Mini-batch stochastic approximation methods for nonconvex stochastic composite optimization,” _Math. Program._ , vol. 155, no. 1-2, pp. 267–305, 2016.
  * [22] M. Frank, P. Wolfe _et al._ , “An algorithm for quadratic programming,” _Nav. Res. Logistics Quart._ , vol. 3, no. 1-2, pp. 95–110, 1956.
  * [23] A. Mokhtari, H. Hassani, and A. Karbasi, “Stochastic conditional gradient methods: From convex minimization to submodular maximization,” _J. Mach. Learn. Res._ , vol. 21, no. 105, pp. 1–49, 2020.
  * [24] J. Mairal, F. Bach, J. Ponce, and G. Sapiro, “Online dictionary learning for sparse coding,” in _Int. Conf. on Mach. Learn._ , 2009, pp. 689–696.
  * [25] N. Akhtar and A. Mian, “Nonparametric coupled bayesian dictionary and classifier learning for hyperspectral classification,” _IEEE Trans. on Neural Netw. and Learning Syst._ , vol. 29, no. 9, pp. 4038–4050, 2017.
  * [26] Y. Xue, V. K. N. Lau, and S. Cai, “Efficient sparse coding using hierarchical riemannian pursuit,” _IEEE Trans. on Signal Process._ , pp. 1–1, 2021.
  * [27] M. Journée, Y. Nesterov, P. Richtárik, and R. Sepulchre, “Generalized power method for sparse principal component analysis,” _J. Mach. Learn. Res._ , vol. 11, no. Feb, pp. 517–553, 2010.
  * [28] E. Hazan and H. Luo, “Variance-reduced and projection-free stochastic optimization,” in _Int. Conf. on Mach. Learn._ , vol. 48. New York, New York, USA: PMLR, 20–22 Jun 2016, pp. 1263–1271.
  * [29] S. J. Reddi, S. Sra, B. Póczos, and A. Smola, “Stochastic Frank-Wolfe methods for nonconvex optimization,” in _54th Annu. Allert. Conf. Commun. Control Comput._ IEEE, 2016, pp. 1244–1251.
  * [30] M. Jaggi, “Revisiting Frank-Wolfe: Projection-free sparse convex optimization,” in _Int. Conf. on Mach. Learn._ , 2013, pp. 427–435.
  * [31] J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra, “Efficient projections onto the l1-ball for learning in high dimensions,” in _Int. Conf. on Mach. Learn._ , ser. ICML ’08. New York, NY, USA: Association for Computing Machinery, 2008, p. 272–279.
  * [32] N. J. Higham and P. Papadimitriou, “A parallel algorithm for computing the polar decomposition,” _Parallel Comput._ , vol. 20, no. 8, pp. 1161–1173, 1994.
  * [33] P.-A. Absil and J. Malick, “Projection-like retractions on matrix manifolds,” _SIAM J. on Optim._ , vol. 22, no. 1, pp. 135–158, 2012.
  * [34] T. Lu, W. Xia, X. Zou, and Q. Xia, “Adaptively compressing iot data on the resource-constrained edge,” in _3rd $\\{$USENIX$\\}$ Workshop on Hot Topics in Edge Comput. (HotEdge 20)_, 2020.
  * [35] S. Kasiviswanathan, H. Wang, A. Banerjee, and P. Melville, “Online l1-dictionary learning with application to novel document detection,” _Adv. in Neural Inf. Process. Syst._ , vol. 25, pp. 2258–2266, 2012.
  * [36] C. D. Manning, P. Raghavan, and H. Schütze, _Introduction to Information Retrieval_. USA: Cambridge University Press, 2008.
  * [37] K. Żbikowski, “Using volume weighted support vector machines with walk forward testing and feature selection for the purpose of creating stock trading strategy,” _Expert Systems with Applications_ , vol. 42, no. 4, pp. 1797–1805, 2015.
  * [38] Airly, “Air quality data from extensive network of sensors: PM1, PM2.5, PM10, temp, pres and hum data for 2017 year from Krakow, Poland.” [Online]. Available: https://www.kaggle.com/datascienceairly/air-quality-data-from-extensive-network-of-sensors
  * [39] D. Harvey, S. Leybourne, and P. Newbold, “Testing the equality of prediction mean squared errors,” _Int. J. Forecasting_ , vol. 13, no. 2, pp. 281–291, 1997.
  * [40] H. Shen and J. Z. Huang, “Sparse principal component analysis via regularized low rank matrix approximation,” _J. Multivariate Anal._ , vol. 99, no. 6, pp. 1015–1034, 2008.
  * [41] W. H. Young, “On classes of summable functions and their Fourier series,” _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ , vol. 87, no. 594, pp. 225–229, 1912\.
  * [42] A. Mokhtari, H. Hassani, and A. Karbasi, “Conditional gradient method for stochastic submodular maximization: Closing the gap,” in _Int. Conf. on Artificial Intell. and Statistics_ , ser. Proceedings of Machine Learning Research, A. Storkey and F. Perez-Cruz, Eds., vol. 84. PMLR, 09–11 Apr 2018, pp. 1886–1895.