As mentioned above, the functional tensor-train decomposition was proposed as an ansatz for representing functions by Oseledets (2013), and computational routines for compressing a function into FT format have been developed before (Gorodetsky et al., 2016). In that setting, an approximation of a black-box function is sought to a prescribed accuracy. A sampling procedure and associated algorithm was designed to optimally evaluate the function in order to obtain an FT representation. In this work, we consider the setting of fixed data.

There has also been some recent work on regression in low-rank format (Doostan et al., 2013; Mathelin, 2014; Chevreuil et al., 2015). These approaches rely on linearity between the parameters of the low-rank format and the output of the function. Utilizing this relationship, they convert the low-rank function approximation to one of low-rank tensor decomposition for the coefficients of a tensor-product basis. In Section 3.2 we show how the representation presented in those works can be obtained as a particular instantiation of the formulation we present here.

In spirit, our approach is also similar to the recent use of the tensor-train decomposition within fully connected neural networks (Novikov et al., 2015). There, the layers of a neural network are assumed to be a linear transformation of an input vector, and the weight matrix of the transformation is estimated from data. Their contribution is representing the weight matrix in a compressed TT-format. In this work, our low-rank format can be thought of as an alternative to the linear transformation layer of a neural network. Indeed, future work can focus on utilizing our proposed methodology within a hierarchy of low-rank approximations.

Let $\mathbb{R}$ be the set of real numbers and $\mathbb{Z}_{+}$ be the set of positive integers. Let $n\in\mathbb{Z}_{+}$, $d\in\mathbb{Z}_{+}$ and suppose that we are given i.i.d. data $\left(x^{(i)},y^{(i)}\right)_{i=1}^{n}$ such that each datum is sampled from a distribution $\mu_{xy}$ on a compact space $(x^{(i)},y^{(i)})\in{\cal X}\times{\cal Y}\subset\mathbb{R}^{d}\times\mathbb{R}.$ Let ${\cal X}={\cal X}_{1}\times\cdots\times{\cal X}_{d}$ be a tensor product space with ${\cal X}_{k}\subset\mathbb{R}$. Let the marginal distribution of $x=(x_{1},\ldots,x_{d})\in{\cal X}$ be the tensor product measure $\mu_{x}=\mu_{x_{1}}\times\cdots\times\mu_{x_{d}}$, and assume that all integrals appearing in this paper are with respect to this measure. For example let $f:{\cal X}\to{\cal Y}$ and $g:{\cal X}\to{\cal Y}$, then the inner product is defined as $\left\langle f,g\right\rangle=\int f(x)g(x)d\mu(x).$ Similarly, the $L_{2}$ norm is defined as $\lVert f\rVert^{2}_{2}=\left\langle f,f\right\rangle.$

In this paper scalars and scalar-valued functions are denoted by lowercase letters, vectors are denoted by bold lower-case letters such as $\mathbf{x},\mathbf{y}$; matrices are denoted with upper boldface such as $\mathbf{X},\mathbf{Y}$; tensors are denoted with upper boldface caligraphic letters such as $\bm{\mathcal{X}},\bm{\mathcal{Y}}$; and matrix-valued functions are denoted by upper caligraphic letters such as ${\cal F},{\cal G}.$ An ordered sequence of elements of the same set are distinguished by parenthesized superscripts such as $x^{(i)}$ or $y^{(i)}$.

The supervised learning task seeks a function $f:{\cal X}\to{\cal Y}$ that minimizes the average of a cost function $g_{f}:{\cal X}\times{\cal Y}{}\to\mathbb{R}$

For example, the standard least squares objective is specified with $g_{f}(x,y)=\left(y-f(x)\right)^{2}$. In this work, the cost functional cannot be exactly calculated; instead, data $\left(x^{(i)},y^{(i)}\right)_{i=1}^{n}$ is used to estimate its value. The cost functional then becomes a sum over the data instead of an integral

