# Large-Scale Learning with Fourier Features and Tensor Decompositions

Frederiek Wesel
Delft Center for Systems and Control
Delft University of Technology
<EMAIL_ADDRESS>
&Kim Batselier
Delft Center for Systems and Control
Delft University of Technology
<EMAIL_ADDRESS>

###### Abstract

Random Fourier features provide a way to tackle large-scale machine learning
problems with kernel methods. Their slow Monte Carlo convergence rate has
motivated the research of deterministic Fourier features whose approximation
error can decrease exponentially in the number of basis functions. However,
due to their tensor product extension to multiple dimensions, these methods
suffer heavily from the curse of dimensionality, limiting their applicability
to one, two or three-dimensional scenarios. In our approach we overcome said
curse of dimensionality by exploiting the tensor product structure of
deterministic Fourier features, which enables us to represent the model
parameters as a low-rank tensor decomposition. We derive a monotonically
converging block coordinate descent algorithm with linear complexity in both
the sample size and the dimensionality of the inputs for a regularized squared
loss function, allowing to learn a parsimonious model in decomposed form using
deterministic Fourier features. We demonstrate by means of numerical
experiments how our low-rank tensor approach obtains the same performance of
the corresponding nonparametric model, consistently outperforming random
Fourier features.

##  1 Introduction

Kernel methods such as Support Vector Machines and Gaussian Processes are
commonly used to tackle problems such as classification, regression and
dimensionality reduction. Since they can be universal function approximators
[13], kernel methods have received renewed attention in the last few years and
have shown equivalent or superior performance to Neural Networks [22, 26, 11].
The main idea behind kernel methods is to lift the data into a higher-
dimensional (or even infinite-dimensional) Reproducing Kernel Hilbert Space
$\mathcal{H}$ by means of a feature map
$\bm{\phi}\left(\cdot\right):\mathcal{X}\rightarrow\mathcal{H}$. Considering
then the pairwise similarities between the mapped data allows to tackle
problems which are highly nonlinear in the original sample space
$\mathcal{X}$, but linear in $\mathcal{H}$. This can be done equivalently by
considering a kernel function
$k\left(\cdot,\cdot\right):\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$
such that
$\left\langle\bm{\phi}\left(\bm{x}\right),\bm{\phi}\left(\bm{x^{\prime}}\right)\right\rangle=k\left(\bm{x},\bm{x^{\prime}}\right)$
and performing thus said mapping implicitly. Although effective at learning
nonlinear patterns in the data, kernel methods are known to scale poorly as
the number of data points $N$ increases. For example, when considering Kernel
Ridge Regression (KRR), Gaussian Process Regression (GPR) [33] or Least-
Squares Support Vector Machine (LS-SVM) [41, 42] training usually consist in
inverting the Gram matrix $k_{ij}=k\left(\bm{x}_{i},\bm{x}_{j}\right)$, which
encodes the pairwise relation between all data. As a consequence, the
associated storage complexity is $\mathcal{O}\left(N^{2}\right)$ and the
computational complexity is $\mathcal{O}\left(N^{3}\right)$, rendering these
methods unfeasible for large data. In order to lower the computational cost,
data dependent methods approximate the kernel function by means of $M$ data
dependent basis functions. Due to their reduced-rank formulation, the
computational complexity is reduced to $\mathcal{O}\left(NM^{2}\right)$ for
$N\gg M$. However, e.g. for the Nyström method [46], the convergence rate is
only of $\mathcal{O}\left(\nicefrac{{1}}{{\sqrt{M}}}\right)$ [7] limiting its
effectiveness. As the name suggests, data independent methods approximate the
kernel function by $M$ data independent basis functions. A good example is the
celebrated Random Fourier Features approach by Rahimi and Recht [31], where
the authors propose for stationary kernels a low-dimensional random mapping
$\bm{z}\left(\cdot\right):\mathbb{R}^{D}\rightarrow\mathbb{R}^{M}$ such that

$k\left(\bm{x},\bm{x^{\prime}}\right)=\left\langle\bm{\phi}\left(\bm{x}\right),\bm{\phi}\left(\bm{x^{\prime}}\right)\right\rangle\approx\left\langle\bm{z}\left(\bm{x}\right),\bm{z}\left(\bm{x^{\prime}}\right)\right\rangle.$
(1)

As in the case of data dependent methods, the reduced-rank formulation allows
for a computational complexity of $\mathcal{O}\left(NM^{2}\right)$ for $N\gg
M$. Probabilistic error bounds on the approximation are provided which result
in a convergence rate of $\mathcal{O}\left(\nicefrac{{1}}{{\sqrt{M}}}\right)$,
which is again the Monte Carlo rate.

Improvements in this sense were achieved by considering deterministic features
resulting from dense quadrature methods [5, 25, 36] and kernel eigenfunctions
[16, 39]. These methods are able to achieve exponentially decreasing upper
bounds on uniform convergence guarantees when certain conditions are met.
However, for a $D$-dimensional input space these methods take the tensor
product of $D$ vectors, resulting in an exponential increase in the number of
basis functions and thus of model weights, effectively limiting the
applicability of deterministic features to low-dimensional data. In this work
we consider deterministic features. In order to take advantage of the tensor
product structure which arises when mapping the inputs to $\mathcal{H}$, we
represent the weights as a low-rank tensor. This allows us to learn the inter-
modal relations in the tensor product of (low-dimensional) Hilbert spaces,
avoiding the exponential computational and storage complexities in $D$. In
this way we are able to obtain a linear computational complexity in both the
number of samples and in the input dimension during training, without having
to resort to the use of sparse grids or additive modeling of kernels.

The main contribution of this work is in lifting the curse of dimensionality
affecting deterministic Fourier features by modeling the weights as a low-rank
tensor. This enables the efficient solution of large-scale and high-
dimensional kernel learning problems. We derive an iterative algorithm under
the exemplifying case of regularized squared loss and test it on regression
and classification problems.

##  2 Related work

Fourier Features (FF) are a collection of data independent methods that
leverage Bochner’s theorem [35] from harmonic analysis to approximate
stationary kernels by numerical integration of their spectral density
$p\left(\cdot\right)$:

$\displaystyle
k\left(\bm{x,x^{\prime}}\right)\overset{\mathrm{Stationarity}}{\coloneqq}k\left(\bm{x}-\bm{x^{\prime}}\right)\overset{\mathrm{Bochner}}{=}\int
p\left(\bm{\omega}\right)\exp\left(\left\langle
i\bm{\omega},\left(\bm{x-x^{\prime}}\right)\right\rangle\right)d\bm{\omega}\overset{\mathrm{FF}}{\approx}\left\langle
z\left(\bm{x}\right),z\left(\bm{x^{\prime}}\right)\right\rangle.$ (2)

Rahimi and Recht [31] proposed to approximate the integral by Monte Carlo
integration i.e. by drawing $M$ random frequencies $\bm{\omega}\sim
p\left(\cdot\right)$. In their work they show how the method converges
uniformly at the Monte Carlo rate. See [23] for an overview of random FF. In
order to achieve faster convergence and a lower sample complexity, a multitude
of approaches that rely on deterministic numerical quadrature of the Fourier
integral were developed. These methods generally consider product kernels
whose spectral density factors in the frequency domain, which enables in turn
to factor the Fourier integral. The resulting deterministic Fourier features
are then the tensor product of $D$ one-dimensional deterministic features.

