A weighted subspace exponential kernel for support tensor machines]A weighted subspace exponential kernel for support tensor machines

[1]Kirandeep<EMAIL_ADDRESS>These authors contributed equally to this work.

2]Sergey<EMAIL_ADDRESS>
1,3]Peter<EMAIL_ADDRESS>
3]Martin<EMAIL_ADDRESS>
3]Max<EMAIL_ADDRESS>These authors contributed equally to this work.

*[1]Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstraße -1, Magdeburg, 39106, Germany

[2]Department of Mathematical Sciences, University of Bath, Bath, BA2 7AY, United Kingdom

[3]Department of Mathematics, TU Chemnitz,  Reichenhainer Str. 41, Chemnitz, 09126, Germany

High-dimensional data in the form of tensors are challenging for kernel classification methods.
To both reduce the computational complexity and extract informative features, kernels based on low-rank tensor decompositions have been proposed.
However, what decisive features of the tensors are exploited by these kernels is often unclear.
In this paper we propose a novel kernel that is based on the Tucker decomposition. For this kernel the Tucker factors are computed based on re-weighting of the Tucker matrices with tuneable powers of singular values from the HOSVD decomposition. This provides a mechanism to balance the contribution of the Tucker core and factors of the data.
We benchmark support tensor machines with this new kernel on several datasets. First we generate synthetic data where two classes differ in either Tucker factors or core, and compare our novel and previously existing kernels. We show robustness of the new kernel with respect to both classification scenarios. We further test the new method on real-world datasets.
The proposed kernel has demonstrated a higher test accuracy than the state-of-the-art tensor train multi-way multi-level kernel, and a significantly lower computational time.

§ INTRODUCTION

Support Vector Machines (SVM) <cit.>, also known as Support Vector Network, maximum margin classifier, are popular machine learning methods, which allow for soft margins and high-dimensional feature embedding by using kernels. However, the standard SVM <cit.> is based on vectors, and may struggle (in terms of both computational complexity and overfitting) for multi-dimensional (tensor) data. Many applications (e.g. in healthcare or signal processing) contain multidimensional data, hence studying kernel methods for tensorial data is an important topic that is addressed in this paper.

The approximation of tensors based on low-rank decompositions has received a lot of attention in scientific computing over recent years <cit.>.
The tensor-based SVM was introduced as Supervised Tensor Learning (STL) in <cit.>.
Using low-rank tensor approximations such as the Canonical Polyadic (CP) <cit.>, Tucker <cit.>, and Tensor Train formats <cit.> within STL alleviates the curse of dimensionality, and allows one to reduce both computational complexity, by computing existing kernels faster,
and overfitting, by designing new dedicated kernels using directly the components of the low-rank decomposition <cit.>.

In the context of kernel methods, the Dual Structure-preserving Kernel (DuSK) for STL, which is particularly tailored to SVM and tensor data, was introduced in <cit.>. This kernel is defined using the CP format.
Later, further kernelization in factors, specifically the Kernelized-CP (KCP) factorization, have been introduced <cit.>, and the entire technique has been called the Multi-way Multi-level Kernel (MMK) method.
Once an accurate CP approximation is available, DuSK and MMK typically deliver an accurate and efficient classification.
However, the CP approximation of arbitrary data can be numerically unstable and difficult to compute <cit.>. In general, any optimization method (Newton, Steepest Descent or Alternating Least Squares) to obtain the CP decomposition might return only a locally optimal solution, and it is difficult to assess whether this is a local or global optimum.

In contrast, the Tucker approximation problem is well-posed, and a quasi-optimal Tucker approximation can be computed reliably by a sequence of singular value decompositions (SVD) <cit.>.
Therefore, the Tucker format is also used often in data science.
In <cit.> the authors have adopted the Tucker decomposition of the weight parameter to retain more structural information, and <cit.> extended this by using a Genetic Algorithm (GA) prior to the Support Tucker Machine (STuM) for the contraction of the input feature tensor. In <cit.>, the authors proposed to minimize the rank of the weight parameter with the orthogonality constraints on the columns of the weight parameter instead of the classical maximum-margin criterion, and in <cit.> the orthogonality constraints are relaxed to further improve Wolf's method.

Further understanding of the KCP approach <cit.> is provided by a kernelized Tucker model, inspired by <cit.>.

The Tensor Train (TT) decomposition offers a stable approximation similarly to the Tucker format, whereas scaling to higher dimensions like the CP format.	A straightforward generalization of DuSK (MMK) to the TT format was proposed in <cit.>.

However, why exactly the kernels based on low-rank decompositions are good for classification remains unclear.
Moreover, since any tensor decomposition is a nonlinear parametrization of the tensor, its representation may be not unique.
For example, Tucker and TT decompositions are invariant to rotation and scaling of the factors.
These formats can also be converted from one to another, albeit with a change of ranks.
Eliminating redundancy in rotation, scaling, and TT to CP conversion in the TT-MMK method has significantly improved the classification accuracy <cit.>.

Further attempts to understand the key features of tensorial data and design the kernel accordingly include the Tucker subspace kernel <cit.>.
Here, the kernel compares projectors onto the subspaces spanned by Tucker factors.
The latter are known to effectively capture the multilinear structure of the data. For example, <cit.> used HOSVD for unsupervised feature extraction. Multiway analysis enables one to effectively capture the multilinear structure of the data, which is usually available as a priori information about the data. In <cit.> a subspace learning technique for Face Recognition was introduced.
The factor match score, a consistent way of comparing the feature vectors of tensor decompositions, has been introduced in <cit.>.

However, the data may contain decisive features not only in the Tucker subspaces, but also in the Tucker core. How to capture both in a computationally efficient way remained largely an open problem. This paper aims to fill this gap by introducing a novel kernel that show high robustness with respect to where the classification information are contained within the tensor.

§.§ Novel contributions

The main aim of this paper is to introduce a novel kernel that shows high robustness with respect to where the classification information is contained within the tensor. The main contributions with this novel kernel are summarized as follows:

* We propose a new form of writing the Tucker (HOSVD) decomposition with weighted factors.

* Based on this form, we propose a new kernel for support tensor machines, which admits fast computation,
whereas the weighting takes into account both Tucker factors and core in building the nonlinear decision boundary.

* Using synthetic data with class assignment based on either Tucker factors or core, we confirm that the new kernel provides an accurate classification in all cases, in contrast to existing Tucker-based kernels.

* Finally, we test that the new kernel provides higher classification accuracy than state-of-the-art methods also on real datasets.

§ NOTATION AND BACKGROUND

This section sets up and extends notations for tensors and the binary classification problem to multi-dimensional data.

§.§ Tensor Algebra

In context of numerical multi-linear algebra, a multidimensional array is called a tensor. Tensors are a generalization of matrices (2-modes; rows and columns) with a higher number of dimensions / modes. We denote all tensors by a calligraphic letter $\tensx$. We assume that all tensors are real-valued. For a general introduction to tensors and their properties we refer to <cit.> and the references mentioned therein. We summarize the common notations encountered in this paper in Table <ref>.

Symbol     Description     Definition

$\overline{i_1,\ldots,i_M}$     multi-index     $\overline{i_1,\ldots,i_M} = 1 + \sum_{k=1}^{M}(i_k-1) \prod_{m=1}^{k-1} I_m$

$x$         scalar value     $x \in \R$

vector      $\bx \in \R^{I}$

matrix      $\bX \in \R^{I \times J}$

tensor     $\tensx \in \dims R I M$

$\tenem x$     $(i_1,i_2,\ldots,i_M)$-entry of a tensor

m     $m$-mode matricization     $(x_{(m)})_{i_m,\overline{i_1,\ldots,i_{m-1},i_{m+1},\ldots,i_M}}= x_{i_1,\ldots,i_M}$

$\tensx \nmod m \bA$      $m$-mode product    $\nten {(\tensx \nmod m  \bA)} m = \bA \nten \tensx m,$  $\bA \in$ R

$\tenz = \tensx \circ \tensy$      outer product      $z_{i_1, \ldots, i_{M}, j_1, \ldots, j_N} =  x_{i_1,\ldots,i_M} y_{j_1, \ldots, j_N}.$

$\bA \otimes \bB$          Kronecker product     $
\begin{bmatrix}
a_{i,1} \bB & \cdots & a_{i,J} \bB\\
\end{bmatrix} \in \RIJKL$

$\bA \in \RIJ, \bB \in \RKL$

$\langle M \rangle$     integer values from 1 to $M$     $\{1,2, \cdots, M\}$

$\langle \tensx, \tensy \rangle$          inner product of tensors  and     $ \teninner.$

$\| \tensx \|$     Frobenius norm of the tensor $\tensx$     $\sqrt{\langle \tensx, \tensx \rangle}$

Tensor notation used in this paper.

§.§ Support Tensor Machines for supervised learning