where we have reused the notation $J[f]$, and further references to $J$ will use this definition. To obtain an optimization problem over a finite set of variables, the search space of functions must be restricted. Nonparametric representations generally allow this space to vary depending upon the data, and parametric representations typically fix the representation. In either case, we denote this space as $\mathscr{F}$ to seek an optimal function $f^{*}\in\mathscr{F}$ such that

One example of a common function space is the reproducing kernel Hilbert space, and resulting algorithms using this space include Gaussian process regression or radial basis function interpolation. Other examples include linear functions (resulting in linear regression) or polynomial functions. In this work, we consider a function space that incorporates all functions with a prescribed rank, which will be defined in the next section.

When solving the supervised learning problem (4) it is often useful to introduce a regularization term to minimize overfitting. In this paper we will consider the following regularized learning problem

where $\lambda\in\mathbb{R}_{+}$ denotes a scaling factor and $\Omega:\mathscr{F}\to\mathbb{R}_{+}$ is a functional that penalizes certain functions in $\mathscr{F}$. For example, the function $\Omega(f)=\lVert f\rVert_{2}$ penalizes functions that have large $L_{2}$ norms, and such a penalty has been used for a certain type of low-rank functional approximation before (Doostan et al., 2013). In this work, we impose a group sparsity type regularization that has been shown to improve the prediction performance in other approximation settings (Turlach et al., 2005; Yuan and Lin, 2006). Such an approach seeks to increase parsimony by reducing the number of non-zero elements in an approximation. In Section 4.1 we describe what this type of constrains means in the context of low-rank functional decompositions.

The function space that we use to constrain our supervised learning problem is related to the concept of low-rank decompositions of multiway arrays, or tensors. Tensor decompositions are multidimensional analogues of the matrix SVD and are used to mitigate the curse of dimensionality associated with storing multivariate arrays (Kolda and Bader, 2009). Consider an array with $\bm{\mathcal{A}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{d}}$, this array contains an exponentially growing number of elements as the dimension increases. Tensor decompositions can provide a compact representation of the tensor and have found wide spread usage for data compression and learning tasks (Cichocki et al., 2009; Novikov et al., 2015; Yu and Lie, 2016).

In this work, we consider the tensor-train (TT) decomposition (Oseledets, 2011). Specifically, we use the TT as an ansatz for representating low-rank functions. A TT representation of a tensor $\bm{\mathcal{A}}$ is defined by a list of $3$-way arrays $\left(\bm{\mathcal{A}}_{k}\in\mathbb{R}^{r_{k-1}\times p_{k}\times r_{k}}\right)_{k=1}^{d}$ called TT-cores where $r_{0}=r_{d}=1$ and $r_{k}\in\mathbb{Z}_{+}$ for $k=2,\ldots,d-1$. The sizes of the cores $\left(r_{k}\right)_{k=0}^{d}$ are called the TT-ranks. In this format a tensor element is obtained according to the following multiplication

where $\bm{\mathcal{A}}_{k}[i_{k}]\in\mathbb{R}^{r_{k-1}\times r_{k}}$ are matrices and the above equation describes a sequence of matrix multiplication.

In this format, storage complexity scales linearly with dimension and quadratically with TT-rank. In the next section we discuss how it can be used as a framework for function approximation.

An FT core is comprised of $d$ sets of univariate functions as shown in Figure 1. Each set of univariate functions, contains all of the parameters of the associated univariate functions. As a result the full FT is parameterized through the parameterization of its univariate functions. Let $p_{kij}\in\mathbb{Z}_{+}$ denote the number of parameters describing ${f}_{k}^{(ij)}.$ Let $\bm{\theta}$ denote the vector of parameters of all the univariate functions. Then, there are a total of $\sum_{k=1}^{d}\sum_{i=1}^{r_{k-1}}\sum_{j=1}^{r_{k}}p_{kij}$ parameters describing the FT, i.e., $\bm{\theta}\in\mathbb{R}^{p_{t}}$.