For example, in [5] the authors give an analysis of the sample complexity of
features resulting from dense Gaussian Quadrature (GQ). In [25] the authors
present a similar quadrature approach for kernels whose spectral density
factors over the dimensions of the inputs. They provide an explicit
construction for their Quadrature Fourier Features relying on a dense
Cartesian grid, and note that their method, as well as GQ, can attain
exponentially decreasing uniform convergence bounds in the total number of
basis functions per dimensions $\hat{M}$ [25, Theorem 1]. To avoid the curse
of dimensionality, they make use of additive modeling of kernels. Variational
Fourier Features (VFF) [16] are derived probabilistically in a one-dimensional
setting for Matérn kernels by projecting a Gaussian Process onto a set of
Fourier basis. An extension to multiple dimensions when considering product
kernels is then achieved by taking the tensor product of the one-dimensional
features, incurring however again in exponentially rising computational costs
in $D$. In what they call Hilbert-GP [39] the authors diagonalize stationary
kernels in terms of the eigenvalues and eigenfunctions of the Laplace operator
with Dirichlet boundary conditions. Due to the multiplicative structure of the
eigenfunctions of the Laplace operator, the complexity of the basis function
increases exponentially in $D$ when one considers product kernels. Like with
other deterministic approximations, bounds on the uniform convergence error
which decrease exponentially in $\hat{M}$ can be achieved [39, Theorem 8].

Tensor decompositions have been used extensively in machine learning in order
to benefit from structure in data [38], in particular to obtain an efficient
representation of the model parameters. In the context of Gaussian Process
Regression, the tensor product structure which arises when the inputs are
located on a Cartesian grid have been exploited to speedup inference [12]. In
their variational inference framework, [19] proposed to parameterize the
posterior mean of a Gaussian Process by means of a tensor decomposition in
order to exploit the tensor product structure which arises when interpolating
the kernel function. In the context of neural networks, [27] proposed to
represent the weights in deep neural networks as a tensor decomposition to
speed up training. A similar approach was carried out in case of recurrent
neural networks [48] [43].

Tensor decompositions have also been used to learn from simple features. In
[45] multivariate functions are learned from a Fourier basis, which
corresponds to assuming a uniform spectral density in (2). Motivated by spin
vectors in quantum mechanics, [40] consider trigonometric feature maps which
they argue induce uninformative kernels. On a similar note, [28] and [1]
leverage the tensor product structure of polynomial mappings to induce a
polynomial kernel and consider model parameters in decomposed form. While
these existing approaches are able to overcome the curse of dimensionality
affecting tensor product feature maps by modeling the parameter tensor as a
tensor network, they consider simple and empirical feature maps which induce
uninformative kernels. In the following section we show how deterministic
Fourier features can be used for supervised learning in both large and high-
dimensional scenarios by assuming that the model weights are a low-rank
tensor, thereby linking tensor decompositions with stationary product kernels.

##  3 Learning with deterministic Fourier features

###  3.1 Notation

All tensor operations and mappings can be formulated over both the complex and
the real field. For the remainder of this article real-valued mappings will be
considered. Vectors and matrices are denoted by boldface lowercase and
uppercase letters, respectively. Higher-order tensors are denoted by boldface
calligraphic letters. Tensor entries are always written as lightface lowercase
letters with a subscript index notation. For example, element
$i_{1},i_{2},i_{3}$ of the third-order tensor
$\bm{\mathcal{A}}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is written as
$a_{i_{1}i_{2}i_{3}}$. The vectorization
$\text{vec}\left(\bm{\mathcal{A}}\right)$ of a tensor
$\bm{\mathcal{A}}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{D}}$ is
the vector such that

$\displaystyle\text{vec}\left(\bm{\mathcal{A}}\right)_{i}=a_{i_{1}i_{2}\cdots
i_{D}},$

where the relationship between the linear index $i$ and $i_{1}i_{2}\ldots
i_{D}$ is given by

$\displaystyle i=i_{1}+\sum_{d=2}^{D}(i_{d}-1)\prod_{k=1}^{d-1}I_{k}.$

The Frobenius inner product between tensors
$\bm{\mathcal{A}}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{D}}$
and $\bm{\mathcal{B}}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times
I_{D}}$ is defined as

$\displaystyle\left\langle\bm{\mathcal{A}},\bm{\mathcal{B}}\right\rangle_{\text{F}}\coloneqq\sum_{i_{1}=1}^{I_{1}}\sum_{i_{2}=1}^{I_{2}}\cdots\sum_{i_{D}=1}^{I_{D}}a_{i_{1}i_{2}\cdots
i_{D}}b_{i_{1}i_{2}\cdots
i_{D}}={\text{vec}\left(\bm{\mathcal{A}}\right)}^{\text{T}}\text{vec}\left(\bm{\mathcal{B}}\right).$

Depending on the context, the symbol $\otimes$ either denotes the tensor outer
product or the matrix Kronecker product. The Khatri-Rao product
$\bm{A}\circledast\bm{B}$ of the matrices $\bm{A}\in\mathbb{R}^{N\times R}$
and $\bm{B}\in\mathbb{R}^{M\times R}$ is the $NM\times R$ matrix that is
obtained by taking the column-wise Kronecker product. The Hadamard (element-
wise) matrix product of $\bm{A}\in\mathbb{R}^{N\times R}$ and
$\bm{B}\in\mathbb{R}^{N\times R}$ is denoted by $\bm{A}\odot\bm{B}$.

###  3.2 The model

In this article we assume models of the form

$\displaystyle
f\left(\bm{x}\right)=\left\langle\bm{w},\bm{\phi}\left(\bm{x}\right)\right\rangle,\quad
y=f\left(\bm{x}\right)+\epsilon.$ (3)

Here $\bm{x}\in\mathbb{R}^{D}$ is the input vector,
$\bm{\phi}\left(\cdot\right):\mathbb{R}^{D}\rightarrow\mathcal{H}$ is a
feature map, $\bm{w}\in\mathcal{H}$ are the model weights, $y\in\mathbb{R}$ is
the corresponding available observation corrupted by a zero mean i.i.d. noise
$\epsilon$. Learning such a kernelized model consists in finding as set of
weights $\bm{w}$ such that:

$\displaystyle\sum_{n=1}^{N}L\left(y_{n},f\left(\bm{x}_{n}\right)\right)+\lambda\left\langle\bm{w},\bm{w}\right\rangle_{\text{p}},$
(4)

is minimized. Here,
$L\left(\cdot,\cdot\right):\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}_{+}$
is a symmetric and positive loss function and
$\lambda\left\langle\bm{w},\bm{w}\right\rangle_{\text{p}}$ is a $p$-th norm
regularization term. A variety of primal machine learning formulations arise
when considering different combinations of loss functions and regularization
terms. For example, considering $p=2$ and squared loss leads to Kernel Ridge
Regression [42], while considering hinge loss leads to SVM [4], and so on.
Furthermore, applying the kernel trick when considering
$\bm{\bm{\phi}}\left(\cdot\right)$ recovers the nonparametric formulation
which, depending on the choice of kernel, enables to perform inference using
infinitely many basis functions with a computational complexity of
$\mathcal{O}(N^{3})$. Considering instead a low-dimensional finite mapping
$\bm{z}\left(\cdot\right):\mathcal{X}\rightarrow\mathbb{R}^{M}$ such as Random
Fourier Features or the Nyström method leads a primal approach and to
computational savings with a computational complexity of
$\mathcal{O}\left(NM^{2}\right)$ for $N\gg M$. However, as already discussed,
these low-dimensional random mappings converge at the slow Monte Carlo rate,
motivating our approach.

