# SGD$\\_$Tucker: A Novel Stochastic Optimization Strategy for Parallel Sparse
Tucker Decomposition

Hao Li, _Student Member_ , _IEEE_ , Zixuan Li, Kenli Li, _Senior Member_ ,
_IEEE_ , Jan S. Rellermeyer, _Member_ , _IEEE_ , Lydia Chen, _Senior Member_ ,
_IEEE_ , Keqin Li, _Fellow_ , _IEEE_ Hao Li, Zixuan Li, Kenli Li and Keqin Li
are with the College of Computer Science and Electronic Engineering, Hunan
University, and National Supercomputing Center in Changsha, Hunan, China,
410082. Corresponding author: Kenli Li.
E-mail<EMAIL_ADDRESS>(H.Li-9@tudelft.nl<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Hao Li is also with TU Delft, Netherlands. J.S.Rellermeyer,
Lydia Chen are with TU Delft, Netherlands. Keqin Li is also with the
Department of Computer Science, State University of New York, New Paltz, New
York 12561, USA.

###### Abstract

Sparse Tucker Decomposition (STD) algorithms learn a core tensor and a group
of factor matrices to obtain an optimal low-rank representation feature for
the High-Order, High-Dimension, and Sparse Tensor (HOHDST). However, existing
STD algorithms face the problem of intermediate variables explosion which
results from the fact that the formation of those variables, i.e., matrices
Khatri-Rao product, Kronecker product, and matrix-matrix multiplication,
follows the whole elements in sparse tensor. The above problems prevent deep
fusion of efficient computation and big data platforms. To overcome the
bottleneck, a novel stochastic optimization strategy (SGD$\\_$Tucker) is
proposed for STD which can automatically divide the high-dimension
intermediate variables into small batches of intermediate matrices.
Specifically, SGD$\\_$Tucker only follows the randomly selected small samples
rather than the whole elements, while maintaining the overall accuracy and
convergence rate. In practice, SGD$\\_$Tucker features the two distinct
advancements over the state of the art. First, SGD$\\_$Tucker can prune the
communication overhead for the core tensor in distributed settings. Second,
the low data-dependence of SGD$\\_$Tucker enables fine-grained
parallelization, which makes SGD$\\_$Tucker obtaining lower computational
overheads with the same accuracy. Experimental results show that
SGD$\\_$Tucker runs at least 2$X$ faster than the state of the art.

###### Index Terms:

High-order, High-dimension and Sparse Tensor; Low-rank Representation
Learning; Machine Learning Algorithm; Sparse Tucker Decomposition; Stochastic
Optimization; Parallel Strategy.

##  1 Introduction

Tensors are a widely used data representation style for interaction data in
the Machine Learning (ML) application community [1], e.g, in Recommendation
Systems [2], Quality of Service (QoS) [3], Network Flow [4], Cyber-Physical-
Social (CPS) [5], or Social Networks [6]. In addition to applications in which
the data is naturally represented in the form of tensors, another common used
case is the fusion in multi-view or multi-modality problems [7]. Here, during
the learning process, each modality corresponds to a feature and the feature
alignment involves fusion. Tensors are a common form of feature fusion for
multi-modal learning [7, 8, 9, 10]. Unfortunately, tensors can be difficult to
process in practice. For instance, an $N$-order tensor comprises of the
interaction relationship between $N$ kinds of attributes and if each attribute
has millions of items, this results in a substantially large size of data
[11]. As a remedy, _dimensionality reduction_ can be used to represent the
original state using much fewer parameters [12].

Specifically in the ML community, Tensor Factorization (TF), as a classic
dimensionality reduction technique, plays a key role for low-rank
representation. Xu et al., [13] proposed a Spatio-temporal multi-task learning
model via TF and in this work, tensor data is of $5$-order, i.e., weather,
traffic volume, crime rate, disease incidents, and time. Meanwhile, this model
made predictions through the time-order for the multi-task in weather, traffic
volume, crime rate, and disease incidents orders and the relationship
construction between those orders is via TF. In the community of Natural
Language Processing (NLP), Liu et al., [14] organized a mass of texts into a
tensor and each slice is modeled as a sparse symmetric matrix. Furthermore,
the tensor representation is a widely-used form for Convolutional Neural
Networks (CNNs), e.g., in the popular TensorFlow [15] framework, and Kossaifi
et al., [16] took tensorial parametrization of a CNNs and pruned the
parametrization by Tensor Network Decomposition (TND). Meanwhile, Ju, et al.,
[17] pruned and then accelerated the Restricted Boltzmann Machine (RBM)
coupling with TF. In the Computer Vision (CV) community, Wang et al., [18]
modeled various factors, i.e., pose and illumination, as an unified tensor and
make disentangled representation by adversarial autoencoder via TF. Zhang et
al., [19] constructed multi subspaces of multi-view data and then abstract
factor matrices via TF for the unified latent of each view.

High-Order, High-Dimension, and Sparse Tensor (HOHDST) is a common situation
in the big-data processing and ML application community [20, 21].
Dimensionality reduction can also be used to find the low-rank space of HOHDST
in ML applications [22], which can help to make prediction from existing data.
Therefore, it is non-trivial to learn the pattern from the existed information
in a HOHDST and then make the corresponding prediction. Sparse Tucker
Decomposition (STD) is one of the most popular TF models for HOHDST, which can
find the $N$-coordinate systems and those systems are tangled by a core tensor
between each other [23]. Liu et al., [24] proposed to accomplish the visual
tensor completion via STD. The decomposition process of STD involves the
entanglement of $N$-factor matrices and core tensor and the algorithms follow
one of the following two directions: 1) search for optimal orthogonal
coordinate system of factor matrices, e.g., High-order Orthogonal Iteration
(HOOI) [25]; 2) designing optimization solving algorithms [26]. HOOI is a
common solution for STD [27] and able to find the $N$ orthogonal coordinate
systems which are similar to Singular Value Decomposition (SVD), but requires
frequent intermediate variables of Khatri-Rao and Kronecker products [28].

An interesting topic is that stochastic optimization [29], e.g., Count Sketch
and Singleshot [30, 31], etc. can alleviate this bottleneck to a certain
extent, depending on the size of dataset. However, those methods depend on the
count sketch matrix, which cannot be easily implemented in a distributed
environment and is notoriously difficult to parallelize. The Stochastic
Gradient Descent (SGD) method can approximate the gradient from randomly
selected subset and it forms the basis of the state of art methods, e.g.,
variance SGD [32], average SGD [33], Stochastic Recursive Gradient [34], and
momentum SGD [35]. SGD is adopted to approximate the eventual optimal solver
with lower space and computational complexities; meanwhile, the low data-
dependence makes the SGD method amenable to parallelization [36]. The idea of
construction for the computational graph of the mainstream platforms, e.g.,
Tensorflow [15] and Pytorch [37], is based on the SGD [38] and practitioners
have already demonstrated its powerful capability on large-scale optimization
problems. The computational process of SGD only needs a batch of training
samples rather than the full set which gives the ML algorithm the low-
dependence between each data block and low communication overhead [39].

There are three challenges to process HOHDST in a fast and accurate way: 1)
how to define a suitable optimization function to find the optimal factor
matrices and core tensor? 2) how to find an appropriate optimization strategy
in a low-overhead way and then reduce the entanglement of the factor matrices
with core tensor which may produce massive intermediate matrices? 3) how to
parallelize STD algorithms and make distributed computation with low
communication cost? In order to solve these problems, we present the main
contributions of this work which are listed as follows:

  1. 1.

A novel optimization objective for STD is presented. This proposed objective
function not only has a low number of parameters via coupling the Kruskal
product (Section 4.1) but also is approximated as a convex function;

  2. 2.

A novel stochastic optimization strategy is proposed for STD, SGD$\\_$Tucker,
which can automatically divide the high-dimension intermediate variables into
small batches of intermediate matrices that only follows the index of the
randomly selected small samples; meanwhile, the overall accuracy and
convergence are kept (Section 4.3);

  3. 3.

The low data-dependence of SGD$\\_$Tucker creates opportunities for fine-
grained parallelization, which makes SGD$\\_$Tucker obtaining lower
computational overhead with the same accuracy. Meanwhile, SGD$\\_$Tucker does
not rely on the specific compressive structure of a sparse tensor (Section
4.4).

To our best knowledge, SGD$\\_$Tucker is the first work that can divide the
high-dimension intermediate matrices of STD into small batches of intermediate
variables, a critical step for fine-grained parallelization with low
communication overhead. In this work, the related work is presented in Section
2. The notations and preliminaries are introduced in Section 3. The
SGD$\\_$Tucker model as well as parallel and communication overhead on
distributed environment for STD are showed in Section 4. Experimental results
are illustrated in Section 5.

##  2 Related Studies

For HOHDST, there are many studies to accelerate STD on the state of the art
parallel and distributed platforms, i.e., OpenMP, MPI, CUDA, Hadoop, Spark,
and OpenCL. Ge et al., [40] proposed distributed CANDECOMP/PARAFAC
Decomposition (CPD) which is a special STD for HOHDST. Shaden et al., [41]
used a Compressed Sparse Tensors (CSF) structure which can optimize the access
efficiency for HOHDST. Tensor-Time-Matrix-chain (TTMc)[42] is a key part for
Tucker Decomposition (TD) and TTMc is a data intensive task. Ma et al., [42]
optimized the TTMc operation on GPU which can take advantage of intensive and
partitioned computational resource of GPU, i.e., a warp threads (32) are
automatically synchronized and this mechanism is apt to matrices block-block
multiplication. Non-negative Tucker Decomposition (NTD) can extract the non-
negative latent factor of a HOHDST, which is widely used in ML community.
However, NTD need to construct the numerator and denominator and both of them
involve TTMc. Chakaravarthy et al., [43] designed a mechanism which can divide
the TTMc into multiple independent blocks and then those tasks are allocated
to distributed nodes.

HOOI is a common used TF algorithm which comprises of a series of TTMc
matrices and SVD for the TTMc. [44] focused on dividing the computational task
of TTMc into a list of independent parts and a distributed HOOI for HOHDST is
presented in [45]. However, the intermediate cache memory for TTMc of HOOI
increased explosively. Shaden et al., [41] presented a parallel HOOI algorithm
and this work used a Compressed Sparse Tensors (CSF) structure which can
improve the access efficiency and save the memory for HOHDST. Oh and Park [46,
47] presented a parallelization strategy of ALS and CD for STD on OpenMP. A
heterogeneous OpenCL parallel version of ALS for STD is proposed on [48].
However, the above works are still focused on the naive algorithm
parallelization which cannot solve the problem of fine-grained
parallelization.

##  3 Notation and Preliminaries

The main notations include the tensors, matrices and vectors, along with their
basic elements and operations (in Table I). Fig. 1 illustrates the details of
TD including the tanglement of the core tensor $\mathcal{G}$ with the $N$
factor matrices $\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$. The data styles for
TF include sparse and dense tensors and STD is devoted to HOHDST. Here, basic
definitions for STD and models are rendered in Section 3.1. Finally, the STD
process for the HOHDST is illustrated in Section 3.2.

TABLE I: Table of symbols. Symbol | Definition
---|---
$I_{n}$ | The size of row in the $n$th factor matrix;
$J_{n}$ | The size of column in the $n$th factor matrix;
$\mathcal{X}$ | Input $N$ order tensor $\in$ $\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}_{+}$;
$x_{i_{1},i_{2},\cdots,i_{n}}$ |  $i_{1},i_{2},\cdots,i_{n}$th element of tensor $\mathcal{X}$;
$\mathcal{G}$ | Core $N$ order tensor $\in$ $\mathbb{R}^{J_{1}\times J_{2}\times\cdots\times J_{N}}$;
X | Input matrix $\in$ $\mathbb{R}^{I_{1}\times I_{2}}_{+}$;
$\textbf{X}^{(n)}$ |  $n$th unfolding matrix for tensor $\mathcal{X}$;
Vecn($\bm{\mathcal{X}})$ |  $n$th vectorization of a tensor $\mathcal{X}$;
$\Omega$ | Index $(i_{1},\cdots,i_{n},\cdots,i_{N})$ of a tensor $\mathcal{X}$;
$\Omega^{(n)}_{M}$ | Index $(i_{n},j)$ of $n$th unfolding matrix $\textbf{X}^{(n)}$;
$(\Omega^{(n)}_{M})_{i}$ | Column index set in $i$th row of $\Omega^{(n)}_{M}$;
$(\Omega^{(n)}_{M})^{j}$ | Row index set in $j$th column of $\Omega^{(n)}_{M}$;
$\Omega^{(n)}_{V}$ | Index $i$ of $n$th unfolding vector Vecn($\bm{\mathcal{X}}$);
$\\{N\\}$ | The ordered set $\\{1,2,\cdots,N-1,N\\}$;
$\textbf{A}^{(n)}$ |  $n$th feature matrix $\in$ $\mathbb{R}^{I_{n}\times J_{n}}$;
$a_{i_{n},:}^{(n)}$ |  $i_{n}$th row vector $\in$ $\mathbb{R}^{K_{n}}$ of $\textbf{A}^{(n)}$;
$a_{:,j}^{(n)}$ |  $j$th column vector $\in$ $\mathbb{R}^{K_{n}}$ of $\textbf{A}^{(n)}$;
$a_{i_{n},k_{n}}^{(n)}$ |  $k_{n}$th element of feature vector $a_{i_{n},:}^{(n)}$;
$\cdot\big{/}(-,/)$ | Element-wise multiplication$\big{/}$ division;
$\circ$ | Outer production of vectors;
$\odot$ | Khatri-Rao (columnwise Kronecker) product;
$\times$ | Matrix product;
$\times_{(n)}$ |  $n$-Mode Tensor-Matrix product;
$\otimes$ | Kronecker product.