Although tensor objects can be reshaped into vectors, the structural information encoded in the tensorial data are lost. For example, in an fMRI image, the values of adjacent voxels are typically close to each other <cit.>.
It is natural to replace the vector-valued SVM by a tensor-valued SVM  (cf. Supervised Tensor Learning (STL) introduced in <cit.>. An extension of the STL using kernelized tensor factorization with maximum-margin criterion (SVM) is given in <cit.>. This preserves the nonlinear structure and enhances the overall performance of the STL model, called Kernelized Support Tensor Machines (KSTM).

§.§.§ Kernelized Support Tensor Machine

The KSTM is a binary classification model for $N$ tensor input data points $\lbrace \left( \tensx_i, y_i\right) \rbrace _ {i = 1}^N$
where each tensor is of the form $\tensx_i \in \dims R I M$ with labels $y_i \in \{0,1\}$ leading to a nonlinear decision boundary.
The method follows a maximum margin approach to get the separation hyperplane. Hence, the objective function for a nonlinear boundary can be written as follows <cit.>:

\begin{align}\label {eq: tenPrimOpt}
&\underset{w,b}{\text{min}} \quad \frac{1}{2} \norm{w}^2 + C \sum_{i=1}^{N} \xi_i \\  \nonumber
\text{subject to} \quad & y_i (\langle \Psi(\tensx_i), w \rangle + b) \geq 1-\xi_{i} ~~ \quad \xi_{i}~\geq 0 ~~~ \forall i.
\end{align}

The classification setup given in (<ref>) is known as Support Tensor Machine (STM) <cit.>. The dual formulation of the corresponding primal problem is given as follows:

\begin{align}\label {eq: tenDualOpt}
&\underset{\alpha_1,\ldots,\alpha_{N}}{\text{max}} \quad  \sum_{i = 1}^{N} \balpha_i - \frac{1}{2} \sum_{i = 1}^{N}\sum_{j = 1}^{N} \balpha_i \balpha_j y_i y_j \langle \Psi(\tensx_i),\Psi(\tensx_j) \rangle \nonumber\\
&  \text{subject to} \quad 0 \leq \balpha_i \leq C,  \quad \sum_{i = 1}^{N} \balpha_i y_i = 0 ~~ \forall i.
\end{align}

The nonlinear feature embedding for tensor inputs in a tensor space to a feature space is analogous to working with vector inputs in a vector space. We can define an embedding from low-dimensional tensor-product space to the tensor-product Reproducing Kernel Hilbert Space <cit.>. And by using Mercer's Theorem, i.e. having a kernel that is positive semidefinite, we can construct a feature embedding $\Psi$ such that,

\begin{align*}
\Psi : \dims R I M \text{ (input space)}  &\rightarrow \mathbb{F}  \text{ (feature space)}\\
\exists~ K \colon \dims R I M \times \dims R I M &\mapsto \R ~~~\text{s.t.}~~~ K\left(\tensx,\tensx'\right) = \langle \Psi(\tensx), \Psi(\tensx') \rangle_{\mathbb{F}}.
\end{align*}

This is not only computationally tractable but also avoids the explicit computation of the function $\Psi$. The kernel matrix that results from the continued evaluation of the kernel function on the data points is then positive semidefinite. With the help of the kernel, a linear learning algorithm can learn a  nonlinear boundary, without explicitly knowing the nonlinear function $\Psi$.

Therefore, by using the kernel trick, KSTM is defined as follows:

\begin{align} \label{eq: kerSTM}
&\underset{\alpha_1,\ldots,\alpha_{N}}{\text{max}} \quad  \sum_{i = 1}^{N} \balpha_i - \frac{1}{2} \sum_{i = 1}^{N}\sum_{j = 1}^{N} \balpha_i \balpha_j y_i y_j K(\tensx_i,\tensx_j) \nonumber \\
&  \text{subject to} \quad 0 \leq \balpha_i \leq C, \quad  \sum_{i = 1}^{N} \balpha_i y_i = 0 ~~ \forall i.
\end{align}

Once we have the real-valued function (kernel) value for each pair of tensors, we can use the state-of-the-art LIBSVM implementation <cit.> , which relies on the Sequential Minimal Optimization algorithm to optimize the weights $\alpha_i$. Hence, the preeminent part is the kernel function $K(\tensx_i,\tensx_j)$.

The STM classifier for predicting labels for unseen test data in tensor form is then given by

\begin{align}\label{eq:classifier}
G (\tensx) & = \text{sign} \left( \sum_{i = 1}^{N} \balpha_i y_i K(\tensx_i,\tensx) + b_{0} \right),
\end{align}

\begin{align}\label{eq:compbval}
b_0 &=\frac{1}{N_0} \sum_{i: \balpha_i\in(0,C)} \left(y_i-\sum_{j=1}^{N} \balpha_j K(\tensx_j, \tensx_i) \right), ~~~ \text{with} ~~~ N_0= \sum_{i: \alpha_i \in (0,C)} 1 .
\end{align}

The only task needed for the KSTM is thus to choose a well-suited kernel function. This way, we can work with the input data in a high-dimensional space, while all computations are performed in the original low-dimensional space. We will discuss possible choices for tensor kernels in Sec. <ref>, where we will also introduce a novel tensor kernel. These kernels are based on low-rank formats for tensors, which we introduce now.

§.§ Low-rank Tensor Decompositions

Given the complexity of storing the full tensor $\tensx$, it is often desirable to have a different potentially more economic representation.  As such, tensor decomposition methods have seen much progress over the last two decades, and they are applied to solve problems of varying computational complexity. The main goal is the linear (or at most polynomial) scaling of the computational complexity in the dimension (order) of a tensor. The key ingredient is the separation of variables via approximate low-rank factorizations.

Canonical Polyadic decomposition
The Canonical Polyadic (CP) decomposition of an $M^{th}-$order tensor $\tensx \in \dims R I M$ is a factorization into a sum of rank-one components <cit.>, which is given element-wise as

\begin{align}\label{eq:cpkrush}
\tenem x & \cong \sum_{r=1}^{R} a^{(1)}_{i_1,r} a^{(2)}_{i_2,r} \cdots a^{(M)}_{i_M,r},  \nonumber\\
\mbox{or shortly,} ~~\qquad \quad \quad \tensx & \cong \llbracket \cpd \rrbracket,
\end{align}

where  $\mathbf{A}^{(m)} = \left[ a^{(m)}_{i_m,r} \right] \in \R^{I_m \times R}$, $m=1,\ldots,M$, are called factor matrices of the CP decomposition, see Fig. <ref>, and $R$ is called the CP-rank. The notation $\llbracket \cpd \rrbracket$ is also called the Kruskal representation of the tensor.
Despite the simplicity of the CP format, the problem of the best CP approximation is often ill-posed <cit.>.
A practical CP approximation can be computed via the Alternating Least Squares (ALS) method <cit.>, but convergence may be slow.
It may also be difficult to choose the rank $R$.

CP decomposition of a 3-way tensor of rank $R$.

Tucker Decomposition

The Tucker decomposition consists of a decomposition of the tensor into matrices and a core tensor, where the core tensor has smaller dimension compared to the original tensor.  For a given tensor $ \tensx \in \dims R I M$ the Tucker decomposition is as then given by,

\begin{align}\label{eq:tucker}
\tensx & \cong \sum_{r_1 = 1}^{R_1} \sum_{r_2 = 1}^{R_2} \ldots \sum_{r_M= 1}^{R_M} g_{r_1 r_2 \ldots r_M} \left( \bu_{r_1}^{(1)} \circ \bu_{r_2}^{(2)} \circ \ldots \circ \bu_{r_M}^{(M)}\right) \nonumber \\
& = \llbracket \teng ;  \bU^{(1)}, \bU^{(2)}, \ldots, \bU^{(M)} \rrbracket . \\
& \nonumber
\end{align}

Tucker decomposition of a  3-way tensor.

Here, $\teng$ is a tensor of size $\R^{R_{1} \times R_{2} \times \cdots \times R_{M}}$ and $R_{m}$ is the Tucker rank in each mode matricization of the tensor $\tensx$. A crucial advantage of the Tucker format (and all tree tensor networks) is the ability to perform algebraic operations directly on the component tensors, avoiding full tensors. Moreover, we can compute a quasi-optimal Tucker approximation of any given tensor using the SVD. This builds on the fact that the Tucker decomposition constitutes a successive matrix factorization, where each Tucker rank is the matrix rank of the appropriate unfolding of the tensor, and hence the Tucker approximation problem is well-posed <cit.>. One unique advantage of the Tucker format is the interpretability of the leaf-components $\bU^{(m)}$: Since they result directly from an SVD of the $m$-th matricization, their columns constitute an orthonormal basis of the subspace of $\mathbb R^{I_m}$ that the data lies in.

Tensor Train decomposition
The Tensor Train (TT) <cit.> decomposition of an $M^{th}-$order tensor  $\tensx \in \dims R I M$ is defined element-wise as

\begin{align}\label{eq:tt}
\tenem x & \cong \sum_{r_0,\ldots,r_{M}}c_{r_0,i_1,r_1}^{(1)} c_{r_1,i_2,r_2}^{(2)} \cdots c_{r_{M-1},i_M,r_M}^{(M)}, \nonumber\\
\tensx & \cong\llangle \tendtt C 1, \tendtt C 2, \ldots, \tendtt C M \rrangle,
\end{align}

where C m $\in$ m, $m=1,\ldots,M,$ are 3rd-order tensors called TT-cores (see Fig. <ref>), and $R_0,\ldots,R_M$ with $R_0 = R_M = 1$ are called TT-ranks.

TT decomposition of a 3-way tensor.

Since the TT decomposition is also a tree-based tensor format, all above-mentioned advantages of the Tucker format also translate to the TT format. Furthermore, complexity of the TT format is quadratic in the ranks, whereas the Tucker format scales exponentially with the Tucker ranks. This reduction of complexity however comes at the cost of interpretability: There is no straightforward way to interpret the meaning of the TT components.

§.§.§ Converting Tucker and TT into CP

It is easy to see from (<ref>), (<ref>), and (<ref>) that one can convert a tensor in Tucker or TT format into the CP format without too much effort: Summing over all Tucker ranks $R_1,\ldots,R_M$ or all TT-ranks $r_0,\ldots,r_M$ will yield a sum of rank-one tensors as required in the CP format. We emphasize here that this does not result in a minimal CP decomposition but only in a CP representation of the tensor. This can however still be useful, as long as the obtained CP representation retains the interpretable qualities of the Tucker or TT decomposition.

The main difficulty with obtaining a meaningful conversion is that none of the decomposition formats is really unique: CP allows for a rescaling of columns of the factor matrices, and for Tucker and TT, one can insert identity matrices $I = Q Q^{-1}$ between the modes without changing the tensor.

In a previous work <cit.>, the problem of meaningfully converting TT into CP was overcome by enforcing uniqueness in the TT-SVD, then converting into CP, and then equilibrating the column norms of the factor matrices in order to avoid ambiguity in the CP representation.

This technique can be used analogously for the conversion of Tucker into CP: In the Higher Order SVD (HOSVD), we enforce uniqueness by fixing the sign of each singular vector. This is done by finding the element of maximum absolute value and making it positive. The Tucker tensor is then converted to CP by simply summing over all ranks $R_1,\ldots,R_M$ and a norm equilibration is performed in order to distribute the scalar $g_{r_1r_2\ldots r_M}$ in (<ref>) across all factors.

§ KERNELS FOR TENSOR DATA IN SVM

In this section, we discuss possible choices for tensor kernels. First, we briefly recapitulate existing tensor kernels before we introduce our new kernel. This, together with the numerical study of tensor kernel in Sec. <ref>, is the main result of this article. At the end of this section, we compare the complexity of computing the different kernels.

§.§ Existing kernels

The Gaussian kernel
The natural idea of defining a kernel for tensorial data would be to extend the classical Gaussian kernel directly from vector to tensor format. That is, the computation can be given directly as follows,

\begin{align}
K(\tensx, \tensy) &= \exp \Big(\frac{- \norm{\tensx - \tensy}_F^2}{2 g^2} \Big),
\end{align}

with $g$ being the length scale of the kernel. We compute the distance between the two input tensors using the Frobenius norm. This norm can be computed efficiently in each of the tensor formats introduced above  (Table <ref> shows the leading terms of the complexity estimates). However, by treating the tensor as a simple vector, we lose valuable information about the different tensor modes. It has been observed (e.g. in <cit.>) that this straightforward idea yields suboptimal classification results and that it can be improved by introducing more sophisticated tensor kernels.

Dual Structure-preserving Kernel
The Dual Structure-preserving Kernel (DuSK) was introduced first in <cit.> for a rank-one tensor factorization and was later extended for the Kernelized CP decomposition in <cit.>. For given tensors $\tensx \in \dims R I M$ and $\tensy \in \dims R I M$ and their corresponding CP decomposition given by $\llbracket \bA^{(1)}, \bA^{(2)}, \ldots, \bA^{(M)} \rrbracket$ and $\llbracket \bB^{(1)}, \bB^{(2)}, \ldots, \bB^{(M)} \rrbracket$, the formulation of the kernel approximation by DuSK is given as follows:

\begin{align}\label{eq:kerAPRX}
\langle \Psi(\tensx), \Psi(\tensy) \rangle & =  K(\tensx,\tensy)\nonumber\\
& = K \left( \sum_{i=1}^{R} \ba_i^{(1)}\otimes \ba_i^{(2)}\otimes \cdots \otimes \ba_i^{(M)}, \sum_{j=1}^{R} \bb_j^{(1)}\otimes \bb_j^{(2)}\otimes \cdots \otimes \bb_j^{(M)}  \right)\nonumber\\
& = \sum_{i,j = 1}^{R}\nk(\ba_{i}^{(1)},\bb_{j}^{(1)}) \nk(\ba_{i}^{(2)},\bb_{j}^{(2)}) \cdots \nk(\ba_{i}^{(M)},\bb_{j}^{(M)}),
\end{align}

\begin{align}
\nk(\ba,\bb) = \exp \left(\frac{-\norm{\ba-\bb}^2}{2 g^2} \right).
\end{align}

In short, we evaluate the kernel function $\nk(\cdot,\cdot)$ on the individual factors of the CP decomposition. The motivation of DuSK is simple: Since the CP decomposition is often unique (up to norm equilibration), comparing the feature vectors in each mode directly will most likely improve classification. However, this kernel is inherently designed for CP tensors, which is why we have to convert other tensor formats into the CP format first (see Sec. <ref> and <cit.>).

The subspace kernel
Instead of comparing the feature vectors in the CP decomposition, one can use a similar approach for the Tucker format. Here, the feature vectors are stored in the leaf-components $\bU^{(m)}$ and they span the column spaces of the $m$-th matricizations. Thus, it makes sense to compare the projections onto these subspaces. This kernel was introduced in <cit.>.

Let $\tensx \in \dims R I M$ denote an $M^{th}$-order tensor, when the SVD is applied on the mode-$m$ unfolding as  $\tensx_{(m)} = \bU_{\tensx}^{(m)} \Sigma_{\tensx}^{(m)} \bV_{\tensx}^{(m)^{\mathsf T}}$ and similarly for $\tensy \in \dims R I M$,  $\tensy_{(m)} = \bU_{\tensy}^{(m)} \Sigma_{\tensy}^{(m)} \bV_{\tensy}^{(m)^{\mathsf T}}$, then the Chordal distance-based kernel is defined as,

\begin{align}
K(\tensx,\tensy) &= \prod_{m = 1}^{M} \exp \Big(-{\frac{1}{2 g^2} \norm{ \bU_{\tensx}^{(m)} \bU_{\tensx}^{(m)^{\mathsf T}} - \bU_{\tensy}^{(m)} \bU_{\tensy}^{(m)^{\mathsf T}}}}_F^2 \Big)
\end{align}

This kernel provides us with rotation and reflection invariance for elements on the Grassmann manifold. Furthermore, the kernel does not see the core tensor $\teng$ and is invariant under rescaling of the feature vectors (i.e., the columns of the leaf matrices).

§.§ The weighted subspace exponential kernel

The subspace kernel performs well when the information about the classification is stored in the subspaces $\bU_{\tensx}^{(m)}$ and it performs poorly when the information is mostly contained in the core tensor $\teng$. This is confirmed in our synthetic numerical experiments below. The DuSK uses a similar strategy as the subspace kernel: Here, we compare all the feature vectors in the CP decomposition. DuSK therefore also performs better if most of the information is in the subspaces, i.e., in the feature vectors.

The main contribution of this article is an improved tensor kernel for Tucker tensors that also includes information of the core tensor $\teng$,  combining the strengths of the subspace kernel and DuSK, and that outperforms both of them in common scenarios, and can be computed more efficiently than DuSK.

We first observe that in the computation of the SVD of a matricized tensor $\tensx_{(m)}$, we can shift any power of the singular values into either the left or the right singular matrices:

\begin{equation*}
\tensx_{(m)} = \bU^{(m)} \boldsymbol\Sigma^{(m)} (\bV^{(m)})^T = \bU^{(m)} \left( \boldsymbol\Sigma^{(m)} \right)^p \left( \boldsymbol\Sigma^{(m)} \right)^{1-p} (\bV^{(m)})^T
\end{equation*}

for any $p\in \mathbb{R}$ (assuming no zero singular values, or defining $0^0=1$).
Using this, we can distribute singular values over the Tucker factors in the HOSVD (see Algorithm <ref>) and we use the resulting rescaled features $\bar \bU^{(m)} = \bU^{(m)} \left( \boldsymbol\Sigma^{(m)} \right)^p$ for the computation of the kernel. We note that by default, we choose $p = \frac{1}{M}$ as this corresponds to distributing the norm $\| \tensx \| = \| \boldsymbol \Sigma^{(m)} \|_F$ equally over the $M$ feature matrices, which provides accurate classification in practice.

Since the subspace kernel is invariant under rescaling of the features, we compute the Euclidean distances instead and sum over the exponential kernels of all combinations, noting that if the distances are large, these terms will be negligible. The result is similar to DuSK, but it uses the feature vectors from the Tucker decomposition and the order of the sum over the ranks and the product over the tensor order is reversed:

\begin{align}\label{eq:rev-wFMS}
K(\tensx,\tensy) &=   \prod_{m = 1}^{M} \sum_{i,j=1}^{R_m} \nk((\bu_{\tensx}^{(m)})_i, (\bu_{\tensy}^{(m)})_j) = \SD{\prod_{m = 1}^{M}  \sum_{i,j = 1}^{R_m}\exp \left(- \frac{1}{2 g^2}\norm{\bar \bu_{\tensx,i}^{(m)} - \bar \bu_{\tensy,j}^{(m)}}_F^2\right),}
\end{align}

where, as in the subspace kernel, we distinguish the SVD of the two tensors, writing $\tensx_{(m)} = \bU_{\tensx}^{(m)} \Sigma_{\tensx}^{(m)} \bV_{\tensx}^{(m)^T}$ and $\tensy_{(m)} = \bU_{\tensy}^{(m)} \Sigma_{\tensy}^{(m)} \bV_{\tensy}^{(m)^T}$.
We call this kernel the weighted subspace exponential kernel or WSEK for short. We also considered leaving the order of the sum and the product the same as in DuSK, but this would need further considerations if the Tucker ranks across the modes are not all the same. Furthermore, our experiments have not shown any significant difference in classification accuracy with respect to the order of $\prod_{m = 1}^{M}$ and $\sum_{i,j=1}^{R_m}$.

Weighted HOSVD

Given tensor  $\tensx \in \dims R I M$, Tucker ranks $R_1,\ldots,R_M$, weighting power $p$ (default $p=1/M$).
Tucker factors $\bar \bU^{(1)},\ldots,\bar \bU^{(M)}$ and core $\teng$.
$m = 1$ to $M$
Step 1: Computing uniqueness-enforced HOSVD
Compute SVD $\left[ \bU^{(m)}, \boldsymbol\Sigma^{(m)}, (\bV^{(m)})^{\mathsf T} \right] = \text{svd}(\tensx_{(m)})$,
$\boldsymbol\Sigma^{(m)} = \text{diag}(\sigma_1^{(m)}, \sigma_2^{(m)}, \ldots, \sigma_{I_m}^{(m)})$
$r_m = 1$ to $R_m$
$i_{r_m}^* = \arg\max_{i=1,\ldots,I_m} |u_{i,r_m}^{(m)}|$
$\mathbf{\hat u}_{r_m}^{(m)}:= \mathbf{u}_{r_m}^{(m)} / \mathrm{sign}(u_{i^*_{r_m},r_m}^{(m)})$
$\hat \bU^{(m)} = [\mathbf{\hat u}_1^{(m)}, \mathbf{\hat u}_2^{(m)}, \ldots, \mathbf{\hat u}_{R_m}^{(m)}]$

Step 2: Computing norm weighted factors
$\bar \bU^{(m)} = \hat \bU^{(m)} \boldsymbol\Sigma_p^{(m)}$
where $\boldsymbol\Sigma_p^{(m)} =  \text{diag}((\sigma_1^{(m)})^{p}, (\sigma_2^{(m)})^{p}, \ldots, (\sigma_{R_m}^{(m)})^{p}) $
$\teng \leftarrow \texttt{ttm}(\tensx, \hat \bU^{(m)} (\boldsymbol\Sigma_p^{(m)})^{-1}, m)$      <cit.>

§.§ Computational complexity of the different kernels

For a large number $N$ of data points, for large tensor order $M$ or large mode dimensions $I_m$, computing the different tensor kernels can be very time consuming. If the data input is already given in CP format, computing the DuSK is not too expensive. But most often, the data is given as a full tensor. In these cases, we prefer to compute first a TT or Tucker decomposition of the tensor and then convert it into CP, in order to circumvent the aforementioned numerical issues with the computation of the CP decomposition. Here, we will however only note the complexity of the kernel computation with respect to the given ranks. We denote the maximal dimensions or ranks by $I = \max_{m = 1,\ldots,M} I_m$, $R_{Tucker} = \max_{m = 1,\ldots,M} R_m$, and $R_{TT} = \max_{m = 1,\ldots,M-1} r_m$. The ranks are also to be understood as the maximal respective rank of all data inputs.

Theoretical complexity of computing a single entry of the different kernel matrices from data given in different formats. Note that some kernel-format combinations are not defined.

Full     CP     Tucker     TT

Gaussian     $\mathcal{O}(I^M)$     $\mathcal{O}(MIR^2)$     $\mathcal{O}(MR_{Tucker}^{(M+1)} + MIR_{Tucker}^2)$     $\mathcal{O}(MIR_{TT}^3)$

DuSK     –     $\mathcal{O}(MIR^2)$     $\mathcal{O} (MIR_{Tucker}^{2M})$     $\mathcal{O} (MIR_{TT}^{2(M-1)})$

Subspace     –     $\mathcal{O}(MIR(I+R))$     $\mathcal{O}(MI R_{Tucker} (I+R_{Tucker}))$     –

WSEK     –     $\mathcal{O}(MIR^2)$     $\mathcal{O}(MIR_{Tucker}^2)$     –

Table <ref> summarizes the computational complexity for a single entry of the kernel matrix. We note that all kernels except the Gaussian kernel can only be computed if the tensor is in a low-rank format. Also, we interpret a tensor in CP format to be a Tucker tensor with diagonal core tensor $\teng$ and thus the subspace kernel and WSEK are computed using the factor matrices. Furthermore, in the cases of Tucker and TT tensors, we report on the complexity using naive matrix multiplication ($\mathcal O(n^3)$).

We observe that for large tensor order $M$, computation of the Gauss kernel is prohibitive if it is not done in a low-rank format. If the CP rank $R$ is small, all kernels can be computed efficiently. However, if the CP decomposition has to be obtained by conversion from Tucker or TT, these ranks can be large and DuSK suffers from the curse of dimensionality. WSEK is even more efficient than the subspace kernel. We will report on the CPU times for our numerical experiments in Sec. <ref>.

§ A NUMERICAL STUDY ON SYNTHETIC DATA

As mentioned above, knowledge about the tensor structure can and should be exploited when choosing the tensor kernel $K$. In this section, we explore why the DuSK performs well in many cases by comparing it to the Gaussian kernel and the Subspace kernel in a synthetic experimental setting. Furthermore, we show that our proposed WSEK retains the advantages of DuSK, while outperforming it in cases where DuSK is less suitable.

§.§ Interpreting CP and Tucker

The CP decomposition of a tensor can be seen as a special case of the Tucker decomposition with a diagonal core tensor. More precisely, if the matrix ranks of the factor matrices $\mathbf A^{(1)},\ldots,\mathbf A^{(M)}$ are smaller than $R$, we can find a Tucker decomposition of the tensor with Tucker ranks $R_m < R$ for $m = 1,\ldots,M$. The core is then no longer diagonal but it will have many zero entries if $R < R_1 + \cdots + R_M$.

In any case, for $m = 1,\ldots,M$, the factor matrix $\mathbf A^{(m)}$ in the CP format, or the leaf $\bU^{(m)}$ in the Tucker format, spans the subspace of $\mathbb R^{I_m}$ that the data of the $m$-th mode lies in (the columns of the $m$-mode matricization of the full tensor also span this subspace). The information that is stored in the Tucker tensor (and by inclusion also in the CP tensor) is therefore twofold: What are the subspaces that our data lies in? This information is stored in the orthogonalized leafs of the Tucker tensor. And what combination of feature vectors is present (and to what degree) in the data? This information is stored in the core tensor $\teng$.

This observation can be exploited when designing a tensor kernel for KSTM. The subspace kernel introduced in Section <ref> can only see the subspaces, i.e. the leafs of the Tucker tensor. DuSK includes the core data via the norm equilibration. The Gaussian kernel uses the whole tensor but it is less sensible in spotting the subspaces, because the tensor is simply vectorized and this information is hidden. Our proposed WSEK includes information of the core tensor because the feature vectors are weighted with the $p$-th power of the corresponding singular values.

§.§ A Synthetic Experiment

We substantiate our considerations by creating two artificial experimental scenarios: In one (the leaf-scenario), all the information necessary for classification is stored in the leafs (i.e. the subspaces) and in the other (the core-scenario), all the information is in the core of the Tucker tensor. We then test the performance of the aforementioned kernels on this data for different noise levels and tensor ranks.

The detailed experimental setup is as follows: Let $M = 3$ and $I_1 = I_2 = I_3 = 100$. For different Tucker ranks $R_1 = R_2 = R_3 = r_{approx} = 1,\ldots,10$, we simulate the approximation of a tensor with Tucker ranks $R_1' = R_2' = R_3' = r_{exact} = 3$ plus some noise. That is, the core $\teng$ of the simulated tensor will have size $r_{approx} \times r_{approx} \times r_{approx}$ and the leafs will have sizes $I_m \times r_{approx}$ for $m = 1,2,3$. The core consists of random noise that is normally distributed with mean 0 and variance $\vartheta^2$ (the noise level to be chosen later). To the small upper-left-and-foremost cube of size $\min(r_{approx},3) \times \min(r_{approx},3) \times \min(r_{approx},3)$, we add the information tensor. In the core-scenario, this is the same tensor for all samples in the same class, generated by drawing the entries from a standard normal distribution. In the leaf-scenario, the information tensor is different for all samples (also drawn from the standard normal distribution).

The leafs of size $I_m \times r_{approx}$ also consist of random noise (normally distributed with mean 0 and variance $\vartheta^2$) and to the first $\min(r_{approx},3)$ columns we add vectors $\cos(\pi*\nu*{\mathbf v})$, where $\nu \in \mathbb R^{100}$ is a uniform discretization of $[-1,1]$ and the frequency $\nu$ is chosen uniformly at random (with mean 0 and variance 1). In the leaf-scenario, these frequencies are the same for all samples in the same class, and in the core-scenario, these frequencies are different for all samples. After the construction, the leafs are orthogonalized (using the QR-decomposition and discarding the $R$-matrix), so that the resulting Tucker tensor is already in the form of a HOSVD.

The reasoning is that these two scenarios yield Tucker tensors of rank $r_{approx}$ that are approximations of noisy tensors of rank $r_{exact}$, and the cluster information is stored exclusively in either the core or the leafs. We generate 100 samples in two classes with 50 samples each for each noise level $\vartheta^2 = 0.01,0.02,0.05,0.1,0.2,0.5,1$ and Tucker ranks $r_{approx} = 1,3,5,10$. We then perform the SVM 20 times with 5-fold cross validation in order to determine the hyperparameters $C \in \{2^{-8 :1:8} \}$ (the soft-margin parameter) and $g \in \{2^{-4 :1:12} \}$ (the variance parameter in the Gaussian kernel).

Classification Accuracy with 95% confidence interval for different noise levels in the core-scenario of the generated synthetic data with different rank trucation $r_{approx} = 1,3,5,10$ (DuSK not included for $r_{approx} = 10$)

§.§ Results and Interpretation

In Fig. <ref> and Fig. <ref>, we can see the results for the two experimental scenarios. In both cases, we plotted the test error (i.e. the classification accuracy on the test set) for the ranks $r_{approx} = 1,3,5,10$. For each rank, we plotted the accuracies of the different kernels for all noise levels in one picture. The computation of DuSK was too expensive in the case $r_{approx} = 10$ as this requires the computation of norms of tensors with CP rank 1000 many times.

In the core-scenario, the Gaussian kernel outperforms both DuSK and the Subspace kernel. The WSEK performs similarly to the Gaussian kernel. In the leaf-scenario, the Subspace kernel gives 100% accuracy in all cases, and all but the Gaussian kernel performed very well on this data.

These results are not surprising following the considerations in Sec. <ref>: The Subspace kernel sees only the information in the leafs and it can therefore not perform well in the core-scenario. The DuSK kernel includes extra information so that it performs better in the core-scenario (especially when we guessed the rank correctly, $r_{approx} = r_{exact} = 3$) but still not as good as the Gaussian kernel. In the leaf-scenario, all the information is in the subspaces and therefore both DuSK and the Subspace kernel perform very well. Here, the Gaussian kernel performed much worse than all other kernels. The WSEK was designed to do well in both scenarios and this is shown also in these experiments.

Classification Accuracy with 95% confidence interval for different noise levels in the leaf-scenario of the generated synthetic data with different rank trucation $r_{approx} = 1,3,5,10$ (DuSK not included for $r_{approx} = 10$)

We will see in the next chapter that especially the leaf-scenario is realistic: In the ADNI dataset explored below, the subspace kernel performs better than the Gaussian kernel and it even outperforms DuSK. We conclude that including information on the $m$-mode subspaces is crucial for the design of a tensor kernel. This is also a possible explanation why DuSK has performed well in many settings. Our new kernel however outperforms DuSK in all our experiments.

§ CLASSIFICATION OF REAL DATASETS

In this section, we test the performance of the discussed tensor kernels on two real world datasets: The ADNI dataset contains fMRI images of patients with and without Alzheimer's disease and the ADHD dataset contains fMRI images of ADHD patients and healthy subjects. We use the KSTM with the different kernels to distinguish the two classes of subjects in each dataset.

All numerical experiments have been done in . Low-rank tensor approximations are computed using  <cit.> and  <cit.>. We have run all experiments on a compute cluster which is equipped with  Harddisk, $2\times$ CPUs with 10 cores per CPU, and  of RAM.The hyperparameters in KSTM (<ref>) are tuned similarly to Synthetic experiments (see Section <ref>), the only difference is that we report results for real data experiments for each of the rank  $R \in \{1,2, \cdots, 10\}$, where $R_1 = R_2 = R_3 = R$. The SVM problem (<ref>) is solved using the  function from LIBSVM <cit.> library, and  computes the classification accuracy using  (<ref>) and (<ref>).

(left) Comparison of mean classification accuracy with variance for the different kernels with different rank truncation for ADNI dataset, (right) Comparison of mean classification accuracy with variance for the different kernels with different rank truncation for ADHD dataset.

§.§ Resting-state fMRI data collection

* Alzheimer Disease (ADNI): The ADNI[<http://adni.loni.usc.edu/>] stands for Alzheimer Disease Neuroimaging Initiative. It contains the resting state fMRI images of 33 subjects. The dataset was collected from the authors of <cit.>. The images belong to either Mild Cognitive Impairment (MCI) with Alzheimer Disease (AD), or normal controls. Each image is a tensor of size $61 \times 73 \times 61$, containing $271633$ elements in total. The AD+MCI images are labeled with $-1$, and the normal control images are labeled with $1$. Preprocessing of the data sets is explained in <cit.>.

* Attention Deficit Hyperactivity Disorder (ADHD): The ADHD dataset is collected from the ADHD-200 global competition dataset[<http://neurobureau.projects.nitrc.org/ADHD200/Data.html>]. It is a publicly available preprocessed fMRI dataset from eight different institutes. The original dataset is unbalanced, so we have chosen 200 subjects randomly, ensuring that 100 of them are ADHD patients (assigned the classification label $-1$) and the 100 other subjects are healthy (denoted with label $1$). Each of the 200 resting state fMRI samples contains $49  \times 58 \times 47  = 133574$ voxels.

The dataset taken here is exactly the same as that used in <cit.>, so the TT-MMK results can be compared one to one. However, the particular indices of the collected data are not similar to those selected in <cit.>, so the accuracy of DuSK reported below is not directly comparable to that in <cit.>.

§.§ Numerical results

In this section, we summarize the results for the two fMRI datasets:

* Classification accuracy: In Fig. <ref>, we show the average classification accuracy resulting from the cross validation. In Table <ref>, we show the best classification accuracy between rank $\left[1,10\right]$. On both datasets the proposed WSEK gives the best average classification accuracy
($79\%$ for ADNI and $64\%$ for ADHD) compared to other state of the art tensor kernels. We note that the accuracy of the subspace kernel improves for higher ranks in the ADHD dataset and is then similar to that of WSEK but this choice of rank is very high for a low-rank truncation method. On the other hand, our proposed kernel gives good classification accuracy already at rank 2. The Gaussian kernel was computed for different Tucker approximations of the full tensor. Using the full tensor in the computation of the kernel did not improve the accuracy. As in <cit.>, the DuSK kernel was computed using a CP approximation of the full tensor (CP-DuSK), because computing the Tucker decomposition and then converting to CP yielded high CP ranks and DuSK was too slow.

Maximum average classification accuracy in percentage $\pm$ standard deviation for different methods, data sets, and rank $R \in \left[1,10\right]$. The values for TTCP-DuSK are taken from <cit.> for comparison.

Methods     ADNI     ADHD

Gauss      53        50

CP-DuSK     64 $\pm$ 0.05 (R = 5)     58 $\pm$ 0.02 (R = 6)

TTCP-DuSK     73 $\pm$ 0.03 (R = 4)     63 $\pm$ 0.01 (R = 5)

Tucker Subspace     $73 \pm 0.04$ (R = 7)     61 $\pm 0.02$ (R = 10),

WSEK     79 $\pm 0.03$ (R = 8)      64 $\pm \mathbf{0.01}$ (R = 2)

* Running Time: Table <ref> shows the running times for the computation of the different kernels on the ADNI dataset. For small tensors, computing the Gauss kernel is fast. Both the subspace kernel and WSEK can be computed efficiently. Computation of DuSK however quickly becomes prohibitive and takes a long time in these experiments.

Comparison of Kernel computation time for $ R\in \left[ 1,10 \right]$. The values for TTCP-DuSK are taken from <cit.> for comparison.

Kernel     Format     Parameters      CPUtime for ADNI (# run = 1)

Gaussian      Ktensor      $C, g$        $20$ seconds

DuSK     CP (<ref>)       $C, g$,  $R$       $17$ minutes

DuSK     TTCP (<ref>)        $C, g, ~R_{TT}$        $ 3.5$ hours

Subspace     Tucker (<ref>)      $C, g,~ R_{Tucker}$       $25$ seconds

WSEK     SqrtmHOSVD (<ref>)      $C, g,~ R_{Tucker}$       $30$ seconds

* Statistic comparison: In Tabel <ref> and Figure <ref>, the variance for corresponding mean accuracy is shown. The WSEK STM shows a good trade-off between classification accuracy and variance. For the ADHD dataset, it even gives the best classfication accuracy with the lowest variance value for small ranks.

§.§ Conclusion

Our real world experiments show superior performance of the WSEK both in terms of classification accuracy and running time. We conclude that the classification information of the datasets is hidden mostly in the subspaces of the Tucker decomposition (hence the previously observed good performance of DuSK) but that classification can be improved by taking the singular values into account (as done in WSEK). Furthermore, computing the Tucker decomposition of tensor inputs is straightforward and efficient, resulting in an all-around very robust tensor kernel.

§ AUTHOR CONTRIBUTIONS

All authors contributed to the conception and design of the tensor kernel method. Implementation and analysis of the new kernel were performed by Kirandeep Kour and Max Pfeffer. The first draft of the manuscript was written by Max Pfeffer and Kirandeep Kour. All authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

§ ACKNOWLEDGMENTS

K.K. is supported by the International Max Planck Research School for Advanced Methods in Process and System Engineering-IMPRS ProEng, Magdeburg. M.P. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 448293816.

§.§ Statements and Declaration

Conflict of interest The authors declare that they have no conflict of interests. Also, there are no financial and non-financial interest to disclose.

, Kolda

Acar, E.
, Kolda, T.G.
Dunlavy, D.M.

All-at-once Optimization for Coupled Matrix and Tensor
Factorizations All-at-once optimization for coupled matrix and tensor

tensortoolbox Version 3.4
tensortoolbox Version 3.4

Bader, B.W.
, Kolda, T.G.
Tensor Toolbox for MATLAB Version 3.4. Tensor toolbox for
matlab version 3.4.

, He
, Wen
, Han

Cai, D.
, He, X.
, Wen, JR.
, Han, J.
Ma, WY.

Support tensor machines for text categorization
Support tensor machines for text categorization.
The University of Illinois at Urbana-Champaign Computer
Science Department.

Chang, CC.
Lin, CJ.

LIBSVM: A Library for Support Vector Machines Libsvm:
A library for support vector machines.
ACM Transactions on Intelligent Systems and

, Batselier
, Ko

Chen, C.
, Batselier, K.
, Ko, CY.
Wong, N.

A support tensor train machine A support tensor train
2019 International Joint Conference on Neural Networks
(IJCNN) 2019 International Joint Conference on Neural Networks (IJCNN)	( 1–8).

Cichocki, A.

Tensor decompositions: new concepts in brain data
analysis? Tensor decompositions: new concepts in brain data
Journal of the Society of Instrument and Control

Cichocki, A.
, Lee, N.
, Oseledets, I.
, Phan, AH.
, Zhao, Q.
Mandic, D.P.

Tensor Networks for Dimensionality Reduction and
Large-scale Optimization: Part 1 Low-Rank Tensor Decompositions Tensor
networks for dimensionality reduction and large-scale optimization: Part 1
low-rank tensor decompositions.
FNT in Machine Learning94-5249-429.

Cortes, C.
Vapnik, V.

Support-vector networks Support-vector
Machine Learning20273-297.

de Silva

de Silva

de Silva, V.
Lim, LH.

Tensor rank and the ill-posedness of the best low-rank
approximation problem Tensor rank and the ill-posedness of the best
low-rank approximation problem.
SIAM Journal on Matrix Analysis and

, Kotsia

Guo, W.
, Kotsia, I.
Patras, I.

Tensor Learning for Regression Tensor learning for
IEEE Transactions on Image Processing212816-827.

, He
, Chen

Hao, Z.
, He, L.
, Chen, B.
Yang, X.

A Linear Support Higher-Order Tensor Machine for
Classification A linear support higher-order tensor machine for
IEEE Transactions on Image

He, L.
, Kong, X.
, Yu, P.S.
, Yang, X.
, Ragin, A.B.
Hao, Z.

DuSK: A Dual Structure-preserving Kernel for Supervised
Tensor Learning with Applications to Neuroimages Dusk: A dual
structure-preserving kernel for supervised tensor learning with applications
to neuroimages.
Proceedings of the 2014 SIAM International Conference
on Data Mining (SDM)127-135.

, Lu
, Ding
, Lu
, Ding

He, L.
, Lu, CT.
, Ding, H.
, Wang, S.
, Shen, L.
, Yu, P.S.
Ragin, A.B.

Multi-Way Multi-Level Kernel Modeling for Neuroimaging
Classification Multi-way multi-level kernel modeling for neuroimaging
The IEEE Conference on Computer Vision and Pattern
Recognition (CVPR)6846-6854.

, Lu
, Ma
, Lu
, Ma

He, L.
, Lu, CT.
, Ma, G.
, Wang, S.
, Shen, L.
, Yu, P.S.
Ragin, A.B.

Kernelized support tensor machines Kernelized support
tensor machines.
Proceedings of the 34th International Conference on
Machine Learning-Volume 701442–1451.

Hitchcock, F.L.

The Expression of a Tensor or a Polyadic as a Sum of
Products The expression of a tensor or a polyadic as a sum of
Journal of Mathematics and Physics61-4164-189.

Kolda, T.G.
Bader, B.W.

Tensor Decompositions and Applications Tensor
decompositions and applications.
SIAM Review513455–500.

Kotsia, I.
Patras, I.

Support Tucker Machines Support Tucker
CVPR 2011633-640.

, Dolgov
, Stoll

Kour, K.
, Dolgov, S.
, Stoll, M.
Benner, P.

Efficient Structure-preserving Support Tensor Train
Machine Efficient structure-preserving support tensor train
Journal of Machine Learning Research2441–22.

, Moor

Lathauwer, L.D.
, Moor, B.D.
Vandewalle, J.

A Multilinear Singular Value Decomposition A
multilinear singular value decomposition.
SIAM Journal on Matrix Analysis and

, Guo
, He

Liu, X.
, Guo, T.
, He, L.
Yang, X.

A Low-Rank Approximation-Based Transductive Support
Tensor Machine for Semisupervised Classification A low-rank
approximation-based transductive support tensor machine for semisupervised
IEEE Transactions on Image

Nion, D.
Lathauwer, L.D.

Fast Communication: An Enhanced Line Search Scheme for
Complex-valued Tensor Decompositions Fast communication: An enhanced line
search scheme for complex-valued tensor decompositions.
Signal Processing883749–755.

Oseledets, I.
Dolgov, S.

TT-Toolbox. Tt-toolbox.

I. Oseledets

I. Oseledets

Oseledets, I.
Tyrtyshnikov, E.

TT-cross approximation for multidimensional arrays
Tt-cross approximation for multidimensional arrays.
Linear Algebra and its Applications432170-88.

I.V. Oseledets

I.V. Oseledets

Oseledets, I.V.

Tensor-train decomposition Tensor-train
SIAM Journal on Scientific

, Ramanan

Pirsiavash, H.
, Ramanan, D.
Fowlkes, C.C.

Bilinear classifiers for visual recognition Bilinear
classifiers for visual recognition.
Advances in Neural Information Processing Systems

, Olivetti
, Lathauwer

Signoretto, M.
, Olivetti, E.
, Lathauwer, L.D.
Suykens, J.A.K.

A kernel-based framework to tensorial data analysis A
kernel-based framework to tensorial data analysis.
Neural Networks248861 - 874.

, Olivetti
, Lathauwer

Signoretto, M.
, Olivetti, E.
, Lathauwer, L.D.
Suykens, J.A.K.

Classification of Multichannel Signals With
Cumulant-Based Kernels Classification of multichannel signals with
cumulant-based kernels.
IEEE Transactions on Signal

, Tran Dinh
, De Lathauwer

Signoretto, M.
, Tran Dinh, Q.
, De Lathauwer, L.
Suykens, J.A.K.

Learning with tensors: a framework based on convex
optimization and spectral regularization Learning with tensors: a framework
based on convex optimization and spectral regularization.
Machine Learning943303–351.

Taguchi, YH.
Turki, T.

Mathematical formulation and application of kernel
tensor decomposition based unsupervised feature extraction Mathematical
formulation and application of kernel tensor decomposition based unsupervised
feature extraction.
Knowledge-Based Systems217106834.

, Li
, Hu
, Maybank

Tao, D.
, Li, X.
, Hu, W.
, Maybank, S.
Wu, X.

Supervised Tensor Learning Supervised tensor
Proceedings of the Fifth IEEE International Conference
on Data Mining450-457.

, Li
, Hu
, Maybank

Tao, D.
, Li, X.
, Hu, W.
, Maybank, S.
Wu, X.

Supervised tensor learning Supervised tensor
Knowledge and Information Systems1-42.

Tucker, L.R.

Some mathematical notes on three-mode factor analysis
Some mathematical notes on three-mode factor analysis.

, Vandebril

Vannieuwenhoven, N.
, Vandebril, R.
Meerbergen, K.

A New Truncation Strategy for the Higher-Order Singular
Value Decomposition A new truncation strategy for the higher-order singular
value decomposition.
SIAM Journal on Scientific

Vapnik, V.

The nature of statistical learning theory The nature of
statistical learning theory.
New YorkSpringer-Verlag.

Vapnik, V.

Statistical Learning Theory Statistical learning theory.

, Jhuang

Wolf, L.
, Jhuang, H.
Hazan, T.

Modeling Appearances with Low-Rank SVM Modeling
appearances with low-rank SVM.
IEEE Conference on Computer Vision and Pattern

Yan, S.
, Xu, D.
, Yang, Q.
, Zhang, L.
, Tang, X.
Zhang, HJ.

Multilinear Discriminant Analysis for Face Recognition
Multilinear discriminant analysis for face recognition.
IEEE Transactions on Image Processing161212-220.

, Wang
, Shen

Zeng, D.
, Wang, S.
, Shen, Y.
Shi, C.

A GA-based feature selection and parameter
optimization for support tucker machine A GA-based feature selection and
parameter optimization for support tucker machine.
Procedia Computer Science11117 - 23.

, Zhou
, Adali
, Zhang

Zhao, Q.
, Zhou, G.
, Adali, T.
, Zhang, L.
Cichocki, A.

Kernelization of Tensor-Based Models for Multiway Data
Analysis: Processing of Multidimensional Structured Data Kernelization of
tensor-based models for multiway data analysis: Processing of
multidimensional structured data.
IEEE Signal Processing Magazine304137-148.

, Zhou
, Adali
, Zhang

Zhao, Q.
, Zhou, G.
, Adali, T.
, Zhang, L.
Cichocki, A.

Kernelization of Tensor-Based Models for Multiway Data
Analysis: Processing of Multidimensional Structured Data Kernelization of
tensor-based models for multiway data analysis: Processing of
multidimensional structured data.
IEEE Signal Processing Magazine304137-148.

, Li

Zhou, H.
, Li, L.
Zhu, H.

Tensor Regression with Applications in Neuroimaging Data
Analysis Tensor regression with applications in neuroimaging data
Journal of the American Statistical


# Federated Double Deep Q-learning for Joint Delay and Energy Minimization in
IoT networks

Sheyda Zarandi and Hina Tabassum, Senior Member, IEEE S. Zarandi and H.
Tabassum are with the Lassonde School of Engineering at York University, ON,
Canada (e-mail<EMAIL_ADDRESS>hinat@yorku.ca). This work is supported by the
MITACS research trainee award and Discovery Grant from the Natural Sciences
and Engineering Research Council of Canada.

###### Abstract

In this paper, we propose a federated deep reinforcement learning framework to
solve a multi-objective optimization problem, where we consider minimizing the
expected long-term task completion delay and energy consumption of IoT
devices. This is done by optimizing offloading decisions, computation resource
allocation, and transmit power allocation. Since the formulated problem is a
mixed-integer non-linear programming (MINLP), we first cast our problem as a
multi-agent distributed deep reinforcement learning (DRL) problem and address
it using double deep Q-network (DDQN), where the actions are offloading
decisions. The immediate cost of each agent is calculated through solving
either the transmit power optimization or local computation resource
optimization, based on the selected offloading decisions (actions). Then, to
enhance the learning speed of IoT devices (agents), we incorporate federated
learning (FDL) at the end of each episode. FDL enhances the scalability of the
proposed DRL framework, creates a context for cooperation between agents, and
minimizes their privacy concerns. Our numerical results demonstrate the
efficacy of our proposed federated DDQN framework in terms of learning speed
compared to federated deep Q network (DQN) and non-federated DDQN algorithms.
In addition, we investigate the impact of batch size, network layers, DDQN
target network update frequency on the learning speed of the FDL.

##  I Introduction

Massive connectivity is among one of the most challenging requirements of
Internet-of-Things (IoT) networks which necessitates efficient, scalable, and
low-complexity network resource management. Furthermore, due to limited
computation and battery capacity of the IoT devices, it is often impossible
for them to process their resource-intensive tasks within a predefined
deadline. In the sequel, mobile cloud computing (MCC) and mobile edge
computing (MEC) enable IoT devices to offload their tasks to the cloud or edge
servers to access their substantial processing capabilities at the expense of
having to transmit the tasks over dynamic wireless channels. Subsequently, to
take full advantage of the MCC and MEC paradigms, it becomes essential to
carefully optimize offloading decisions, communication, and computation
resources.

Most of the existing research works solved the joint offloading decision,
communication, and computation resource allocation problem leveraging on tools
from optimization theory [1, 2]. However, the algorithms were typically non-
scalable, time-consuming, and computationally expensive. Unlike optimization
frameworks, deep reinforcement learning (DRL) enables agents to learn by
interacting with the environment. This unique approach to learning, turns DRL
into an ideal problem-solving tool in dynamic environments. Yet, most DRL
algorithms are centralized and thus suffer from lack of scalability when the
number of devices grow. Also, the computational complexity of finding an
optimal policy may increase exponentially as the state space and action space
grow. Furthermore, the centralized learning requires IoT devices to share
their information in order to train the global model which may violate their
privacy and create unnecessary communication overhead on the already scarce
frequency spectrum.

Recently, federated learning (FDL) has emerged as a new paradigm for
cooperative learning, where multiple nodes contribute in training a single
global model. The devices use their local datasets to train and then offload
their local models to the central unit for global aggregation. FDL enhances
the cooperation between agents and scalability of the network resource
management algorithms. Furthermore, FDL does not require local agents to share
their data with any external entity, thereby preserves the privacy of each
agent [3].

To date, several research works have explored the problem of minimizing delay
and energy consumption of the IoT devices considering a FDL system [4, 5, 6].
Specifically, in the aforementioned research works, FDL features were
incorporated in the problem formulation and the problem was then solved using
traditional optimization theory. Similarly, [7, 8, 9] focused on optimizing
different aspects of FDL, such as compression of weights, convergence
analysis, reduction in the number of iterations, and incentive mechanisms. In
these research works, the capabilities of FDL as a part of resource allocation
solution approach were not investigated. In [10], the problem of computation
resource allocation was addressed considering an FDL system. However, FDL is
only considered to formulate an optimization problem which is later on solved
by using a centralized actor-critic agent and without using FDL in the
solution approach.

None of the aforementioned research works applied FDL to enhance the efficacy
of solving a realistic wireless resource allocation problem. Very recently,
[11, 12] adopted FDL to facilitate the learning process in DRL, i.e., local
DRL models were trained and then integrated together to cooperatively develop
a comprehensive global DRL model. However, in [11], a cooperative caching
scheme was proposed and offloading decisions were not considered. In [12],
computational offloading was considered; however, the network was modeled as a
queuing system, transmit power was modeled as an integer variable whose
maximum value is equal to the maximum length of the energy queue. Also, in
[12], no explicit quality-of-service (QoS) was guaranteed for users’ tasks and
computation resource allocation was overlooked.

In this paper, we employ federated DRL (FedRL) to address the problem of
minimizing joint expected task completion delay and energy consumption of IoT
devices with offloading decisions, computation resources, and transmit powers
as variables. Considering the mixed-integer non-linear programming programming
(MINLP) nature of the problem, we first reformulate our problem as a multi-
agent DRL problem and solve it using double deep Q-network (DDQN). In this
problem, offloading decisions would be the actions and the immediate cost is
calculated through solving either the transmit power or local computation
resource optimization. To improve learning quality and speed of DRL, we
incorporate FDL at the end of each episode. Using FDL results in a privacy-
preserving and scalable framework and creates a context for cooperation
between agents. Our numerical results demonstrate the efficacy of our
federated DDQN framework in terms of learning speed compared to federated deep
Q-network (DQN) and non-federated DDQN algorithms.

The rest of this paper is organized as follows: Section II describes the
system model and the problem is formulated in Section III. Then, in Section
IV, our proposed algorithm is presented. Finally, simulation results and
related discussions are provided in Section V followed by conclusions in
Section VI.

##  II System model and Assumptions

We consider a network containing one MEC server, one cloud server, and a set
of $\mathcal{N}=\\{1,...,N\\}$ IoT devices with limited computation and energy
resources. We consider a given time horizon $\mathcal{T}$ which is divided
into $T$ time steps. At each time $t$, device $i$ needs to process one of the
tasks in its queue, defined with the tuple $(L_{i,t},C_{i,t},\hat{T}_{i,t})$,
where $L_{i,t}$ is the size of the task (in bits), $C_{i,t}$ denotes the CPU
cycle requirement of the task, and $\hat{T}_{i,t}$ denotes the maximum delay
threshold of the task. At any time $t$, devices can either execute their task
locally or offload it to edge or cloud server.

Let us denote local offloading decision of device $i$ at time $t$ as
$x_{i,t}$, where $x_{i,t}=1$ means the task would be performed locally, and
$x_{i,t}=0$, otherwise. Similarly, we define MCC and MEC offloading variable
of device $i$ by $z_{i,t}$ and $y_{i,t}$, respectively. If device $i$ offloads
its task to the cloud $z_{i,t}=1$ and if it offloads the task to the edge
server $y_{i,t}=1$. As a binary offloading decision is considered, we have:

$x_{i,t}+y_{i,t}+z_{i,t}=1.~{}\forall{t}\in\mathcal{T}.$ (1)

When device $i$ offloads its task (whether to MEC server or to the cloud), the
delay and energy consumption would depend on the channel condition, the size
of the task, and the power with which the device transmits its task and, in
the case of local computation they depend on computation resource utilization.
In what follows, we model the delay and energy consumption that an IoT device
would experience, given its offloading decision.

If the device $i$ decides to offload its task, it should first transmit it to
the MEC-enabled base station (BS) through wireless channels. At time step $t$,
the transmission data rate of this user, denoted by $r_{i,t}$ is calculated as
follows:

$r_{i,t}=B\log_{2}(1+\frac{p_{i,t}h_{i,t}}{\sigma^{2}}),$ (2)

where $B$ and $p_{i,t}$ denote the bandwidth and transmit power of device $i$
at time step $t$, respectively. Also, $h_{i,t}$ and $\sigma^{2}$ represent the
path-gain of device $i$ at time $t$ and the receiver noise. Thus, the
communication delay and energy consumption of device $i$, while offloading is
given, respectively, as follows:

$T_{i,t}^{\mathrm{comm}}=\frac{L_{i,t}}{r_{i,t}},$ (3)
$E_{i,t}^{\mathrm{comm}}=p_{i,t}T_{i,t}^{\mathrm{comm}}=\frac{p_{i,t}L_{i,t}}{B\log_{2}(1+\frac{p_{i,t}h_{i,t}}{\sigma^{2}})}$
(4)

If device $i$ offloads its task to the edge server the computation delay would
be $T^{e}_{i,t}=\frac{C_{i,t}}{F^{e}},$ where $F^{e}$ denotes the average
computation capacity of edge server. Also, if device $i$ offloads the task to
cloud server, the computation delay would be
$T^{c}_{i,t}=\frac{C_{i,t}}{F^{c}},$ where $F^{c}$ represents the average
computation capacity of the cloud server.

From the perspective of IoT device, the energy that is consumed for processing
a task when it is offloaded to either of the servers, is the energy spent on
the transfer of the task. Therefore, both cloud and edge computing energy
utilization at step $t$, denoted by $E^{c}_{i,t}$ and $E^{e}_{i,t}$, would be
equal to $E_{i,t}^{\mathrm{comm}}$.

If device $i$ chooses to perform its task locally, the local computation delay
and energy consumption would depend on the amount of computation resource
allocated to process the task at time $t$, which we denote by $f_{i,t}^{L}$.
Thus, the local delay and energy consumption of the device is modeled as
follows:

$T_{i,t}^{L}=\frac{C_{i,t}}{f_{i,t}^{L}},~{}~{}E_{i,t}^{L}=\kappa_{i}(f_{i,t}^{L})^{2},$
(5)

where $\kappa$ is a constant coefficient that depends on the chip architecture
in devices. Note that higher resource utilization (transmit power or
computation capacity) , decreases the task completion delay at the expense of
increased energy consumption. Therefore, this trade-off must be carefully
managed through efficient offloading decision making and precise optimization
of $f_{i,t}$ and $p_{i,t}$ in the case of local computation and offloading,
respectively.

##  III Multi-objective Problem Statement

In this section, we formulate the multi-objective problem of jointly
minimizing the long-term delay and energy consumption of an IoT device in a
decentralized manner over a specified time horizon $\mathcal{T}$. The long-
term expected cost (weighted sum of delay and energy consumption) for each IoT
device $i$ is formulated, respectively, as follows:

$\displaystyle T_{i}(\mathbf{p}_{i},$
$\displaystyle\mathbf{f}_{i},\mathbf{x}_{i},\mathbf{y}_{i},\mathbf{z}_{i})=\mathbb{E}[\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{j=0}^{t}\left(x_{i,j}T_{i,j}^{L}+\right.$
$\displaystyle\left.y_{i,j}(T^{e}_{i,j}+T_{i,j}^{\mathrm{comm}})+z_{i,j}(T_{i,j}^{c}+T_{i,j}^{\mathrm{comm}}+\Psi)\right)],$
(6) $\displaystyle
E_{i}(\mathbf{p}_{i},\mathbf{f}_{i},\mathbf{x}_{i},\mathbf{y}_{i},\mathbf{z}_{i})$
$\displaystyle=\mathbb{E}[\lim_{t\rightarrow\infty}\frac{1}{t}\sum_{j=0}^{t}\left(x_{i,j}E_{i,j}^{L}+y_{i,j}E^{e}_{i,j}+z_{i,j}E^{c}_{i,j}\right)],$
(7)

where $\mathbf{p}_{i}=[p_{i,1},p_{i,2},\cdots,p_{i,t}]$,
$\mathbf{f}_{i}=[f_{i,1},f_{i,2},\cdots,f_{i,t}]$,
$\mathbf{x}_{i}=[x_{i,1},x_{i,2},\cdots,x_{i,t}]$,
$\mathbf{y}_{i}=[y_{i,1},y_{i,2},\cdots,y_{i,t}]$, and
$\mathbf{z}_{i}=[z_{i,1},z_{i,2},\cdots,z_{i,t}]$, represent the vectors of
transmit powers, computation resource allocation, local computing, edge
offloading, and cloud offloading decision of device $i$, respectively. As
cloud server is generally located far from the IoT devices, the delay of
accessing cloud is commonly more than the delay of offloading to the edge
server which is located at the edge of the network. In the equation (III),
$\Psi$ denotes the delay of accessing the cloud server, including the time
necessary to transfer the task from BS to the cloud, the possible routing in
the path, and the response delay.

For any given device $i$ at time $t$, we model our problem as:

$\displaystyle\min_{\mathbf{p}_{i},\mathbf{f}_{i},\mathbf{x}_{i},\mathbf{y}_{i},\mathbf{z}_{i}}~{}T_{i}+\lambda_{i}E_{i}$
(8) subject to:
$\displaystyle\mathrm{C1}:f_{i,t}^{L}<=F_{\mathrm{max},i},~{}~{}\forall
t\in\mathcal{T},$
$\displaystyle\mathrm{C2}:x_{i,t}E^{L}_{i,t}+y_{i,t}E^{e}_{i,t}+z_{i,t}E^{c}_{i,t}<=E_{\mathrm{max},i},~{}~{}\forall
t\in\mathcal{T},$ $\displaystyle\mathrm{C3}:T_{i,t}<=\hat{T}_{i,t},~{}\forall
t\in\mathcal{T},$
$\displaystyle\mathrm{C4}:x_{i,t}+y_{i,t}+z_{i,t}=1,~{}~{}\forall
i\in\mathcal{N},\forall t\in\mathcal{T},$
$\displaystyle\mathrm{C5}:x_{i,t},y_{i,t},z_{i,t}\in\\{0,1\\},~{}~{}\forall
i\in\mathcal{N},\forall t\in\mathcal{T}.$

In the above optimization problem, $\lambda_{i}$ is a weighting factor whose
value should be carefully selected based on the heterogeneity of resources
available at each individual IoT device. If device $i$ is more restricted in
the energy resource compared to computation resource, the value of
$\lambda_{i}$ should be set to a larger number. Otherwise, $\lambda_{i}$
should be a small number. Furthermore, constraint C1 indicates the local
computation capacity with the maximum threshold $F_{\mathrm{max},i}$.
Constraint C2 represents the restriction on the energy resource of the device
and that energy utilization should not exceed $E_{\mathrm{max},i}$.
Furthermore, constraints C4 and C5 define the binary offloading scheme adopted
in this paper. It can be proven that both equations (III) and (III) are convex
with respect to the variables $p_{i,t}$ and $f_{i,t}$, respectively. However,
with binary offloading variables ($x_{i}$, $y_{i}$, and $z_{i}$) included, (8)
turns into a MINLP that cannot be solved in an acceptable time span.

##  IV Proposed Federated DDQN Algorithm

To solve (8) at each IoT device, we propose a DDQN algorithm, and solve the
problem in the following two phases:

$\bullet$ Offloading Decision Optimization: Since each IoT device has three
options to process a task (namely local, edge server, or cloud server
computing), there are almost $2^{3N}$ possible offloading options (from the
perspective of a centralized controller) at each given time step. As the
number of devices increases, this complexity would also surge exponentially.
To address this problem, we apply a multi-agent DDQN framework where each IoT
device would train their local DDQN models using their local data.

$\bullet$ Computing and Communication Resource allocation: Given the
offloading decision, we optimize computation capacity or transmit power of the
devices to minimize the weighted sum of energy consumption and delay. We use
optimization theory to address this part of the problem and then feed the
results into the DDQN framework as the immediate cost function. In this way,
we provide the learning agent with a real sense of the quality of the adopted
offloading policy that reflects many important aspects of the system model
(such as limitation of resources in each device and their QoS demands).

After the DDQN agent is trained through the above mentioned process for one
training round, we apply a federated learning framework where each IoT device
will train its DDQN models, share their models with the centralized
controller, and update their models to the central aggregating unit. This
mechanism is detailed in the flowchart provided in Fig. 1.

In what follows, we first focus on developing local models through DDQN
algorithm and then explain how FDL would be deployed.

###  IV-A Double deep Q-network for offloading decision making

In the first step of our algorithm, we model our problem as a multi-agent DDQN
problem. For each device (DRL agent), we have following components:

  1. 1.

State space: the state space for each agent $i$, denoted by $s_{i}$, consists
of the following components: the length of the task queue of device $i$ (tasks
that are not yet processed or are not successfully processed, would be kept in
this queue) which is denoted by $\mathcal{L}_{i,t}$, the path gain of the IoT
device $h_{i,t}$, the size of the task currently being processed $L_{i,t}$,
its CPU cycle requirement, $C_{i,t}$, and available resources. Thus,
$s_{i}=\\{\mathcal{L}_{i,t},h_{i,t},L_{i,t},C_{i,t},E_{\mathrm{max},i},F_{\mathrm{max},i}\\}$.
If a task is successfully processed under a given offloading decision policy
(its QoS requirement is satisfied), it would be removed from the task queue of
the device. Otherwise, it will remain at the top of the queue to be processed
under another offloading decision.

  2. 2.

Action space: The action space of agents, denoted by $\mathcal{A}$, contains
possible offloading decisions, i.e., whether to process the task locally or
offload.

  3. 3.

Cost: (8) suggests that the cost of an agent is equal to the weighted
summation of the delay and energy consumption given in the objective function.
The value of this objective function and thus the cost depends on the value of
$p_{i,t}$ in the case of offloading and $f_{i,t}$ if local computation is
selected. Therefore, to ensure that the cost function accurately reflects the
benefit of a given offloading decision, these variables should be carefully
optimized. To this end, when local computation is selected ($x_{i,t}=1$), the
cost would be calculated by solving the instantaneous optimization problem
below:

$\displaystyle~{}\min_{{f}_{i,t}}~{}T_{i,t}^{L}+\lambda E_{i,t}^{L}$ (9)
subject to: C1, C2, C3.

in the case offloading is selected ($y_{i,t}=1$ or $z_{i,t}=1$), the transmit
power would be optimized by solving the following optimization problem:

$\displaystyle~{}\min_{{p}_{i,t}}~{}y_{i,t}(T_{i,t}^{e}+T_{i,t}^{\mathrm{comm}}+\lambda
E_{i,t}^{\mathrm{comm}})$ (10)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+z_{i,t}(T_{i,t}^{c}+T_{i,t}^{\mathrm{comm}}+\Psi+\lambda
E_{i,t}^{\mathrm{comm}})$ subject to: C2, C3.

The design of state and cost function has a significant impact on the success
of DDQN in finding the optimal offloading policy, $\pi^{*}$. By using a multi-
agent approach, we are in fact limiting the state and action space and focus
on each device separately. Also, by modeling the cost function as an
optimization problem not only can we optimize the local resource utilization
and enforce system constraints, but also we can provide the agent with an
accurate estimation of the quality of an offloading decision. Note that it can
be easily proved that both (9) and (10) are convex single variable
optimization problems that can be solved using standard softwares.

Algorithm 1 Proposed Federated DDQN Algorithm

1: Initialize the global model $\theta^{\mathrm{global}}$, and set maximum FDL
iterations to $K$.

2: For each agent initialize online and target networks as:
$\theta_{i}^{\mathrm{online}}=\theta_{i}^{\mathrm{target}}=\theta^{\mathrm{global}}$.

3:While ($K\geq 0$):

4: Set
$\theta_{i}^{\mathrm{online}}=\theta_{i}^{\mathrm{target}}=\theta^{\mathrm{global}}$,
$\forall i\in\mathcal{N}$,

5: Select the set of participating devices, $\mathcal{I}$, based on (13),

6: For each device $i$ in $\mathcal{I}$ do:

7: For each time step $t$ if $|\mathcal{L}_{i,t}|>0$ do:

8: Interact with environment and calculate the cost using (9) or (10),

9: Save the experience in replay memory $\mathcal{M}_{i,t}$.

10: Train the local model on $\mathcal{M}_{i,t}$,

11: Transmit $\theta_{i}^{\mathrm{online}}$ to the central aggregation unit,

12: End For

13: End For

14: update $\theta^{\mathrm{global}}$ using (14).

Let us denote the immediate cost of each device $i$ obtained from the solution
of the above mentioned optimization process as $u_{i}(s,a)$. Using Bellman
equation, the action-state value is:

$Q_{i}({s,a})=u_{i}(s,a)+\gamma\sum_{s^{\prime}\in\mathcal{S}}P_{ss^{\prime}}(a)\max_{a^{\prime}\in\mathcal{A}}Q_{i}(s^{\prime},a^{\prime}),$
(11)

where $\mathcal{S}$, $P_{ss^{\prime}}(a)$, and $\gamma$ are the set of states,
the transition probability function, and the discount factor, respectively. To
overcome the need for having a full model of environment, calculating the
transition probability function, and to acquiring a more stable learning
process, DDQN is employed in this work. Each agent $i$ has two neural networks
working alongside each other, one called online network with parameters
$\theta_{i}^{\mathrm{online}}$ and the other called target network with
parameters $\theta_{i}^{\mathrm{target}}$. At each training iteration the
target value for training the online network in device $i$ is calculated as:

$L_{i}=u_{i}(s,a)+\gamma
Q_{i}(s^{\prime},\arg\max_{a^{\prime}\in\mathcal{A}}Q_{i}(s^{\prime},a^{\prime};\theta_{i}^{\mathrm{online}}),\theta_{i}^{\mathrm{target}})$
(12)

While $\theta_{i}^{\mathrm{online}}$ is updated at every iteration, the
frequency of change in $\theta_{i}^{\mathrm{target}}$ is typically much lower
and only once in every $f_{\mathrm{update}}$ rounds,
$\theta_{i}^{\mathrm{target}}$ would be set equal to
$\theta_{i}^{\mathrm{online}}$.

As discussed before, training a DRL agent in a centralized manner can lead to
critical issues related to scalability, agents’ privacy, and additional
communication overheads. On the other hand, training a DRL agent in a
distributed manner can impact the overall performance gains (e.g., an agent
might consume a longer time to train its model). As such, we consider FDL to
combine the benefits of both centralized and distributed learning. FDL enables
each agent to train its own local model, using its own local data. Then these
local models are sent to a central aggregation unit to be combined together.
This process continues until a criterion is met.

Figure 1: The federated reinforcement learning process

###  IV-B Federated DRL Approach

The steps to train the FedRL agents are presented in the following:

####  IV-B1 Device selection strategy

At the beginning of each iteration of FDL, a set of IoT agents are selected to
participate in the FDL. Thus, of all $N$ devices in the network, only a small
subset, denoted by $\mathcal{I}=\\{1,...,I\\}$ is selected to contribute in
FDL. In this paper, the device selection is done based on the following
criterion:

$\arg\max_{i\in\mathcal{N}}~{}\text{Var}(\frac{d_{i}P_{\mathrm{max},i}}{F_{\mathrm{max},i}}),$
(13)

where $d_{i}$ represents the distance of device from BS and the function Var
stands for variance. This criterion helps in identifying devices whose
experiences are more heterogeneous and thus can contribute more in the the
learning process.

####  IV-B2 Training local models

As explained previously , all IoT devices use DDQN to train their local
models. After this local training is finished (no more unprocessed task
remains in the queue), the weights of online network,
$\theta_{i}^{\mathrm{online}}$, is extracted in each agent and is then sent to
the central aggregating unit.

####  IV-B3 Model Aggregation

When central unit receives the models of participating IoT devices, it would
aggregate the models which results in a single global model that would be then
transmitted to all agents. For the purpose of aggregation, we utilize FedAvg
[13], and perform model aggregation as:

$\theta^{\mathrm{global}}=\frac{\sum_{i\in\mathcal{I}}\theta_{i}^{\mathrm{online}}}{|\mathcal{I}|}.$
(14)

This global model, which has integrated the experiences of all devices, is
then transmitted back to IoT devices and the three steps above would be
repeated. The details of our proposed framework is provided in Algorithm 1 as
well as the flowchart given in Fig. 1.

##  V Simulation results

Here, we present our simulation results and extract useful insights related to
the performance of our proposed federated DDQN framework in comparison to
federated DDQN and distributed DDQN algorithms. In addition, we investigate
the impact of batch size, network layers, target network update frequency on
the convergence of the FDL. In what follows, we first focus on the impact of
parameters of DRL on the learning speed of our proposed algorithm and then the
comparison of the proposed algorithm with benchmarks would be presented.

To simulate our system, we consider a network of 100 IoT devices among which
only 20 devices are selected in each round to contribute in the FDL process.
Without loss of generality and for the sake of fair comparison, we assume the
maximum computation capacity and energy consumption limit of the IoT devices
are 1 Gbps and 23 dBm, respectively.

Fig. 2 demonstrates the effect of network architecture on the convergence of
our proposed FedRL algorithm. Here we have some shallow networks with up to
five layers and deeper networks that are obtained by stacking multiple layers
with [16,32,32] neurons on each other. We note that by increasing the number
of layers, faster model training can be achieved. The reason behind this
observation is that, by exploiting deeper neural networks, we can better find
the patterns in data (here devices’ experiences), which subsequently improves
the quality of our local models. Thus, the global model is trained much faster
as its underlying components, local models, are more accurate.

However, since our algorithm will be executed on IoT devices that often lack
necessary resources to train a deep network, it may be infeasible to implement
deeper neural networks. Therefore, in the next figure, we select a rather
simple network architecture with [30,64,16,32,32] neurons in each layer and
instead look for other parameters that may facilitate the learning process.

Figure 2: The impact of neural network architecture on the convergence of the
proposed FedRL algorithm.

The other parameter we focus on in Fig. 3, is the batch size. We can observe
from this figure that as the batch size increases the convergence of our
proposed FedRL algorithm becomes faster. When batch size is equal to 10 and is
extremely small it takes up to 200 iterations to finally converge to a
relative global model, whereas in case batch size is 30, convergence is
achieved almost around iteration number 40. By increasing the size of batches,
we are basically training our model using more data instances. which results
in enhancing the quality of local models and faster training process.

Figure 3: The impact of batch size on the convergence of the proposed FedRL
algorithm.

Similar to network architecture and the stated concern regarding the limited
computation capacity, memory is another bottleneck in learning process of IoT
devices. Larger batch size means higher memory consumption. If the device is
limited both in CPU and memory capacity, neither very deep neural networks nor
increased batch size can be a proper solution to facilitate deployment of
FedRL on IoT devices.

To this end, in Fig. 4, we illustrate the effect of one of the parameters of
DDQN, namely frequency of updating target network with online network. We can
observe here that while the effect of this parameter on performance of DDQN
algorithm is well investigated, this parameter is also considerably effective
in the performance of federated DDQNs. Since many of the components in our
state space, such as path-gain, QoS of tasks, and the length of the task
queue, are constantly changing, efficient choice of the frequency of updating
target network can stabilize the environment enough for the agent to track it
better and obtain a better solution. This effect on local models is quite
notable on the FDL as well.

Figure 4: The impact of target network update frequency on the convergence of
the proposed FedRL algorithm.

In Fig. 5, we compare the performance of our proposed federated DDQN approach
with those of federated DQN and simple distributed DDQN without any
aggregation. As can be seen, the performance of federated DDQN is superior to
federated DQN in terms of learning speed. As previously explained, the main
advantage of DDQN over DQN is the capability to keep target network
stationary, helping with the tractability of states’ values and subsequently a
faster convergence to the correct estimation of them. The impact of this
approach is even more notable when DRL is combined with FDL, since if the
local models are not correctly trained, their errors would be propagated to
other devices’ local model through aggregation. Therefore, aggregation can in
fact negatively effect the result.

The comparison between distributed DDQN and federated DDQN underlines that the
benefits of federated DRL are not limited to its scalability and privacy
preservation. The aggregation incorporated in FDL provides IoT devices a great
context to cooperatively train their models and merge their intelligence
together while preserving privacy of their information. Exploiting federated
learning, at every training round is almost as if devices’ models are trained
with $I$ times more data than their local information. The significance of
this share of knowledge and not data is quite notable in Fig. 5, where as the
result of this aggregation step, federated DDQN is working much better and
faster than simple distributed DDQN. Note that, averaging is used as the
smoothing function.

Figure 5: Performance of federated DDQN compared to federated DQN and
distributed non-federated DDQN.

##  VI Conclusion

In this paper, we investigated the problem of joint delay and energy
minimization in an IoT network with a three-tier offloading scheme. To solve
this problem we combined FDL, DDQN, and optimization theory. Combination of
these tools helped us to achieve a scalable, privacy-preserving, and
computationally efficient framework for joint power and computation resource
allocation and offloading decision optimization. In simulation results, we
compared our work with those of 1) federated DQN to demonstrate the
superiority of DDQN, especially in dynamic environments and 2) with
distributed DDQN to signify the impact of aggregation step incorporated in FDL
on the performance of the framework.

## References

  * [1] S. Zarandi and H. Tabassum, “Delay minimization in sliced multi-cell mobile edge computing (mec) systems,” _IEEE Communications Letters_ , pp. 1–1, 2021.
  * [2] A. Khalili, S. Zarandi, and M. Rasti, “Joint resource allocation and offloading decision in mobile edge computing,” _IEEE Communications Letters_ , vol. 23, no. 4, pp. 684–687, 2019.
  * [3] W. Y. B. Lim, N. C. Luong, D. T. Hoang, Y. Jiao, Y. C. Liang, Q. Yang, D. Niyato, and C. Miao, “Federated learning in mobile edge networks: A comprehensive survey,” _IEEE Communications Surveys Tutorials_ , vol. 22, no. 3, pp. 2031–2063, 2020.
  * [4] S. Luo, X. Chen, Q. Wu, Z. Zhou, and S. Yu, “Hfel: Joint edge association and resource allocation for cost-efficient hierarchical federated edge learning,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 10, pp. 6535–6548, Oct 2020.
  * [5] J. X. X. Mo, “Energy-efficient federated edge learning with joint communication and computation design,” _arXiv_ , Feb 2020.
  * [6] J. Yao and N. Ansari, “Enhancing federated learning in fog-aided iot by cpu frequency and wireless power control,” _IEEE Internet of Things Journal_ , pp. 1–1, 2020.
  * [7] Y. Chen, X. Sun, and Y. Jin, “Communication-efficient federated deep learning with layerwise asynchronous model update and temporally weighted aggregation,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 31, no. 10, pp. 4229–4238, 2020.
  * [8] J. Mills, J. Hu, and G. Min, “Communication-efficient federated learning for wireless edge intelligence in iot,” _IEEE Internet of Things Journal_ , vol. 7, no. 7, pp. 5986–5994, 2020.
  * [9] Y. Zhan, P. Li, Z. Qu, D. Zeng, and S. Guo, “A learning-based incentive mechanism for federated learning,” _IEEE Internet of Things Journal_ , vol. 7, no. 7, pp. 6360–6368, 2020.
  * [10] Y. Zhan, P. Li, and S. Guo, “Experience-driven computational resource allocation of federated learning by deep reinforcement learning,” in _2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS)_ , 2020, pp. 234–243.
  * [11] X. Wang, C. Wang, X. Li, V. C. M. Leung, and T. Taleb, “Federated deep reinforcement learning for internet of things with decentralized cooperative edge caching,” _IEEE Internet of Things Journal_ , vol. 7, no. 10, pp. 9441–9455, 2020.
  * [12] J. Ren, H. Wang, T. Hou, S. Zheng, and C. Tang, “Federated learning-based computation offloading optimization in edge computing-supported internet of things,” _IEEE Access_ , vol. 7, pp. 69 194–69 201, 2019.
  * [13] D. R. H. B. McMahan, E. Moore and B. A. Y. Arcas, “Federated learning of deep networks using model averaging,” _arXiv preprint arXiv:1602.05629, 2016._ , 2016.