In order to approximate the kernel function with faster and possibly
exponential convergence, we consider deterministic, finite-dimensional
features which are the tensor product of $D$ vectors. For a given
$D$-dimensional input point $\bm{x}$ we define the deterministic feature
mapping
$\bm{\mathcal{Z}}\left(\cdot\right):\mathbb{R}^{D}\rightarrow\mathbb{R}^{\hat{M}_{1}}\otimes\mathbb{R}^{\hat{M}_{2}}\otimes\cdots\otimes\mathbb{R}^{\hat{M}_{D}}$
as

$\bm{\mathcal{Z}}\left(\bm{x}\right)=\bm{z}^{(1)}\left(x_{1}\right)\otimes\bm{z}^{(2)}\left(x_{2}\right)\cdots\otimes\bm{z}^{(D)}\left(x_{D}\right),$
(5)

where $\bm{z}^{\left(d\right)}$ is a deterministic mapping applied to the
$d$-th dimension. The dimension of the feature space is denoted
$M=\prod_{d=1}^{D}\hat{M}_{d}$. Without loss of generality it is from here on
assumed that $\hat{M}_{1}=\cdots=\hat{M}_{D}=\hat{M}$ such that
$M=\hat{M}^{D}$. It is easy to verify by applying the kernel trick that the
features of (5) yield product kernels. The mapping $\bm{\mathcal{Z}}(\cdot)$
encompasses for instance the mappings derived by quadrature [5, 25, 36] and by
projection [16, 39]. Note that these mappings cover many popular kernels such
as the Gaussian and Matérn kernels. To give a concrete example, in the
framework of A. Solin and S. Särkkä [39, Equation 60], for input data centered
in a hyperbox
$\mathcal{X}=\left[-U_{1},U_{1}\right]\otimes\cdots\otimes\left[-U_{D},U_{D}\right]$,
the Gaussian kernel is approximated by means of $D$ tensor products of
$\hat{M}$ weighted sinusoidal basis functions with frequencies lying on a
harmonic scale such that:

$\left({\bm{z}^{\left(d\right)}\left(x_{d}\right)}\right)_{i_{d}}=\frac{1}{\sqrt{U_{d}}}\
p\left(\frac{\pi i_{d}}{2U_{d}}\right)\ \sin\left(\frac{\pi
i_{d}\left(x_{d}+U_{d}\right)}{2U_{d}}\right),\ i_{d}=1,\ldots,\hat{M}.$ (6)

Here $p\left(\cdot\right)$ is the spectral density of the Gaussian kernel with
one-dimensional inputs, which is known in closed-form [33, page 83]. This
deterministic mapping then approximates the Gaussian kernel function [39,
Equation 59] such that

$\displaystyle
k\left(\bm{x},\bm{x^{\prime}}\right)\approx\left\langle\text{vec}\left(\bm{\mathcal{Z}}\left(\bm{x}\right)\right),\text{vec}\left(\bm{\mathcal{Z}}\left(\bm{x^{\prime}}\right)\right)\right\rangle=\langle\bm{\mathcal{Z}}\left(\bm{x}\right),\bm{\mathcal{Z}}\left(\bm{x^{\prime}}\right)\rangle_{\text{F}},$

and converges uniformly with exponentially decreasing bounds [39, Theorem 8].
We therefore use $\bm{\mathcal{Z}}(\cdot)$ instead of $\bm{\phi}(\cdot)$ in
(3) to obtain

$f\left(\bm{x}\right)=\langle\bm{w},\text{vec}(\bm{\mathcal{Z}}(\bm{x}))\rangle\\\
=\langle\bm{\mathcal{W}},\bm{\mathcal{Z}}(\bm{x})\rangle_{\text{F}},$ (7)

where the weight vector $\bm{w}$ has been reshaped into a $D$-dimensional
tensor
$\bm{\mathcal{W}}\in\mathbb{R}^{\hat{M}_{1}\times\hat{M}_{2}\times\cdots\times\hat{M}_{D}}$.
Learning the exponential number of model parameters $\bm{\mathcal{W}}$ in (7)
under a hinge loss leads to Support Vector Machines, while considering a
squared loss leads to the primal formulation of Kernel Ridge Regression, which
we will consider as exemplifying case from here on:

$\min_{\bm{\mathcal{W}}}\sum_{n=1}^{N}\left(y_{n}-\langle\bm{\mathcal{W}},\bm{\mathcal{Z}}\left(\bm{x}_{n}\right)\rangle_{\text{F}}\right)^{2}+\lambda\langle\bm{\mathcal{W}},\bm{\mathcal{W}}\rangle_{\text{F}}.$
(8)

Since the number of elements $M$ in $\bm{\mathcal{W}}$ and $\bm{\mathcal{Z}}$
grows exponentially in $D$, this primal approach is advantageous compared to
the nonparametric dual approach only if $\hat{M}^{D}\ll N$, limiting it to
low-dimensional inputs. In order to lift this curse of dimensionality, we
propose to represent and to learn $\bm{\mathcal{W}}$ directly as a low-rank
tensor decomposition. A low-rank tensor decomposition allows us to exploit
redundancies in $\bm{\mathcal{W}}$ in order to obtain a parsimonious model
with a storage complexity that scales linearly in both $\hat{M}$ and $D$. The
low-rank structure will, as explicitly shown in the experiments, act as a form
of regularization by limiting the total number of degrees of freedom of the
model.

###  3.3 Low-rank tensor decompositions

Tensor decompositions (also called tensor networks) can be seen as
generalizations of the Singular Value Decomposition (SVD) of a matrix to
tensors [20]. Three common tensor decompositions are the Tucker decomposition
[6, 44], the Tensor Train decomposition [29] and the Canonical Polyadic
Decomposition (CPD) [17, 21], where each of them encompasses different
properties of the matrix SVD. In this subsection we briefly discuss these
three decompositions and list both their advantages and disadvantages for
modeling $\bm{\mathcal{W}}$ in (8). For a detailed exposition on tensor
decompositions we refer the reader to [2] and the references therein. An
important property of tensor decompositions is that they can always be written
linearly in their components, which implies that applying a block coordinate
descent algorithm to solve (8) results in a series of linear least-squares
problems.

A rank-$R$ CPD of $\bm{\mathcal{W}}$ consists of $D$ factor matrices
$\bm{W}^{(d)}\in\mathbb{R}^{\hat{M}\times R}$ and a vector
$\bm{s}\in\mathbb{R}^{R}$ such that

$\displaystyle\bm{w}$
$\displaystyle=\left(\bm{W}^{(1)}\circledast\bm{W}^{(2)}\circledast\cdots\circledast\bm{W}^{(D)}\right)\bm{s}=\sum_{r=1}^{R}\bm{w}^{\left(1\right)}_{r}\otimes\bm{w}^{\left(2\right)}_{r}\otimes\cdots\otimes\bm{w}^{\left(D\right)}_{r}.$

Note that the entries of the $\bm{s}$ vector were absorbed in one of the
factor matrices when writing the CPD as a sum of $R$ terms. The CPD has been
shown, in contrast to matrix factorizations, to be unique under mild
conditions [37]. The storage complexity comes mainly from the $D$ factor
matrices and is therefore $\mathcal{O}(R\hat{M}D)$. The Tucker decomposition
generalizes the CPD in two ways. First, the Khatri-Rao product is replaced
with a Kronecker product. The $\bm{s}$ vector has to grow in length
accordingly from $R$ to $R^{D}$. Second, the $R$-dimension of each of the
factor matrices is allowed to vary, resulting in a multi-linear rank
$(R_{1},R_{2},\ldots,R_{D})$:

$\displaystyle\bm{w}$
$\displaystyle=\left(\bm{W}^{(1)}\otimes\bm{W}^{(2)}\otimes\cdots\otimes\bm{W}^{(D)}\right)\bm{s}.$

The Tucker decomposition is inherently non-unique and its storage complexity
$\mathcal{O}(R^{D})$ is dominated by the $\bm{s}$ vector, which limits its use
to low-dimensional input data. For this reason we do not consider the Tucker
decomposition any further in this article.

The Tensor Train (TT) decomposition consists of $D$ third-order tensors
$\bm{\mathcal{W}}^{(d)}\in\mathbb{R}^{R_{d}\times\hat{M}\times R_{d+1}}$ such
that

$\displaystyle w_{i_{1}i_{2}\cdots
i_{D}}=\sum_{r_{1}=1}^{R_{1}}\cdots\sum_{r_{D+1}=1}^{R_{D+1}}w^{(1)}_{r_{1}i_{1}r_{2}}\cdots
w^{(D)}_{r_{D}i_{D}r_{D+1}}.$ (9)

The auxiliary dimensions $R_{1},R_{2},\ldots,R_{D+1}$ are called the TT-ranks.
In order to ensure that the right-hand side of (9) is a scalar, the boundary
condition $R_{1}=R_{D+1}=1$ is enforced. The TT decomposition is, just like
the Tucker decomposition, non-unique and its storage complexity
$R^{2}\hat{M}D$ is due to the $D$ tensor components $\bm{\mathcal{W}}^{(d)}$.
Considering their storage complexity, both the CPD and TT decomposition are
viable candidates to replace $\bm{\mathcal{W}}$ in (8). The CP-rank $R$ and
TT-ranks $R_{2},\ldots,R_{D}$ are additional hyperparameters, which favors the
CPD in practice. For the TT, one could choose $R_{2}=R_{3}=\cdots=R_{D}$ to
reduce the number of additional hyperparameters but this constraint turns out
in practice to lead to suboptimal results. For these reasons we limit the
discussion of the learning algorithm to the CPD case.

###  3.4 Tensor learning with deterministic Fourier features

We now wish to minimize the standard regularized squared loss function as in
(8) with the additional constraint that the weight tensor has a rank-$R$ CPD
structure:

$\displaystyle\min_{\bm{\mathcal{W}}}\quad$
$\displaystyle\sum_{n=1}^{N}\left(y_{n}-\left\langle\bm{\mathcal{W}},\bm{\mathcal{Z}}\left(\bm{x}_{n}\right)\right\rangle_{\text{F}}\right)^{2}+\lambda\left\langle\bm{\mathcal{W}},\bm{\mathcal{W}}\right\rangle_{\text{F}},$
(10) subject to $\displaystyle\text{CP-rank}\left(\bm{\mathcal{W}}\right)=R.$
(11)

Note that if $R$ equals the true CP-rank of the underlying weight tensor then
the exact solution of (8) would be obtained. In practice a low-rank solution
for $\bm{\mathcal{W}}$ achieves a sufficiently complex decision boundary that
is practically indistinguishable from the full-rank solution, as is
demonstrated in the experiments. Imposing the rank-$R$ reduces the number of
unknowns from $\hat{M}^{D}$ to $R\hat{M}D$ and allows the application of a
block coordinate descent algorithm (also known as alternating linear scheme),
which is shown to converge monotonically [3, 18, 34]. Each of the factor
matrices $\bm{W}^{\left(d\right)}$ is optimized in an iterative fashion while
keeping the others fixed. Such a factor matrix update is obtained by solving a
linear least-squares problem with $\hat{M}R$ unknowns. In what follows we
derive the linear problem for the $d$-th factor matrix $\bm{W}^{(d)}$. The
data-fitting term
$\langle\bm{\mathcal{W}},\bm{\mathcal{Z}}(\bm{x})\rangle_{\text{F}}$ can be
rewritten linearly in terms of the unknown factor matrix as

$\displaystyle\langle\bm{\mathcal{W}},\bm{\mathcal{Z}}(\bm{x})\rangle_{\text{F}}=\left\langle\text{vec}\left(\bm{W}^{\left(d\right)}\right),\bm{g}^{\left(d\right)}\left(\bm{x}\right)\right\rangle,$
(12)

where
$\bm{g}^{\left(d\right)}\left(\bm{x}\right)\coloneqq\bm{z}^{\left(d\right)}\otimes\left({\bm{z}^{\left(1\right)}}^{\text{T}}{\bm{W}^{\left(1\right)}}^{\text{T}}\odot\cdots\odot{\bm{z}^{\left(D\right)}}^{\text{T}}{\bm{W}^{\left(D\right)}}^{\text{T}}\right)$.
Similarly, for the regularization term:

$\displaystyle\left\langle\bm{\mathcal{W}},\bm{\mathcal{W}}\right\rangle_{\text{F}}=\left\langle\text{vec}\left(\bm{W}^{{\left(d\right)}^{\text{T}}}\bm{W}^{\left(d\right)}\right),\text{vec}\left(\bm{H}^{\left(d\right)}\right)\right\rangle.$
(13)

Here
$\bm{H}^{\left(d\right)}\coloneqq\left(\bm{W}^{{\left(1\right)}^{\text{T}}}\bm{W}^{\left(1\right)}\odot\cdots\odot\bm{W}^{{\left(D\right)}^{\text{T}}}\bm{W}^{\left(D\right)}\right)$.
Derivations of (12) and (13) can be found in Appendix. Substitution of the
expressions (12) and (13) into (10) and (11) leads to a linear least-squares
problem for $\text{vec}\left(\bm{W}^{\left(d\right)}\right)$:

$\displaystyle\min_{\text{vec}\left(\bm{W}^{\left(d\right)}\right)}\sum_{n=1}^{N}\left(y_{n}-\left\langle\text{vec}\left(\bm{W}^{\left(d\right)}\right),\bm{g}^{\left(d\right)}\left(\bm{x}_{n}\right)\right\rangle\right)^{2}+\lambda\left\langle\text{vec}\left(\bm{W}^{{\left(d\right)}^{\text{T}}}\bm{W}^{\left(d\right)}\right),\text{vec}\left(\bm{H}^{\left(d\right)}\right)\right\rangle.$
(14)

Its unique solution can be computed exactly in an efficient manner by solving
the normal equations, requiring $N\hat{M}^{2}R^{2}+\hat{M}^{3}R^{3}$
operations. Since we are solving a non-convex optimization problem that
consists of a series of convex and exactly solvable sub-problems, our
algorithm is indeed monotonically decreasing, and, although it is not
guaranteed to converge to the global optimum, it is the standard choice when
dealing with tensors [3, 20]. The total computational complexity of the
algorithm when $N\gg\hat{M}R$ is then $\mathcal{O}(ND\hat{M}^{2}R^{2})$,
rendering it suitable for problems which are large in both $N$ and $D$,
provided that $R$ and $\hat{M}$ are small. The necessary memory is equal to
$R\hat{M}D+2R^{2}\hat{M}^{2}+2R\hat{M}$ (storing the weight tensor
$\bm{\mathcal{W}}$ in decomposed form, the rank-$R\hat{M}$ Gram matrix and
regularization matrix, the transformed responses and the solution of the
linear system), leading to a storage complexity of
$\mathcal{O}\left(R^{2}\hat{M}^{2}\right)$ for $R\hat{M}\gg D$. Notably the
cost is independent of $N$, which allows processing large data on modest
hardware.