###  3.1 Basic Definitions

######  Definition 1 (Tensor Unfolding (Matricization)).

$n$th tensor unfolding (Matricization) refers to that a low order matrix
$\textbf{X}^{(n)}\in\mathbb{R}^{I_{n}\times I_{1}\cdots I_{n-1}\cdot
I_{n+1}\cdot\cdots\cdot I_{N}}$ stores all information of a tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots I_{n}\cdots\times
I_{N}}$ and the matrix element $x^{(n)}_{i_{n},j}$ of $\textbf{X}^{(n)}$ at
the position $j=1+\sum_{k=1,n\neq
k}^{N}\left[(i_{k}-1)\mathop{\prod}_{m=1,m\neq n}^{k-1}I_{m}\right]$ contains
the tensor element $x_{i_{1},i_{2},\cdots,i_{n}}$ of a tensor
$\bm{\mathcal{X}}$.

######  Definition 2 (Vectorization of a tensor $\bm{\mathcal{X}}$).

$n$th tensor vectorization refers to that a vector $x^{(n)}$
(Vecn($\bm{\mathcal{X}}$) and Vec($\textbf{X}^{(n)}$)) stores all elements in
the $n$th matricization $\textbf{X}^{(n)}$ of a tensor $\mathcal{X}$ and
$x^{(n)}_{k}$ $=$ $\textbf{X}^{(n)}_{i,j}$, where $k=(j-1)I_{n}+i$.

######  Definition 3 (Tensor Approximation).

A $N$-order tensor $\bm{\mathcal{X}}$ $\in$
$\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ can be approximated by
$\widehat{\bm{\mathcal{X}}}$ $\in$ $\mathbb{R}^{I_{1}\times\cdots\times
I_{N}}$, as well as a $N$-order residual or noisy tensor $\bm{\mathcal{E}}$
$\in$ $\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$. The low-rank approximation
problem is defined as
$\bm{\mathcal{X}}=\widehat{\bm{\mathcal{X}}}+\bm{\mathcal{E}}$, where
$\widehat{\bm{\mathcal{X}}}$ is denoted by a low-rank tensor.

Figure 1: Illustration of TD.

######  Definition 4 ($n$-Mode Tensor-Matrix product).

$n$-Mode Tensor-Matrix product is an operation which can reflect coordinate
projection of a tensor $\bm{\mathcal{X}}$ $\in$
$\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ with projection matrix U $\in$
$\mathbb{R}^{I_{n}\times J_{n}}$ into a tensor
$(\bm{\mathcal{X}}\times_{(n)}\textbf{U})$ $\in$
$\mathbb{R}^{I_{1}\times\cdots\times I_{n-1}\times J_{n}\times\cdots I_{N}}$
where $(\bm{\mathcal{X}}\times_{(n)}\textbf{U})_{i_{1}\times\cdots\times
i_{n-1}\times j_{n}\times i_{n+1}\times\cdots\times i_{N}}$ $=$
$\sum\limits_{i_{n}=1}^{I_{n}}$ $x_{i_{1}\times\cdots\times
i_{n}\times\cdots\times i_{N}}u_{j_{n},i_{n}}$.

######  Definition 5 ($R$ Kruskal Product).

For an $N$-order tensor $\widehat{\bm{\mathcal{X}}}$ $\in$
$\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$, the $R$ Kruskal Product of
$\widehat{\bm{\mathcal{X}}}$ is given by $R$ Kruskal product as:
$\widehat{\bm{\mathcal{X}}}=\sum_{r=1}^{R}a^{(1)}_{:,r}\circ\cdots\circ
a^{(n)}_{:,r}\circ\cdots\circ a^{(N)}_{:,r}$.

######  Definition 6 (Sparse Tucker Decomposition (STD)).

For a $N$-order sparse tensor $\bm{\mathcal{X}}$ $\in$
$\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$, the STD of the optimal
approximated $\widehat{\bm{\mathcal{X}}}$ is given by
$\widehat{\bm{\mathcal{X}}}=\mathcal{G}\times_{(1)}\textbf{A}^{(1)}\times_{(2)}\cdots\times_{(n)}\textbf{A}^{(n)}\times_{(n+1)}\cdots\times_{(N)}\textbf{A}^{(N)},$
where $\mathcal{G}$ is the core tensor and $\textbf{A}^{(n)}$, $n\in\\{N\\}$
are the low-rank factor matrices. The rank of TF for a tensor is
$[J_{1},\cdots,J_{n},\cdots,J_{N}]$. The determination process for the core
tensor $\mathcal{G}$ and factor matrices $\textbf{A}^{(n)}$, $n\in\\{N\\}$
follows the sparsity model of the sparse tensor $\bm{\mathcal{X}}$.

In the convex optimization community, the literature [49] gives the definition
of Lipschitz-continuity with constant $L$ and strong convexity with constant
$\mu$.

######  Definition 7 ($L$-Lipschitz continuity).

A continuously differentiable function $f(\textbf{x})$ is called $L$-smooth on
$\mathbb{R}^{r}$ if the gradient $\nabla f(\textbf{x})$ is $L$-Lipschitz
continuous for any x, y $\in$ $\mathbb{R}^{r}$, that is $\|$ $\nabla
f(\textbf{x})$ $-$ $\nabla f(\textbf{y})$ $\|_{2}$ $\leq$ $L$ $\|$ x $-$ y
$\|_{2}$, where $\|\bullet\|_{2}$ is $L_{2}$-norm $\|\textbf{x}\|_{2}$ $=$
$(\mathbb{\sum}_{k=1}^{r}x_{k}^{2})^{1/2}$ for a vector x.

######  Definition 8 ($\mu$-Convex).

A continuously differentiable function $f(\textbf{x})$ is called strongly-
convex on $\mathbb{R}^{r}$ if there exists a constant $\mu$ $>$ $0$ for any x,
y $\in$ $\mathbb{R}^{r}$, that is $f(\textbf{x})$ $\geq$ $f(\textbf{y})$ $+$
$\nabla$ $f(\textbf{y})$ $(\textbf{x}-\textbf{y})^{T}$ $+$
$\frac{1}{2}\mu\|\textbf{x}-\textbf{y}\|_{2}^{2}$.

Due to limited spaces, we provide the description of basic optimization as the
supplementary material.

###  3.2 Problems for Large-scale Optimization Problems

Many ML tasks are transformed into the solvent of optimization problem [32,
35, 33, 34] and the basis optimization problem is organized as:

$\displaystyle\mathop{\arg\min}_{w\in\mathbb{R}^{R}}f(w)$
$\displaystyle=\underbrace{L\bigg{(}w\bigg{|}y_{i},x_{i},w\bigg{)}}_{Loss\
Function}+\underbrace{\lambda_{w}R(w)}_{Regularization\ Item}$ (1)
$\displaystyle=\sum\limits_{i=1}^{N}L_{i}\bigg{(}w\bigg{|}y_{i},x_{i},w\bigg{)}+\lambda_{w}R_{i}(w),$

where $y_{i}$ $\in$ $\mathbb{R}^{1}$, $x_{i}$ $\in$ $\mathbb{R}^{R}$,
$i\in\\{N\\}$, $w\in\mathbb{R}^{R}$ and the original optimization model needs
gradient which should select all the samples $\\{x_{i}|i\in\\{N\\}\\}$ from
the dataset $\Omega$ and the GD is presented as:

$\displaystyle w$ $\displaystyle\leftarrow w-\gamma\frac{\partial
f_{\Omega}(w)}{\partial w}$ (2)
$\displaystyle=w-\gamma\frac{1}{N}\sum\limits_{i=1}^{N}\frac{\partial\bigg{(}L_{i}(w)+\lambda_{w}R_{i}(w)\bigg{)}}{\partial
w}.$

The second-order solution, i.e., ALS and CD, etc, are deduced from the
construction of GD from the whole dataset. In large-scale optimization
scenarios, SGD is a common strategy [32, 35, 33, 34] and promises to obtain
the optimal accuracy via a certain number of training epoches [32, 35, 33,
34]. An $M$ entries set $\Psi$ is randomly selected from the set $\Omega$, and
the SGD is presented as:

$\displaystyle w$ $\displaystyle\leftarrow w-\gamma\frac{\partial
f_{\Psi}(w)}{\partial w}$ (3) $\displaystyle\approx
w-\gamma\frac{1}{M}\sum\limits_{i\in\Psi}\frac{\partial\bigg{(}L_{i}(w)+\lambda_{w}R_{i}(w)\bigg{)}}{\partial
w}.$

Equ. (3) is an average SGD [33], and the averaged SG can be applied to build
the basic tool of the modern stochastic optimization strategies, e.g.,
Stochastic Recursive Gradient [34], variance SGD [32], or momentum SGD [35],
which can retain robustness and fast convergence. The optimization function
can be packaged in the form of $SGD\bigg{(}M,\lambda,\gamma,w,\frac{\partial
f_{\Psi}(w)}{\partial w}\bigg{)}$.

Figure 2: Illustration of Optimization Process for SGD$\\_$Tucker: 1) Original
problem for core tensor, 2) Kruskal product for approximating core tensor, 3)
Optimization process for factor matrices.

##  4 SGD$\\_$Tucker

Overview: In this section, SGD$\\_$Tucker is proposed to decompose the
optimization objectives which involves frequent operations of intermediate
matrices into a problem which only needs the manipulations of small batches of
intermediate matrices. Fig. 2 illustrates the framework for our work. As Fig.
2 shows, SGD$\\_$Tucker comprises of an approximation process of the core
tensor and an optimization for factor matrices. Table II records the
intermediate variables which follow the problem deduction process. Section 4.1
presents the deduction process for the core tensor (Lines 1-16, Algorithm 1),
and Section 4.2 presents the proceeding process for the factor matrices of
SGD$\\_$Tucker (Lines 17-26, Algorithm 1). Section 4.3 shows the stochastic
optimization process for the proposed novel optimization objectives, which
illustrates the details of the process that how can the SGD$\\_$Tucker divide
the high-dimension intermediate matrices
$\bigg{\\{}\textbf{H}^{(n)},\textbf{S}^{(n)},\textbf{E}^{(n)}\big{|}n\in\\{N\\}\bigg{\\}}$
into a small batches of intermediate matrices. Section 4.4 shows the parallel
and distributed strategies for SGD$\\_$Tucker. Finally, section 4.5
illustrates the analysis of space and time complexities. Because the
optimization process for approximating the core tensor and factor matrices are
non-convex. We fix the core tensor and optimize a factor matrix and fix factor
matrices to optimize Kruskal product for core tensor. This alternative
minimization strategy is a convex solution [46, 47, 48, 50].

TABLE II: Table of Intermediate Varibles for SGD$\\_$Tucker. Symbol | Description | Emerging
---|---|---
|  | Equations
$\widehat{\bm{\mathcal{G}}}$ $\in$ $\mathbb{R}^{J_{1}\times J_{2}\times\cdots\times J_{N}}_{+}$ |  $\mathbb{R}_{core}$ Kruskal | Equ. (4);
| product for $\bm{\mathcal{G}}$ |
$\textbf{B}^{(n)}\in\mathbb{R}^{J_{n}\times R_{core}}$ |  $n$th Kruskal matrix for $\widehat{\bm{\mathcal{G}}}$ | Equ. (4);
$\widehat{g}^{(n)}$ $\in$ $\mathbb{R}^{\prod\limits_{n=1}^{N}J_{n}}$ |  $n$th vectorization | Equ. (5);
| of tensor $\widehat{\bm{\mathcal{G}}}$ |
$\textbf{H}^{(n)}$ $\in$ $\mathbb{R}^{\prod\limits_{n=1}^{N}I_{n}\times\prod\limits_{n=1}^{N}J_{n}}$ | Coefficient for $\widehat{g}^{(n)}$ | Equ. (5);
$x^{(n)}$ $\in$ $\mathbb{R}^{\prod\limits_{n=1}^{N}I_{n}}$ |  $n$th vectorization of | Equ. (5);
| Tensor $\mathcal{X}$ |
$\textbf{Q}^{(n)}$ $\in$ $\mathbb{R}^{\prod\limits_{r=1,r\neq n}^{N}J_{r}\times R_{core}}$ | Coefficient of $\textbf{B}^{(n)}$ | Equ. (6);
| for constructing $\widehat{g}^{(n)}$ |
$\textbf{U}^{(n)}$ $\in$ $\mathbb{R}^{J_{n}\times J_{n}}$ | Unity matrix | Equ. (7);
$\textbf{O}^{(n)}_{r}$ $\in$ $\mathbb{R}^{\prod\limits_{r=1}^{N}J_{r}\times J_{n}}$ | Coefficient of $b^{(n)}_{:,r}$ | Equ. (8);
| for constructing $\widehat{g}^{(n)}$ |
$x^{(n)}_{r_{core}}$ $\in$ $\mathbb{R}^{\prod\limits_{r=1}^{N}I_{r}}$ | Intermediate vector of $x^{(n)}$ | Equ. (9);
$\textbf{H}^{(n)}_{r_{core}}$ $\in$ $\mathbb{R}^{\prod\limits_{n=1}^{N}I_{n}\times J_{n}}$ | Coefficient of $b^{(n)}_{:,r}$ | Equ. (9);
| for constructing $\widehat{g}^{(n)}$ |
$\textbf{S}^{(n)}$ $\in$ $\mathbb{R}^{\prod\limits_{k=1,k\neq n}^{N}J_{n}\times\prod\limits_{k=1,k\neq n}^{N}I_{k}}$ | Coefficient of $\textbf{A}^{(n)}\widehat{\textbf{G}}^{(n)}$ | Equ. (10);
| for constructing $\widehat{\bm{X}}^{(n)}$ |
$\textbf{E}^{(n)}$ $\in$ $\mathbb{R}^{J_{n}\times\prod\limits_{k=1,k\neq n}^{N}I_{k}}$ | Coefficient of $\textbf{A}^{(n)}$ | Equ. (11);
| for constructing $\widehat{\bm{X}}^{(n)}$ |
$x^{(n)}_{\Psi^{(n)}_{V}}$ $\in$ $\mathbb{R}^{|\Psi^{(n)}_{V}|}$ | Partial vector $x^{(n)}$ from set $\Psi$ | Equ. (14);
$\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)}$ $\in$ $\mathbb{R}^{|\Psi^{(n)}_{V}|\times\prod\limits_{n=1}^{N}J_{n}}$ | Partial matrix of $\textbf{H}^{(n)}$ | Equ. (14);
$\bm{X}_{:,(\Psi^{(n)}_{M})_{i_{n}}}^{(n)}$ $\in$ $\mathbb{R}^{I_{n}\times|(\Psi^{(n)}_{M})_{i_{n}}|}$ | Partial matrix of $\bm{X}^{(n)}$ | Equ. (17);
$\textbf{E}^{(n)}_{:,(\Psi^{(n)}_{M})_{i_{n}}}$ $\in$ $\mathbb{R}^{J_{n}\times|(\Psi^{(n)}_{M})_{i_{n}}|}$ | Partial matrix of $\textbf{E}^{(n)}$ | Equ. (17).

