Tensor-based algorithms for image classification

The interest in machine learning with tensor networks has been growing rapidly in recent years. We show that tensor-based methods developed for learning the governing equations of dynamical systems from data can, in the same way, be used for supervised learning problems and propose two novel approaches for image classification. One is a kernel-based reformulation of the previously introduced MANDy (multidimensional approximation of nonlinear dynamics), the other an alternating ridge regression in the tensor-train format. We apply both methods to the MNIST and fashion MNIST data set and show that the approaches are competitive with state-of-the-art neural network-based classifiers.


Introduction
Tensor-based methods have become a powerful tool for scientific computing over the last years. In addition to many application areas such as quantum mechanics and computational dynamics where low-rank tensor approximations have been successfully applied, using tensor networks for supervised learning has gained a lot of attention recently. In particular the canonical format and the tensor-train format have been considered for quantum machine learning 1 problems, see, e.g., [1,2,3]. A tensor-based algorithm for image classification using sweeping techniques inspired by the density matrix renormalization group (DMRG) [4] was proposed in [5,6] and further discussed in [7,8]. Interestingly, also researchers at Google are currently developing a tensor-based machine learning framework called TensorNetwork 2 [9,10]. The goal is to expedite the adoption of such methods by the machine learning community.
Our goal is to show that recently developed methods for recovering the governing equations of dynamical systems can be generalized in such a way that they can also be used for supervised learning tasks, e.g., classification problems. In order to learn the governing equations from simulation or measurement data, regression methods such as SINDy (sparse identification of nonlinear dynamics) [11,12] and its tensor-based reformulation MANDY (multidimensional approximation of nonlinear dynamics) [13] can be applied. The main challenge is often to choose the right function space from which the system representation is learned. While SINDy and MANDy essentially select functions from a potentially large set of basis functions by applying regularized regression methods, other approaches allow nested functions and typically result in nonlinear optimization problems, which are then frequently solved using (stochastic) gradient descent. By constructing a basis comprising tensor products of simple functions (e.g., functions depending only on one variable), extremely high-dimensional feature spaces can be generated.
In this work, we explain how to compute the pseudoinverse required for solving the minimization problem directly in the tensor-train (TT) format, i.e., we replace the iterative approach from [5,6] by a direct computation of the least-squares solution and point out similarities with the aforementioned system identification methods. The reformulated algorithm can be regarded as a kernelized variant of MANDy, where the kernel is based on tensor products. This is also related to quantum machine learning ideas: As pointed out in [14], the basic idea of quantum computing is similar to kernel methods in that computations are performed implicitly in otherwise intractably large Hilbert spaces. While kernel methods were popular in the 1990s, the focus of the machine learning community has shifted to deep neural networks in recent years [14]. We will show that for simple image classification tasks, kernels based on tensor products are competitive with neural networks.
In addition to the kernel-based approach, we propose another DMRG-inspired method for the construction of TT decompositions of weight matrices containing the coefficients for the selected basis functions. Instead of computing pseudoinverses, a core-wise ridge regression [15] is applied to solve the minimization problem. While the approach introduced in [5,6] only involves tensor contractions corresponding to single images of the training data set, we use TT representations of transformed data tensors, see [13,16], in order to include the entire training data set at once for constructing low-dimensional systems of linear equations. Combining an efficient computational scheme for the corresponding subproblems and truncated singular value decompositions [17], we call the resulting algorithm alternating ridge regression (ARR) and discuss connections to MANDy and other regularized regression techniques.
While we describe the classification problems using the example of the iconic MNIST data set [18] and the fashion MNIST data set [19], the derived algorithms can be easily applied to other classification problems. There is a plethora of kernel and deep learning methods for image classification, a list of the most successful methods for the MNIST and fashion MNIST data sets including nearest-neighbor heuristics, support vector machines, and convolutional neural networks can be found on the respective website. 3,4 We will not review these methods in detail, but instead focus on relationships with data-driven methods for analyzing dynamical systems. The main contributions of this paper are: • Extension of MANDy : We show that the efficacy of the pseudoinverse computation in the tensor-train format can be improved by eliminating the need to left-and right-orthonormalize the tensor. While this is a straightforward modification of the original algorithm, it enables us to consider large data sets. The resulting method is closely related to kernel ridge regression.
• Alternating ridge regression: We introduce a modified TT representation of transformed data tensors for the development of a tensor-based regression technique which computes low-rank representations of coefficient tensors. We show that it is possible to obtain results which are competitive with those computed by MANDy, and, at the same time, reduce the computational costs and the memory consumption significantly.
• Classification of image data: While originally designed for system identification, we apply these methods to classification problems and visualize the learned classifier, which allows us to interpret features detected in the images.
The remainder is structured as follows: In Section 2, we will describe methods to learn governing equations of dynamical systems from data as well as a tensor-based iterative scheme for image classification and highlight their relationships. In Section 3, we describe how to apply MANDy to classification problems and introduce the ARR approach based on the alternating optimization of TT cores. Numerical results are presented in Section 4, followed by a brief summary and conclusion in Section 5.