The parameter vector $\bm{\theta}$ is indexed by a multi-index $\bm{\alpha}=(k,i,j,\ell)$ where $k=1,\ldots,d$ corresponds to an input variable, $i=1,\ldots,r_{k-1}$ and $j=1,\ldots,r_{k}$ correspond to a univariate function within the $k$th core, and $\ell=1,\ldots,p_{kij}$ corresponds to a specific parameter within that univariate function. In other words, we adopt the convention that $\bm{\theta}_{\bm{\alpha}}=\bm{\theta}_{kij\ell}$ refers to the $\ell$th parameter of the univariate function in the $i$th row and $j$th column of the $k$th core.

The additional flexibility of the FT representation allows both linear and nonlinear parameterizations of univariate functions. A linear parameterization represents a univariate function as an expansion of basis functions $\left({\phi}_{k\ell}^{(ij)}:{\cal X}_{k}\to\mathbb{R}\right)_{\ell=1}^{p_{kij}}$ according to

Linear parameterizations are often convenient within the context of gradient-based optimization methods since the gradient with respect to a parameter is independent of the parameter values. Thus, one only needs to compute and store the evaluation of the basis functions a single time. Nonlinear parameterizations are more flexible and general, but often incur a greater computational cost within optimization. One example of a nonlinear parameterization is that of a Gaussian kernel, which can be written as

where $\sigma>0$ Here, the first $p_{kij}/2$ parameters refer to the coefficients of radial basis functions and that second half of the parameters refer to the centers.

The functional tensor train can be used by independently parameterizing the univariate functions of each core, and both linear and nonlinear parameterizations are possible. As described below, the advantage of this representation includes a naturally sparse storage scheme for the cores. Another advantage is the availablitiy of efficient computational algorithms for multilinear algebra that can adapt the representation of each univariate function individually as needed (Gorodetsky et al., 2016; Gorodetsky, 2017) in the spirit of continuous computation pioneered by Chebfun (Platte and Trefethen, 2010).

Although the functional form of the TT is general, most of the literature makes two simplifying assumptions (Doostan et al., 2013; Mathelin, 2014; Chevreuil et al., 2015):

1. 1.

   a linear parameterization of each ${f}_{k\ell}^{(ij)}$, and

2. 2.

   an identitical basis for the functions within each FT core, i.e., $p_{kij}=p_{k}$ and ${\phi}_{k\ell}^{(ij)}=\phi_{k\ell}$ for all $i=1,\ldots r_{k-1}$, $j=1,\ldots,r_{k},$ and $\ell=1,\ldots,p_{k}$.

These assumptions transform the problem of learning low-rank functions to the problem of learning low-rank coefficients. This transformation of the problem allows one to facilitate the use of discrete TT algorithms and theory. We will refer to representations using these assumptions as a functional tensor-train with TT coefficients or FT-c, where “c” stands for coefficients.

The FT-c representation stores the coefficients of a tensor-product basis $\phi_{k\ell}$ for all $k$ and $\ell$ in TT format. Let $\bm{\mathcal{F}}_{k}\in\mathbb{R}^{r_{k-1}\times p_{k}\times r_{k}}$ be a tensor of the following form

for $\ell=1,\ldots,p_{k}$. Function evaluations can be obtained from the from coefficients stored in TT format by performing the the following summation

From (8), (12), and (13) we can see that the relationship between the TT cores $\bm{\mathcal{F}}_{k}$ and the FT cores ${\cal F}_{k}$ is

where the basis function multiplies every element of the tensor. In other words, the TT cores represent a TT decomposition of the $p_{1}\times p_{2}\times\cdots\times p_{d}$ coefficient tensor of a tensor-product basis.