###  4.1 Optimization Process for the Core Tensor

Due to the non-trivial coupling of the core tensor $\bm{\mathcal{G}}$, the
efficient and effective way to infer it is through an approximation. A tensor
$\bm{\mathcal{G}}$ can be approximated by $R_{core}$ $\leq$
${J_{n},n\in\\{N\\}}$ Kruskal product of low-rank matrices
$\\{\textbf{B}^{(n)}\in\mathbb{R}^{J_{n}\times R_{core}}|n\in\\{N\\}\\}$ to
form $\widehat{\bm{\mathcal{G}}}$:

$\displaystyle\widehat{\bm{\mathcal{G}}}=\sum_{r_{core}=1}^{R_{core}}b^{(1)}_{:,r_{core}}\circ\cdots\circ
b^{(n)}_{:,r_{core}}\circ\cdots\circ b^{(N)}_{:,r_{core}}.$ (4)

As the direct approximation for the core tensor $\bm{\mathcal{G}}$ may result
in instability, we propose to apply the coupling process to approximate STD
and tensor approximation. Specifically, we use Kruskal product of low-rank
matrices$\\{\textbf{B}^{(n)}\in\mathbb{R}^{J_{n}\times
R_{core}}|n\in\\{N\\}\\}$ as follows:

$\displaystyle\mathop{\arg\min}_{\widehat{\bm{\mathcal{G}}}}f$
$\displaystyle\bigg{(}\widehat{g}^{(n)}\bigg{|}x^{(n)},\\{\textbf{A}^{(n)}\\}\bigg{)}$
(5) $\displaystyle=$
$\displaystyle\bigg{\|}x^{(n)}-\textbf{H}^{(n)}\widehat{g}^{(n)}\bigg{\|}_{2}^{2}+\lambda_{\widehat{g}^{(n)}}\bigg{\|}\widehat{g}^{(n)}\bigg{\|}_{2}^{2},$

where $\widehat{g}^{(n)}$ $=$ $Vec(\textbf{B}^{(n)}\textbf{Q}^{(n)^{T}})$ and
$\textbf{Q}^{(n)}$ $=$
$\textbf{B}^{(N)}\odot\cdots\odot\textbf{B}^{(n+1)}\odot\textbf{B}^{(n-1)}\odot\cdots\odot\textbf{B}^{(1)}$.
The tanglement problem resulted from $R_{core}$ Kruskal product of low-rank
matrices $\\{\textbf{B}^{(n)}\in\mathbb{R}^{J_{n}\times
R_{core}}|n\in\\{N\\}\\}$ leads to a non-convex optimization problem for the
optimization objective (5). The alternative optimization strategy [50] is
adopted to update the parameters
$\\{\textbf{B}^{(n)}\in\mathbb{R}^{J_{n}\times R_{core}}|n\in\\{N\\}\\}$ and
then the non-convex optimization objective (5) is turned into the objective
as:

$\displaystyle\mathop{\arg\min}_{\textbf{B}^{(n)},n\in\\{N\\}}$ $\displaystyle
f\bigg{(}\textbf{B}^{(n)}\bigg{|}x^{(n)},\\{\textbf{A}^{(n)}\\},\\{\textbf{B}^{(n)}\\}\bigg{)}$
(6) $\displaystyle=$
$\displaystyle\bigg{\|}x^{(n)}-\textbf{H}^{(n)}Vec(\textbf{B}^{(n)}\textbf{Q}^{(n)^{T}})\bigg{\|}_{2}^{2}+\lambda_{\textbf{B}}\bigg{\|}\textbf{B}\bigg{\|}_{2}^{2},$

where $Vec(\textbf{B}^{(n)}\textbf{Q}^{(n)^{T}})$ $=$
$Vec(\sum\limits_{r=1}^{R_{core}}b^{(n)}_{:,r}q_{:,r}^{(n)^{T}})$.

The optimization objective (6) is not explicit for the variable
$\textbf{B}^{(n)},n\in\\{N\\}$ and it is hard to find the gradient of the
variables $\textbf{B}^{(n)},n\in\\{N\\}$ under the current formulation. Thus,
we borrow the intermediate variables with a unity matrix as
$\textbf{O}^{(n)}_{r}$ $\in$ $\mathbb{R}^{\prod_{k=1}^{N}J_{k}\times J_{n}}$,
$r$ $\in$ $\\{R_{core}\\}$ and the unity matrix $\textbf{U}^{(n)}$ $\in$
$\mathbb{R}^{J_{n}\times J_{n}}$, $n\in\\{N\\}$ which can present the
variables $\textbf{B}^{(n)},n\in\\{N\\}$ in a gradient-explicit formation. The
matrix $\textbf{O}^{(n)}_{r}$ is defined as

$\displaystyle\textbf{O}^{(n)}_{r}=\bigg{[}$ $\displaystyle
q_{1,r}\textbf{U}^{(n)},\cdots,q_{m,r}\textbf{U}^{(n)},\cdots,$ (7)
$\displaystyle q_{\prod\limits_{k=1,k\neq
n}^{N}J_{k},r}\textbf{U}^{(n)}\bigg{]}^{T}.$

The key transformation of the gradient-explicit formation for the variables
$\textbf{B}^{(n)},n\in\\{N\\}$ is as
$Vec(\sum\limits_{r=1}^{R_{core}}b^{(n)}_{:,r}q_{:,r}^{(n)^{T}})=\sum\limits_{r=1}^{R_{core}}\textbf{O}^{(n)}_{r}b^{(n)}_{:,r}$.
Then the optimization objective (6) is reformulated into:

$\displaystyle\mathop{\arg\min}_{\textbf{B}^{(n)},n\in\\{N\\}}$ $\displaystyle
f\bigg{(}\textbf{B}^{(n)}\bigg{|}x^{(n)},\\{\textbf{A}^{(n)}\\},\\{\textbf{B}^{(n)}\\}\bigg{)}$
(8)
$\displaystyle=\bigg{\|}x^{(n)}-\textbf{H}^{(n)}\sum\limits_{r=1}^{R_{core}}\textbf{O}^{(n)}_{r}b^{(n)}_{:,r}\bigg{\|}_{2}^{2}+\lambda_{\textbf{B}}\bigg{\|}\textbf{B}\bigg{\|}_{2}^{2}.$

The cyclic block optimization strategy [51] is adopted to update the variables
$\big{\\{}b^{(n)}_{:,r_{core}}|r_{core}\in\\{R_{core}\\}\big{\\}}$ within a
low-rank matrix $\textbf{B}^{(n)},n\in\\{N\\}$ and eventually the optimization
objective is reformulated into:

$\displaystyle\mathop{\arg\min}_{b^{(n)}_{:,r_{core}},n\in\\{N\\}}$
$\displaystyle
f\bigg{(}b^{(n)}_{:,r_{core}}\bigg{|}x^{(n)},\\{\textbf{A}^{(n)}\\},\\{\textbf{B}^{(n)}\\}\bigg{)}$
(9) $\displaystyle=$
$\displaystyle\bigg{\|}x^{(n)}_{r_{core}}-\textbf{H}^{(n)}_{r_{core}}b^{(n)}_{:,r_{core}}\bigg{\|}_{2}^{2}+\lambda_{\textbf{B}}\bigg{\|}b^{(n)}_{:,r_{core}}\bigg{\|}_{2}^{2},$

where $x^{(n)}_{r_{core}}$ $=$ $x^{(n)}-\textbf{H}^{(n)}\sum\limits_{r=1,r\neq
r_{core}}^{R_{core}}\textbf{O}^{(n)}_{r}b^{(n)}_{:,r}$ and
$\textbf{H}^{(n)}_{r_{core}}$ $=$
$\textbf{H}^{(n)}\textbf{O}^{(n)}_{r_{core}}\in\mathbb{R}^{\prod\limits_{n=1}^{N}I_{n}\times
J_{n}}$.

######  Theorem 1.

From the function form of (9), the optimization objective for
$b^{(n)}_{:,r_{core}}$ is a $u$-convex and $L$-smooth function.

###### Proof.

The proof is divided into the following two parts: 1) The non-convex
optimization problem is transformed into fixing factor matrices
$\textbf{A}^{(n)},n\in\\{N\\}$ and updating core tensor part, and this
transformation can keep promise that convex optimization for each part and
appropriate accuracy [46, 47, 48, 50]. 2) The distance function
$\big{\|}x^{(n)}_{r_{core}}-\textbf{H}^{(n)}_{r_{core}}b^{(n)}_{:,r_{core}}\big{\|}_{2}^{2}$
is an Euclidean distance with $L_{2}$ norm regularization
$\lambda_{\textbf{B}}\big{\|}b^{(n)}_{:,r_{core}}\big{\|}_{2}^{2}$ [52]. Thus,
the optimization objective (9) is a $u$-convex and $L$-smooth function
obviously [50], [53], [49]. The proof details are omitted which can be found
on [32], [54]. ∎

###  4.2 Optimization Process for Factor Matrices

After the optimization step is conducted for
$\\{\textbf{B}^{(n)}\in\mathbb{R}^{J_{n}\times R_{core}}|n\in\\{N\\}\\}$, the
Kruskal product for approximated core tensor $\hat{g}^{(n)}$ is constructed by
Equ. (4). Thus, we should consider the optimization problem for factor
matrices $\textbf{A}^{(n)},\\{n|n\in\\{N\\}\\}$ as:

$\displaystyle\mathop{\arg\min}_{\textbf{A}^{(n)},n\in\\{N\\}}$ $\displaystyle
f\bigg{(}\textbf{A}^{(n)}\bigg{|}\bm{X}^{(n)},\big{\\{}\textbf{A}^{(n)}\big{\\}},\widehat{\textbf{G}}^{(n)}\bigg{)}$
(10)
$\displaystyle=\bigg{\|}\bm{X}^{(n)}-\widehat{\bm{X}}^{(n)}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}\textbf{A}^{(n)}\bigg{\|}_{2}^{2},$

where $\widehat{\bm{X}}^{(n)}$ $=$
$\textbf{A}^{(n)}\widehat{\textbf{G}}^{(n)}\textbf{S}^{(n)^{T}}$ and
$\textbf{S}^{(n)}$ $=$
$\textbf{A}^{(N)}\otimes\cdots\otimes\textbf{A}^{(n+1)}\otimes\textbf{A}^{(n-1)}\otimes\cdots\otimes\textbf{A}^{(1)}$.

The optimization process for the whole factor matrix set
$\\{\textbf{A}^{(1)},\cdots,\textbf{A}^{(N)}\\}$ is non-convex. The
alternative optimization strategy [50] is adopted to transform the non-convex
optimization problem into $N$ convex optimization problems under a fine-tuned
initialization as:

$\displaystyle\left\\{\begin{aligned}
\mathop{\arg\min}_{\textbf{A}^{(1)}}f\bigg{(}\textbf{A}^{(1)}\bigg{)}=&\bigg{\|}\bm{X}^{(1)}-\textbf{A}^{(1)}\textbf{E}^{(1)}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}\textbf{A}^{(1)}\bigg{\|}_{2}^{2};\\\
\vdots\\\
\mathop{\arg\min}_{\textbf{A}^{(n)}}f\bigg{(}\textbf{A}^{(n)}\bigg{)}=&\bigg{\|}\bm{X}^{(n)}-\textbf{A}^{(n)}\textbf{E}^{(n)}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}\textbf{A}^{(n)}\bigg{\|}_{2}^{2};\\\
\vdots\\\
\mathop{\arg\min}_{\textbf{A}^{(N)}}f\bigg{(}\textbf{A}^{(N)}\bigg{)}=&\bigg{\|}\bm{X}^{(N)}-\textbf{A}^{(N)}\textbf{E}^{(N)}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}\textbf{A}^{(N)}\bigg{\|}_{2}^{2},\end{aligned}\right.$
(11)

where $\textbf{E}^{(n)}$ $=$ $\widehat{\textbf{G}}^{(n)}\textbf{S}^{(n)^{T}}$
$\in$ $\mathbb{R}^{J_{n}\times\prod\limits_{k=1,k\neq n}^{N}I_{k}}$. The
optimization objectives of the $n$th variables are presented as an independent
form as:

$\displaystyle\left\\{\begin{aligned}
\mathop{\arg\min}_{a_{1,:}^{(n)}}f(a_{1,:}^{(n)})=&\bigg{\|}\textbf{X}_{1,:}^{(n)}-a_{1,:}^{(n)}\textbf{E}^{(n)}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}a_{1,:}^{(n)}\bigg{\|}_{2}^{2};\\\
\vdots\\\
\mathop{\arg\min}_{a_{i_{n},:}^{(n)}}f(a_{i_{n},:}^{(n)})=&\bigg{\|}\textbf{X}_{i_{n},:}^{(n)}-a_{i_{n},:}^{(n)}\textbf{E}^{(n)}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}a_{i_{n},:}^{(n)}\bigg{\|}_{2}^{2};\\\
\vdots\\\
\mathop{\arg\min}_{a_{I_{n},:}^{(n)}}f(a_{I_{n},:}^{(n)})=&\bigg{\|}\textbf{X}_{I_{n},:}^{(n)}-a_{I_{n},:}^{(n)}\textbf{E}^{(n)}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}a_{I_{n},:}^{(n)}\bigg{\|}_{2}^{2}.\end{aligned}\right.$
(12)

######  Theorem 2.

From the function form of (12), the optimization objective for
$\textbf{A}^{(n)}$ is a $u$-convex and $L$-smooth function.

###### Proof.

By adopting the alternative strategy [46, 47, 48, 50], we fix
$\widehat{\bm{\mathcal{G}}}$ and $\textbf{A}^{(k)},k\neq n,k\in\\{N\\}$. Then,
we update $\textbf{A}^{(n)}$. The distance function
$\bigg{\|}\textbf{X}_{i_{n},:}^{(n)}-a_{i_{n},:}^{(n)}\textbf{E}^{(n)}\bigg{\|}_{2}^{2}$
is an Euclidean distance with $L_{2}$ norm regularization
$\lambda_{\textbf{A}}\bigg{\|}a_{i_{n},:}^{(n)}\bigg{\|}_{2}^{2}$,
$i_{n}\in{I_{n}}$ [52]. Thus, the optimization objective (9) is a $u$-convex
and $L$-smooth function obviously [50], [53], [49]. Due to the limited space,
the proof details are omitted which can be found on [32], [54]. ∎

###  4.3 Stochastic Optimization

The previous Sections 4.1 and 4.2 presented the transformation of the
optimization problem. In this section, the solvent is introduced. The average
SGD method is a basic part of state of the art stochastic optimization models.
Thus, in this section, we present the average SGD for our optimization
objectives. The optimization objectives for the core tensor are presented in
Equs. (8) and (9). The optimization objectives for the factor matrix are
presented in Equ. (12) which are in the form of a basic optimization model. In
the industrial and big-data communities, the HOHDST is very common. Thus, SGD
is proposed to replace the original optimization strategy.

The solution for $b^{(n)}_{:,r_{core}}$ is presented as:

$\displaystyle\mathop{\arg\min}_{b^{(n)}_{:,r_{core}},n\in\\{N\\}}$
$\displaystyle
f\bigg{(}b^{(n)}_{:,r_{core}}\bigg{|}x^{(n)},\\{\textbf{A}^{(n)}\\},\\{\textbf{B}^{(n)}\\}\bigg{)}$
(13) $\displaystyle=$
$\displaystyle\sum\limits_{i\in\Omega^{(n)}_{V}}L_{i}\bigg{(}b^{(n)}_{:,r_{core}}\bigg{|}x^{(n)}\bigg{)}+\lambda_{\textbf{B}}\bigg{\|}b^{(n)}_{:,r_{core}}\bigg{\|}_{2}^{2}.$

If a set $\Psi$ including $M$ randomly entries is selected, the approximated
SGD solution for $b^{(n)}_{:,r_{core}}$ is presented as:

$\displaystyle\mathop{\arg\min}_{b^{(n)}_{:,r_{core}},n\in\\{N\\}}f_{\Psi^{(n)}_{V}}\bigg{(}b^{(n)}_{:,r_{core}}\bigg{|}x^{(n)}_{\Psi^{(n)}_{V}},\\{\textbf{A}^{(n)}\\},\\{\textbf{B}^{(n)}\\}\bigg{)}$
(14)
$\displaystyle=\sum\limits_{i\in\Psi^{(n)}_{V}}L_{i}\bigg{(}b^{(n)}_{:,r_{core}}\bigg{|}x^{(n)}\bigg{)}+\lambda_{\textbf{B}}\bigg{\|}b^{(n)}_{:,r_{core}}\bigg{\|}_{2}^{2}$
$\displaystyle=\bigg{\|}x^{(n)}_{\Psi^{(n)}_{V}}-\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)}\sum\limits_{r=1}^{R_{core}}\textbf{O}^{(n)}_{r}b^{(n)}_{:,r}\bigg{\|}_{2}^{2}+\lambda_{\textbf{B}}\bigg{\|}b^{(n)}_{:,r_{core}}\bigg{\|}_{2}^{2}$
$\displaystyle=\bigg{\|}(x^{(n)}_{r_{core}})_{\Psi^{(n)}_{V}}-\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)}\textbf{O}^{(n)}_{r_{core}}b^{(n)}_{:,r_{core}}\bigg{\|}_{2}^{2}+\lambda_{\textbf{B}}\bigg{\|}b^{(n)}_{:,r_{core}}\bigg{\|}_{2}^{2}.$

The SGD for the approximated function
$f_{\Psi^{(n)}_{V}}\bigg{(}b^{(n)}_{:,r_{core}}\bigg{|}x^{(n)}_{\Psi^{(n)}_{V}},\\{\textbf{A}^{(n)}\\},\\{\textbf{B}^{(n)}\\}\bigg{)}$
is deduced as:

$\displaystyle\frac{\partial
f_{\Psi^{(n)}_{V}}\bigg{(}b^{(n)}_{:,r_{core}}\bigg{|}x^{(n)}_{\Psi^{(n)}_{V}},\\{\textbf{A}^{(n)}\\},\\{\textbf{B}^{(n)}\\}\bigg{)}}{\partial
b^{(n)}_{:,r_{core}}}$ (15) $\displaystyle=$
$\displaystyle\frac{1}{M}\bigg{(}-\textbf{O}^{(n)T}_{r_{core}}\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)T}(x^{(n)}_{r_{core}})_{\Psi^{(n)}_{V}}$
$\displaystyle+\textbf{O}^{(n)T}_{r_{core}}\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)T}\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)}\textbf{O}^{(n)}_{r_{core}}b^{(n)}_{:,r_{core}}\bigg{)}+\lambda_{\textbf{B}}b^{(n)}_{:,r_{core}}.$

The solution for factor matrices $a_{i_{n},:}^{(n)}$, $i_{n}$ $\in$
$\\{I_{N}\\}$, $n$ $\in$ $\\{N\\}$ is presented as:

$\displaystyle\mathop{\arg\min}_{a_{i_{n},:}^{(n)},n\in\\{N\\}}$
$\displaystyle
f\bigg{(}a_{i_{n},:}^{(n)}\bigg{|}\bm{X}_{i_{n},:}^{(n)},\big{\\{}\textbf{A}^{(n)}\big{\\}},\widehat{\textbf{G}}^{(n)}\bigg{)}$
(16) $\displaystyle=$
$\displaystyle\sum\limits_{j\in(\Omega^{(n)}_{M})_{i_{n}}}L_{j}\bigg{(}a_{i_{n},:}^{(n)}\bigg{|}\bm{X}_{i_{n},j}^{(n)}\bigg{)}+\lambda_{\textbf{A}}\bigg{\|}a_{i_{n},:}^{(n)}\bigg{\|}_{2}^{2}.$

If a set $\Psi$ including $M$ randomly chosen entries is selected, the SGD
solution for $a_{i_{n},:}^{(n)}$ can be expressed as:

$\displaystyle\mathop{\arg\min}_{a_{i_{n},:}^{(n)},n\in\\{N\\}}f_{\Psi^{(n)}_{M}}\bigg{(}a_{i_{n},:}^{(n)}\bigg{|}\bm{X}_{i_{n},(\Psi^{(n)}_{M})_{i_{n}}}^{(n)},\big{\\{}\textbf{A}^{(n)}\big{\\}},\widehat{\textbf{G}}^{(n)}\bigg{)}$
(17)
$\displaystyle=\sum\limits_{j\in(\Psi^{(n)}_{M})_{i_{n}}}L_{j}\bigg{(}a_{i_{n},:}^{(n)}\bigg{|}\bm{X}_{i_{n},j}^{(n)}\bigg{)}+\lambda_{\textbf{A}}\bigg{\|}a_{i_{n},:}^{(n)}\bigg{\|}_{2}^{2}$
$\displaystyle=\bigg{\|}\textbf{X}_{i_{n},(\Psi^{(n)}_{M})_{i_{n}}}^{(n)}-a_{i_{n},:}^{(n)}\textbf{E}^{(n)}_{:,(\Psi^{(n)}_{M})_{i_{n}}}\bigg{\|}_{2}^{2}+\lambda_{\textbf{A}}\bigg{\|}a_{i_{n},:}^{(n)}\bigg{\|}_{2}^{2}.$

The stochastic gradient for the approximated function
$f_{\Psi^{(n)}_{V}}\bigg{(}a_{i_{n},:}^{(n)}\bigg{|}\bm{X}_{i_{n},(\Psi^{(n)}_{M})_{i_{n}}}^{(n)},\big{\\{}\textbf{A}^{(n)}\big{\\}},\textbf{G}^{(n)}\bigg{)}$
is deduced as:

$\displaystyle\frac{\partial
f_{\Psi^{(n)}_{V}}\bigg{(}a_{i_{n},:}^{(n)}\bigg{|}\bm{X}_{i_{n},(\Psi^{(n)}_{M})_{i_{n}}}^{(n)},\big{\\{}\textbf{A}^{(n)}\big{\\}},\widehat{\textbf{G}}^{(n)}\bigg{)}}{\partial
a_{i_{n},:}^{(n)}}$ (18) $\displaystyle=$
$\displaystyle\frac{1}{M}\bigg{(}-\textbf{X}_{i_{n},(\Psi^{(n)}_{M})_{i_{n}}}^{(n)}\textbf{E}^{(n)T}_{:,(\Psi^{(n)}_{M})_{i_{n}}}$
$\displaystyle+a_{i_{n},:}^{(n)}\textbf{E}^{(n)}_{:,(\Psi^{(n)}_{M})_{i_{n}}}\textbf{E}^{(n)T}_{:,(\Psi^{(n)}_{M})_{i_{n}}}\bigg{)}+\lambda_{\textbf{A}}a_{i_{n},:}^{(n)}.$

The computational details are presented in Algorithm 1, which shows that
SGD$\\_$Tucker is able to divide the high-dimension intermediate matrices
$\bigg{\\{}\textbf{H}_{\Omega^{(n)}_{V},:}^{(n)},\textbf{S}^{(n)}_{(\Omega^{(n)}_{M})_{i_{n}},:},\textbf{E}^{(n)}_{:,(\Omega^{(n)}_{M})_{i_{n}}}\big{|}i_{n}\in\\{I_{n}\\},n\in\\{N\\}\bigg{\\}}$
into small batches of intermediate matrices
$\bigg{\\{}\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)},\textbf{S}^{(n)}_{(\Psi^{(n)}_{M})_{i_{n}},:},\textbf{E}^{(n)}_{:,(\Psi^{(n)}_{M})_{i_{n}}}\big{|}i_{n}\in\\{I_{n}\\},n\in\\{N\\}\bigg{\\}}$.
We summarized all steps of SGD$\\_$Tucker in Algorithm 1.

######  Theorem 3.

From the function forms of (9) and (12), the optimization objectives for the
core tensor and factor matrics are both $u$-convex and $L$-smooth functions.
The stochastic update rules of (15) and (18) can ensure the convergency of
alternative optimization objectives (9) and (12), respectively.

###### Proof.

The two function forms of (9) and (12) can be conclude as
$f(x)=\frac{1}{n}\mathop{\sum}_{i=1}^{n}f_{i}(\textbf{x})$, where
$f(\textbf{x})$ is a strongly-convex with constant $\mu$, and each
$f_{i}(\textbf{x})$ is smooth and Lipschitz-continuous with constant $L$. At
$t$th iteration, for chosen $f_{i_{t}}$ randomly, and a learning rate sequence
$\gamma_{t}>0$, the expectation $\mathbb{E}[\nabla
f_{i}(\textbf{x}_{t})|\textbf{x}_{t}]$ of $\nabla f_{i_{t}}(\textbf{x}_{t})$
is equivalent to $\nabla f(\textbf{x}_{t})$ [53, 49, 54, 55]. ∎

Algorithm 1 Stochastic Optimization Strategies on a Training Epoch for
SGD$\\_$Tucker.

Input: Sparse tensor $\mathcal{X}$, randomly selected set $\Psi$ with $M$
entries, Learning step $\gamma_{\textbf{A}}$, learning step
$\gamma_{\textbf{B}}$;
Initializing $\textbf{B}^{(n)}$, $n$ $\in$ $\\{N\\}$, core tensor
$\widehat{\mathcal{G}}$, $\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$,
$\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)}$ $\in$
$\mathbb{R}^{M\times\prod_{n=1}^{N}J_{n}}$, $r$ $\in$ $\\{R_{core}\\}$, $n$
$\in$ $\\{N\\}$, $\textbf{O}^{(n)}_{r}$ $\in$
$\mathbb{R}^{\prod_{k=1}^{N}J_{k}\times J_{n}}$, $r$ $\in$ $\\{R_{core}\\}$,
$n$ $\in$ $\\{N\\}$, $\widehat{x}^{(n)}_{\Psi^{(n)}_{V}}\in\mathbb{R}^{M}$.
Output: $\textbf{B}^{(n)}$, $n$ $\in$ $\\{N\\}$, $\widehat{\mathcal{G}}$,
$\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$.