Prerequisites
We will introduce the original MNIST and the fashion MNIST data set, which will serve as guiding examples. Afterwards SINDy and MANDy as well as tensor-based methods for image classification problems will be briefly discussed. In what follows, we will use the notation summarized in Table 1.

MNIST and fashion MNIST
The MNIST data set [18], see Figure 1 (a), contains grayscale 5 images of hand-written digits and the associated labels. The data set is split into 60 000 images for training and 10 000 images for testing. Each image is of size 28 × 28. Let d = 784 be the number of pixels of one image and let the images, reshaped as vectors, be denoted by x (j) ∈ R d and the corresponding labels by y (j) ∈ R d , where d = 10 is the number of different classes. Each label encodes a number in {0, . . . , 9} and the entries y (j) i of the vector y (j) are given by 3 http://yann.lecun.com/exdb/mnist/ 4 http://github.com/zalandoresearch/fashion-mnist 5 The methods described below can be easily extended to color images by defining basis functions for each primary color.  i.e., y (j) = [1, 0, 0, . . . , 0] represents 0, y (j) = [0, 1, 0, . . . , 0] represents 1, etc. This is also called one-hot encoding in machine learning. The fashion MNIST data set [19] can be regarded as a shoo-in replacement for the original data set. There are again 60 000 training and 10 000 test images of size 28 × 28. Some samples are shown in Figure 1 (b) and the corresponding labels in Figure 1 (c). Given a picture of a clothing item, the goal now is to identify the correct category, which is encoded as described above.

SINDy
SINDy [11] was originally developed to learn the governing equations of dynamical systems from data. We will show how it can, in the same way, be used for classification problems. Consider an autonomous ordinary differential equation of the formẋ = f (x), with f : R d → R d . Given m measurements of the state of the system, denoted by x (j) , j = 1, . . . , m, and the corresponding time derivatives y (j) :=ẋ (j) , the goal is to reconstruct the function f from the measurement data. Let X = [x (1) , . . . , x (m) ] ∈ R d×m and Y = [y (1) , . . . , y (m) ] ∈ R d×m . That is, d = d in this case. In order to represent f , we select a vector-valued basis function Ψ : R d → R n and define the transformed data matrix Ψ X ∈ R n×m by Omitting sparsity constraints, SINDy then boils down to solving where is the coefficient matrix. Each column vector ξ i then represents a function f i , i.e., We thus obtain a model of the formẋ = Ξ Ψ(x), which approximates the possibly unknown dynamics. The solution of the minimization problem (3) with minimal Frobenius norm is given by where + denotes the pseudoinverse, see [20].

Tensor-based learning
We will now briefly introduce the basic concepts of tensor decompositions and tensor formats as well as the tensor-based reformulation of SINDy, called MANDy, proposed in [13]. Additionally, recently introduced methods for supervised learning with tensor networks will be discussed.

