Tensor-based algorithms for image classification

Tensor-based algorithms for image

Stefan Klus    Patrick Gelß stefan.klus@fu-berlin.de

The interest in machine learning with tensor networks has been growing rapidly in recent years. The goal is to exploit tensor-structured basis functions in order to generate exponentially large feature spaces which are then used for supervised learning. We will propose two different tensor approaches for quantum-inspired machine learning. One is a kernel-based reformulation of the previously introduced MANDy, the other an alternating ridge regression in the tensor-train format. We will apply both methods to the MNIST and fashion MNIST data set and compare the results with state-of-the-art neural network-based classifiers.

quantum machine learning, image classification, tensor-train format, kernel-based methods, ridge regression

15A69, 62J07, 65D18, 68Q32


Department of
Mathematics and
Computer Science,
Freie Universität Berlin,
Berlin 14195, Germany Department of
Mathematics and
Computer Science,
Freie Universität Berlin,
Berlin 14195, Germany

1 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 learning111There are different research directions in the field of quantum machine learning, we here understand it as using quantum computing capabilities for machine learning problems. 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 TensorNetwork222http://github.com/google/TensorNetwork [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.333http://yann.lecun.com/exdb/mnist/,444http://github.com/zalandoresearch/fashion-mnist 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:

  • [leftmargin=*]

  • 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.

2 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.

Symbol Description
data matrix in
label matrix in
mode dimensions of tensors
ranks of tensor trains
basis functions
/ transformed data matrices/tensors
/ coefficient matrices/tensors
Table 1: Notation used in this work.

2.1 MNIST and fashion MNIST

The MNIST data set [18], see Figure 1 (a), contains grayscale555The methods described below can be easily extended to color images by defining basis functions for each primary color. images of hand-written digits and the associated labels. Let be the number of pixels of one image and let the images, reshaped as vectors, be denoted by and the corresponding labels by , where is the number of different classes. Each label encodes a number in and the entries of the vector are given by

i.e., represents , represents , etc. This is also called one-hot encoding in machine learning. The data set is split into images for training and images for testing.

0 T-shirt/top
1 Trousers
2 Pullover
3 Dress
4 Coat
5 Sandal
6 Shirt
7 Sneaker
8 Bag
9 Ankle boot
Figure 1: (a) Samples of the MNIST data set. (b) Samples of the fashion MNIST data set. Each row represents a different item type. (c) Corresponding labels for the fashion MNIST data set.

The fashion MNIST data set [19] can be regarded as a shoo-in replacement for the original data set. There are again training and test images of size . 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.

2.2 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 , with . Given measurements of the state of the system, denoted by , , and the corresponding time derivatives , the goal is to reconstruct the function from the measurement data. Let and . That is, in this case. In order to represent , we select a vector-valued basis function and define the transformed data matrix by


Omitting sparsity constraints, SINDy then boils down to solving



is the coefficient matrix. Each column vector then represents a function , i.e.,

We thus obtain a model of the form , which approximates the possibly unknown dynamics. The solution of the minimization problem (2) with minimal Frobenius norm is given by


where denotes the pseudoinverse, see [20].

2.3 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.

2.3.1 Tensor decompositions

In order to mitigate the curse of dimensionality when working with tensors , where , we will exploit low-rank tensor approximations. The simplest approximation of a tensor of order is a rank-one tensor, i.e., a tensor product of vectors given by

where , , are vectors in . If a tensor is written as the sum of rank-one tensors, i.e.,

with , 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 is said to be in the TT format if

The tensors of order 3 are called TT cores. The numbers are called TT ranks and have a strong influence on the expressivity of a tensor train. It holds that and for . 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 are given by the matrices

respectively. Here, the indices of two modes of 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 left-orthonormal if its left-unfolding is orthonormal with respect to the rows, i.e., . Correspondingly, a core is called right-orthonormal if its right-unfolding is orthonormal with respect to the columns, i.e., . 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.


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 . (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 given TT core can be left- or right-orthonormalized, respectively, by computing an SVD of its unfolding. For instance, the components of an SVD of the form can be interpreted as a left-orthonormalized version of coupled with the matrices and . 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 builds the updated core, while the non-orthonormal part 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] to propose 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 is written as


Then, the expression is used for representing tensor trains , cf. [26, 27, 13].

2.3.2 MANDy

MANDy [13] is a tensorized version of SINDy and constructs counterparts of the transformed data matrices (1) directly in the TT format. Two different types of decompositions, namely the coordinate- and 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 and basis functions , , the tensor-based representation of the corresponding transformed data tensors enables us to solve the reformulated minimization problem


so that the coefficients are given in form of a tensor train , 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 now contains the transformed images and the corresponding labels. Since the matrix may have different dimensions than , i.e., , the aim is to find the optimal solution of (5) in form of a tensor train . We will discuss the explicit representation of transformed data tensors and their pseudoinversion in Section 3.