1: for $n$ from $1$ to $N$ do

2: Obtain $\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)}$;

3: for $r_{core}$ from $1$ to $R_{core}$ do

4: Obtain $\textbf{O}^{(n)}_{r_{core}}$ by Equ. (7);

5: Obtain $\textbf{W}^{(n)}_{r_{core}}$ $\leftarrow$
$\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)}\textbf{O}^{(n)}_{r_{core}}$;

6: end for

7: for $r_{core}$ from $1$ to $R_{core}$ do

8: $\widehat{x}^{(n)}_{\Psi^{(n)}_{V}}$ $\leftarrow$
$x^{(n)}_{\Psi^{(n)}_{V}}$;

9: for $r$ from $1$ to $R_{core}$ ($r$ $\neq$ $r_{core}$) do

10: $\widehat{x}^{(n)}_{\Psi^{(n)}_{V}}$ $\leftarrow$
$\widehat{x}^{(n)}_{\Psi^{(n)}_{V}}-\textbf{W}^{(n)}_{r}b^{(n)}_{:,r}$;

11: end for

12: $\textbf{C}_{r_{core}}$ $\leftarrow$
$\textbf{W}^{(n)^{T}}_{r_{core}}\textbf{W}^{(n)}_{r_{core}}$;

13: $\textbf{V}^{(n)}_{r_{core}}$ $\leftarrow$
$-\textbf{W}^{(n)^{T}}_{r_{core}}\widehat{x}^{(n)}_{\Psi^{(n)}_{V}}+\textbf{C}_{r_{core}}b^{(n)}_{:,r_{core}}$;

14: $b^{(n)}_{:,r_{core}}$ $\leftarrow$
$SGD\bigg{(}M,\lambda_{\textbf{B}},\gamma_{\textbf{B}},b^{(n)}_{:,r_{core}},\textbf{V}^{(n)}_{r_{core}}\bigg{)}$;

15: end for

16: end for

17: Obtain $\widehat{\bm{\mathcal{G}}}$ by Equ. (4);

18: for $n$ from $1$ to $N$ do

19: for $i_{n}$ from $1$ to $I_{n}$ do

20: Obtain $\textbf{S}^{(n)}_{(\Psi^{(n)}_{M})_{i_{n}},:}(cache_{S})$;

21: $\textbf{E}^{(n)}_{:,(\Psi^{(n)}_{M})_{i_{n}}}$ $\leftarrow$
$\widehat{\textbf{G}}^{(n)}\textbf{S}^{(n)^{T}}_{(\Psi^{(n)}_{M})_{i_{n}},:}(cache_{E})$;

22: $\textbf{C}^{(n)}$ $\leftarrow$
$\textbf{E}^{(n)}_{:,(\Psi^{(n)}_{M})_{i_{n}}}\textbf{E}^{(n)T}_{:,(\Psi^{(n)}_{M})_{i_{n}}}$;

23: $\textbf{F}^{(n)}$ $\leftarrow$
$\underbrace{-\textbf{X}_{i_{n},(\Psi^{(n)}_{M})_{i_{n}}}^{(n)}\textbf{E}^{(n)T}_{:,(\Psi^{(n)}_{M})_{i_{n}}}}_{cache_{Fact1}}+\overbrace{\underbrace{a_{i_{n},:}^{(n)}\textbf{C}^{(n)}}_{cache_{Fact2}}}^{cache_{Factp},cache_{Factvec}}$;

24: $a_{i_{n},:}^{(n)}$ $\leftarrow$
$SGD\bigg{(}|(\Psi^{(n)}_{M})_{i_{n}}|,\lambda_{\textbf{A}},\gamma_{\textbf{A}},a_{i_{n},:}^{(n)},\textbf{F}^{(n)}\bigg{)}$.

25: end for

26: end for

27: Return: $\textbf{B}^{(n)}$, $n$ $\in$ $\\{N\\}$, $\widehat{\mathcal{G}}$,
$\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$.

###  4.4 Parallel and Distributed Strategy

We first explain the naive parallel strategy (Section 4.4.1), which relies on
the coordinate continuity and the specific style of matricization unfolding of
the sparse tensor to keep the entire continuous accessing in the process of
updating $\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$. Then, we present the
improved parallel strategy (Section 4.4.2), which can save the $N$ compressive
styles of the sparse tensor $\mathcal{X}$ to just a compressive format. At
last, the analysis of the communication cost is reported on Section 4.4.3.

####  4.4.1 Naive Parallel Strategy

To cater to the need of processing big data, the algorithm design shall
leverage the increasing number of cores while minimizing the communication
overehead among parallel threads. If there are $L$ threads, the randomly
selected set $\Psi$ is divided into $L$ subsets
$\big{\\{}\Psi_{1},\cdots,\Psi_{l},\cdots,\Psi_{L}\big{|}l\in\\{L\\}\big{\\}}$
and the parallel computation analysis follows each step in Algorithm 1 as:

(i) Updating the core tensor:

Line 2: The computational tasks of the $L$ Intermediate matrices
$\bigg{\\{}\textbf{H}_{\Psi_{l_{V},:}^{(n)}}^{(n)}$ $\in$
$\mathbb{R}^{|\Psi_{l}|\times\prod\limits_{n=1}^{N}J_{n}}\bigg{|}l\in\\{L\\},n\in\\{N\\}\bigg{\\}}$
be allocated to $L$ threads;

Line 4: Intermediate matrix $\textbf{O}^{(n)}_{r}$ $\in$
$\mathbb{R}^{\prod\limits_{k=1}^{N}J_{k}\times J_{n}}$, $r$ $\in$
$\\{R_{core}\\}$, $n$ $\in$ $\\{N\\}$; thus, the independent computational
tasks of the $\prod\limits_{k=1,k\neq n}^{N}J_{k}$ diagonal sub-matrices can
be allocated to the $L$ threads;

Line $5$: The computational tasks of the $L$ intermediate matrices
$\bigg{\\{}\textbf{H}_{\Psi_{l_{V},:}^{(n)}}^{(n)}\textbf{O}^{(n)}_{r}$ $\in$
$\mathbb{R}^{|\Psi_{l}|\times J_{n}}\bigg{|}l\in\\{L\\},n\in\\{N\\}\bigg{\\}}$
can be allocated to $L$ threads;

Line $8$: The $L$ assignment tasks
$\bigg{\\{}\widehat{x}^{(n)}_{\Psi^{(n)}_{l_{V}}}\bigg{|}l\in\\{L\\}\bigg{\\}}$
can be allocated to the $L$ threads;

Figure 3: Naive Parallel Strategy for SGD$\\_$Tucker.

Line $10$: The computational tasks of the $L$ intermediate matrices
$\bigg{\\{}(\textbf{W}_{r}^{(n)})_{\Psi_{l_{V},:}^{(n)}}b^{(n)}_{:,r}\in\mathbb{R}^{|\Psi_{l}|}\bigg{|}r\in\\{R_{core}\\},n\in\\{N\\}\bigg{\\}}$
can be allocated to $L$ threads;

Line $12$: The computational tasks of the $L$ intermediate matrices
$\bigg{\\{}\textbf{C}_{r_{core}}^{l}=(\textbf{W}_{r}^{(n)})^{T}_{\Psi_{l_{V},:}^{(n)}}(\textbf{W}_{r}^{(n)})_{\Psi_{l_{V},:}^{(n)}}\in\mathbb{R}^{J_{n}\times
J_{n}}\bigg{|}l\in\\{L\\},r\in\\{R_{core}\\},n\in\\{N\\}\bigg{\\}}$ can be
allocated to $L$ threads and the $l$ thread does the intra-thread summation
$\textbf{C}_{r_{core}}^{l}$ $=$
$(\textbf{W}_{r}^{(n)})^{T}_{\Psi_{l_{V},:}^{(n)}}(\textbf{W}_{r}^{(n)})_{\Psi_{l_{V},:}^{(n)}}$.
Then, the main thread sums
$\textbf{C}_{r_{core}}=\sum\limits_{l=1}^{L}\textbf{C}_{r_{core}}^{l}$;

Line $13$: The computational tasks of the $L$ intermediate matrices
$\bigg{\\{}(\textbf{W}_{r}^{(n)})^{T}_{\Psi_{l_{V},:}^{(n)}}\widehat{x}^{(n)}_{\Psi_{l_{V}}^{(n)}}\in\mathbb{R}^{J_{n}}$,
$\textbf{C}_{r}b^{(n)}_{:,r}\in\mathbb{R}^{J_{n}}\bigg{|}r\in\\{R_{core}\\},n\in\\{N\\}\bigg{\\}}$
can be allocated to $L$ threads.

Line $14$: This step can be processed by the main thread.

(ii) Updating factor matrices: $I_{n}$ loops including the Lines $20-24$, $n$
$\in$ $\\{N\\}$ are independent. Thus, the $I_{n}$, $n$ $\in$ $\\{N\\}$ loops
can be allocated to $L$ threads.

The parallel training process for $\textbf{B}^{(n)}$, $n$ $\in$ $\\{N\\}$ does
not need the load balance. The computational complexity of
$\textbf{S}^{(n)}_{(\Psi^{(n)}_{M})_{i_{n}},:}$ for each thread is
proportional to $|(\Psi^{(n)}_{M})_{i_{n}}|$, and the $I_{n}$, $n$ $\in$
$\\{N\\}$ independent tasks are allocated to the $L$ threads. Then, there is
load imbalance problem for updating the factor matrices $\textbf{A}^{(n)}$,
$n$ $\in$ $\\{N\\}$. Load balance can fix the problem of non-even distribution
for non-zero elements.

A toy example is shown in Fig. 3 which comprises of:

(i) Each thread $l$ first selects the even number of non-zero elements: The
$3$ threads $\\{1,2,3\\}$ select $\big{\\{}$ $\\{$
$\bm{x}^{(1)}_{1},\bm{x}^{(1)}_{6}$ $\\}$, $\\{$
$\bm{x}^{(1)}_{13},\bm{x}^{(1)}_{16}$ $\\}$, $\\{$
$\bm{x}^{(1)}_{22},\bm{x}^{(1)}_{27}$ $\\}$ $\big{\\}}$, respectively. Step 1:
The $3$ threads $\\{1,2,3\\}$ construct $\big{\\{}$ $\bm{H}^{(1)}_{1,:}$,
$\bm{H}^{(1)}_{13,:}$, $\bm{H}^{(1)}_{22,:}$ $\big{\\}}$, respectively. Then,
the $3$ threads $\\{1,2,3\\}$ compute $\big{\\{}$
$\textbf{W}^{(1)}_{1}=\bm{H}^{(1)}_{1,:}\textbf{O}^{(1)}$,
$\textbf{W}^{(1)}_{13}=\bm{H}^{(1)}_{13,:}\textbf{O}^{(1)}$,
$\textbf{W}^{(1)}_{22}=\bm{H}^{(1)}_{22,:}\textbf{O}^{(1)}$ $\big{\\}}$,
respectively. Step 2: The $3$ threads $\\{1,2,3\\}$ construct $\big{\\{}$
$\bm{H}^{(1)}_{6,:}$, $\bm{H}^{(1)}_{16,:}$, $\bm{H}^{(1)}_{27,:}$
$\big{\\}}$, respectively. Then, the $3$ threads $\\{1,2,3\\}$ compute
$\big{\\{}$ $\textbf{W}^{(1)}_{6}=\bm{H}^{(1)}_{6,:}\textbf{O}^{(1)}$,
$\textbf{W}^{(1)}_{16}=\bm{H}^{(1)}_{16,:}\textbf{O}^{(1)}$,
$\textbf{W}^{(1)}_{27}=\bm{H}^{(1)}_{27,:}\textbf{O}^{(1)}$ $\big{\\}}$,
respectively. Each thread does the summation within the thread and the $3$
threads do the entire summation by the code $\\#\ pragma\ omp\ parallel\ for\
reduction\ (+:sum)$ for multi-thread summation. We observe that the
computational process of
$\textbf{W}^{(1)^{T}}_{k}\textbf{W}^{(1)}_{k}b^{(1)}_{1,:}$ , $k$ $\in$
$\\{1,6,13,16,22,27\\}$ is divided into the vectors reduction operation
$p=\textbf{W}^{(1)}_{k}b^{(1)}_{1,:}$ and $vec=\textbf{W}^{(1)^{T}}_{k}p$.
Step 3: The $b^{(1)}_{:,1}$ is updated. The process of updating
$b^{(3)}_{:,1}$ is similar the process of updating $b^{(1)}_{:,1}$. The
description is omitted. We observe that each thread selects $2$ elements.
Thus, the load for the $3$ threads is balanced.