Tensor decompositions
In order to mitigate the curse of dimensionality when working with tensors T ∈ R n1×···×np , where n µ ∈ N, we will exploit low-rank tensor approximations. The simplest approximation of a tensor of order p is a rank-one tensor, i.e., a tensor product of p vectors given by where T (µ) , µ = 1, . . . , p, are vectors in R nµ . If a tensor is written as the sum of r rank-one tensors, i.e., :,k ⊗ T :,k ⊗ · · · ⊗ T (p) with T (µ) ∈ R nµ×r , this results in the so-called canonical format. In fact, any tensor can be expressed in this format, but we are particularly interested in low-rank representations of tensors in order to reduce the storage consumption as well as the computational costs. The same requirement applies to tensors expressed in the tensor-train format (TT format), where a high-dimensional tensor is represented by a network of multiple low-dimensional tensors [21,22]. A tensor T ∈ R n1×···×np is said to be in the TT format if The tensors T (µ) ∈ R rµ−1×nµ×rµ of order 3 are called TT cores. The numbers r µ are called TT ranks and have a strong influence on the expressivity of a tensor train. It holds that r 0 = r p = 1 and r µ ≥ 1 for µ = 1, . . . , p − 1. Figure 2 (a) shows the graphical representation of a tensor train, which is also called Penrose notation, see [23]. The left-and right-unfoldings of a TT core T (µ) are given by the matrices respectively. Here, the indices of two modes of T (µ) are lumped into a single row or column index while the remaining mode forms the other dimension of the unfolding matrix. We call the TT core T (µ) left-orthonormal if its left-unfolding is orthonormal with respect to the rows, i.e., L µ · L µ = Id ∈ R rµ×rµ . Correspondingly, a core is called right-orthonormal if its right-unfolding is orthonormal with respect to the columns, i.e., R µ · R µ = Id ∈ R rµ−1×rµ−1 . In Penrose notation, orthonormal components are depicted by half-filled circles, cf. Figure 2 (b), where a tensor train with left-orthonormal cores is shown. A given TT core can be left-or right-orthonormalized, respectively, by computing a singular value decomposition (SVD) of its unfolding. For instance, the components of an SVD of the form L µ = U ·Σ·V can be interpreted as a left-orthonormalized version of T (µ) coupled with the matrices Σ and V . When we talk about, e.g., left-orthonormalization of the cores of a tensor train, we mean the application of sequential SVDs from left to right (also called HOSVD, cf. [24]) where U builds the updated core, while the non-orthonormal part Σ · V is contracted with the subsequent TT core. As described in [25,13,16], left-and right-orthonormalization can be used to construct pseudoinverses of tensors. The general idea is to construct a global SVD of a given tensor train by left-and right-orthonormalizing its cores. However, in Section 3.2 we will exploit the structure of transformed data tensors as introduced in [13 Figure 2: Graphical representation of tensor trains: (a) A core is depicted by a circle with different arms indicating the modes of the tensor and the rank indices. The first and the last TT core are regarded as matrices due to the fact that r 0 = r p = 1. (b) Left-orthonormalized tensor train obtained by, e.g., sequential SVDs. Note that the TT ranks may change due to orthonormalization, e.g., when using (reduced/truncated) SVDs.
a different method for the construction of pseudoinverses which significantly reduces the computational effort.
We also represent TT cores as two-dimensional arrays containing vectors as elements. In this notation, a single core of a tensor train T ∈ R n1×···×np is written as Then, the expression T = T (1) ⊗ · · · ⊗ T (p) is used for representing tensor trains T, cf. [26,27,13].

MANDy
MANDy [13] is a tensorized version of SINDy and constructs counterparts of the transformed data matrices (2) directly in the TT format. Two different types of decompositions, namely the coordinateand the function-major decomposition, were introduced in [13]. In [16], the technique for the construction of the transformed data tensors was generalized to arbitrary lists of basis functions. This will be explained in more detail in Section 3.1. Given data matrices X, Y ∈ R d×m and basis functions ψ µ : R d → R nµ , µ = 1, . . . , p, the tensor-based representation of the corresponding transformed data tensors Ψ X ∈ R n1×···×np×m enables us to solve the reformulated minimization problem so that the coefficients are given in form of a tensor train Ξ ∈ R n1×···×np×d , cf. Section 2.2. Instead of identifying the governing equations of dynamical systems from data, see [13], we seek to classify images using MANDy. The only difference is that Ψ X now contains the transformed images and Y the corresponding labels. Since the matrix Y may have different dimensions than X, i.e., Y ∈ R d ×m , the aim is to find the optimal solution of (12) in form of a tensor train Ξ ∈ R n1×···×np×d . We will discuss the explicit representation of transformed data tensors and their pseudoinversion in Section 3.

Supervised learning with tensor networks
It has been shown in [5,6] that tensor-based optimization schemes can be adapted to supervised learning problems. A given input vector x is mapped into a higher-dimensional space using a feature map Ψ before being classified by a decision function f : where Ξ is a coefficient tensor in TT format. The ith entry of the vector f (x) then represents the likelihood that the image x belongs to the class with label i − 1. The transformation defined in [5,6] reads as follows: where α is a parameter. It turns out that the originally proposed choice of α = π 2 is often not optimal. This will be discussed in more detail below. The function Ψ assigns each pixel of the image a two-dimensional vector, inspired by the spin-vectors encountered in quantum mechanics [6]. It was illustrated in [14] how such a transformation can be implemented as a quantum feature map, where the information is encoded in the amplitudes of qubits. Embedding data into quantum Hilbert spaces might be interesting in cases where the quantum device evaluates kernels faster or where kernels cannot be simulated by classical computers anymore [14].
Due to the tensor structure, Ψ(x) is a tensor with 2 d entries, which, for the original MNIST image size amounts to n ≈ 10 236 basis functions. In [5,6], the image size is first reduced to 14 × 14 pixels by averaging groups of four pixels, which then results in "only" n ≈ 10 59 basis functions. Thus, storing the full coefficient matrix is clearly infeasible since Ξ ∈ R 2×···×2×d ∼ = R n×d . Here, d appears as an additional tensor index since the decision function is computed for all d labels simultaneously.
In order to learn the tensor Ξ from training data, a DMRG/ALS-related algorithm (cf. [4,28]) that sweeps back and forth along the cores and iteratively minimizes the cost function is devised. The suggested algorithm varies two neighboring cores at the same time, which allows for adapting the tensor ranks, and computes an update using a gradient descent step. The tensor ranks are reduced by truncated SVDs to control the computational costs. The truncation of the TT ranks can also be interpreted as a form of regularization. For more details, we refer to [5,6].
Different techniques to improve the original algorithm presented in [5] were proposed. In [29], the image data is preprocessed using a discrete cosine transformation and the ordering of the pixels is optimized in order to reduce the ranks. In [10], the DMRG-based sweeping method was replaced by a stochastic gradient descent approach, where the gradient is computed with the aid of automatic differentiation. Furthermore, it was shown that GPUs allow for an efficient solution of such problems.

Tensor-based classification algorithms
We will now describe two different tensor-based classification approaches. First, we show how to combine MANDy with kernel-based regression techniques in order to derive an efficient method for the computation of the pseudoinverse of the transformed data tensor. Then, a classification algorithm based on the alternating optimization of the TT cores of the coefficient tensor is proposed.
For m different vectors stored in a data matrix X = [x (1) , . . . , x (m) ] ∈ R d×m , we want to construct transformed data tensors Ψ X ∈ R n1×···×np×m with (Ψ X ) :,...,:,j = Ψ(x (j) ). In [13,16], this was achieved by multiplying (with the aid of the tensor product) the rank-one decompositions given in (16) for all vectors x (1) , . . . , x (m) by additional unit vectors and subsequently summing them up. The transformed data tensor can then be represented using the following canonical/TT decompositions: where e j , j = 1, . . . , m, denote the unit vectors of the standard basis in the m-dimensional Euclidean space. An entry of Ψ X is given by for 1 ≤ i k ≤ n k and 1 ≤ j ≤ m. Thus, the matrix-based counterpart of Ψ X , see (2), would be given by the mode-p unfolding That is, the modes n 1 , . . . , n p represent row indices of the unfolding, and mode m is the column index. However, for the purpose of this paper, we modify the representation of our transformed data tensors. First, realize that the last core of the TT representation in (17) can be neglected since it is only a reshaped identity matrix. The result is then a tensor network with an "open arm" which can be regarded as a tensor train with an additional column mode located at the last core, see Figure 3 (a). Second, this additional mode can be shifted to any TT core of the decomposition. This is shown in Figure 3 (b). We will benefit from these modifications in Section 3.3 when constructing the subproblems for the ALS-inspired approach. Consider the TT decomposition Ψ X given by Note that this tensor is an element of the tensor space R n1×···×np , i.e., Ψ X has no additional column dimension and it holds that Now, we define Ψ X,µ ∈ R n1×···×np×m to be the tensor derived from Ψ X by replacing the µth core by where the outer modes correspond to the rank dimensions while the inner modes represent the dimensions of the matrices. Analogously, for the first and the last core of Ψ X,µ the non-diagonal core structure  Figure 3: TT representation of transformed data tensors: (a) As in [13], the first p cores (blue circles) are given by (17). The direct contraction of the two last TT cores in (17) can be regarded as an operator-like TT core with a row and column mode (green circle). (b) The additional column mode can be shifted to any of the p TT cores.
has to be used. The 4-dimensional TT core (22) naturally represents a component of a TT operator.
In what follows, we will not need to store the whole TT core given in (22). Otherwise this would mean that we have to save m 3 · n scalar entries (not using a sparse format). However, from a theoretical point of view, Ψ X in Figure 3 (a) and Ψ X,µ in Figure 3 (b) represent the same tensor in R n1×···×np×m , see Appendix A.1.

Kernel-based MANDy
Given a training set X ∈ R d×m , the corresponding label matrix Y ∈ R d ×m , and a set of basis functions ψ µ : R d → R nµ , µ = 1, . . . , p, we exploit the canonical representation of Ψ X given in (17) for kernelbased MANDy. The aim is to solve the optimization problem (12), i.e., we try to find a coefficient tensor Ξ ∈ R n1×···×np×d such that Ξ Ψ X is as close as possible to the corresponding label matrix Y ∈ R d ×m . The solution of (12) with minimal Frobenius norm is given by Ξ = Y Ψ + X , cf. (6). Note that, compared to standard SINDy/MANDy, the matrix Y here does not necessarily have the same dimensions as X. Due to potentially large ranks of the transformed data tensor Ψ X , the direct computation of the pseudoinverse using left-and right-orthonormalization as proposed in [13] would be computationally expensive. However, using the identity Ψ + X = (Ψ X Ψ X ) + Ψ X , we can rewrite the coefficient tensor as The contraction of Ψ X and Ψ X yields a Gram matrix G ∈ R m×m whose entries are given by the resulting kernel function k(x, x ) = Ψ(x), Ψ(x ) , i.e., Note that due to the tensor structure of Ψ X , we obtain i.e., a product of p local kernels. (14), this can be simplified to

Remark 1: For the basis functions defined in
which is a product of cosine kernels, cf. [5].
The product structure of the kernel allows us to compute the Gram matrix G as a Hadamard product (denoted by ) of p matrices, that is, where Θ µ ∈ R m×m is given by We now define Z := Y G + ∈ R d ×m , which can be obtained by solving the system Z G = Y (in the least-squares sense if G is singular). The decision function f , cf. (13), is then given by and again only requires kernel evaluations. As above, we can use a sequence of Hadamard products to compute Ψ X Ψ(x). The classification problem can thus be solved as summarized in Algorithm 1.

Input:
Training set X and label matrix Y , test setX, basis functions. Output: Label matrixỸ .
1: Compute G using (27) and (28). 2: Solve Z G = Y . 3: Define the decision function f using (29). 4: Apply f to every vectorx in the test set, and store the resulting vectorsỹ in the matrixỸ . 5: The index i of the largest entry ofỹ determines the detected label, i.e., setỹ = e i .
We could also replace the pseudoinverse G + by the regularized inverse (G + ε Id) −1 , where ε is the regularization parameter, which would lead to a slightly different system of linear equations. However, for the numerical experiments in Section 4, we do not use regularization. Algorithm 1 is equivalent to kernel ridge regression (see, e.g., [15]) with a tensor-product kernel. This is not surprising since we are solving simple least-squares problems.

Remark 2:
Note that the kernel does not necessarily have to be based on tensor products of basis functions for this method to work, we could also simply use, e.g., a Gaussian kernel, which for the MNIST data set leads to slightly lower but similar classification rates. Tensor-based kernels, however, have an exponentially large yet explicit feature space representation and additional structure that could be exploited to speed up computations. Moreover, the kernel-based algorithm outlined above can in the same way be applied to time-series data to learn governing equations in potentially infinite-dimensional feature spaces.
Compared to the method proposed in [5,6], the advantage of our approach, which can be regarded as a kernel-based formulation of MANDy (or SINDy), is that we can compute a closed-form solution without necessitating any iterations or sweeps. However, even though this approach for classification problems computes an optimal solution of the minimization problem (12), the runtime as well as the memory consumption of the algorithm depend crucially on the size of the training data set (and also the number of labels) and the resulting coefficient tensor Ξ has no guaranteed low-rank structure. We will now propose an alternating optimization method which circumvents this problem.

Alternating ridge regression (ARR)
In what follows, we will use the TT representation illustrated in Figure 3 (b) for the transformed data tensor Ψ X ∈ R n1×···×np×m . Even though we do not consider a TT operator, the proposed approach is closely related to the DMRG method [4], also called alternating linear scheme (ALS) [28]. As in [5,6], the idea here is to compute a low-rank TT approximation of the coefficient tensor Ξ by an alternating scheme. That is, a low-dimensional system of linear equations has to be solved for each TT core. Our approach is outlined in Algorithm 2. ].