2.3.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 is mapped into a higher-dimensional space using a feature map before being classified by a decision function of the form


where is a coefficient tensor in TT format. The th entry of the vector then represents the likelihood that the image belongs to the class with label . The transformation defined in [5, 6] reads as follows:


where is a parameter. It turns out that the originally proposed choice of 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, is a tensor with entries, which, for the original MNIST image size amounts to basis functions. In [5, 6], the image size is first reduced to pixels by averaging groups of four pixels, which then results in “only” basis functions. Thus, storing the full coefficient matrix is clearly infeasible since . Here, appears as an additional tensor index since the decision function is computed for all 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.

3 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.

3.1 Basis decomposition

As above, let be a vector and , , basis functions. We consider the rank-one tensors


For different vectors stored in a data matrix , we want to construct transformed data tensors with . In [13, 16], this was achieved by multiplying (with the aid of the tensor product) the rank-one decompositions given in (9) for all vectors by additional unit vectors and subsequently summing them up. The transformed data tensor can then be represented using the following canonical/TT decompositions:


where , , denote the unit vectors of the standard basis in the -dimensional Euclidean space. An entry of is given by

for and . Thus, the matrix-based counterpart of , see (1), would be given by the mode- unfolding

That is, the modes represent row indices of the unfolding, and mode 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 (10) 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 given by

Note that this tensor is an element of the tensor space , i.e., has no additional column dimension and it holds that


Figure 3: TT representation of transformed data tensors: (a) As in [13], the first cores (blue circles) are given by (10). The direct contraction of the two last TT cores in (10) 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 TT cores.

Now, we define to be the tensor derived from 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 the non-diagonal core structure has to be used. The -dimensional TT core (11) naturally represents a component of a TT operator. In what follows, we will not need to store the whole TT core given in (11). Otherwise this would mean that we have to save scalar entries (not using a sparse format). However, from a theoretical point of view, in Figure 3 (a) and in Figure 3 (b) represent the same tensor in , see Appendix A.1.

3.2 Kernel-based MANDy

Given a training set , the corresponding label matrix , and a set of basis functions , , we exploit the canonical representation of given in (10) for kernel-based MANDy. The aim is to solve the optimization problem (5), i.e., we try to find a coefficient tensor such that is as close as possible to the corresponding label matrix . The solution of (5) with minimal Frobenius norm is given by , cf. (3). Note that, compared to standard SINDy/MANDy, the matrix here does not necessarily have the same dimensions as . Due to potentially large ranks of the transformed data tensor , the direct computation of the pseudoinverse using left- and right-orthonormalization as proposed in [13] would be computationally expensive. However, using the identity , we can rewrite the coefficient tensor as

The contraction of and yields a Gram matrix whose entries are given by the resulting kernel function , i.e.,

Note that due to the tensor structure of , we obtain

i.e., a product of local kernels. {remark} For the basis functions defined in (7), this can be simplified to

which is a product of cosine kernels, cf. [5]. The product structure of the kernel allows us to compute the Gram matrix as a Hadamard product (denoted by ) of matrices, that is,


where is given by


We now define , which can be obtained by solving the system (in the least-squares sense if is singular). The decision function , cf. (6), is then given by


and again only requires kernel evaluations. As above, we can use a sequence of Hadamard products to compute . The classification problem can thus be solved as summarized in Algorithm 1.



1:Compute using (12) and (13).
2:Solve .
3:Set using (14).
4:The index of the largest entry of determines the detected label.
Algorithm 1 Kernel-based MANDy for classification.

We could also replace the pseudoinverse by the regularized inverse , 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.


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 (5), 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.

3.3 Alternating ridge regression (ARR)

In what follows, we will use the TT representation illustrated in Figure 3 (b) for the transformed data tensor . 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.



1:for  do(parallelizable)
2:     Define .
3:     Define initial guess and right-orthonormalize.
4:     Compute right stack .
5:     for  do(first half sweep)
6:         Compute .
7:         Construct micro-matrix from .
8:         Determine truncated SVD solution of .
9:         Apply QR decomposition to extract updated core.      
10:     for  do(second half sweep)
11:         Compute .
12:         Construct micro-matrix from .
13:         Determine truncated SVD solution of .
14:         if  then
15:              Apply QR decomposition to extract updated core.
16:         else
17:              Set the updated core to a reshape of .               
18:     Repeat lines 517 to increase accuracy (if needed).
19:Define using (16) and set using (6).
20:The index of the largest entry of determines the detected label.
Algorithm 2 Alternating ridge regression (ARR) for classification.

First, note that instead of solving the minimization problem (5) we can also find separate solutions of