The two assumptions required for converting the problem from one of low-rank function approximation to low-rank coefficient approximation often result in a computational burden in practice. The burden of the first assumption is clear, some functions are more compressible if nonlinear parameterizations are allowed. The second assumption can also limit compressibility, but also results in a larger computational and storage burden. One example in which this burden becomes obvious is when attempting to represent additively separable functions, e.g., Equation (1), in the FT format. As mentioned in the introduction, this representation is common for high-dimensional regression, and an FT can represent this function using cores that fave the following rank 2 structure.

Suppose that each of the univariate functions can be represented with $p_{k}$ parameters and the constants $0$ and $1$ can be stored with a single parameter. Then, the storage requirement is $p_{k}+3$ parameters for the $d-2$ middle cores and $p_{k}+1$ for the outer cores, totaling $(p_{kij}+3)d-4$ floating point values. In the TT case, the $0$ and $1$ terms must be stored with the same complexity as the other terms, since the TT is a three-way array. Thus, the total number of parameters becomes $4p_{k}(d-2)$. Almost four times less storage is required for the FT format than the TT format, in the limit of large $p_{k}$ and $d$. Essentially, the TT format does take into account the sparsity of the cores. In this case, one can think of the FT format as storing the TT cores in a sparse format. This burden is exacerbated if we add interaction terms to Equation 1.

The function space $\mathscr{F}$ in Equation (4) constrains the search space. The space constrained by low-rank functions $\mathscr{F}_{\mathbf{r}}$ is defined as follows. Let $\mathbf{r}=[1,r_{1},\ldots,r_{d-1},1]$ such that $r_{i}\in\mathbb{Z}_{+}$ for $i=1,\ldots,d-1$. Then

denotes the space of functions with FT ranks $\mathbf{r}.$ Note that this function space also includes functions with smaller ranks. For example, a function that differs in ranks $\hat{r}_{1},\hat{r}_{2}$ such that $\hat{r}_{1}<r_{1}$ and $\hat{r}_{2}<r_{2}$ can be obtained by setting all univariate functions, aside from the top left $\hat{r}_{1}\times\hat{r}_{2}$ block, to zero; analogous modifications must be made to all other cores.

In this paper we also add a regularization to attempt to minimize the number of nonzero functions in each core. The structure of the FT-cores in Equation (9) admits a natural grouping in terms of the univariate functions. Setting a univariate function to zero inside of a core essentially lowers the number of interactions allowable between univariate functions in neighbouring dimensions, with the overall effect beign a smaller number of sums in the representation of the function in Equation (7). Minimizes the number of nonzero univariate functions not only produces a more parsimonious model, but also minimizes the number of interactions between function variables.

Specifically, we penalize the regression problem using the sum of the norms of the univariate functions

Note that this regularization method is different from that first proposed in (Mathelin, 2014) where the FT-c representation was used and the TT-cores of the coefficient tensor themselves were forced to be sparse. In our case we do not look for sparse coefficients, instead our goal is to increase the number of functions that are identically zero. As a surrogate for this goal we follow the standard practice of replacing the non-differential “$L_{0}$” norm by the another norm, in our case a sum of the norms of the functions. This replacement also seeks to directly minimize the number of interacting univariate functions and also provides a differentiable objective for optimization.

A careful selection of the number of parameters $p_{kij}$ in each univariate function, the rank $r_{k}$ of each core, and the Lagrangian multiplier $\lambda$ in the regularized learning problem is required to avoid over fitting. For example if a polynomial basis is chosen, setting the degree too high can result in spurious oscillations. Cross validation provides a mechanism to estimate hyper-parameters. In this paper we use 20 fold cross validation to select the hyper-parameters of an FT that minimizes the expected prediction error.

Alternatively, we can design a more effective rank adaptation scheme by combining cross validation with a rounding procedure. The problem with the simplest cross validation option is that a scheme that optimizes separately over each $r_{k}$ imposes a computational burden, and a scheme that optimizes for a single rank across all cores, i.e., $r_{k}=r$ does not allow enough flexibility. Instead, we propose to combine FT rounding (Gorodetsky et al., 2016) and cross validation to avoid overfitting.