13:
Determine truncated SVD solution of min v w − vM µ 2 . 14: if µ > 1 then 15: Apply QR decomposition to extract updated core. 16: else 17: Set the updated core to a reshape of v.

18:
Repeat lines 5-17 to increase accuracy (if needed). 19: Define Ξ using (31) and set y = f (x) using (13). 20: The index of the largest entry of y determines the detected label, see Algorithm 1.
First, note that instead of solving the minimization problem (12) we can also find separate solutions of for each row of Y . Since these systems can be solved independently, Algorithm 2 can be easily parallelized. We then use a DMRG/ALS-inspired scheme to split the optimization problem (30) into p subproblems. The micro-matrix M µ of such a subproblem can be built from three different parts, namely Ψ (µ) X,µ , P µ , and Q µ . The latter are both collected in a left and right stack in order to avoid repetitive computations. Note that P µ is determined by contracting P µ−1 with the (µ − 1)th cores of Ξ i and Ψ X . Analogously, Q µ is build from Q µ+1 and the (µ + 1)th cores of Ξ i and Ψ X . During the first half sweep of Algorithm 2, we only have to compute the matrices P µ since the used matrices Q µ are not based on any updated cores. Afterwards, the matrices Q µ are (re-)computed during the second half. See [28] for further details and Figure 4 for a graphical illustration of the construction of the subproblems and the extraction of the optimized core. Note that it is not necessary to store the (sparse) core Ψ (µ) X,µ in its full representation as a 4-dimensional array in order to construct the matrix M µ . By using, e.g., NumPy's einsum the TT core can be replaced a (dense) matrix containing the corresponding function evaluations.
By orthonormalizing the fixed cores of Ξ and using truncated SVDs [17] for solving the subsystems we can interpret our approach as a core-wise ridge regression approximating the solution obtained by kernel-based MANDy, see Appendix A.2. After approximating the coefficient tensor the decision function f is given by (13). The main difference between our approach and the method introduced in [5, 6] is that we do not update the TT cores of Ξ using gradient descent steps. Instead we solve a low-dimensional system of linear equations corresponding to the entire training data set whose solution yields the updated core. Moreover, we solve a minimization problem for each row of The TT core (red circle) obtained by solving the low-dimensional minimization problem is decomposed (e.g., using QR factorization) into a orthonormal tensor and a triangular matrix. The orthonormal tensor then yields the updated core. the label matrix Y . Using the modified basis decomposition introduced in Section 3.1, it is possible to significantly reduce the storage consumption of the stack, see Algorithm 2 Lines 4 and 11. If we only use a fixed representation for Ψ X as given in (17), the additional mode would lead to a much higher storage consumption of the right stack. Thus, our method provides an efficient construction of the subproblems.