(ii) Each thread selects the independent rows and $L$ threads realize naive
parallelization for $\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$ by the $n$th
matricization $\textbf{X}^{(n)}$, $n$ $\in$ $\\{N\\}$: As show in Fig. 3, the
$3$ threads $\\{1,2,3\\}$ update
$\\{a_{1,:}^{(1)},a_{2,:}^{(1)},a_{3,:}^{(1)}\\}$ explicitly. Thus, the
description is omitted. The $2$ threads $\\{1,2\\}$ update
$\bigg{\\{}a_{1,:}^{(1)},\big{\\{}a_{2,:}^{(1)},a_{3,:}^{(1)}\big{\\}}\bigg{\\}}$,
respectively and independently by $\bigg{\\{}$ $\big{\\{}$
$\\{\bm{X}^{(3)}_{1,1}$, $\bm{S}^{(3)}_{:,1}\\}$, $\\{\bm{X}^{(3)}_{1,5}$,
$\bm{S}^{(3)}_{:,5}\\}$, $\\{\bm{X}^{(3)}_{1,8}$, $\bm{S}^{(3)}_{:,8}$,
$\\{\bm{X}^{(3)}_{1,9}$, $\bm{S}^{(3)}_{:,9}\\}$ $\big{\\}}$, $\big{\\{}$
$\\{\bm{X}^{(3)}_{2,2}$, $\bm{S}^{(3)}_{:,2}\\}$, $\\{\bm{X}^{(3)}_{2,7}$,
$\bm{S}^{(3)}_{:,7}\\}$ $\big{\\}}$, $\big{\\{}$ $\\{\bm{X}^{(3)}_{3,5}$,
$\bm{S}^{(3)}_{:,5}\\}$ $\big{\\}}$ $\bigg{\\}}$, respectively, with the
shared matrix $\widehat{\textbf{G}}^{(3)}$. Thread $1$ selects $4$ elements
for updating $a_{1,:}^{(3)}$. Thread $2$ selects $2$ elements and $1$ element
for updating $a_{2,:}^{(3)}$ and $a_{3,:}^{(3)}$, respectively, which can
dynamically balance the load. In this condition, the load for the $3$ threads
is slightly more balanced.

As shown in Fig. 3, the naive parallel strategy for
$\\{\textbf{A}^{(1)},\cdots,\textbf{A}^{(N)}\\}$ relies on the matricization
format $\\{\bm{X}^{(1)},\cdots,\bm{X}^{(N)}\\}$ of a sparse tensor
$\mathcal{X}$, respectively, which is used to avoid the read-after-write and
write-after-read conflicts. Meanwhile, the transformation process for the
compressive format to another one is time consuming. Thus, the space and
computational overheads are not scalable.

Figure 4: Improved Parallel Strategy for SGD$\\_$Tucker.

####  4.4.2 Improved Parallel Strategy

The improved parallel strategy is developed to use only one compressive format
for the entire updating process and hence save the memory requirements. At the
same time, it can avoid the read-after-write or write-after-read conflicts. In
Fig. 4, a toy example of improved parallel strategy is illustrated. As show in
Fig. 4, the $3$ threads $\\{1,2,3\\}$ select $6$ non-zeros elements and the
Coordinate Format (COO) is $\bigg{\\{}$ $\big{\\{}$ $(1,1,1,2.0)$,
$(3,2,1,1.0)$ $\big{\\}}$, $\big{\\{}$ $(1,2,2,2.0)$, $(1,3,2,4.0)$
$\big{\\}}$, $\big{\\{}$ $(1,2,3,5.0)$, $(3,3,3,2.0)$ $\big{\\}}$
$\bigg{\\}}$, respectively. By the structure of COO, the training process of
$B^{(1)}$ and $B^{(3)}$ does not need a specific shuffling order of a sparse
tensor. Thus, the description of updating $b_{:,1}^{(1)}$ and $b_{:,1}^{(3)}$
is omitted.

As for $\bm{X}^{(3)}$ to update $\textbf{A}^{(3)}$, we neglect this condition
because the selected elements of $3$ threads lie in independent rows and it is
trivial to parallelize. In the style of $\bm{X}^{(1)}$ for updating
$\textbf{A}^{(1)}$, updating $a_{1,:}^{(1)}$ relies on
$\big{\\{}(1,1,1,2.0),\textbf{S}_{:,1}^{(1)}\big{\\}}$ (selected by thread
$1$),
$\bigg{\\{}\big{\\{}(1,2,2,2.0),\textbf{S}_{:,5}^{(1)}\big{\\}},\big{\\{}(1,3,2,4.0),\textbf{S}_{:,6}^{(1)}\big{\\}}\bigg{\\}}$
(selected by thread $2$), and
$\big{\\{}(1,2,3,5.0),\textbf{S}_{:,8}^{(1)}\big{\\}}$ (selected by thread
$3$) with the shared matrix $\widehat{\textbf{G}}^{(1)}$. It has following
three steps.

(i) (Lines 20-21 in Algorithm 1) The $3$ threads $\\{1,2,3\\}$ compute
$\bigg{\\{}$
$\textbf{E}_{:,1}^{(1)}=\widehat{\textbf{G}}^{(1)}\textbf{S}_{:,1}^{(1)^{T}}$,
$\big{\\{}\textbf{E}_{:,5}^{(1)}=\widehat{\textbf{G}}^{(1)}\textbf{S}_{:,5}^{(1)^{T}},\textbf{E}_{:,6}^{(1)}=\widehat{\textbf{G}}^{(1)}\textbf{S}_{:,6}^{(1)^{T}}\\}$,
$\textbf{E}_{:,8}^{(1)}=\widehat{\textbf{G}}^{(1)}\textbf{S}_{:,8}^{(1)^{T}}$
$\bigg{\\}}$ independently, by the private cache matrix
$\\{cache_{S},cache_{E}\\}$ of each thread;

(ii) (Lines 22-23 in Algorithm 1) The computational process of
$a_{i_{n},:}^{(n)}\textbf{E}_{:,k}^{(1)}\textbf{E}_{:,k}^{(1)^{T}}$, $k$ $\in$
$\\{1,5,6,8\\}$ is divided into the vectors reduction operation
$cache_{Factp}=a_{i_{n},:}^{(n)}\textbf{E}_{:,k}^{(1)}$ and
$cache_{Factvec}=cache_{Factp}\textbf{E}_{:,k}^{(1)^{T}}$. The $3$ threads
$\\{1,2,3\\}$ compute $\bigg{\\{}$
$a_{i_{n},:}^{(n)}\textbf{E}_{:,1}^{(1)}\textbf{E}_{:,1}^{(1)^{T}}$,
$\big{\\{}a_{i_{n},:}^{(n)}\textbf{E}_{:,5}^{(1)}\textbf{E}_{:,5}^{(1)^{T}},a_{i_{n},:}^{(n)}\textbf{E}_{:,6}^{(1)}\textbf{E}_{:,6}^{(1)^{T}}\\}$,
$a_{i_{n},:}^{(n)}\textbf{E}_{:,8}^{(1)}\textbf{E}_{:,8}^{(1)^{T}}$
$\bigg{\\}}$, respectively and independently, by the private cache
$cache_{Factp}$ and $cache_{Factvec}$ of each thread. Then, the $3$ threads
$\\{1,2,3\\}$ can use the synchronization mechanisms, i.e., $atomic$,
$cirtical$ or $mutex$, of OpenMP to accomplish
$\prod\limits_{k=1,5,6,8}a_{i_{n},:}^{(n)}\textbf{E}_{:,k}^{(1)}\textbf{E}_{:,k}^{(1)^{T}}$.
Then, the results are returned to global shared cache $cache_{Fact2}$;

(iii) (Line 23 in Algorithm 1). The $3$ threads $\\{1,2,3\\}$ compute
$\bigg{\\{}$ $x_{1,1,1}\textbf{E}_{:,1}^{(1)}$,
$\big{\\{}x_{1,2,2}\textbf{E}_{:,5}^{(1)},x_{1,3,2}\textbf{E}_{:,6}^{(1)}\\}$,
$x_{1,2,3}\textbf{E}_{:,8}^{(1)}$ $\bigg{\\}}$, respectively and
independently. Then, the $3$ threads $\\{1,2,3\\}$ can use the synchronization
mechanisms, i.e., $atomic$, $cirtical$ or $mutex$, to accomplish
$\textbf{F}^{(1)}_{1}$. Then, the results are returned to the global shared
cache $cache_{Fact1}$. Eventually, the $I_{1}$ tasks
$SGD\bigg{(}6,\lambda_{\textbf{A}},\gamma_{\textbf{A}},a^{(1)}_{1,:},\textbf{F}^{(1)}_{1}\bigg{)}$
be allocated to the the $3$ threads $\\{1,2,3\\}$ in a parallel and load
balance way. Due to the same condition of updating $a_{3,:}^{(1)}$ and limited
space, the description of updating $a_{3,:}^{(1)}$ is omitted.

By the global shared caches and private caches, $SGD$_$Tucker$ can handle the
parallelization on OpenMP by just a compressive format and the space overhead
is much less than the compressive data structure of a sparse tensor
$\mathcal{X}$; meanwhile, this strategy does not increase extra computational
overhead.

####  4.4.3 Communication in Distributed Platform

In distributed platform, the communication overhead for a core tensor is
$O\bigg{(}\prod\limits_{n=1}^{N}J_{n}\bigg{)}$, which is non-scalable in
HOHDST scenarios. SGD$\\_$Tucker can prune the number of the parameters for
constructing an updated core tensor from
$O\bigg{(}\prod\limits_{n=1}^{N}J_{n}\bigg{)}$ to
$O\bigg{(}\sum\limits_{n=1}^{N}J_{n}R_{core}\bigg{)}$ where $R_{core}$ $\ll$
$J_{n},n\in\\{N\\}$ while maintaining the same overall accuracy and lower
computational complexity. Hence, nodes only need to communicate the Kruckal
product matrices $\bigg{\\{}\textbf{B}_{r}^{(n)}\in\mathbb{R}^{J_{n}\times
R_{core}}\bigg{|}r\in\\{R_{core}\\},n\in\\{N\\}\bigg{\\}}$ rather than the
entire core tensor $\mathcal{G}$ $\in$ $\mathbb{R}^{J_{1}\times
J_{2}\times\cdots\times J_{N}}_{+}$. Hence, SGD$\\_$Tucker features that 1) a
stochastic optimization strategy is adopted to core tensor and factor matrices
which can keep the low computation overhead and storage space while not
compromising the accuracy; 2) each minor node can select sub-samples from
allocated sub-tensor and then compute the partial gradient and the major node
gathers the partial gradient from all minor nodes and then obtains and returns
the full gradient; 3) the communication overhead for information exchange of
the core tensor $\bm{\mathcal{G}}$ is
$O\bigg{(}\sum\limits_{n=1}^{N}J_{n}R_{core}\bigg{)}$.

###  4.5 Complexity Analysis

The space and complexity analyses follow the steps which are presented in
Algorithm 1 as:

######  Theorem 4.

The space overhead for updating $\textbf{B}^{(n)}$, $n$ $\in$ $\\{N\\}$ is
$O\bigg{(}\big{(}MN+R_{core}NJ_{n}\big{)}\prod\limits_{k=1}^{N}J_{k}+R_{core}MJ_{n}+MN+J_{n}\bigg{)}$.

######  Theorem 5.

The computational complexity for updating $\textbf{B}^{(n)}$, $n$ $\in$
$\\{N\\}$ is
$O\bigg{(}(MN+R_{core}NJ_{n})\prod\limits_{k=1}^{N}J_{k}+R_{core}MJ_{n}+R_{core}(R_{core}-1)MN+R_{core}NJ_{n}\bigg{)}$.

###### Proof.

In the process updating $\textbf{B}^{(n)}$, $n$ $\in$ $\\{N\\}$, the space
overhead and computational complexity of intermediate matrices
$\bigg{\\{}\textbf{H}_{\Psi^{(n)}_{V},:}^{(n)},\textbf{O}^{(n)}_{r_{core}},\textbf{W}^{(n)}_{r_{core}},\widehat{x}^{(n)}_{\Psi^{(n)}_{V}},\textbf{V}^{(n)}_{r_{core}}\bigg{|}n\in\\{N\\},r_{core}\in\\{R_{core}\\}\bigg{\\}}$
are
$\bigg{\\{}M\prod\limits_{k=1}^{N}J_{k},J_{n}\prod\limits_{k=1}^{N}J_{k},MJ_{n},M,J_{n}\bigg{|}n\in\\{N\\},r_{core}\in\\{R_{core}\\}\bigg{\\}}$
and
$\bigg{\\{}M\prod\limits_{k=1}^{N}J_{k},\prod\limits_{k=1}^{N}J_{k},MJ_{n}\prod\limits_{k=1}^{N}J_{k}$,
$(R_{core}-1)MJ_{n},MJ_{n}+MJ_{n}^{2}\bigg{|}n\in\\{N\\},r_{core}\in\\{R_{core}\\}\bigg{\\}}$,
respectively. ∎

######  Theorem 6.

The space overhead for updating $\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$ is
$O\bigg{(}(max(|(\Psi^{(n)}_{M})_{i_{n}}|)+1)\prod\limits_{k=1}^{N}J_{k}+max(|(\Psi^{(n)}_{M})_{i_{n}}|)J_{n}+J_{n}\bigg{)}$.

######  Theorem 7.

The computational complexity for updating $\textbf{A}^{(n)}$, $n$ $\in$
$\\{N\\}$ is
$O\bigg{(}(N+MN+M)\prod\limits_{n=1}^{N}J_{n}+\sum\limits_{n=1}^{N}I_{n}J_{n}\bigg{)}$.

###### Proof.

In the process updating $\textbf{A}^{(n)}$, $n$ $\in$ $\\{N\\}$, the space
overhead and computational complexity of intermediate matrices
$\bigg{\\{}\widehat{\textbf{G}}^{(n)},\textbf{S}^{(n)}_{(\Psi^{(n)}_{M})_{i_{n}},:},\textbf{E}^{(n)}_{:,(\Psi^{(n)}_{M})_{i_{n}}},\textbf{F}^{(n)}\bigg{|}i_{n}\in\\{I_{n}\\},n\in\\{N\\}\bigg{\\}}$
are
$\bigg{\\{}\prod\limits_{k=1}^{N}J_{k},max(|(\Psi^{(n)}_{M})_{i_{n}}|)\prod\limits_{k=1}^{N}J_{k},max(|(\Psi^{(n)}_{M})_{i_{n}}|)J_{n},J_{n}\bigg{|}i_{n}\in\\{I_{n}\\},n\in\\{N\\}\bigg{\\}}$
and
$\bigg{\\{}\prod\limits_{k=1}^{N}J_{k},|(\Psi^{(n)}_{M})_{i_{n}}|\prod\limits_{k=1}^{N}J_{k},|(\Psi^{(n)}_{M})_{i_{n}}|J_{n}$,
$J_{n}\bigg{|}i_{n}\in\\{I_{n}\\},n\in\\{N\\}\bigg{\\}}$, respectively. ∎