Rounding is a procedure to reapproximate an FT by one with lower ranks to a given tolerance. The tolerance criterion can be thought of as a regularization term, since it limits the complexity of the representation and will be investigated in furture work. The full rank adaptation scheme is provided by Algorithm 1. In that algorithm $\texttt{cv}(\mathbf{r})$ refers to a function that provides a cross-validation error for an optimization over the space $\mathscr{F}_{\mathbf{r}}$, and the function $\texttt{rounding-rank}(f,\delta)$ provides the ranks of a rounded function $f$ with a particular tolerance $\delta$.

The scheme increases ranks until either the cross-validation error increases or until the rounding procedure decreases every rank. The first termination criterion targets overfitting, and the second termination criterion seeks to limit rank growth when data is no longer informative enough to necessitate the increase in expressivity.

In this section, we overview three gradient-based optimization algorithms for solving supervised learning problems in low-rank format: alternating minimization, gradient decent and stochastic gradient decent. The computational complexity of these algorithms in analyzed in Section 5.

In this section we discuss how to obtain the partial derivatives with respect to parameters of a specific core. Without loss of generality , we will make the following assumption to ease the notational burden

Under this assumption the FT cores have the following structure.

Now if we let $G(p)$ denote the number of operations required to compute the gradient of a univariate function with respect to all of its parameters then the cost of computing $\partial_{\bm{\alpha}}{\cal F}_{k}(x_{k})$ is $\mathcal{O}\left(G(p)r_{k-1}r_{k}\right)$ operations.

Computation of the partial derivative of the FT (20) requires the evaluation of the left and right product cores ${\cal F}_{<k}(x_{<k})$ and ${\cal F}_{>k}(x_{>k})$. A single forward and backward sweep the cores can be used to obtain these values using the following recursion identities

Thus we can obtain the gradient of the FT with respect to parameters of the $k$th core using the Algorithms 3 and 4.

These two algorithms consist of a triple nested loop. The innermost loop of coregrad-left requires computing the gradient of a univariate function with respect to its parameters, and multipling each element of the gradient by the left multiplier. Thus if $r_{k}<r$ then $\mathcal{O}(G(p)r^{2})$ operations are needed for the gradient and $\mathcal{O}(pr^{2})$ products between floating point operations are needed. Finally, the partial derivative with respect to each parameter of the core is stored with a complexity of $\mathcal{O}(pr^{2})$ floating point numbers. Each of these partial derivatives is updated in coregrad-right for a computational cost of $\mathcal{O}(pr^{2})$ and no additional storage. Since each of these functions is called $d$ times, the additional computational cost they incur is $\mathcal{O}(dr^{2}(G(p)+p))$ and the additional storage complexity they incur is $\mathcal{O}(dpr^{2})$.

Th entire gradient of the FT can be obtained with a single forward and backward sweep as presented in Algorithm 5. In Algorithm 5 we use $\partial f_{k\bm{\alpha}_{k}}$ to denote the partial derivative of $f$ with respect to the parameters of the $k$th core. Each step of the forward sweep requires: (i) creating an intermediate result for the gradient of the FT with respect to each parameter of the core in line 4 of Algorithm 3, (ii) evaluating and storing the core of the FT in line 6 of Algorithm 5, and (iii) updating the product of the cores through the current step in line 7 of Algorithm 5). The backward sweep involves involves updating the the intermediate gradient result obtained from the forward sweep (line 4 in Algorithm 4), and then updating the product of the cores in line 14 of Algorithm 5.