Numerical results
We apply the tensor-based classification algorithms described in Sections 3.2 and 3.3 to both the MNIST and fashion MNIST data set, choosing the basis defined in (14) and setting α ≈ 0.59. This value was determined empirically for the MNIST data set, but also leads to better classification rates for the fashion MNIST set. Kernel-based MANDy as well as ARR are available in Scikit-TT. 6 The numerical experiments were performed on a Linux machine with 128 GB RAM and an Intel Xeon processor with a clock speed of 3 GHz and 8 cores.
For the first approach, using kernel-based MANDy, we do not apply any regularization techniques. For the ARR approach, we set the TT rank for each solution Ξ i , see Algorithm 2, to 10 and repeat the scheme 5 times. Here, we use regularization, i.e., truncated SVDs with a relative threshold of 10 −2 are applied to the minimization problems given in Algorithm 2 (Lines 8 and 13). The obtained classification rates for the reduced and full MNIST and fashion MNIST data are shown in Figure 5.
Similarly to [5,6], we first apply the classifiers to the reduced data sets, see Figure 5 (a). Using MANDy, we obtain classification rates of up to 98.75 % for the MNIST and 88.82 % for the fashion MNIST data set. Using the ARR approach, the classification rates are not monotonically increasing, which simply may be an effect of the alternating optimization scheme. The highest classification rates we obtain are 98.16 % for the MNIST data and 87.55 % for the fashion MNIST data. We typically obtain a 100 % classification rate for the training data (as a consequence of the richness of the feature space). This is not necessarily a desired property since the learned model might not generalize well to new data, but seems to have no detrimental effects for the simple MNIST classification problem. As shown in Figure 5 (b), kernel-based MANDy can still be applied when considering the full data sets without reducing the image size. Here, we obtain classification rates of up to 97.24 % for the MNIST and 88.37 % for the fashion MNIST data set. That we obtain lower classification rates for the full images as compared to the reduced ones might be due to the fact that pixel-by-pixel comparisons of images are not expedient. The averaging effect caused by downscaling the images helps to detect coarser features. This is similar to the effect of convolutional kernels and pooling layers. In principle, ARR can also be used for the classification of the full data sets. So far, however, our numerical experiments produced only classification rates significantly lower than those obtained by applying MANDy (95.94 % for the MNIST and 82.18 % for fashion MNIST data set). This might be due to convergence issues caused by the kernel. The application to higher-order transformed data tensors and potential improvements of ARR will be part of our future research. Figure 5 also shows a comparison with tensorflow. We run the code provided as a classification tutorial 7 ten times and compute the average classification rate. The input layer of the network comprises 784 nodes (one for each pixel; for the reduced data sets, we thus have only 196 input nodes), followed by two dense layers with 128 and 10 nodes. The layer with 10 nodes is the output layer containing probabilities that a given image belongs to the class represented by the respective neuron. Note that although more sophisticated methods and architectures for these problems exists-see the (fashion) MNIST website for a ranking-, the results show that our tensor-based approaches are competitive with state-of-the-art deep-learning techniques.
In order to understand the numerical results for the MNIST data set (obtained by applying kernelbased MANDy to all 60 000 training images), we analyze the misclassified images, examples of which are displayed in Figure 6 (a). For misclassified images x, the entries of f (x), see (29), are often numerically zero, which implies that there is no other image in the training set that is similar enough so that the kernel can pick up the resemblance. Some of the remaining misclassified digits are hard to recognize even for humans. Histograms demonstrating which categories are misclassified most often are shown in Figure 6 (b). Here, we simply count the instances where an image with label i was assigned the wrong label j. The digits 2 and 7 as well as 4 and 9 are confused most frequently. Additionally, we wish to visualize what the algorithm detects in the images. To this end, we perform a sensitivity analysis as follows: Starting with an image whose pixel values are constant everywhere (zero or any other value smaller than one, we choose 0.5), we set pixel (i, j) to one and compute y = f (x) for this image. The process is repeated for all pixels. For each label, we then plot a heat map of the values of y. This tells us which pixels contribute most to the classification of the images. The resulting maps are shown in Figure 6 (c). Except for the digit 1, the results are highly similar to the images obtained by averaging over all images containing a certain digit. Figure 7 shows examples of misclassified images and the corresponding histogram as well as the results of the sensitivity analysis for the fashion MNIST data set. We see that the images of shirts (6) are most difficult to classify (due to the ambiguity in the category definitions), whereas trousers (1) and bags (8) have the lowest misclassification rates (probably due to their distinctive shapes). In contrast to the MNIST data set, the results of the sensitivity analysis differ widely from the average images. The classifier for coats (4), for instance, "looks for" a zipper and coat pockets, which are not visible in the average coat, and the classifier for dresses (3) seems to base the decision on the presence of creases, which are also not distinguishable in the average dress. The interpretation of other classifiers is less clear, e.g., the ones for sandals (5) and sneakers (7) seem to be contaminated by other classes.
Comparing the runtimes of both approaches applied to the reduced data sets with 60 000 training images, kernel-based MANDy needs approximately one hour for the construction of the decision function (29). On the other hand, ARR needs less than 10 minutes to compute the coefficient tensor assuming we parallelize Algorithm 2.

