{'fraction_non_alphanumeric': 0.02600973788440677, 'fraction_numerical': 0.012210197028179319, 'mean_word_length': 4.547935619314206, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Relevant Features of Data with Multi-scale Tensor Networks

Inspired by coarse-graining approaches used in physics, we show how similar algorithms can be adapted for data. The resulting algorithms are based on layered tree tensor networks and scale linearly with both the dimension of the input and the training set size. Computing most of the layers with an unsupervised algorithm, then optimizing just the top layer for supervised classification of the MNIST and fashion-MNIST data sets gives very good results. We also discuss mixing a prior guess for supervised weights together with an unsupervised representation of the data, yielding a smaller number of features nevertheless able to give good performance.

A particularly interesting class of tensor decompositions are tensor networks. These are factorizations of a very high-order tensor into a contracted network of low-order tensors-see Figs. 1,2. Tensor networks break the curse of dimensionality by allowing operations such as contracting very high-order tensors or retrieving their components to be accomplished with polynomial cost by manipulating the low-order factor tensors [12].
A key motivation for the development of tensor networks was the idea of coarse graining. Known in physics as the renormalization group, coarse graining is a way to gain insight into a complicated statistical system by marginalizing over its smallest length scales while preserving properties at larger scales [13,14]. After each round of coarse graining, one identifies a new smallest scale and repeats the process in a hierarchical fashion. A related idea is a wavelet multiresolution analysis in applied mathematics [15]. For the case of deep neural networks, there is evidence that certain networks implement coarse graining, with each layer of neurons learning progressively coarser features [16]. Such an analogy between the renormalization group and neural nets has been made precise for the case of deep belief networks [17]. Renormalization group ideas could yield insights into data when applied in the framework of PCA [18].
Numerical implementations of the renormalization group motivated the development of tensor networks [19][20][21][22][23] which not only underpin powerful algorithms for solving quantum [24][25][26][27][28] and classical [29,30] systems in physics, but lead to the insight that the solution has a structure imparted by the coarse-graining procedure itself. This structure is most striking in the case of the MERA family of tensor networks Fig. 2(d) [23,31], where the quantum wavefunction acquires an emergent extraf`( 1. The models f (x) we will construct are defined by a tree tensor network that coarse grains the feature vector Φ(x) via the isometric tree tensor layer U. In Section III, the layers U are determined using a purely unsupervised algorithm generalizable to various tasks. In Section V the tensors in U are instead computed by a mixture of the unsupervised approach and an initial guess of a lower-accuracy solution for a supervised task. After the coarse graining, the resulting tensors are contracted with the top tensor w to compute f (x) (in Sec. VI the top tensor is generalized to be a second type of tensor network). Only the top tensor is optimized for the specific task, such as supervised learning. In the multi-task setting, the top tensor has an extra uncontracted index running over labels as shown in the figure.
Recently a framework has been discussed by several groups [4][5][6] that offers a particularly straightforward way to apply tensor network algorithms to machine learning tasks, such as supervised [5,6,[36][37][38] and unsupervised [39] learning. This framework can be viewed equivalently as a neural network architecture with linear activation and product pooling [4,10,36,40] or an approach to kernel learning with the weights represented as a tensor network [5,6]. We will take the kernel learning perspective here, but nevertheless construct models resembling deep neural networks.
When parameterizing a kernel learning model with a tensor network, all operations other than the initial feature map are linear when viewed as transformations on the entire feature space, which is exponentially higher dimensional than the input space. The fact that tensor networks are composed from linear maps is what makes them amenable to theoretical analysis [32,41,42] and useful for devising algorithms [26,43,44]. In the context of parameterizing distributions, many interesting tensor networks are tractable and can be directly sampled [39,45]. Sums and products of tensor networks can be controllably approximated as a tensor network of the same type [26,46]. Tensor network representations of distributions can be proven to have exponentially or power-law decaying correlations depending on their geometry [47]. Algorithms to optimize tensor networks are often adaptive, allowing dimensions of internal indices to adjust as needed [28].
In what follows we will use tensor network coarse graining to compress data originally represented in a very high dimensional space. The approach is unsupervised, based purely on statistical properties of the data. Each step has a controlled accuracy, and the method is scalable to large data set sizes and input dimensions. The resulting reduced description takes the form of a layered tensor network with a tree structure-a tree tensor network Fig. 2(c)-and can be used to obtain good results on learning tasks by optimizing tensors only at the top layer. One could further specialize all of the tensors in the network for a specific task, but we will not do so here.
The cost of training each of the models discussed below is linear in both training set size and input dimension, assuming a fixed number of parameters. The cost of evaluating the model on a test input is independent of training set size. The scaling of the optimization algorithm with training set size could be reduced to sub-linear with stochastic optimization techniques.
Throughout we will use tensor diagram notation. For a brief introduction to this notation, see the Appendix. The experiments were implemented using the ITensor software [48].

II. MOTIVATION AND BACKGROUND
The algorithm we will develop is motivated by the fact that within kernel learning, the optimal weights belong to the span of the training data within feature space. More specifically, consider a model defined in terms of a high-dimensional feature map Φ and weights W . Given a set of training inputs {x j } N T j=1 , it can be shown that for a broad set of learning tasks the optimal weights have the form where only the α j parameters remain to be optimized for the specific task. The well-known fact that W can be expressed this way is called the representer theorem, which applies to many common supervised tasks, as well as certain unsupervised tasks such as kernel PCA [49].

A. Alternate Basis for Weight Parameters
Parameterizing the weights W by a set of numbers α j of the size of the training set is a major improvement compared to representing W within the possibly infinite dimensional feature space defined by Φ. But standard approaches to optimizing W in terms of the α j typically exhibit quadratic or worse dependence on training set size, which can be prohibitive for state-of-the-art tasks with training sets reaching millions in size.
To make progress, notice that the main content of Eq.
(2) is that W resides in the span of the {Φ † (x j )}. Thus the optimal weights W can be expanded in any basis of vectors U † n W = n β n U † n (3) as long as the U † n span the same space as the Φ † (x j ). If cost was no concern, a natural way to obtain a basis U would be to think of the training set feature vectors as a matrix Φ s j = Φ s (x j ) and obtain U from singular value decomposition of Φ s j as where S is the diagonal matrix of singular values s n = S n n . By inserting this decomposition into Eq.
(2), the optimal weights W can indeed be written which explicitly relates the α j and β n parameters. One advantage of expressing W in terms of U is that the columns of U are orthonormal. But more importantly, if the dimension of feature space is much larger than the minimum statistically significant training set size, then many singular values s n will be very small or zero and the corresponding rows of U † can be discarded. Following such a truncation, Eq. (5) says that to a good approximation, the optimal weights can be parameterized within a significantly reduced space of parameters β n and U † is the transformation from the entire feature space to the reduced parameter space.
Computing the singular value decomposition of Φ s j directly would not scale well for large training sets or highdimensional feature maps, yet as we will show it is nevertheless possible to efficiently determine the transformation U in truncated form. Observe that U diagonalizes the feature space covariance matrix [50] defined as where P n = (S n n ) 2 are the eigenvalues of the Hermitian matrix ρ. As we demonstrate in Sec. III below, the feature space covariance matrix ρ is amenable to decomposition as a layered tensor network. Computing every layer of this network can provide an efficient expression for the elements of the basis U corresponding to the largest eigenvalues of ρ. Computing only some of the layers still has the beneficial effect of projecting out directions in feature space along which ρ has small or zero eigenvalues. By carrying out an iterative procedure to truncate directions in feature space along which ρ has a very small projection, one can rapidly reduce the size of the space needed to carry out learning tasks.
We will also see that ρ is not the only choice of matrix for determining a tensor network basis for features. As demonstrated in Sec. V, other choices result in a network more adapted for a specific task, and can have fewer latent parameters without reducing model performance.

B. Tensor Product Feature Maps
Before describing the algorithm to partially or fully diagonalize ρ as a tensor network, we briefly review the class of feature maps which lead to a natural representation of model parameters as a tensor network, as discussed in Refs. 4-6. These are feature maps Φ(x) which map inputs x from a space of dimension N into a space of dimension d N with a tensor product structure. The simplest case of such a map begins by defining a local feature map φ sj (x j ) where s j = 1, 2, . . . , d. These local feature maps define the full feature map as: as shown in Fig. 3(a), where placement of tensors next to each other implies an outer product. This choice of feature map leads to models of the form which are depicted in Fig. 3(b). As evident from the above expression, the weight parameters are indexed by N indices of dimension d. Thus there are d N weight parameters and W is a tensor of order N . We will be interested in the case where d is small (of order one or ten) and N is many hundreds or thousands in size. Of course, manipulating or even storing d N parameters quickly becomes impossible as N increases. A solution that is both practical and interesting is to assume that the optimal weights W can be efficiently approximated by a tensor network [12,28], an idea proposed recently by several groups [4][5][6]36].
A tensor network is a factorization of an order N tensor into the contracted product of low-order tensors. Key examples of well-understood tensor networks for which efficient algorithms are known are depicted in Fig. 2 and include: • the matrix product state (MPS) [19][20][21]51] or tensor train decomposition [52], Fig. 2
Each of these networks makes various tradeoffs in terms of how complicated they are to manipulate versus their ability to represent statistical systems with higherdimensional interactions or more slowly decaying correlations. A good introduction to tensor networks in the physics context is given by Orús in Ref. 12 and in a mathematics context by Cichocki in Ref. 56. Other detailed reviews include Refs. 28, 32, 57, and 58.

III. UNSUPERVISED COARSE GRAINING
As discussed in the previous section, if one can compute the eigenvectors U s n of the feature space covariance matrix ρ, defined as then the optimal weights for a wide variety of learning tasks can be represented as Furthermore, for a specific task the eigenvectors of ρ with small enough eigenvalues can be discarded without reducing performance. As an example, for the task of supervised learning with a quadratic cost, if the cost function includes a quadratic weight penalty λ|W | 2 , eigenvectors whose eigenvalues are much smaller than λ do not contribute significantly to W and can be projected out. Now we will outline a strategy to find a controlled approximation for the dominant eigenvectors of ρ using a layered tensor network. For the purposes of this section, it will be useful to think of as a tensor of order (N + 1) as shown in Fig. 4. In this view, ρ is formed by contracting Φ and Φ † over the training data index j as shown in Fig. 5.
Because it is not feasible to diagonalize ρ directly, the strategy we will pursue is to compute local isometries which combine two indices into one and project out subspaces of the feature space spanned by eigenvectors of ρ with small eigenvalues, as defined by some pre-defined cutoff or threshold . The term isometry here refers to a third-order tensor U s1s2 where δ t t is the Kronecker delta tensor (or identity matrix). This isometric constraint is depicted in Fig. 6(b).
It is sometimes convenient to view the collection of feature vectors Φ(xj) indexed over each training input xj as a single tensor.
The feature space covariance matrix can be formed by contracting two copies of Φ over the training set index j. In practice one does not form this entire object but uses efficient algorithms to compute reduced covariance matrices.
which is an isometry mapping two vector spaces whose bases are labeled by s1 and s2 into a single vector space labeled by t obeys the condition that (b) contracting U and U † over the s1 and s2 indices yields the identity matrix (represented diagrammatically as a line).
The dimension of t can be less than or equal to the product of the dimensions of s 1 and s 2 . The isometric property of U means it can be interpreted as a unitary rotation followed by a projection. This straightforward interpretation could be useful for interpreting and analyzing a learned model after training.
To compute the first isometry U 1 in the network we want to construct, define U 1 such that when it acts on the first two feature space indices s 1 and s 2 of the tensor Φ s j it maximizes the fidelity F , defined as After coarse graining the feature map using the isometry U 1 the fidelity of the resulting approximation to ρ is as shown in Fig. 7. Because U 1 is an isometry, the fidelity F 1 of the coarse-grained feature vectors U † 1 Φ j is always less than or equal to F . The fidelity will be maximized if U 1 projects out a subspace of the full feature space within which the components of the feature vectors Φ s j are on average very small.
To solve for the optimal isometry U 1 , it is convenient to introduce the reduced covariance matrix ρ 12 , defined by tracing over all of the indices s 3 , s 4 , . . . , s N of ρ as shown in Fig. 8(a). The motivation for introducing the reduced covariance matrix can be seen in the last two panels of Fig. 7. The manipulations there show the fidelity F 1 can be written in terms of ρ 12 as It follows that the optimal isometry U 1 can be computed by diagonalizing ρ 12 as here viewing ρ 12 as a matrix with row index (s 1 s 2 ) and column index (s 1 s 2 ) as shown in Fig. 8. The matrix P 12 is a diagonal matrix whose diagonal elements are the eigenvalues of ρ 12 . After the diagonalization, the columns of U 1 are chosen to be the eigenvectors corresponding to the D largest eigenvalues of ρ 12 . Let the rank of the matrix ρ 12 be r and call its eigenvalues {p i } r i=1 . One way to determine the number D of eigenvalues to keep is to choose a predetermined threshold and define D such that the truncation error E is less than , where the truncation error is defined as This is the same procedure proposed by White in the context of the density matrix renormalization group The fidelity F is defined as the average inner product of the training set feature vectors, or equivalently the trace of the covariance matrix ρ. The isometry U1 is chosen to maximize the fidelity following coarse graining (second panel above) which is equivalent to maximizing the trace of ρ12 after coarse graining (last panel above). (DMRG) algorithm used in quantum mechanics, where the Φ s j is analogous to an ensemble of wavefunctions enumerated by j; ρ is the full density matrix; and ρ 12 is a reduced density matrix [24].
To compute the remaining isometries which will form the first layer, the procedure continues in an analogous fashion by next computing the reduced covariance matrix ρ 34 as in Fig. 8(c) and diagonalizing it to obtain the isometry U 34 . Note that the calculation of the reduced covariance matrices as well as the individual summations over the training data used to produce them can be performed in parallel. What is more, we find that when summing over the training data in a random order, the reduced covariance matrices typically converge before summing over the entire training set, and this convergence can be monitored to ensure a controlled approximation. After diagonalizing the reduced covariance matrices for every pair of local indices (s 2i−1 , s 2i ), one obtains the first layer of isometries depicted in Fig. 9.
The isometry layer can now be used to coarse grain each of the training set feature vectors Φ(x j ). After the coarse graining, one can repeat the process again to make a second layer by making a new covariance matrix in the coarse grained feature space and diagonalizing its reduced covariance matrices.
If the coarse-graining procedure is repeated log 2 (N ) times (assuming N is a power of two), the end result is that one has approximately diagonalized ρ in the form shown in Fig. 10, with U expressed in tree tensor network form. With the particular choice of ρ made above, the resulting algorithm is an approximate implementation of kernel principal component analysis (kernel PCA) [49], but computed directly in feature space, not in the dual formulation. If one made no truncations throughout the algorithm, the equivalence to kernel PCA would be exact, but in practice the dominant eigenvectors of ρ are approximated by a tensor network with smaller internal indices than needed to capture U exactly. The fact that the calculation can be performed efficiently results from both the choice of a feature map with low-rank structure Eq. 8 and on the iterative algorithm based on tree tensor networks. But whether the approximation is accurate as well as efficient depends on the particular data set.
Similar to how a principal component analysis (PCA) can be used as a preprocessing step for other learning tasks, in the next section we will use the tree tensor network algorithm to compute a reduced set of features for supervised learning. However, unlike typical kernel PCA approaches, the tree tensor network approach discussed above imparts an explicit layered structure to the resulting model. In future work it would be very interesting to analyze the layers to see if they learn a hierarchy of features as observed in works on deep neural networks [16,17]. It may also be the case that the learned tree tensor network representation of the data acts as a form of regularization when used within other tasks.

IV. SUPERVISED OPTIMIZATION OF THE TOP TENSOR
Having determined a tree like network U of isometry tensors, we can now use this tree network as a starting point for optimizing a model for a supervised task. The specific model we explore in this section is shown in Fig. 1, and consists of the isometry layers U and a top tensor w. We emphasize that the isometry layers forming the tree tensor network U are obtained in an unsupervised manner, just by using the criterion of projecting out directions in feature space not spanned by the most significant eigenvectors of ρ. Only the top tensor w will be optimized for the supervised task we are interested in; the layers U will be left fixed. Of course we could optimize all of the layers for the supervised task, but here we are interested in the question of whether layers computed using a purely unsupervised algorithm are a good enough representation of the data for other tasks.
Recall that when producing the layers U, one progressively coarse grains the training set feature vectors as shown in Fig. 9. Thus for the purposes of optimizing the top tensor w, one only needs to use the reduced representation of the training inputs x j which are order-two tensorsΦ t1t2 (x j ) defined as These coarse-grained feature tensors reside in the same space as w t1t2 , so we can conveniently write the supervised model as In the experiments below, we choose to define the supervised task using a quadratic cost function.
To extend the class of supervised models proposed above to the multi-task setting, one generalizes the model f (x) to a family of models f (x) by training a collection of top tensors w , one for each label . In this multi-task setting, we found it sufficient to compute a single tree network U which is shared between the different models (different top tensors) specialized for each task.
To test the proposal, we consider the MNIST dataset of grayscale images of handwritten digits, consisting of 60,000 training images and 10,000 test images of size 28 × 28 with ten labels 0-9 [59]. For the local feature maps φ sn (x n ) we choose as proposed in Ref. 5. This is certainly not the only choice of feature map one could make, but has the advantage of providing a simple interpretation of the resulting model as a high-order polynomial.
To format the input, we scale each pixel value x n to be in the interval [0, 1]. To arrange the pixels of the two-dimensional images into components x n of a vector, we rasterize the image, simply proceeding along the first row, then the second row, etc. Other arrangements of the pixel data might lead to better outcomes and would be interesting to explore, such as choosing pixels which are spatially close in two dimensions to be grouped together as components the vector x.
The main hyper parameter controlling the experiments is the truncation error cutoff used when making the tree isometry tensors. To optimize the top tensor, we performed conjugate gradient optimization until the cost function reached convergence, typically requiring a few hundred iterations. The bond dimensions of the isometry layers (dimensions of the internal indices between layers) was very sensitive to the truncation error cutoff Eq. (19) that was used, yet fortunately the performance of the model also increased rapidly with lower cutoffs.
For a truncation cutoff = 10 −3 we found the bond dimensions connecting layers of isometries grew adaptively from about 3 between the first and second layers, to a maximum of 15 between the sixth and seventh layers, finally leading to the top two indices t 1 and t 2 having dimensions 107 and 151 respectively. Optimizing the top tensor of this model for the supervised task gave a cost function value of C = 0.078 yielding 98.75% accuracy on the training set and 97.44% accuracy on the test set.
For a cutoff = 6 × 10 −4 we found the bond dimensions connecting layers of isometries became gradually larger in intermediate layers, then rapidly increased in the topmost layers resulting in top indices t 1 and t 2 of size 328 and 444 respectively. Though this is perhaps a large coarse-grained feature space in absolute terms, it represents a significant reduction (coarse graining) of the original feature space size of 2 784 . Optimizing the top tensor of this network for the supervised task gave a cost function value of C = 0.042 yielding 99.68% accuracy on the training set and 98.08% accuracy on the test set.
Neither of the above experiments used an explicit form of regularization, such as a weight penalty. Nevertheless the test set generalization was good. This suggests that the form of the model based on training only the top layer for the specific task while computing the lower layers just from statistical properties of the data could aid generalization. However, more work and experiments on more challenging data sets is needed to verify this intuition.

V. MIXED TASK-SPECIFIC / UNSUPERVISED ALGORITHM
An interesting extension of the above approach is that one can mix the feature space covariance matrix ρ defined in Eq. (10) together with another matrix based on a specific task when computing the tree isometry tensors. In this way one can provide a prior guess for a supervised task, for example, by adapting or biasing the tree tensors toward the prior guess for the task, but still resulting in models with better performance than the guess due to the partially unsupervised nature of the algorithm. (Of course the specific task could itself be unsupervised as well, such as fitting the distribution of the data by minimizing the negative log likelihood.) To provide an initial guess for the supervised problem, one can use a low-cost approach such as a linear classifier f lin (x) = V · x. Then a model of the form Eq. (9) can be written in terms of an MPS representation of the weights W whose output is guaranteed to be equal to that of the linear classifier. See Appendix B for the details of this construction. Of course the prior weight MPS does not have to be initialized by a linear classifier only, but could be further optimized for the supervised task using the techniques of Refs. 5 and 6. By representing the initial guess as a tensor network-such as a matrix product state (MPS)-one can retain the ability to efficiently compute the reduced covariance matrices needed for computing the tree tensors.
Assuming we have obtained a prior guess for supervised task weights in the form of an MPS, begin the coarse graining procedure by defining a covariance matrix ρ W from the provided weights as or in the multi-task case, Again define ρ as a sum of outer products of the training data feature vectors as before in Eq. (10). Normalize both ρ and ρ W to have unit trace, defininĝ Now choose an empirical mixing parameter µ ∈ [0, 1] and define ρ µ as so that for µ = 0 the covariance matrix ρ µ will be the same as for the purely unsupervised case (up to normalization) but for µ > 0 the tree tensors making up U will be adapted to represent both the training data and the weights W . For the case µ = 1, if one takes a small enough cutoff when computing the tree tensors then the network U will be adapted to exactly reproduce the provided weights W regardless of the data.
Having defined ρ µ , the procedure to determine each layer of tree tensors is quite similar to the purely unsupervised algorithm of Section III, just with ρ µ substituted for ρ. Figure 11 shows the first step of the algorithm which is to compute the reduced covariance matrix ρ µ 12 from ρ µ , then diagonalize it to compute the first tree tensor. The key algorithmic difference from the unsupervised case is that after determining each layer, one must also coarse grain the provided weights W along with the training data so one can compute the reduced covariance matrices from ρ µ at the next scale. Although the weights W in MPS form have additional internal indices as shown in Fig. 2(a), it is straightforward to coarse grain an MPS with a tree tensor network layer: one simply contracts each isometry with pairs of MPS tensors.
For the case of a multi-class supervised task there will be multiple prior weight MPS W , one for each label (or one can equivalently provide a single MPS with an external or uncontracted label index). To generalize the above algorithm to the multi-task setting, one defines the covariance matrix ρ W as the sum over the covariance matrices of each of the prior supervised weights W To test whether the strategy of mixing in a prior estimate of the supervised task weights results in an improved model, we experiment again on the MNIST handwritten digits data set. Using a mixing parameter µ = 0.5 and a truncation error cutoff = 4×10 −4 results in a tree tensor network with top index sizes 279 and 393, where after making the tree layers in a single pass, only the top tensor is optimized further for the supervised task. Despite the top index sizes being significantly smaller than those for the best experiment in Sec. IV (where the sizes were 328 and 444), the results are slightly better: the cost function value is C = 0.0325, training set accuracy is 99.798%, and test set accuracy is 98.110%. This experiment strongly suggests that mixing weights trained for the supervised task with the covariance matrix based purely on the data leads to a representation of the data more suited for the specific task, which can be compressed further without diminishing performance.

VI. PARTIAL COARSE GRAINING: TREE CURTAIN MODEL
While the approaches in the previous sections involved computing tree tensor networks with the maximum number of layers, computing fewer layers can balance the benefits of a compressed data representation against the loss of expressiveness from accumulated truncations when computing more layers.
One interesting aspect of computing fewer tree layers is that after coarse graining, the data are still represented as high-order tensors, similar to Fig. 9. Specifically, the order of the data tensors after R rescalings will be N top = N/2 R . Therefore to complete the model, the top tensor must also be a tensor of order N top if the output is to be a scalar, or order N top + 1 for vector-valued output in the multi-task case. So to complete the model one can use another type of tensor network to represent the top tensor, such as a matrix product state.
Choosing a matrix product state (MPS) form of the top tensor results in the architecture shown in Fig. 12. After coarse graining the training data through the tree layers, one can optimize the top MPS using previously developed methods for supervised [5,6] or unsupervised [39] learning tasks. The resulting model resembles an MPS with a tree "curtain" attached.
To test the effectiveness of the partial coarse-grained approach, and the resulting tree-curtain model, we study the fashion MNIST dataset [60]. Similarly to MNIST, the data set consists of 28 × 28 grayscale images with ten labels, with 60,000 training and 10,000 test images. However, the supervised learning task is significantly more challenging than MNIST because the images consist of photographs of a variety of clothing (shirts, shoes, etc.). Example images from the data set are shown in Fig. 13.
To compute a tree curtain model for fashion MNIST, we first training a linear classifier resulting in only 83% test accuracy. We then represent each of the linear classifier vectors V as MPS and used the mixed covariance matrix approach (Sec. V) with mixing parameter µ = 0.9 to optimize four tree tensor layers. Using a truncatioǹ Computing only a few tree layers results in a model with a high-order top tensor. In the model above the top tensor w is represented as a matrix product state with a label index appropriate for the multi-task case.
FIG. 13. Examples from the fashion MNIST data set [60]. error cutoff = 2 × 10 −9 , the final coarse-grained feature indices attached to the top MPS reach a maximum of about 30. We then fix the internal bond dimension of the top MPS to be 300 and minimize the quadratic cost using thirty sweeps, or passes of alternating least squares (single-site DMRG) optimization.
The optimized model reaches 95.38% accuracy on the training set, and 88.97% accuracy on the testing set. While ∼ 89% test accuracy is significantly less than achieved on the much easier MNIST handwriting dataset, many of the available benchmarks using state-of-theart approaches for fashion MNIST without preprocessing are in fact comparable to the results here, for example XGBoost (89.8%), AlexNet (89.9%), and a two-layer convolutional neural network trained with Keras (87.6%). Better results are attainable; the best we are aware of is a GoogLeNet reaching 93.7% test accuracy. But the fact that the architecture discussed here yields similar performance to other powerful and standard approaches for a challenging data set is an encouraging result. With further hyperparameter optimization and more efficient training algorithms such as stochastic gradient descent, we expect better results can be achieved.

VII. RELATED WORK
The feature space covariance matrix ρ s s = j Φ s j Φ †j s is closely related to the kernel matrix K j j = s Φ †j s Φ s j , which is a central quantity in the theory of kernel learning. That both ρ and K have the same spectrum can be seen from the SVD of Φ s j . Thus the idea of exploiting low-rank approximations to the kernel matrix is closely related to the present work [61][62][63], as well as the idea of feature sampling or random features [64,65]. But in contrast, the the present work constructs a low-rank approximation for ρ directly in feature space, with an algorithm nevertheless scaling linearly in both training set size and input space dimension. The tensor network form also imparts an interesting structure onto the resulting model, allowing us to experiment with training the top layer only. And we considered other choices besides ρ, namely a mixed covariance matrix ρ µ , which produces a more compressed representation of the data without reducing performance on a supervised task.
The algorithm in Sec. III has many precedents in the tensor network literature, such as discussions by Vidal [23,66] of computing tree networks layer-wise to optimize fidelity of density matrices (pure and mixed states), and multigrid DMRG [67]. Similar algorithms have also been employed for quantum state tomography [68,69]. Recent work by Nouy [70] develops a related algorithm to reconstruct a function using a tree tensor network.
Another related line of work is the study of convolutional arithmetic circuits or ConvACs [4,10,40], a type of neural network constructed from a tree tensor network [36]. The present work could be viewed as a proposal for an unsupervised, adaptive algorithm for training ConvAC neural networks, along with variations such as replacing the upper layers with an MPS.

VIII. DISCUSSION
The idea of coarse-graining in statistical physics motivates a similar approach to real-world data. The algorithm presented here resembles kernel PCA, but computed directly in feature space with the diagonalizing unitary approximated as a layered tensor network. Despite, or perhaps because of the approximations made, the resulting coarse-grained representation of the data could be used as a good starting point for other supervised or unsupervised tasks. The idea suggests interesting extensions such as exploring quantities besides the covariance matrix ρ around which to base the algorithm. Using a mixed covariance matrix ρ µ (Sec. V) gave an improved representation, but there could be other good choices. Instead of a tree tensor network, which has certain known deficiencies as a coarse-graining scheme, it would be very interesting to use a MERA tensor network [23,31].
In an engineering sense, there is much room to improve the above algorithms, such as using tree tensors with two-dimensional groupings of indices [37]. Another idea is feeding the optimized model back through the mixed algorithm of Sec. V to further adapt the tree tensors.
That PCA is similar to coarse-graining, or the renormalization group, has been observed recently by Bradde and Bialek [18]. It would be very interesting to use the tools discussed here to analyze various data sets from a renormalization group perspective. An especially welcome outcome would be if the analysis could identify distinct classes of data, and determine which machine learning architectures have an inductive bias [10] suited to each class.

IX. ACKNOWLEDGEMENTS
We thank David J. Schwab for an especially careful reading of the paper and many helpful comments. We also thank Glen Evenbly, William Huggins, Matthew Fishman, Steven White, John Terilla and Yiannis Vlassopoulos for helpful comments and discussions.
Appendix A: Tensor Diagram Notation Tensor diagram notation is a simple, graphical notation for representing networks of contractions of highorder tensors. The main rules of this notation are: • A tensor with n indices (an order n tensor) is represented as a shape with n lines coming out of it. Each line represents a specific index but indices do not have to be labeled or named when they can be distinguished by context.
• Connecting index lines between a pair of tensors indicates that these indices are contracted, or summed over.
Other informal conventions may be used by specific authors, such as choosing special shapes or decorations of tensors to denote special properties. For example, in this work a triangular shape indicates an isometric tensor. Unless otherwise specified, the only information that matters is the connectivity of the network, with no fixed meaning given to the orientation or ordering of the index lines. Furthermore, no special properties of the tensors are assumed such as symmetries or transformation properties, unless specifically stated. We emphasize that tensor diagram notation is completely rigorous, and just a way of notating complicated sums. Unlike similar-looking neural network diagrams, every element of a tensor diagram is a purely linear transformation (of the vector spaces whose basis elements are labeled by the indices).
The main advantage of tensor diagram notation is that it frees one from having to assign names to every index, which can be cumbersome for complicated tensor networks. The notation also makes it easier to compactly express sophisticated algorithms and to visually understand different families of tensor networks. Various operations such as traces, outer products, transposes, or matricization can be expressed implicitly in diagrammatic form without requiring additional notation or symbols. For example, the diagrammatic expression of an outer product is simply the placement of two tensors nearby each other.

Appendix B: Representing a Linear Classifier as an MPS
One can write a model of the form Eq. (9) with the weights in MPS form such that the output of the model is guaranteed to be the same as a linear classifier with weights V , assuming local feature maps of the form Eq. (22) namely φ(x) = [1, x]. This mapping was first discussed in Ref. 5.
Say we want to make our model output equivalent to a linear classifier with parameters V First, assume the weights of our model Eq. (9) are in MPS form, meaning Now define each MPS tensor A sj αj−1αj for each fixed value of s j as follows: For the first A tensor one only takes the second row of each matrix and for the last A tensor one only takes the first column of each matrix. To include a constant shift, either add a fictitious input component x 0 = 1 to the definition of x or suitably adjust the first "A" tensor. This mapping from a linear classifier can also be extended to more general local feature maps φ s (x) of dimension d by suitably defining an "extended linear classifier" model [2] Hanie Sedghi and Anima Anandkumar, "Training inputoutput recurrent neural networks through spectral methods," arxiv:1603.00954 (2016).