(a) Movielen-10M

(b) Movielen-20M

(c) Netflix-100M

(d) Yahoo-250M

Figure 5:  Rank scalability for time overhead on full threads. The
computational scalability for P$-$Tucker, HOOI, CD, and SGD$\\_$Tucker on the
4 datasets with successively increased number of total elements,
i.e.,Movielen-10M, Movielen-20M, Netflix-100M, and Yahoo-250M.

(a) Movielen-10M

(b) Movielen-20M

(c) Netflix-100M

(d) Yahoo-250M

Figure 6: Rank scalability for memory overhead on a thread running. (In this work, GB refers to GigaBytes and MB refers to Megabytes). The space scalability for P$-$Tucker, CD, HOOI, and SGD$\\_$Tucker on the 4 datasets with successively increased total elements, i.e., Movielen-10M, Movielen-20M, Netflix-100M, and Yahoo-250M.  TABLE III: Datasets Parameters Datasets | Movielens-100K | Movielens-1M | Movielens-10M | Movielens-20M | Netflix-100M | Yahoo-250M
---|---|---|---|---|---|---
$N_{1}$ | 943 | 6, 040 | 71, 567 | 138, 493 | 480, 189 | 1, 000, 990
$N_{2}$ | 1, 682 | 3, 706 | 10, 677 | 26, 744 | 17, 770 | 624, 961
$N_{3}$ | 2 | 4 | 15 | 21 | 2, 182 | 133
$N_{4}$ | 24 | 24 | 24 | 24 | - | 24
$|\Omega|$ | 90, 000 | 990, 252 | 9, 900, 655 | 19, 799, 448 | 99, 072, 112 | 227, 520, 273
$|\Gamma|$ | 10, 000 | 9, 956 | 99, 398 | 200, 815 | 1, 408, 395 | 25, 280, 002
Order | 4 | 4 | 4 | 4 | 3 | 4
Max Value | 5.0 | 5.0 | 5.0 | 5.0 | 5.0 | 5.0
Min Value | 0.5 | 0.5 | 0.5 | 0.5 | 1.0 | 1.0

##  5 Evaluation

The performance demonstration of SGD$\\_$Tucker comprises of two parts: 1) SGD
can split the high-dimension intermediate variables into small batches of
intermediate matrices and the scalability and accuracy are presented; 2)
SGD$\\_$Tucker can be parallelized in a straightfoward manner and the speedup
results are illustrated. The experimental settings are presented in Section
5.1. Sections 5.2 and 5.3 show the scalability and the influence of
parameters, respectively. At last, Section 5.4 presents the speedup
performance of SGD$\\_$Tucker and comparison with the state of the art
algorithms for STD, i.e., P$-$Tucker [46], CD [47], and HOOI [41]. Due to the
limited space, the full experimental details for HOOI and 2 small datasets,
i.e., Movielen-100K and Movielen-1M, are presented in the supplementary
material.

(a) Computational Time

(b) RMSE

Figure 7: Computational overhead, RMSE, and MAE VS. various rank values for
SGD$\\_$Tucker on 8 Cores. Due to limited space, each dataset is presented in
$4$-order. As Netflix-100M data set only has $3$-order, its 4th index shall be
neglected.

(a) Movielen-10M

(b) Movielen-20M

(c) Netflix-100M

(d) Yahoo-250M

Figure 8: RMSE and MAE vs time for SGD$\\_$Tucker on training set $\Omega$ and
testing set $\Gamma$

(a) Movielen-10M

(b) Movielen-20M

(c) Netflix-100M

(d) Yahoo-250M

Figure 9: RMSE comparison of SGD$\\_$Tucker, P$-$Tucker, and CD on the 4
datasets. Figure 10: Speedup Performance on the 6 datasets.

###  5.1 Experimental Settings

The CPU server is equipped with 8 Intel(R) Xeon(R) E5-2620 v4 CPUs and each
core has 2 hyperthreads, running on 2.10 GHz, for the state of the art
algorithms for STD, e.g., P$-$Tucker [46], CD [47] and HOOI [41]. The
experiments are conducted 3 public datasets : Netflix
111https://www.netflixprize.com/ , Movielens
222https://grouplens.org/datasets/movielens/ , and Yahoo-music
333https://webscope.sandbox.yahoo.com/catalog.php?datatype . The datasets
which be used in our experiments can be downloaded in this link.
444https://drive.google.com/drive/folders/1tcnAsUSC9jFty7AfWetb5R
7nevN2WUxH?usp=sharing For Movielens, data cleaning is conducted and $0$
values are changed from zero to $0.5$ and we make compression for Yahoo-music
dataset in which the maximum value 100 be compressed to 5.0. The specific
situation of datasets is shown in Table III. In this way, the ratings of the
two data sets are in the same interval, which facilitates the selection of
parameters. Gaussian distribution $N\bigg{(}0.5,0.1^{2}\bigg{)}$ is adopted
for initialization; meanwhile, the regularization parameters
$\\{\lambda_{\textbf{A}},\lambda_{\textbf{B}}\\}$ are set as $0.01$ and the
learning rate $\\{\gamma_{\textbf{A}},\gamma_{\textbf{B}}\\}$ are set as
$\\{0.002,0.001\\}$, respectively. For simplification, we set $M=1$. The
accuracy is measured by $RMSE$ and $MAE$ as:

$\displaystyle RMSE$
$\displaystyle=\sqrt{\bigg{(}\sum_{(i_{1},\cdots,i_{N})\in\Gamma}(x_{i_{1},\cdots,i_{N}}-\widehat{x}_{i_{1},\cdots,i_{N}})^{2}\bigg{)}\bigg{/}|\Gamma|};$
(19) $\displaystyle MAE$
$\displaystyle=\bigg{(}\sum_{(i_{1},\cdots,i_{N})\in\Gamma}\bigg{|}x_{i_{1},\cdots,i_{N}}-\widehat{x}_{i_{1},\cdots,i_{N}}\bigg{|}\bigg{)}\bigg{/}|\Gamma|,$

respectively, where $\Gamma$ denotes the test sets. All the experiments are
conducted on double-precision float numeric.

###  5.2 Scalability of the Approach

When ML algorithms process large-scale data, two factors to influence the
computational overhead are: 1) time spent for data access on memory, i.e.,
reading and writing data; 2) time spent on the computational components for
executing the ML algorithms. Fig. 5 presents the computational overhead per
training epoch. The codes of P$-$Tucker and CD have fixed settings for the
rank value, i.e., $J_{1}=\cdots=J_{n}=\cdots=J_{N}$. The rank value is set as
an increasing order, $\\{3,5,7,9,11\\}$. As shown in Fig. 5, SGD$\\_$Tucker
has the lowest computational time overhead. HOOI needs to construct the
intermediate matrices $\textbf{Y}_{(n)}$ $\in$
$\mathbb{R}^{I_{n}\times\prod\limits_{k\neq n}^{N}J_{k}}$, $n$ $\in$ $\\{N\\}$
and SVD for $\textbf{Y}_{(n)}$ (In supplementary material). P$-$Tucker should
construct Hessian matrix and the computational complexity of each Hessian
matrix inversion is $O\big{(}J_{n}^{3}\big{)}$, $n$ $\in$ $\\{N\\}$.
P$-$Tucker saves the memory overhead a lot at the expense of adding the memory
accessing overhead. Eventually, the accessing and computational overheads for
inverse operations of Hessian matrices make the overall computational time of
P$-$Tucker surpass SGD$\\_$Tucker. The computational structure of CD has
linear scalability. However, the CD makes the update process of each feature
element in a feature vector be discrete. Thus, the addressing process is time-
consuming. Meanwhile, as Fig. 5 shows, the computational overhead satisfies
the constraints of the time complexity analyses, i.e., P$-$Tucker, HOOI, CD
(Supplementary material) and SGD$\\_$Tucker (Section 4.5).

The main space overheads for STD comes from the storage need of the
intermediate matrices. The SVD operation for the intermediate matrices
$\textbf{Y}_{(n)}$ $\in$ $\mathbb{R}^{I_{n}\times\prod\limits_{k\neq
n}^{N}J_{k}}$, $n$ $\in$ $\\{N\\}$ makes the HOOI scale exponentially.
P$-$Tucker constructs the eventual Hessian matrices directly which can avoid
the huge construction of intermediate matrices. However, the overhead for
addressing is huge and the long time-consuming for CD lies in the same
situation. CD is a linear approach from the update equation (Supplementary
material). However, in reality, the linearity comes from the approximation of
the second-order Hessian matrices and, thus, the accessing overhead of
discrete elements of CD overwhelms the accessing overhead of continuous ones
(P$-$Tucker and SGD$\\_$Tucker). As shown in Fig. 6, SGD$\\_$Tucker has linear
scalability for space overhead and HOOI is exponentially scalable. P$-$Tucker
and CD have the near-linear space overhead. However, as Fig. 5 shows, the
computational overheads of P$-$Tucker and CD is significantly higher than
SGD$\\_$Tucker. Meanwhile, as Fig. 6 shows, the space overhead satisfies the
constraints of the space complexity analyses, i.e., HOOI, P$-$Tucker, and CD
(Supplementary material) and SGD$\\_$Tucker (Section 4.5).

###  5.3 The Influence of Parameters

In this section, the influence of the rank on the behavior of the algorithm is
presented. We summarize the influence on the computational time and the RMSE
are presented in Fig. 7. The Netflix-100M dataset has only 3-order. Owing to
the limited space, the performances on Netflix-100M dataset are combined with
other 5 datasets (Netflix-100M dataset has only 3-order and the performances
on the 4th order of Netflix-100M dataset should be negligible). The update for
the Kruskal matrices $\textbf{B}^{(n)}$, $n$ $\in$ $N$ on the steps of
beginning and last epochs can also obtain a comparable results. Thus, the
presentation for the Kruskal matrices $\textbf{B}^{(n)}$, $n$ $\in$ $N$ is
omitted. The influence for computational time is divided into 5 sets, i.e.,
$\big{\\{}J_{1}\in\\{5,10,15,20,25\\},J_{k}=5,k\neq 1,R_{Core}=5\big{\\}}$,
$\big{\\{}J_{2}\in\\{5,10,15,20,25\\},J_{k}=5,k\neq 2,R_{Core}=5\big{\\}}$,
$\big{\\{}J_{3}\in\\{5,10,15,20,25\\},J_{k}=5,k\neq 3,R_{Core}=5\big{\\}}$,
$\big{\\{}J_{4}\in\\{5,10,15,20,25\\},J_{k}=5,k\neq 4,R_{Core}=5\big{\\}}$,
$\big{\\{}J_{n}=5,n\in\\{N\\},R_{Core}\in\\{10,15,20,25\\}\big{\\}}$. As the
Fig. 7 (a) shows, the computation time increases with $J_{n}$, $n$ $\in$
$\\{N\\}$ increasing linearly and the $R_{Core}$ only has slight influence for
computational time.

The codes of P$-$Tucker and CD have fixed settings for the rank value
($J_{1}=\cdots=J_{n}=\cdots=J_{N}$). For a fair comparison, the rank values
that we select should have slight changes with the RMSE performances of other
rank values. The computational overhead of P$-$Tucker and CD is sensitive to
the rank value and RMSE is comparable non-sensitive. Hence, we choose a
slightly small value $J_{n}=5|n\in\\{N\\}$. As Fig. 7 (b) shows, when rank is
set as $J_{n}=5,n\in\\{N\\}$, the RMSE is equivalent to other situations on
average.

###  5.4 Speedup and Comparisons

We test the RMSE performances on the 6 datasets and the RMSE is used to
estimate the missing entries. The speedup is evaluated as $Time_{1}/Time_{T}$
where $Time_{T}$ is the running time with $T$ threads. When $T$ $=$
$\\{1,\cdots,8\\}$, as the Fig. 10 shows, the speedup performance of
SGD$\\_$Tucker has a near-linear speedup. The speedup performance is presented
in Fig. 10 which shows the whole speedup performance of the SGD$\\_$Tucker.
The speedup performance is a bit more slower when the number of threads is
larger than $8$. The reason is that the thread scheduling and synchronization
occupies a large part of time. When we use $16$ threads, there is still a
$11X$ speedup and the efficiency is $11/16=68\%$.

The rank value for $3$-order tensor is set to $[5,5,5]$ and the rank value for
$4$-order tensor is set to $[5,5,5,5]$. To demonstrate the convergence of
SGD$\\_$Tucker, the convergence performances of SGD$\\_$Tucker for the
Movielen and Netflix are presented in Fig. 8. As shown in Fig. 8,
SGD$\\_$Tucker can get more stable RMSE and MAE metrices which mean that
SGD$\\_$Tucker has more robustness in large-scale datasets than in small
datasets. As Fig. 9 shows, SGD$\\_$Tucker can not only run faster than the
state of the art approaches, i.e., P$-$Tucker and CD, but also can obtain
higher RMSE value. SGE$\\_$Tucker runs 2$X$ and 20$X$ faster than P$-$Tucker
and CD, respectively, to obtain the optimal RMSE value.