When learning, the selection of $\lambda$ and of the kernel-related
hyperparameters can be carried out by standard methods such as cross-
validation. The choice of $\hat{M}$ and the additional hyperparameter $R$
which we introduce are directly linked with the available computational
budget. One should in fact choose $\hat{M}$ so that the model has access to a
sufficiently complex set of basis functions to learn from. In practice we
notice that for our choices of kernel function hyperparameters, at most
$\hat{M}=40$ basis functions per dimension suffice. $R$ can then be fixed
accordingly in order to match the computational budget at hand. As we will
show in the next section, learning is possible with small values of $\hat{M}$
and $R$.

##  4 Numerical experiments

We implemented our Tensor-Kernel Ridge Regression (T-KRR) algorithm in
Mathworks MATLAB 2021a (Update 1) [24] and tested it on several regression and
classification problems. Our implementation can be freely downloaded from
https://github.com/fwesel/T-KRR and allows reproduction of all experiments in
this section. In our implementation we avoid constructing
$\bm{g}^{\left(d\right)}\left(\cdot\right)$ and $\bm{H}^{\left(d\right)}$ from
scratch at every iteration by updating their components iteratively. With the
exception of the first experiment 4.1, we further speedup our implementation
by considering only the diagonal of the regularization term
$\bm{H}^{\left(d\right)}$. All experiments were run on a Dell Inc. Latitude
7410 laptop with 16 GB of RAM and an Intel Core i7-10610U CPU running at 1.80
GHz. In all our experiments the Gaussian kernel
$k\left(\bm{x},\bm{x^{\prime}}\right)=\exp\left(-\nicefrac{{|\lvert\bm{x}-\bm{x^{\prime}}\rvert\rvert_{2}^{2}}}{{2l^{2}}}\right)$
was approximated by considering the Hilbert-GP [39] mapping of (6). In all
experiments inputs were scaled to lie in a $D$-dimensional unit hypercube.
When dealing with regression, the responses were standardized around the mean,
while when considering binary classification inference was done by looking at
the sign of the model response (LS-SVM [41]). One sweep of our T-KRR algorithm
is defined as updating factor matrices in the order $1\rightarrow D$ and then
back from $D\rightarrow 1$. All initial factor matrices were initialized with
standard normal numbers and normalized by dividing all entries with their
Frobenius norm. For all experiments the number of sweeps of T-KRR algorithm
are set to $10$. In the following three experiments it is shown how our
proposed T-KRR algorithm is stable and recovers the full KRR estimate in case
of low number of frequencies $\hat{M}$ and low rank $R$, consistently
outperforming RFF. Finally, our model exhibits very competitive performance on
large-scale problems when compared with other kernel methods.

###  4.1 Banana classification

The banana dataset is a two-dimensional classification problem [15] consisting
of $N=5300$ datapoints. Since the data is two-dimensional, we consider
$\hat{M}=12$ frequencies per dimension, which enables us to compare T-KRR with
Hilbert-GP. Furthermore since for tensors of order two (i.e. matrices) the CP-
rank is the matrix rank, we can deduce that our approach should recover the
underlying Hilbert-GP method with $R=\hat{M}^{2}=144$. In this example we fix
the kernel hyperparameters of our method and of the Hilbert-GP to
$l=\nicefrac{{1}}{{2}}$ and $\lambda=10^{-5}$ for visualization purposes. In
Figure 1 we plot the decision boundary of T-KRR (full line) and of the
associated Hilbert-GP (dashed line). We can see that already when $R=6$ the
learned decision boundary is indistinguishable to the one of Hilbert-GP,
meaning that in this example it is possible to obtain a low-rank
representation of the model weights. From a computational point of view, we
solve a series of linear systems with $\hat{M}R$ unknowns instead of
$\hat{M}^{D}$. However due to the low-dimensionality of the dataset, the
computational savings per iteration are very modest, i.e. for $R=6$ and
$\hat{M}=12$, we solve per iteration a linear system with $72$ unknowns as
opposed to the one-time solve of a linear system with $144$ unknowns as in
case of Hilbert-GP. As we show in Subsection 4.2, the linear computational
complexity in $D$ of our algorithm results in ever-increasing performance
benefits as the dimensionality of the problem becomes larger, where it becomes
impossible to consider a full weight tensor.

(a) $R=2$

(b) $R=4$

(c) $R=6$

(d) $R=\hat{M^{2}}=144$

Figure 1: Classification boundary of the two-dimensional banana dataset for
increasing CP-ranks. In the last plot the chosen CP-rank matches the true
matrix rank. The dashed line is the Hilbert-GP, while the full line is T-KRR.

###  4.2 Model performance with baseline

We consider five UCI [8] datasets in order to compare the performance of our
model with RFF and the GPR/KRR baseline. For each dataset, we consider $90\%$
of the data for training and the remaining $10\%$ for testing. In particular,
in all regression datasets, we first obtain an estimate of the lengthscale of
the Gaussian kernel $l$ and regularization term $\lambda$ by log-marginal
likelihood optimization over the whole training set using the GPLM toolbox
[32]. We subsequently train GPR/KRR, RFF and our method with the estimated
hyperparameter. In case of classification, we choose the lengthscale $l$ to be
the sample mean of the sample standard deviation of our data (which is the
default choice in the Matlab Machine Learning Toolbox and scikit-learn [30])
and $\lambda=10^{-5}$. We repeat this procedure ten times over different
random splits and report the sample mean and sample standard deviation of the
predictive mean-squared error (MSE) and the misclassification rate for
regression and classification respectively. In order to test the approximation
capabilities of our approach, we set $\hat{M}$ and $R$ such that $\hat{M}R\ll
N$ in order to benefit from computational gains. Furthermore, we note that due
to the curse of dimensionality it is not possible to compare our approach with
RFF while considering the same total number of basis functions
$M=\hat{M}^{D}$. Hence we select $M_{\text{RFF}}=\hat{M}R$ such that the
computational complexity of both methods is similar. Table 1 shows the MSE for
different UCI datasets. We see that our approach comes close to the
performance of full KRR and outperforms RFF on all datasets when considering
the same number of model parameters. This is because our approach considers
implicitly $\hat{M}^{D}$ frequencies to learn a low-rank parsimonious model
which is fully described by $\hat{M}RD$ parameters, as opposed to RFF where
the number of frequencies is the same as the number of model weights.
Similarly to what is reported in [47, Figure 2] we observe lower performance
of RFF on the Adult dataset than reported in [31]. Figure 2 plots the
monotonically decreasing loss function as well as the corresponding
misclassification rate while training with T-KRR on the Spambase dataset.
After four sweeps the misclassification rate has converged. Since T-KRR is
able to obtain similar performance as the KRR baseline on a range of small
datasets, we proceed to tackle a large-scale regression problem.

Figure 2: Normalized loss and misclassification rate while training on the Spambase dataset. Table 1: Predictive MSE (regression) and misclassification rate (classification) with one standard deviation for RFF, T-KRR and KRR on different UCI datasets. Dataset | $N\downarrow$ | $D$ | $\hat{M}$ | $R$ | RFF | T-KRR | KRR
---|---|---|---|---|---|---|---
Yacht | 308 | 6 | 10 | 25 | $0.0021\pm 0.0018$ | $0.0009\pm 0.0006$ | $\bm{0.0007}\pm 0.0004$
Energy | 768 | 9 | 20 | 10 | $0.0275\pm 0.0075$ | $0.0200\pm 0.0066$ | $0.0200\pm 0.0087$
Airfoil | 1503 | 5 | 20 | 10 | $0.2180\pm 0.0336$ | $0.1679\pm 0.0258$ | $\bm{0.1587}\pm 0.0232$
Spambase | 4601 | 57 | 40 | 10 | $0.3620\pm 0.0188$ | $0.0935\pm 0.0095$ | $\bm{0.0909}\pm 0.0115$
Adult | 45222 | 96 | 40 | 10 | $0.3976\pm 0.0056$ | $\bm{0.1596}\pm 0.0046$ | N/A