The values $G(p)$ and $E(p)$ are dependent upon the types of parameterizations of univariate functions. Consider two examples: one linear and one nonlinear. For both parameterizations we consider kernel basis functions; however, for the nonlinear parameterization we will consider the centers of each kernel as free parameters. The linear parameterization is given by Equation (10) with ${\phi}_{k\ell}^{(ij)}(x_{k})=\exp\left(-\frac{1}{\sigma^{2}}\left(x_{k}-c_{k\ell}\right)^{2}\right)$, where $\sigma\in\mathbb{R}_{+}$ and $c_{k\ell}\in{\cal X}_{k}$ are the centers of the kernel such that each univariate function of the $k$th dimension is represented with as a sum of kernels with the same locations. If the evaluation of the exponential takes a constant number of operations with respect to $x_{k}$ then we have $E(p)=\mathcal{O}(p)$ and $G(p)=\mathcal{O}(p)$ because the the gradient of the univariate function with respect to its $\ell$th parameter is

This gradient is independent of any other parameters. In practice this means that it can be precomputed at each location $x_{k}$ and recalled at runtime. In such a case the storage increases to $\mathcal{O}(ndpr^{2})$ numbers if each univariate function in each core has a different parameterization. If the univariate functions of each core share the same parameterization then the additional storage cost is $\mathcal{O}(npd)$. In either case the online cost becomes a simple lookup, i.e., $G(p)=\mathcal{O}(1)$.

The nonlinear parameterization provided in Equation (11) is different because we are free to optimize over the centers. In this case the gradient of each univariate function becomes The gradient with respect to the second half of the parameters now depends on the parameter value

The parameter dependence of the gradient increases the complexity of optimization algorithms since precomputation of the gradients cannot be performed. In this case we have $E(p)=\mathcal{O}(p)$ and $G(p)=\mathcal{O}(p)$.

The complexity of evaluating the gradient of the least squares objective function

is dominated by the cost of this derivative is evaluating the gradient of $f$ at $n$ points. Consequently, the total complexity is $n$ times that provided by Proposition 2, for a total cost of $\mathcal{O}(ndr^{2}(G(p)+E(p)))$ operations. The storage cost need not increase because one can evaluate the sum by sequentially iterating through each sample.

Now that we have described the computational complexity of the gradient computation, we summarize the computational complexity of the proposed optimization algorithms. Table 1 shows the computational complexity of the all-at-once optimization scheme when using a low-memory BFGS solver. Three parameterizations are considered a linear parameterization with identical parameterizations of each univariate function in a particular core, a more general parameterization with varying linear parmaeterizations within each core and, and the most general case of nonlinear parameterization of each function. We see that using linear parameterizations allows us to precompute certain gradients before running the optimizer. This precomputation reduces the online cost of optimization. The computational complexity of the full solution is dominated by the number of evaluations and gradients of the objective function, and we denote this number as $N_{\textrm{opt,AAO}}$.

The computational cost of stochastic gradient descent is given in Table 2. In this case the cost per epoch (once through all of the training points) is the same as a single gradient evaluation in the all-at-once approach. The total cost of such a scheme is dominated by how many samples are used during the optimization procedures. In the associated table, we represent this number as the number of times one requires iterating through all of the samples $N_{\textrm{epoch}}n$. A fully online algorithm would not have any associated offline cost and its complexity would be equivalent to the nonlinear parameterization case.

Finally, the alternating optimization scheme is of a slightly different nature. In the case of linear parameterization, each sub-optimization problem is quadratic and can be posed as solving a linear system by setting the right hand side of Equation (17) to zero. This linear system has $n$ equations and $pr^{2}$ unknowns and its solution, using a direct method, has an asymptotic complexity of $\mathcal{O}\left(np^{2}r^{4}\right)$. We also need to introduce a new constant called $N_{\textrm{sweep}}$ that represents how many sweeps through all of the dimensions are required. We see that in the most general case, this algorithm is $N_{sweep}$ times more expensive than all-at-once. However, this number is a bit desceptive since each sub problem has only $pr^{2}$ unknowns and therefore we can assume that $N_{\textrm{opt,ALS}}<N_{\textrm{opt,AAO}}$. In Table 3 we summarize the costs of this algorithm.