for each row of . 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 (15) into subproblems. The micro-matrix of such a subproblem can be built from three different parts, namely , , and . The latter are both collected in a left and right stack in order to avoid repetitive computations. Note that is determined by contracting with the th cores of and . Analogously, is build from and the th cores of and . During the first half sweep of Algorithm 2, we only have to compute the matrices since the used matrices are not based on any updated cores. Afterwards, the matrices 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 in its full representation as a -dimensional array in order to construct the matrix . By using, e.g., NumPy’s einsum the TT core can be replaced a (dense) matrix containing the corresponding function evaluations.


Figure 4: Construction and solution of the subproblem for the th core: (a) The -dimensional core of (green circle) is contracted with the matrices and constructed by joining the fixed cores of the coefficient tensor (orange circles) with the corresponding cores of the transformed data tensor. The matricization then defines the matrix . (b) 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.

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 is given by (6). 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 label matrix . 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 as given in (10), 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.

4 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 (7) and setting . 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.666https://github.com/PGelss/scikit_tt The numerical experiments were performed on a Linux machine with GB RAM and an Intel Xeon processor with a clock speed of GHz and 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 , see Algorithm 2, to and repeat the scheme times. Here, we use regularization, i.e., truncated SVDs with a relative threshold of 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.

Figure 5: Results for MNIST and fashion MNIST: (a) Classification rates for the reduced 14x14 images. (b) Classification rates for full 28x28 images. Reducing the image size by averaging over groups of pixels improves the performance of the algorithm.

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 for the MNIST and 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 for the MNIST data and for the fashion MNIST data. We typically obtain a 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 for the MNIST and 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 ( for the MNIST and 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 tutorial777www.tensorflow.org/tutorials/keras/basic_classification 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.

Figure 6: MNIST classification: (a) Images misclassified by kernel-based MANDy described in Section 3.2. The original image is shown in black, the identified label in red, and the correct label in green. (b) Histograms illustrating which categories are misclassified most often. The rows represent the correct labels of the misclassified image and the columns the detected labels. (c) Visualizations of the learned classifiers showing a heat map of the classification function obtained by applying it to images that differ in one pixel.

In order to understand the numerical results for the MNIST data set (obtained by applying kernel-based MANDy to all training images), we analyze the misclassified images, examples of which are displayed in Figure 6 (a). For misclassified images , the entries of , see (14), 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). The digits and as well as and 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 ), we set pixel to one and compute for this image. The process is repeated for all pixels. For each label, we then plot a heat map of the values of . 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 , the results are highly similar to the images obtained by averaging over all images containing a certain digit.

Figure 7: Fashion MNIST classification: (a) Misclassified images. (b) Histogram of misclassified images. (c) Visualizations of the learned classifiers.

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 training images, kernel-based MANDy needs approximately one hour for the construction of the decision function (14). On the other hand, ARR needs less than minutes to compute the coefficient tensor assuming we parallelize Algorithm 2.

5 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 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.


This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 “Scaling Cascades in Complex Systems”. Part of this research was performed while S.K. was visiting the Institute for Pure and Applied Mathematics (IPAM), which is supported by the National Science Foundation (Grant No. DMS-1440415). We would like to thank Michael Götte and Alex Goeßmann from the TU Berlin for interesting discussions related to tensor decompositions and system identification.


Appendix A Appendix

a.1 Representation of transformed data tensors


For all , it holds that

That is, the TT decompositions and represent the same tensor in .


An entry of , , is given by

By definition,

On the other hand, an entry of with and is nonzero if and only if . It follows that

This can be shown in an analogous fashion for and . ∎

a.2 Interpretation of ARR as ALS ridge regression

The following reasoning will elucidate the relation between ARR, ridge regression, and kernel-based MANDy. We only outline the rough idea without concrete proofs. Let denote the retraction operator, see [28], consisting of the fixed TT cores and of the solution at any iteration step of Algorithm 2. Furthermore, assume that are left- and right-orthonormal. In Lines 8 and 13 of Algorithm 2, we consider the system (with a slight abuse of notation)

The application of a truncated SVD to the matricization of (as done in Algorithm 2) is then similar to a regularization in the form of


with appropriate regularization parameter , i.e., for both approaches, see [17, 30]. The formulation (17) is known as Tikhonov’s smoothing functional, ridge regression, or regularization (which, of course, could also directly be applied in Algorithm 2). The solution of (17) is also the solution of the regularized normal equation

see, e.g., [31]. Since , it follows that

In fact, this is a subproblem corresponding to the application of ALS [28] to the tensor-based system


Note that all requirements for the application of ALS are satisfied since is a symmetric positive definite tensor operator and is orthonormal. The system of linear equations given in (18) is then equivalent to the minimization problem

For sufficiently small , it holds that , see [32], meaning Algorithm 2 computes an approximation of the coefficient tensor resulting from the application of kernel-based MANDy, see Section 3.2.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description