###  4.3 Large-scale experiment

In order to highlight the ability of our method to deal with large data we
consider the Airline dataset. The Airline dataset [14, 16] is a large-scale
regression problem originally considered in [14] which is often used to
compare state-of-the-art Gaussian Process Regression approximations due to its
large size and its non-stationary features. The goal is in fact to predict the
delay of a flight given eight features, which are the age of the airplane,
route distance, airtime, departure time, arrival time, day of the week, day of
the month, and month. We follow the same exact preprocessing steps as in the
experiments in [16], [39] and [9], which consider subsets of data of size
$N=10000,100000,1000000,5929413$, each chosen uniformly at random. Training
the model is then accomplished with $\nicefrac{{2}}{{3}}N$ datapoints, with
the remaining portion reserved for testing. The entire procedure is then
repeated $10$ times with random initialization in order to obtain a sample
mean and sample standard deviation estimate of the MSE. Following exactly the
approach in Hilbert-GP [39], we consider $\hat{M}=40$ basis functions per
dimension, with the crucial difference that T-KRR approximates a standard
product Gaussian kernel, as opposed to an additive model. As was the case in
the previous section for classification, we select the lengthscale of the
kernel as the mean of the standard deviations of the eight features and choose
$\lambda_{N}=\nicefrac{{100}}{{N}}$ for the different splits, where the
dependency on $N$ is suggested by standard learning theory. We then train four
distinct models for $R=5,10,15,20$. Table 2 compares T-KRR with the best
performing models found in the literature, which all (except for GPR) rely on
low-rank approximation of the kernel function. The T-KRR method is able to
recover the baseline GPR performance already with $R=5$ and remarkably
outperform all other approaches although we choose the hyperparameters naively
and consider equal lengthscales $l$ for all dimensions. Interestingly, a
higher choice of $R$ does result in better performance in all cases except for
$N=10000$, where model performance very much depends on the random split (note
that we consider different splits for each experiment). The SVIGP [14]
achieves comparable results to T-KRR, indicating that a performance gain is
possible from product kernels. The difficulty with SVIGP is however that the
kernel function is interpolated locally at $M$ nodes in a data-dependent
fashion, requiring an increasing amount of interpolation nodes to cover the
whole domain to allow for good generalization. In contrast, T-KRR considers an
exponentially large amount of basis functions in the frequency domain, and
learns an efficient representation of the model weights. In light of the
performance of the additive Hilbert-GP and additive VFF models, we expect
similar performance when considering other feature maps which induce
stationary kernels. When considering the whole dataset, training our model
with $R=5$ takes $1565\pm 1$ seconds on a laptop, while for $R=20$ it takes
$7141\pm 245$ seconds. Reported training times of SVIGP indicate $18360\pm
360$ seconds [16] on a cluster. Since the computational complexity of our
algorithm is dominated by matrix-matrix multiplications, we expect significant
speedups when relying on GPU computations.

Table 2: Predictive MSE with one standard deviation for T-KRR. The prefix -A indicates that the model in question considers an additive kernel. $N$ | $10000$ | $100000$ | $1000000$ | $5929413$
---|---|---|---|---
T-KRR ($R=5$) | $0.91\pm 0.10$ | $0.82\pm 0.03$ | $0.80\pm 0.02$ | $0.800\pm 0.008$
T-KRR ($R=10$) | $0.89\pm 0.05$ | $0.80\pm 0.05$ | $0.79\pm 0.02$ | $0.785\pm 0.009$
T-KRR ($R=15$) | $0.90\pm 0.07$ | $0.80\pm 0.04$ | $0.78\pm 0.02$ | $0.773\pm 0.007$
T-KRR ($R=20$) | $0.97\pm 0.15$ | $\bm{0.78}\pm 0.04$ | $\bm{0.77}\pm 0.01$ | $\bm{0.763}\pm 0.007$
A-Hilbert-GP [39] | $0.97\pm 0.14$ | $0.80\pm 0.06$ | $0.83\pm 0.02$ | $0.827\pm 0.005$
A-VFF [16] | $0.89\pm 0.15$ | $0.82\pm 0.05$ | $0.83\pm 0.01$ | $0.827\pm 0.004$
SVIGP [14] | $0.89\pm 0.16$ | $0.79\pm 0.05$ | $0.79\pm 0.01$ | $0.791\pm 0.005$
VISH [9] | $0.90\pm 0.16$ | $0.81\pm 0.05$ | $0.83\pm 0.03$ | $0.834\pm 0.055$
GPR [33] | $0.89\pm 0.16$ | N/A | N/A | N/A
A-GPR [10] | $0.89\pm 0.16$ | N/A | N/A | N/A

##  5 Conclusion

In this work a framework to perform large-scale supervised learning with
tensor decompositions which leverages the tensor product structure of
deterministic Fourier features was introduced. Concretely, a monotonically
decreasing learning algorithm with linear complexity in both sample size and
dimensionality was derived. This algorithm leverages the efficient format of
the Canonical Polyadic Decomposition in combination with exponentially fast
converging deterministic Fourier features which allow to implicitly
approximate stationary product kernels up to machine precision. Numerical
experiments show how the performance of the baseline Kernel Ridge Regression
is recovered with a very limited number of parameters. The proposed method can
handle problems which are both large in the number of samples as well as in
their dimensionality, effectively enabling large-scale supervised learning
with stationary product kernels. The biggest limitation of the current
approach is that it is does not allow for uncertainty quantification, which
motivates further work in that direction.

## Funding transparency statement

Frederiek Wesel, and thereby this work, is supported by the Delft University
of Technology AI Labs program. The authors declare no competing interests.

## Acknowledgments

We would like to thank Arno Solin and Vincent Dutordoir for the help received
with the Airline dataset.