##  6 Conclusion and Future Works

STD is widely used in low-rank representation learning for sparse big data
analysis. Due to the entanglement problem of core tensor and factor matrices,
the computational process for STD has the intermediate variables explosion
problem due to following all elements in HOHDST. To solve this problem, we
first derive novel optimization objection function for STD and then propose
SGD$\\_$Tucker to solve it by dividing the high-dimension intermediate
variables into small batches of intermediate matrices. Meanwhile, the low
data-dependence of SGD$\\_$Tucker makes it amenable to fine-grained
parallelization. The experimental results show that SGD$\\_$Tucker has linear
computational time and space overheads and SGD$\\_$Tucker runs at least 2$X$
faster than the state of the art STD solutions. In the future works, we plan
to explore how to accelerate the SGD$\\_$Tucker by the state of the art
stochastic models, e.g., variance SGD [32], Stochastic Recursive Gradient
[34], or momentum SGD [35].

SGD$\\_$Tucker is a linearly scalable method for STD on big-data platforms.
For the future work, we will embed SGD$\\_$Tucker into popular distributed
data processing platforms such as Hadoop, Spark, or Flink. We will also aim to
extend SGD$\\_$Tucker to GPU platforms, i.e., OpenCL and CUDA.

## Acknowledgement

## References

  * [1] A. Cichocki, N. Lee, I. Oseledets, A.-H. Phan, Q. Zhao, and D. P. Mandic, “Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions,” _Foundations and Trends® in Machine Learning_ , vol. 9, no. 4-5, pp. 249–429, 2016\.
  * [2] V. N. Ioannidis, A. S. Zamzam, G. B. Giannakis, and N. D. Sidiropoulos, “Coupled graph and tensor factorization for recommender systems and community detection,” _IEEE Transactions on Knowledge and Data Engineering_ , 2019.
  * [3] X. Luo, H. Wu, H. Yuan, and M. Zhou, “Temporal pattern-aware qos prediction via biased non-negative latent factorization of tensors,” _IEEE transactions on cybernetics_ , 2019.
  * [4] K. Xie, X. Li, X. Wang, G. Xie, J. Wen, and D. Zhang, “Active sparse mobile crowd sensing based on matrix completion,” in _Proceedings of the International Conference on Management of Data_ , 2019, pp. 195–210.
  * [5] X. Wang, L. T. Yang, X. Xie, J. Jin, and M. J. Deen, “A cloud-edge computing framework for cyber-physical-social services,” _IEEE Communications Magazine_ , vol. 55, no. 11, pp. 80–85, 2017.
  * [6] P. Wang, L. T. Yang, G. Qian, J. Li, and Z. Yan, “Ho-otsvd: A novel tensor decomposition and its incremental decomposition for cyber-physical-social networks (cpsn),” _IEEE Transactions on Network Science and Engineering_ , 2019.
  * [7] T. Wang, X. Xu, Y. Yang, A. Hanjalic, H. T. Shen, and J. Song, “Matching images and text with multi-modal tensor fusion and re-ranking,” in _Proceedings of the 27th ACM international conference on multimedia_. ACM, 2019, pp. 12–20.
  * [8] M. Hou, J. Tang, J. Zhang, W. Kong, and Q. Zhao, “Deep multimodal multilinear fusion with high-order polynomial pooling,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 12 113–12 122.
  * [9] P. P. Liang, Z. Liu, Y.-H. H. Tsai, Q. Zhao, R. Salakhutdinov, and L.-P. Morency, “Learning representations from imperfect time series data via tensor rank regularization,” in _Proceedings of the Annual Meeting of the Association for Computational Linguistics_ , 2019.
  * [10] Y. Liu, L. He, B. Cao, S. Y. Philip, A. B. Ragin, and A. D. Leow, “Multi-view multi-graph embedding for brain network clustering analysis,” in _AAAI Conference on Artificial Intelligence_ , 2018.
  * [11] A. Cichocki, A. Phan, I. Oseledets, Q. Zhao, M. Sugiyama, N. Lee, and D. Mandic, “Tensor networks for dimensionality reduction and large-scale optimizations: Part 2 applications and future perspectives,” _Foundations and Trends in Machine Learning_ , vol. 9, no. 6, pp. 431–673, 2017.
  * [12] L. Van Der Maaten, E. Postma, and J. Van den Herik, “Dimensionality reduction: a comparative,” _Journal of Machine Learning Research_ , vol. 10, no. 66-71, p. 13, 2009.
  * [13] J. Xu, J. Zhou, P.-N. Tan, X. Liu, and L. Luo, “Spatio-temporal multi-task learning via tensor decomposition,” _IEEE Transactions on Knowledge and Data Engineering_ , 2019.
  * [14] X. Liu, X. You, X. Zhang, J. Wu, and P. Lv, “Tensor graph convolutional networks for text classification,” _arXiv preprint arXiv:2001.05313_ , 2020\.
  * [15] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard _et al._ , “Tensorflow: A system for large-scale machine learning,” in _$\\{$ USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$)_, 2016, pp. 265–283.
  * [16] J. Kossaifi, A. Bulat, G. Tzimiropoulos, and M. Pantic, “T-net: Parametrizing fully convolutional nets with a single high-order tensor,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 7822–7831.
  * [17] F. Ju, Y. Sun, J. Gao, M. Antolovich, J. Dong, and B. Yin, “Tensorizing restricted boltzmann machine,” _ACM Transactions on Knowledge Discovery from Data_ , vol. 13, no. 3, pp. 1–16, 2019.
  * [18] M. Wang, Z. Shu, S. Cheng, Y. Panagakis, D. Samaras, and S. Zafeiriou, “An adversarial neuro-tensorial approach for learning disentangled representations,” _International Journal of Computer Vision_ , vol. 127, no. 6-7, pp. 743–762, 2019.
  * [19] C. Zhang, H. Fu, J. Wang, W. Li, X. Cao, and Q. Hu, “Tensorized multi-view subspace representation learning,” _International Journal of Computer Vision_ , pp. 1–18, 2020.
  * [20] Y. Luo, D. Tao, K. Ramamohanarao, C. Xu, and Y. Wen, “Tensor canonical correlation analysis for multi-view dimension reduction,” _IEEE transactions on Knowledge and Data Engineering_ , vol. 27, no. 11, pp. 3111–3124, 2015.
  * [21] X. Liu, X. Zhu, M. Li, L. Wang, E. Zhu, T. Liu, M. Kloft, D. Shen, J. Yin, and W. Gao, “Multiple kernel k-means with incomplete kernels,” _IEEE transactions on pattern analysis and machine intelligence_ , 2019.
  * [22] P. Symeonidis, A. Nanopoulos, and Y. Manolopoulos, “Tag recommendations based on tensor dimensionality reduction,” in _Proceedings of the ACM conference on Recommender systems_ , 2008, pp. 43–50.
  * [23] I. Balazevic, C. Allen, and T. Hospedales, “Tucker: Tensor factorization for knowledge graph completion,” in _Proceedings of the Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , 2019, pp. 5188–5197.
  * [24] J. Liu, P. Musialski, P. Wonka, and J. Ye, “Tensor completion for estimating missing values in visual data,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 35, no. 1, pp. 208–220, 2012.
  * [25] S. Zubair and W. Wang, “Tensor dictionary learning with sparse tucker decomposition,” in _International Conference on Digital Signal Processing_. IEEE, 2013, pp. 1–6.
  * [26] S. Oh, N. Park, S. Lee, and U. Kang, “Scalable tucker factorization for sparse tensors-algorithms and discoveries,” in _IEEE International Conference on Data Engineering_. IEEE, 2018, pp. 1120–1131.
  * [27] V. T. Chakaravarthy, J. W. Choi, D. J. Joseph, P. Murali, S. S. Pandian, Y. Sabharwal, and D. Sreedhar, “On optimizing distributed tucker decomposition for sparse tensors,” in _Proceedings of the International Conference on Supercomputing_ , 2018, pp. 374–384.
  * [28] B. N. Sheehan and Y. Saad, “Higher order orthogonal iteration of tensors (hooi) and its relation to pca and glram,” in _Proceedings of the SIAM International Conference on Data Mining_. SIAM, 2007, pp. 355–365.
  * [29] A. Kulunchakov and J. Mairal, “A generic acceleration framework for stochastic composite optimization,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 12 556–12 567.
  * [30] A. Traoré, M. Berar, and A. Rakotomamonjy, “Singleshot: a scalable tucker tensor decomposition,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 6301–6312.
  * [31] O. A. Malik and S. Becker, “Low-rank tucker decomposition of large tensors using tensorsketch,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 10 096–10 106.
  * [32] R. Johnson and T. Zhang, “Accelerating stochastic gradient descent using predictive variance reduction,” in _Advances in neural information processing systems_ , 2013, pp. 315–323.
  * [33] M. Schmidt, N. Le Roux, and F. Bach, “Minimizing finite sums with the stochastic average gradient,” _Mathematical Programming_ , vol. 162, no. 1-2, pp. 83–112, 2017.
  * [34] L. M. Nguyen, J. Liu, K. Scheinberg, and M. Takáč, “Sarah: A novel method for machine learning problems using stochastic recursive gradient,” in _Proceedings of the International Conference on Machine Learning_. JMLR. org, 2017, pp. 2613–2621.
  * [35] H. Yu, R. Jin, and S. Yang, “On the linear speedup analysis of communication efficient momentum sgd for distributed non-convex optimization,” in _International Conference on Machine Learning_ , 2019, pp. 7184–7193.
  * [36] S. Shalev-Shwartz _et al._ , “Online learning and online convex optimization,” _Foundations and Trends in Machine Learning_ , vol. 4, no. 2, pp. 107–194, 2012.
  * [37] N. Ketkar, “Introduction to pytorch,” in _Deep learning with python_. Springer, 2017, pp. 195–208.
  * [38] J. Verbraeken, M. Wolting, J. Katzy, J. Kloppenburg, T. Verbelen, and J. S. Rellermeyer, “A survey on distributed machine learning,” _ACM Computing Survey_ , vol. 53, no. 2, Mar. 2020.
  * [39] S. Dutta, G. Joshi, S. Ghosh, P. Dube, and P. Nagpurkar, “Slow and stale gradients can win the race: Error-runtime trade-offs in distributed sgd,” in _International Conference on Artificial Intelligence and Statistics_ , 2018, pp. 803–812.
  * [40] H. Ge, K. Zhang, M. Alfifi, X. Hu, and J. Caverlee, “Distenc: A distributed algorithm for scalable tensor completion on spark,” in _2018 IEEE International Conference on Data Engineering_. IEEE, 2018, pp. 137–148.
  * [41] S. Smith and G. Karypis, “Accelerating the tucker decomposition with compressed sparse tensors,” in _European Conference on Parallel Processing_. Springer, 2017, pp. 653–668.
  * [42] Y. Ma, J. Li, X. Wu, C. Yan, J. Sun, and R. Vuduc, “Optimizing sparse tensor times matrix on gpus,” _Journal of Parallel and Distributed Computing_ , vol. 129, pp. 99–109, 2019.
  * [43] V. T. Chakaravarthy, S. S. Pandian, S. Raje, and Y. Sabharwal, “On optimizing distributed non-negative tucker decomposition,” in _Proceedings of the ACM International Conference on Supercomputing_ , 2019, pp. 238–249.
  * [44] O. Kaya and B. Uçar, “High performance parallel algorithms for the tucker decomposition of sparse tensors,” in _2016 45th International Conference on Parallel Processing (ICPP)_. IEEE, 2016, pp. 103–112.
  * [45] V. T. Chakaravarthy, J. W. Choi, D. J. Joseph, P. Murali, S. S. Pandian, Y. Sabharwal, and D. Sreedhar, “On optimizing distributed tucker decomposition for sparse tensors,” in _Proceedings of the International Conference on Supercomputing_ , 2018, pp. 374–384.
  * [46] S. Oh, N. Park, S. Lee, and U. Kang, “Scalable tucker factorization for sparse tensors-algorithms and discoveries,” in _IEEE International Conference on Data Engineering (ICDE)_. IEEE, 2018, pp. 1120–1131.
  * [47] M. Park, J.-G. Jang, and S. Lee, “Vest: Very sparse tucker factorization of large-scale tensors,” _arXiv preprint arXiv:1904.02603_ , 2019.
  * [48] S. Oh, N. Park, J.-G. Jang, L. Sael, and U. Kang, “High-performance tucker factorization on heterogeneous platforms,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 30, no. 10, pp. 2237–2248, 2019.
  * [49] Y. Nesterov, _Introductory lectures on convex optimization: A basic course_. Springer Science & Business Media, 2013, vol. 87.
  * [50] Y. Chi, Y. M. Lu, and Y. Chen, “Nonconvex optimization meets low-rank matrix factorization: An overview,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 20, pp. 5239–5269, 2019.
  * [51] Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” _SIAM Journal on imaging sciences_ , vol. 6, no. 3, pp. 1758–1789, 2013.
  * [52] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [53] H. Li, C. Fang, and Z. Lin, “Accelerated first-order optimization algorithms for machine learning,” _Proceedings of the IEEE_ , 2020.
  * [54] S. J. Reddi, A. Hefny, S. Sra, B. Poczos, and A. Smola, “Stochastic variance reduction for nonconvex optimization,” in _International conference on machine learning_ , 2016, pp. 314–323.
  * [55] H. Li, K. Li, J. An, and K. Li, “Msgd: A novel matrix factorization approach for large-scale collaborative filtering recommender systems on gpus,” _IEEE Transactions on Parallel and Distributed Systems_ , vol. 29, no. 7, pp. 1530–1544, 2017.