Conclusion
In this work, we presented two different tensor-based approaches for supervised learning. We showed that a kernel-based extension of MANDy can be utilized for image classification. That is, extending the method to arbitrary least-squares problems (originally, MANDy was developed to learn governing equations of dynamical systems) and using sequences of Hadamard products for the computation of the pseudoinverse, we were able to demonstrate the potential of kernel-based MANDy by applying it to the MNIST and fashion MNIST data sets. Additionally, we proposed the alternating optimization scheme ARR, which approximates the coefficient tensors by low-rank TT decompositions. Here, we used a mutable tensor representation of the transformed data tensors in order to construct low-dimensional regression problems for optimizing the TT cores of the coefficient tensor.
Both approaches use an exponentially large set of basis functions in combination with least-squares regression techniques on a given set of training images. The results are encouraging and show that methods exploiting tensor products of simple basis functions are able to detect characteristic features in image data. The work presented in this paper constitutes a further step towards tensor-based techniques for machine learning.
The reason why we can handle the extremely high-dimensional feature space spanned by basis functions is its tensor-product format. Besides the general questions of the choice of basis functions and the expressivity of these functions, the rank-one tensor products that were used in this work can in principle be replaced by other structures which might result in higher classification rates. For instance, the transformation of an image could be given by a TT representation with higher ranks or hierarchical tensor decompositions (with the aim to detect features on different levels of abstraction). Furthermore, we could define different basis functions for each pixel, vary the number of basis functions per pixel, or define basis functions for groups of pixels.
Even though kernel-based MANDy computes the minimum norm solution of the considered regression problems as an exact TT decomposition, the method is likely to suffer from high ranks of the transformed data tensors and might thus not be competitive for large data sets. At the moment, we are computing the Gram matrix for the entire training data set. However, a possibility to speed up computations and to lower the memory consumption would be to exploit properties of the kernel. That is, if the kernel almost vanishes if two images differ significantly in at least one pixel (as it is the case for the specific kernel used in this work, provided that the originally proposed value α = π 2 is used), the Gram matrix is essentially sparse when setting entries smaller than a given threshold to zero. Using sparse solvers would allow us to handle much larger data sets. Moreover, the construction of the Gram matrix is highly parallelizable and it would be possible to use GPUs to assemble it in a more efficient fashion.
Further modifications of ARR such as different regression methods for the subproblems, an optimized ordering of the TT cores, and specific initial coefficient tensors can help to improve the results. We provided an explanation for the stability of ARR, but the properties of alternating regression schemes have to be analyzed in more detail in the future.