## References

  * Chen et al. [2018] Z. Chen, K. Batselier, J. A. K. Suykens, and N. Wong.  Parallelized Tensor Train Learning of Polynomial Classifiers.  _IEEE Transactions on Neural Networks and Learning Systems_ , 29(10):4621–4632, Oct. 2018.
  * Cichocki et al. [2016] 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_ , 9(4-5):249–429, 2016.
  * Comon et al. [2009] P. Comon, X. Luciani, and A. L. F. de Almeida.  Tensor decompositions, alternating least squares and other tales.  _Journal of Chemometrics_ , 23(7-8):393–405, July 2009.
  * Cortes and Vapnik [1995] C. Cortes and V. Vapnik.  Support-vector networks.  _Machine Learning_ , 20(3):273–297, Sept. 1995\.
  * Dao et al. [2017] T. Dao, C. D. Sa, and C. Ré.  Gaussian quadrature for kernel features.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17, pages 6109–6119, Red Hook, NY, USA, Dec. 2017. Curran Associates Inc.
  * De Lathauwer et al. [2000] L. De Lathauwer, B. De Moor, and J. Vandewalle.  A Multilinear Singular Value Decomposition.  _SIAM Journal on Matrix Analysis and Applications_ , 21(4):1253–1278, Jan. 2000.
  * Drineas and Mahoney [2005] P. Drineas and M. W. Mahoney.  On the Nystrom Method for Approximating a Gram Matrix for Improved Kernel-Based Learning.  _Journal of Machine Learning Research_ , 6(72):2153–2175, 2005.
  * Dua and Graff [2017] D. Dua and C. Graff.  UCI Machine Learning Repository.  http://archive.ics.uci.edu/ml, 2017.
  * Dutordoir et al. [2020] V. Dutordoir, N. Durrande, and J. Hensman.  Sparse Gaussian Processes with Spherical Harmonic Features.  In _International Conference on Machine Learning_ , pages 2793–2802. PMLR, Nov. 2020.
  * Duvenaud et al. [2011] D. Duvenaud, H. Nickisch, and C. E. Rasmussen.  Additive gaussian processes.  In _Proceedings of the 25th International Conference on Neural Information Processing Systems_ , NIPS’11, pages 226–234, Red Hook, NY, USA, 2011. Curran Associates Inc.
  * Garriga-Alonso et al. [2018] A. Garriga-Alonso, C. E. Rasmussen, and L. Aitchison.  Deep Convolutional Networks as shallow Gaussian Processes.  In _International Conference on Learning Representations_ , 2018.
  * Gilboa et al. [2015] E. Gilboa, Y. Saatçi, and J. P. Cunningham.  Scaling Multidimensional Inference for Structured Gaussian Processes.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 37(2):424–436, Feb. 2015.
  * Hammer and Gersmann [2003] B. Hammer and K. Gersmann.  A Note on the Universal Approximation Capability of Support Vector Machines.  _Neural Processing Letters_ , 17(1):43–53, Feb. 2003.
  * Hensman et al. [2013] J. Hensman, N. Fusi, and N. D. Lawrence.  Gaussian processes for Big data.  In _Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence_ , UAI’13, pages 282–290, Arlington, Virginia, USA, Aug. 2013. AUAI Press.
  * Hensman et al. [2015] J. Hensman, A. Matthews, and Z. Ghahramani.  Scalable Variational Gaussian Process Classification.  In _Artificial Intelligence and Statistics_ , pages 351–360. PMLR, Feb. 2015.
  * Hensman et al. [2017] J. Hensman, N. Durrande, and A. Solin.  Variational Fourier features for Gaussian processes.  _The Journal of Machine Learning Research_ , 18(1):5537–5588, Jan. 2017.
  * Hitchcock [1927] F. L. Hitchcock.  The Expression of a Tensor or a Polyadic as a Sum of Products.  _Journal of Mathematics and Physics_ , 6(1-4):164–189, 1927.
  * Holtz et al. [2012] S. Holtz, T. Rohwedder, and R. Schneider.  The alternating linear scheme for tensor optimization in the tensor train format.  _SIAM J. Sci. Comput._ , 34(2):A683–A713, 2012\.
  * Izmailov et al. [2018] P. Izmailov, A. Novikov, and D. Kropotov.  Scalable Gaussian Processes with Billions of Inducing Inputs via Tensor Train Decomposition.  In _International Conference on Artificial Intelligence and Statistics_ , pages 726–735. PMLR, Mar. 2018.
  * Kolda and Bader [2009] T. G. Kolda and B. W. Bader.  Tensor Decompositions and Applications.  _SIAM Review_ , 51(3):455–500, Aug. 2009.
  * Kruskal [1977] J. B. Kruskal.  Three-way arrays: Rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics.  _Linear Algebra and its Applications_ , 18(2):95–138, Jan. 1977.
  * Lee et al. [2018] J. Lee, Y. Bahri, R. Novak, S. S. Schoenholz, J. Pennington, and J. Sohl-Dickstein.  Deep Neural Networks as Gaussian Processes.  _arXiv:1711.00165 [cs, stat]_ , Mar. 2018.
  * Liu et al. [2021] F. Liu, X. Huang, Y. Chen, and J. A. K. Suykens.  Random features for kernel approximation: A survey on algorithms, theory, and beyond.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , PP, July 2021.
  * Matlab [2021] Matlab.  _9.10.0.1649659 (R2021a) Update 1_.  The MathWorks Inc., Natick, Massachussetts, 2021.
  * Mutný and Krause [2018] M. Mutný and A. Krause.  Efficient high dimensional Bayesian optimization with additivity and quadrature fourier features.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS’18, pages 9019–9030, Red Hook, NY, USA, Dec. 2018. Curran Associates Inc.
  * Novak et al. [2018] R. Novak, L. Xiao, Y. Bahri, J. Lee, G. Yang, J. Hron, D. A. Abolafia, J. Pennington, and J. Sohl-dickstein.  Bayesian deep convolutional networks with many channels are gaussian processes.  In _International Conference on Learning Representations_ , 2018.
  * Novikov et al. [2015] A. Novikov, D. Podoprikhin, A. Osokin, and D. P. Vetrov.  Tensorizing neural networks.  In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.
  * Novikov et al. [2017] A. Novikov, M. Trofimov, and I. Oseledets.  Exponential Machines.  _arXiv:1605.03795 [cs, stat]_ , Dec. 2017.
  * Oseledets [2011] I. V. Oseledets.  Tensor-Train Decomposition.  _SIAM Journal on Scientific Computing_ , 33(5):2295–2317, Jan. 2011.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and É. Duchesnay.  Scikit-learn: Machine Learning in Python.  _Journal of Machine Learning Research_ , 12(85):2825–2830, 2011.
  * Rahimi and Recht [2007] A. Rahimi and B. Recht.  Random features for large-scale kernel machines.  In _Proceedings of the 20th International Conference on Neural Information Processing Systems_ , NIPS’07, pages 1177–1184, Red Hook, NY, USA, Dec. 2007. Curran Associates Inc.
  * Rasmussen and Nickisch [2010] C. E. Rasmussen and H. Nickisch.  Gaussian Processes for Machine Learning (GPML) Toolbox.  _Journal of Machine Learning Research_ , 11(100):3011–3015, 2010.
  * Rasmussen and Williams [2006] C. E. Rasmussen and C. K. I. Williams.  _Gaussian Processes for Machine Learning_.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, Mass, 2006.
  * Rohwedder and Uschmajew [2013] T. Rohwedder and A. Uschmajew.  On local convergence of alternating schemes for optimization of convex problems in the tensor train format.  _SIAM Journal on Numerical Analysis_ , 51(2):1134–1162, 2013.
  * Rudin [1962] W. Rudin.  _Fourier Analysis on Groups._ Number 12 in Interscience Tracts in Pure and Applied Mathematics. Interscience Publishers, New York, 1962.
  * Shustin and Avron [2021] P. F. Shustin and H. Avron.  Gauss-Legendre Features for Gaussian Process Regression.  _arXiv:2101.01137 [cs, math]_ , Jan. 2021.
  * Sidiropoulos and Bro [2000] N. D. Sidiropoulos and R. Bro.  On the uniqueness of multilinear decomposition of N-way arrays.  _Journal of Chemometrics_ , 14(3):229–239, 2000\.
  * Sidiropoulos et al. [2017] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos.  Tensor Decomposition for Signal Processing and Machine Learning.  _IEEE Transactions on Signal Processing_ , 65(13):3551–3582, July 2017.
  * Solin and Särkkä [2020] A. Solin and S. Särkkä.  Hilbert space methods for reduced-rank Gaussian process regression.  _Statistics and Computing_ , 30(2):419–446, Mar. 2020.
  * Stoudenmire and Schwab [2016] E. M. Stoudenmire and D. J. Schwab.  Supervised learning with tensor networks.  In _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , NIPS’16, pages 4806–4814, Red Hook, NY, USA, Dec. 2016. Curran Associates Inc.
  * Suykens and Vandewalle [1999] J. Suykens and J. Vandewalle.  Least Squares Support Vector Machine Classifiers.  _Neural Processing Letters_ , 9(3):293–300, June 1999.
  * Suykens et al. [2002] J. A. K. Suykens, T. Van Gestel, J. De Brabanter, B. De Moor, and J. Vandewalle.  _Least Squares Support Vector Machines_.  World Scientific, Nov. 2002.
  * Tjandra et al. [2017] A. Tjandra, S. Sakti, and S. Nakamura.  Compressing recurrent neural network with tensor train.  In _2017 International Joint Conference on Neural Networks (IJCNN)_ , pages 4451–4458, May 2017.
  * Tucker [1966] L. R. Tucker.  Some mathematical notes on three-mode factor analysis.  _Psychometrika_ , 31(3):279–311, Sept. 1966.
  * Wahls et al. [2014] S. Wahls, V. Koivunen, H. V. Poor, and M. Verhaegen.  Learning multidimensional Fourier series with tensor trains.  In _2014 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_ , pages 394–398, Dec. 2014.
  * Williams and Seeger [2001] C. Williams and M. Seeger.  Using the nyström method to speed up kernel machines.  In _Proceedings of the 14th International Conference on Neural Information Processing Systems_ , NIPS’01, pages 682–688. MIT Press, 2001.
  * Yang et al. [2012] T. Yang, Y.-f. Li, M. Mahdavi, R. Jin, and Z.-H. Zhou.  Nyström method vs random fourier features: A theoretical and empirical comparison.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 25. Curran Associates, Inc., 2012.
  * Yang et al. [2017] Y. Yang, D. Krompass, and V. Tresp.  Tensor-Train Recurrent Neural Networks for Video Classification.  In _International Conference on Machine Learning_ , pages 3891–3900. PMLR, July 2017.

##  Appendix A Appendix

Derivation of the data-fitting term of (12). We make use of the multi-
linearity property of the CPD and rely on re-ordering the summations:

$\displaystyle\langle\bm{\mathcal{W}},\bm{\mathcal{Z}}(\bm{x})\rangle_{\text{F}}$
$\displaystyle=\left\langle\sum_{r=1}^{R}\bm{w}^{\left(1\right)}_{r}\otimes\bm{w}^{\left(2\right)}_{r}\otimes\cdots\otimes\bm{w}^{\left(D\right)}_{r},\bm{z}^{(1)}\otimes\bm{z}^{(2)}\cdots\otimes\bm{z}^{(D)}\right\rangle_{\text{F}}$
$\displaystyle=\sum_{i_{1}=1}^{\hat{M}}\cdots\sum_{i_{d}=1}^{\hat{M}}\cdots\sum_{i_{D}=1}^{\hat{M}}\sum_{r=1}^{R}w^{\left(1\right)}_{i_{1}r}\,z^{\left(1\right)}_{i_{1}}\cdots
w^{\left(d\right)}_{i_{d}r}\,z^{\left(d\right)}_{i_{d}}\cdots
w^{\left(D\right)}_{i_{D}r}\,z^{\left(D\right)}_{i_{D}}$
$\displaystyle=\sum_{i_{d}=1}^{\hat{M}}\sum_{r=1}^{R}w^{\left(d\right)}_{i_{d}r}\left(z^{\left(d\right)}_{i_{d}}\sum_{i_{1}=1}^{\hat{M}}w^{\left(1\right)}_{i_{1}r}\,z^{\left(1\right)}_{i_{1}}\cdots\sum_{i_{D}=1}^{\hat{M}}w^{\left(D\right)}_{i_{D}r}\,z^{\left(D\right)}_{i_{D}}\right)$
$\displaystyle=\text{vec}\left(\bm{W}^{\left(d\right)}\right)^{\text{T}}\left(\bm{z}^{\left(d\right)}\otimes\left({\bm{z}^{\left(1\right)}}^{\text{T}}{\bm{W}^{\left(1\right)}}^{\text{T}}\odot\cdots\odot{\bm{z}^{\left(D\right)}}^{\text{T}}{\bm{W}^{\left(D\right)}}^{\text{T}}\right)\right)$
$\displaystyle=\left\langle\text{vec}\left(\bm{W}^{\left(d\right)}\right),\bm{g}^{\left(d\right)}\left(\bm{x}\right)\right\rangle$

The derivation of the regularization term of (13) follows a similar reasoning
as for the data-fitting term:

$\displaystyle\left\langle\bm{\mathcal{W}},\bm{\mathcal{W}}\right\rangle_{\text{F}}$
$\displaystyle=\left\langle\sum_{r=1}^{R}\bm{w}^{\left(1\right)}_{r}\otimes\bm{w}^{\left(2\right)}_{r}\otimes\cdots\otimes\bm{w}^{\left(D\right)}_{r},\sum_{r=1}^{R}\bm{w}^{\left(1\right)}_{r}\otimes\bm{w}^{\left(2\right)}_{r}\otimes\cdots\otimes\bm{w}^{\left(D\right)}_{r}\right\rangle_{\text{F}}$
$\displaystyle=\sum_{i_{1}=1}^{\hat{M}}\cdots\sum_{i_{d}=1}^{\hat{M}}\cdots\sum_{i_{D}=1}^{\hat{M}}\sum_{r=1}^{R}\sum_{p=1}^{R}w^{\left(1\right)}_{i_{1}r}\,w^{\left(1\right)}_{i_{1},p}\cdots
w^{\left(d\right)}_{i_{d}r}\,w^{\left(d\right)}_{i_{d},p}\cdots
w^{\left(D\right)}_{i_{D}r}\,w^{\left(D\right)}_{pi_{D}}$
$\displaystyle=\sum_{r=1}^{R}\sum_{p=1}^{R}\left(\sum_{i_{d}=1}^{\hat{M}}w^{\left(d\right)}_{i_{d}r}\,w^{\left(d\right)}_{i_{d}p}\right)\left(\sum_{i_{1}=1}^{\hat{M}}w^{\left(1\right)}_{i_{1}r}\,w^{\left(1\right)}_{i_{1}p}\cdots\sum_{i_{D}=1}^{\hat{M}}w^{\left(D\right)}_{i_{D}r}\,w^{\left(D\right)}_{i_{D}p}\right)$
$\displaystyle=\sum_{r=1}^{R}\sum_{p=1}^{R}\left(\bm{w}^{{\left(d\right)}^{\text{T}}}_{r}\bm{w}_{p}^{\left(d\right)}\right)\left(\bm{w}^{{\left(1\right)}^{\text{T}}}_{r}\bm{w}_{p}^{\left(1\right)}\odot\cdots\odot\bm{w}^{{\left(D\right)}^{\text{T}}}_{r}\bm{w}_{p}^{\left(D\right)}\right)$
$\displaystyle=\text{vec}\left(\bm{W}^{{\left(d\right)}^{\text{T}}}\bm{W}^{\left(d\right)}\right)^{\text{T}}\text{vec}\left(\bm{W}^{{\left(1\right)}^{\text{T}}}\bm{W}^{\left(1\right)}\odot\cdots\odot\bm{W}^{{\left(D\right)}^{\text{T}}}\bm{W}^{\left(D\right)}\right)$
$\displaystyle=\left\langle\text{vec}\left(\bm{W}^{{\left(d\right)}^{\text{T}}}\bm{W}^{\left(d\right)}\right),\text{vec}\left(\bm{H}^{\left(d\right)}\right)\right\rangle$