{'fraction_non_alphanumeric': 0.03136186354650852, 'fraction_numerical': 0.018085163220601594, 'mean_word_length': 4.583146905294556, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Locally orderless tensor networks for classifying two- and three-dimensional medical images

Tensor networks are factorisations of high rank tensors into networks of lower rank tensors and have primarily been used to analyse quantum many-body problems. Tensor networks have seen a recent surge of interest in relation to supervised learning tasks with a focus on image classification. In this work, we improve upon the matrix product state (MPS) tensor networks that can operate on one-dimensional vectors to be useful for working with 2D and 3D medical images. We treat small image regions as orderless, squeeze their spatial information into feature dimensions and then perform MPS operations on these locally orderless regions. These local representations are then aggregated in a hierarchical manner to retain global structure. The proposed locally orderless tensor network (LoTeNet) is compared with relevant methods on three datasets. The architecture of LoTeNet is fixed in all experiments and we show it requires lesser computational resources to attain performance on par or superior to the compared methods.


Introduction
Support vector machines (SVMs) and other kernel based methods ushered in a new wave of supervised learning methods based on the fundamental insight that challenging tasks in low dimensions may become easier when the data is lifted to higher dimensional spaces (Boser et al., 1992;Cortes and Vapnik, 1995;Hofmann et al., 2008), as illustrated in Figure 1. These methods, however, become prohibitive when dealing with massive datasets in high dimensional feature spaces as their space complexity grows at least quadratically with the number of data points (Bordes et al., 2005;Novikov et al., 2016). Further, SVMs are also known to be sensitive to the specific choice of the kernel parameters and as a consequence the decision boundaries learnt are prone to over-fitting (Burges, 1998;Bordes et al., 2005). These factors have discouraged their successful adoption to tasks involving large datasets comprising high resolution images where deep learning based methods have shown to fare better .
Tensor networks, also known as tensor trains, offer a different and more efficient framework to dealing with such high dimensional spaces. Fundamentally, tensor networks are factorisations of higher order 2 tensors into networks of lower order tensors (Fannes et al., 1992;Oseledets, 2011;Bridgeman and Chubb, 2017). The number of parameters needed to specify an order-N tensor using tensor networks can be drastically reduced, from exponential to polynomial dependence on N (Perez- Garcia et al., 2006). This massive reduction in number of parameters using tensor networks has been predominantly applied to better understand quantum wave functions (Shi et al., 2006). They have also seen applications in data compression (Cichocki et al., 2016), and recently to better understand the expressive power of deep learning models Glasser et al., 2019).
Figure 1: Data that is only non-linearly separable in lower dimensions can become linearly separable in higher dimensions, illustrated here for a simple case of 2-dimensional data that becomes linearly separable when lifted into 3-dimensions. This is the underlying principle in kernel based methods and tensor networks which learn linear decision boundaries in high dimensional spaces.
Recently, there has been an increasing interest in using tensor networks for supervised learning, specifically focused on image classification tasks (Stoudenmire and Schwab, 2016;Klus and Gelß, 2019;Efthymiou et al., 2019;Sun et al., 2020). These methods rely on transforming two dimensional input images (order-2 tensors) into one dimensional vectors (order-1 tensors) to obtain linear decision boundaries in high dimensional spaces. Due to the constraint of flattening to obtain vector inputs these methods are constrained to work with images of small spatial resolution (12x12 px to 28x28 px). Several improved flattening strategies have been attempted to maximize the retained spatial correlation between pixels as the correlation declines exponentially causing loss of information in larger images (Stoudenmire and Schwab, 2016;Efthymiou et al., 2019;Cheng et al., 2019;Sagan, 2012). For small enough images (like in MNIST 3 or Fashion MNIST 4 datasets) there is some residual correlation in the flattened images which can be exploited by lifting the vector input to higher dimensions using tensor networks. This might not be effective in the case of images with higher spatial resolution. Further, the information lost by flattening of images in medical imaging tasks can be critical as the predicted decisions could be dependent on the global structure of the image.
In this work, we present a tensor network based model adapted for classifying two-and three-dimensional medical images which can be of higher spatial resolution. The method presented here relies on lifting small image (slice or volume) patches to high dimensions and performing tensor network operations on them. These small, local regions can be treated as being locally orderless drawing parallels to the classical theory of locally orderless images in Koenderink and Van Doorn (1999). As with locally orderless image analysis, we propose to extract useful representations of small image regions in higher dimensions and aggregate Figure 2: Left: Tensor notation depicting a scalar S, vector V i , matrix M ij and a general order-3 tensor T ijk . Center: Tensor notation for matrix multiplication or tensor contraction, which are used extensively in the matrix product state networks used in this work. We adhere to the convention that the contracted indices are written as subscripts. Right: Tensor notation for trace of product of two matrices. them in a hierarchical manner resulting in our locally orderless tensor network (LoTeNet) model. The proposed LoTeNet model is used to learn decision functions in high dimensional spaces in a supervised learning set-up and is optimized end-to-end by backpropagating the error signal through the tensor network. We conduct experiments for classifying three medical imaging datasets: 2D histopathology images in PCam dataset (Veeling et al., 2018), 2D thoracic computed tomography (CT) slices from LIDC-IDRI dataset (Armato III et al., 2004) and 3D magnetric resonance imaging (MRI) from the OASIS dataset (Marcus et al., 2007). These experiments show that the LoTeNet model fares comparably to relevant stateof-the-art deep learning methods requiring the tuning of a single model hyperparameter while utilising only a fraction of the graphics processing unit (GPU) memory when compared to their convolutional neural network (CNN) counterparts. The work presented here is an extension of an earlier version of the method which operated on two-dimensional data and is published as a conference publication (Selvan and Dam, 2020).
The key contributions in this work are: 1. A novel tensor network based model for classifying 2D and 3D medical image data 2. Extending supervised learning with tensor networks to images of higher resolution 3. Validation of the method on three public datasets and comparison to relevant methods 4. Demonstrating competitive performance using little GPU resources 5. Conceptual connections of tensor networks with deep neural networks In the remainder of this manuscript we present the basics of tensor notation, tensor network fundamentals and formulate the supervised image classification task in Section 2. The components of the proposed LoTeNet model and the final model itself are presented in Section 3. Experiments on the three datasets and results are presented in Section 4, along with the discussions and future research directions in Section 5, and present overall conclusions from this work in Section 6.

Background and Problem Formulation
We introduce key concepts pertaining the use and optimisation of tensor networks in this section which will be put together to describe the proposed method in the following sections.

Tensor Network Notations
Tensor networks and operations on them are described using an intuitive graphical notation, introduced in Penrose (1971). Figure 2 (left) shows the commonly used notations for a scalar S, vector V i , matrix M ij and a general order-3 tensor T ijk . The number of dimensions of a tensor is captured by the number of edges emanating from the nodes denoted by the edge indices. For instance, the vector V i is an order-1 tensor indicated by the edge with index i and an order-3 tensor has three indices (i, j, k) depicted by the three edges, and so on.
Tensor notations can also be used to capture operations on higher-order tensors succinctly as shown in Figure 2 (center) where matrix product is depicted, which is also known as tensor contraction. The edge between the tensor nodes X i j and Y k j is the dimension subsumed in matrix multiplication resulting in the matrix Z ik . More thorough introduction to tensor notations can be found in Bridgeman and Chubb (2017).

Linear Model in Exponentially High Dimensions
A linear model in a sufficiently high dimensional space can be very powerful (Novikov et al., 2016). In SVMs, this is accomplished by the implicit mapping of the input data into an infinite dimensional space using radial basis function kernels (Hofmann et al., 2008). In this section, we describe the procedure followed in recent tensor network based works, including in ours, to map the input data into an exponentially high dimensional space.
Consider an input vector x ∈ [0, 1] N , which can be obtained by flattening a 2D or 3D image with N pixels with intensity values that are normalized in the interval [0, 1]. A commonly used feature map for tensor networks is obtained from the tensor product of pixel-wise feature maps (Stoudenmire and Schwab, 2016): where the local feature map acting on a pixel x j is indicated by φ i j (·). The indices i j run from 1 to d where d is the dimension of the local feature map. The feature maps are usually simple non-linear functions restricted to have unit norm such that the joint feature map in Eq. (1) also has unit norm. A widely used local feature map with d = 2 inspired from quantum wave function analysis is (Stoudenmire and Schwab, 2016) is shown in (2), and a simpler local feature map from Efthymiou et al. (2019) in Eq. (3): The joint feature map Φ(x) in Eq.
(1) is an order-N tensor due to tensor product of the N order-1 local feature maps of dimension d in Eq. (2). The joint feature map Φ(x) can be seen as mapping each image to a vector in the d N dimensional feature space (Stoudenmire and Schwab, 2016). For multi-channel inputs, such as RGB images or other imaging modalities, with C input channels the local feature map can be applied to each channel separately such that the resulting space is of dimension (d · C) N (Efthymiou et al., 2019). Step 3) Series of tensor contractions done with MPS to compute W m · Φ(x) in Eq. (5) Given this high dimensional joint feature map Φ(x) of Eq. (1) for the input data x, a linear decision rule for a multi-class classification task can be formulated as: where m = [0, 1, . . . M − 1] are the M classes, and the weight tensor W m is of order-(N + 1) with the output tensor index m.
The tensor notation representation of the linear model of Eq. (5) is shown in Figure 3 (Step 1) where the first column of gray nodes are the individual pixel feature maps of order-1 with feature dimension d. The local feature maps are connected to the order-(N + 1) weight tensor W m along N edges with one edge for output of dimension M marked with index m.
The order-(N +1) weight tensor W m results in a total of M · d N weight elements. Even for a relatively small gray scale image, say of size 100 × 100, the total number of weight elements in W m can be massive: 2 · 2 10000 ≈ 10 3010 . In the next section we will see how tensor networks can represent such high dimensional tensors with parameters that grow linearly with N instead of growing exponentially with N .

Matrix Product State (MPS)
Consider two vectors, T i and U j with indices i and j respectively. The tensor product 5 of these two vectors yields an order-2 tensor or a matrix X ij . The matrix product state (MPS) (Perez-Garcia et al., 2006) is a type of tensor network that expands on this notion of tensor products allowing the factorization of an order-N tensor (with N edges) into a chain of lower-order tensors. Specifically, the MPS tensor network can factorise an order-N tensor with a chain of order-3 (with three edges) except on the borders where they are of order-2, 5. Note that the tensor product of two order-1 tensors is the same as the vector outer product. as shown in Figure 3 (Step 2). Consider a tensor of order-N with indices i i , i 2 , . . . i N , using MPS it can be approximated using lower-order tensors as where A i j are the lower-order tensors. The subscript indices α j are virtual indices that are contracted and are of dimension β referred as the bond dimension. The components of these intermediate lower-order tensors A i j form the tunable parameters of the MPS approximation. Note that any N dimensional tensor can be represented exactly using an MPS tensor network if β = d N/2 . In most applications, however, β is fixed to a small value or allowed to adapt dynamically when the MPS tensor network is used to approximate higher-order tensors (Perez-Garcia et al., 2006;Miller, 2019).
One of the drawbacks of using MPS tensor networks is that they operate along one dimension (as a chain). This is the primary reason that two dimensional image data has to be first flattened to a vector when working with tensor networks such as the MPS. Tensor networks that can work on arbitrary graphs, which might be more suitable for image data like the projected entangled pair states (PEPS) (Verstraete and Cirac, 2004), are not as well understood and do not yet have efficient algorithms like the MPS. Dot product with MPS: The decision function in Eq. (5) comprising the dot product of the order-(N +1) weight tensor W m and the joint feature map Φ(x) in Eq. (1) can be efficiently computed using the MPS approximation in Eq. (6), depicted in Figure 3 (Step 2). The order in which tensor contractions are performed can yield a computationally efficient algorithm. The original MPS algorithm (Perez-Garcia et al., 2006) starts from one of the ends, contracts a pair of tensors to obtain a new tensor which is then contracted with the next tensor and this process is repeated until the output tensor is reached. The cost of this algorithm scales with {N ·β 3 ·d} when compared to the cost without the MPS approximation which scales with d N . In this work, we use the MPS implementation in Miller (2019) which performs parallel contraction of the horizontal edges and then proceeds to contract them vertically as depicted in Figure 3 (Step 3). The MPS approximation also reduces the number of tunable parameters by an exponential factor, from {M · d N } to {M · d · N · β 2 }.

Methods
The primary contribution in this work is a tensor network based image classification model that can handle high resolution images of both two and three spatial dimensions. Several recent tensor network models for supervised image classification tasks operate on small planar images and flatten them into vectors with different raveling strategies at the expense of global image structure (Stoudenmire and Schwab, 2016;Han et al., 2018;Efthymiou et al., 2019). In contrast to these methods, we only flatten small regions of the images which can be assumed to be locally orderless (Koenderink and Van Doorn, 1999) and obtain expressive representations of local image regions using tensor contractions in high dimensional spaces. We process these locally orderless image regions in a hierarchical fashion using layers of MPS blocks in the final model, which we call the locally orderless tensor network (LoTeNet), shown in Figure 5. We next describe the proposed LoTeNet model in detail.

Squeeze operation
The proposed LoTeNet model operates on small image regions and aggregates them with MPS operations at multiple levels. These small image patches at each level are created using the squeeze operation in two steps, as illustrated in Figure 4 for an order-2 tensor i.e an image with two spatial dimensions.
The input images are converted into vectors with inflated feature dimensions by applying kernels of stride k along each of the spatial dimensions. Consider an input image X hij with S = 3 spatial dimensions (height:H, width:V, depth:D), N pixels and C = d channels as feature dimension: X hij ∈ R H×V ×D×C . In the first step, input image is reshaped into smaller patches controlled by the stride k: X hij ∈ R (H/k)×(V /k)×(D/k)×d such that the feature dimension increases to d = C · k S . The stride of the kernel k decides the extent of reduction in spatial dimensions and the corresponding increase in feature dimension of the squeezed image. Similar strategies have also been used in normalizing flow literature to provide additional spatial context via feature dimensions such as in Dinh et al. (2017). In the second step, the reshaped image patches with inflated feature dimensions are flattened from order-S tensors to order-1 tensors, resulting in X h ≡ x ∈ R (N/k S ) with d = C · k S . This flattening operation provides spatial information in the feature dimension to LoTeNet, thus retaining additional image level structure when compared to flattening the entire image into a single vector x ∈ R N with d = C channels (Efthymiou et al., 2019). Further, the increase in the feature dimension d = C · k S makes the tensor network more expressive as it increases the dimensionality of the feature space (Stoudenmire and Schwab, 2016).
The transformations in dimensions due to the squeeze operation parameterised by the stride k, ψ(·; k), are summarised below: The squeeze operation can also be treated as a local feature map due to the increase in feature dimensions, similar to the more formal pixel level feature maps mentioned in Eq.
(2) and Eq. (3), operating on image patches instead of individual pixels. Finally note that the squeeze kernel stride can be different at each level l of the model denoted k l . In the current formulation, the images are assumed to be rectangular in each plane. As the main purpose of the squeeze operation is to increase the feature dimensions by flattening small local neighbourhoods, this can be achieved by reshaping any small, non-square region of the image into a vector of fixed feature dimension. The simplest strategy to work with non-square (circular images for instance) would be to pad the regions around the borders to make them rectangular. Another approach could be to partition non-square images into cells using Voronoi tesselations (Du et al., 1999) and squeezing these cells into vectors.

Locally orderless tensor network (LoTeNet)
An overview of the proposed LoTeNet model is shown in Figure 5 for an image with two spatial dimensions i.e with S = 2. The number of squeezed vectors and the corresponding MPS blocks at each layer of the proposed model are also indicated in Figure 5.
The LoTeNet model comprises of layers of MPS blocks interleaved with squeeze operations without any non-linear components. The input to the first layer in Figure 5 is the full resolution input image with N pixels, shown with grids marking the k × k patches which are then squeezed to obtain N 1 = N/k S vectors with d = C · k S . Each of these squeezed patches are input into MPS blocks which contract the d = C · k S vectors to output a vector with dimension d = ν. MPS blocks operating on squeezed image patches can be interpreted as summarising them with a vector of size ν using a linear model in a higher dimensional feature space. In LoTeNet, we set ν to be the same as the MPS bond dimension β, so that there is only a single hyperparameter to be tuned.
The output vectors from all N l MPS blocks at a given layer l are reshaped back into the S dimensional image space. However, due to the MPS contractions, these intermediate image space representations will be of lower resolution as indicated by the smaller image with fewer patches in Figure 5. This is analogous to obtaining an average pooled version of the intermediate feature maps in traditional CNNs. The intermediate images of lower resolution formed at layer l are further squeezed and contracted in the subsequent layers of the model. This process is continued for L − 1 layers and the final MPS block acting as layer L contracts the output vector from layer L − 1 to predict the decisions in Eq. (5). The weights of all the MPS blocks in each of the L layers are the model parameters which are tuned in a supervised setting as described next.

Model Optimization
We view the sequence of MPS contractions in successive layers of our model as forward propagation and rely on automatic differentiation to compute the backward computation graph (Efthymiou et al., 2019). Torch MPS package (Miller, 2019) is used to implement MPS blocks and trained in PyTorch (Paszke et al., 2019) to learn the model parameters from training data in an end-to-end fashion. These parameters are analogous to the weights of neural network layers and can be updated in a similar iterative manner by backpropagating a relevant error signal computed between the model predictions and the training labels.
We minimize the cross-entropy loss between the true label y i ∈ [0, . . . , M − 1] for each image x i ∈ D and the predicted label f (y i ) (x i ) in the training set D: where σ(·) is the softmax operation used to obtain normalized scores that can be interpreted as the predicted class probabilities. For binary classes, the loss in Eq. (8) reduces to the binary cross entropy and output dimension M = 1 can be used with the sigmoid nonlinearity to obtain probabilistic predictions in [0, 1].

Data
We report experiments on three public datasets for the task of binary classification from 2D histopathology slides, 2D thoracic computed tomography (CT) and 3D T1-weighted MRI scans.  Figure 6(a). In this work, a modified PCam dataset from the Kaggle challenge 6 is used that removes duplicate image patches included in the original dataset. This dataset consists of about 220, 000 patches for training and an independent test set of about 57, 500 patches is provided for evaluating the 6. https://www.kaggle.com/c/histopathologic-cancer-detection models. Data augmentation is performed during training with 0.5 probability comprising of rotation, horizontal and vertical flips. LIDC Dataset: The LIDC-IDRI datatset comprises of 1018 thoracic CT images with lesions annotated by four radiologists (Armato III et al., 2004). We use the 128 × 128 px 2D slices from Knegt (2018) (shown in Figure 6 (b)) yielding a total of 15, 096 patches. The segmentation masks from the four raters are transformed to binary labels indicating the presence (if two or more radiologists marked a tumour) or absence of tumours (if less than two raters marked tumours). These binary labels capture a majority voting among the radiologists. The dataset is split into 60 : 20 : 20 splits for training, validation and a hold-out test set. OASIS Dataset: The OASIS dataset consists of T1-weighted MRI scans of 416 subjects aged between 18 to 96 years (Marcus et al., 2007). The scans are graded with clinical dementia rating (CDR) into four classes: non-demented=0, very mild dementia=0.5, mild dementia=1.0, moderate dementia=2.0. We create an Alzheimer's disease classification dataset according to Wen et al. (2020) with scans of all subjects above 60 years resulting in a dataset with 155 subjects. Binary labels are extracted from the CDR scores: cognitively normal (CN) if CDR=0 and Alzheimer's disease (AD) if CDR>0, yielding 82:73 split between CN and AD cases. The MRI scans are preprocessed using the extensive pipeline described in Wen et al. (2020) comprising bias field correction, non-linear registration and skull stripping using the ClincaDL package 7 . After preprocessing we obtain volumes of size 128x128x128 voxels (Figure 6-c).

Implementation
The model is implemented in PyTorch (Paszke et al., 2019) based on the Torch MPS package (Miller, 2019). The architecture of LoTeNet model is kept fixed across all the datasets. The proposed LoTeNet model is evaluated with L = 4 layers, squeeze kernel of size k l = 2 except for the input layer where k 1 = 8 in order to reduce the number of MPS blocks in the first layer. The most critical hyperparameter of LoTeNet is its bond dimension β; it was set to β = 5 obtained from the range [2, 3 . . . 20] based on the performance on the PCam validation set. We set the virtual dimension ν to be the same as the bond dimension. We use the Adam optimizer (Kingma and Ba, 2014) with a learning rate of 5 × 10 −4 . We use a batch size of 512 for experiments on the PCam and LIDC datasets, and a batch size of 4 for the OASIS dataset. We incorporate batch normalisation (Ioffe and Szegedy, 2015) after each layer in all the models and we observed this resulted in faster and more robust convergence. All models were assumed to have converged if there was no improvement in validation accuracy over 10 consecutive epochs and the model with the best validation performance was used to predict on the test set. All experiments were run on a single Tesla K80 GPU with 12 GB memory. The development and training of all models (including baselines) in this work is estimated to produce 14.1 kg of CO2eq, equivalent to 120.1 km travelled by car as measured by Carbontracker 8 (Anthony et al., 2020). 7. https://clinicadl.readthedocs.io/en/latest/ 8. https://github.com/lfwa/carbontracker/

Experiments on 2D data
PCam models: Performance of the LoTeNet model is compared to the rotation equivariant CNNs in (Veeling et al., 2018) which is also state-of-the-art for the PCam dataset. Additionally, we compare to the single layer MPS model with local feature map of the form in Eq.
(3) with β = 10 from Efthymiou et al. (2019) reported here as Tensor Net-X and to Densenet  with 4 layers and a growth rate of 12. We compare the binary classification performance of the models using the area under the ROC curve (AUC) as it is not sensitive to arbitrary decision thresholds. The performance metrics are shown in Table 1 (left) along with the maximum GPU memory utilisation for each of the models. We observe the AUC on the test set attained by LoTeNet model (0.943) is comparable to the Rotation Eq. CNN (0.963) and Densenet models (0.962) on the PCam dataset with a drastic reduction in the maximum GPU memory utilisation; a mere 0.8 GB when compared to upto 11 GB for the Rotation Eq.CNN 9 and 10.5 GB for Densenet. We also notice a considerable improvement when compared to Tensor Net-X in the attained AUC (0.908). LIDC models: Similar to the PCam dataset, performance of LoTeNet model is compared using AUC and maximum GPU memory utilisation, and are reported in Table 1 (right). We compare to Densenet and Tensor Net-X with two bond dimensions (β = 5, 10) to highlight the influence of the bond dimension. LoTeNet model fares better than the compared models with an AUC of 0.874 whereas the Densenet model achieves 0.829 and Tensor Net-X achieves 0.847 (β = 10) and 0.823 (β = 5). The GPU memory utlisation follows a similar trend as with the PCam models with LoTeNet requiring only 0.7 GB.

Experiments on 3D data
LoTeNet used in the 3D experiments is identical to the one used in 2D experiments except for the local feature map. For the 3D experiments, we do not apply the local feature map in Eq.
(2) but instead rely on the increase in feature dimension due to the squeeze operation in Eq. (7). This reduces the number of parameters of the LoTeNet model and the risk of overfitting to the limited data in the OASIS dataset. The experimental set-up closely follows from Wen et al. (2020) and we use 5-fold cross validation to report the balanced accuracy on each of the test folds. Balanced accuracy (BA) is binary accuracy normalised by the class skew and helps when the classes are imbalanced. We reimplement the subject-level 3D CNN used in Wen et al. (2020)   We also report two other baseline models based on 3D CNN and a multi-layered perceptron (MLP) model. The CNN baseline model is similar to the subject-level CNN but has an initial feature map size of 32. The MLP baseline has four layers to match the LoTeNet model, and ReLU non-linear activation function. Table 2 summarises the performance of LoTeNet along with the compared methods for two types of inputs. As LoTeNet is agnostic to the dimensionality of the input image due to the squeeze operations, we demonstrate the performance on volume data (3D) and also on 19,840 slices extracted along the z-axis (2D) from the volume data. For the 2D experiments we ensure there is no data leakage by creating folds at subject level. The balanced accuracy across five test folds are shown and we see a clear improvement (0.71 ± 0.09) for the 3D case compared to subject-level CNN model (0.67±0.08). LoTeNet model has a lower score in the 2D case, indicating the model is able to extract more global information from the volume data than from 2D slices. Additionally, we report the number of parameters, maximum GPU utlisation and time per training epoch for all the methods. To show the influence of the number of parameters, we compare to a large MLP baseline model with 78M parameters and observe this increase does not necessarily improve the baseline model's performance.

Discussion
In this section, we focus on discussing high level trends observed with experiments reported in Section 4.2, conceptual connections to other methods, and point out possible directions for future work.

On Performance
The performance of LoTeNet across the three datasets is on par with, or superior, to the other methods as reported in Tables 1 and 2. This robustness across tasks is remarkable because the LoTeNet architecture is fixed across the tasks. This behaviour could be attributed to the single model hyperparameter -bond dimension, β -which is fixed (β = 5) across all experiments. Recollect that the bond dimension controls the quality of MPS approximations of the corresponding operations in high dimensional spaces as described in Section 2.3. Consistent with the reporting in earlier works (Efthymiou et al., 2019), we also find that the value of β after a specific value does not affect the model's performance. This is illustrated in Figure 7, where we show the sensitivity of the validation set performance on LIDC dataset for different bond dimensions. This is a desirable feature as the models can be expected to yield comparable performance when applied to different datasets without entailing dataset specific cost of elaborate architecture search. Additionally, this robustness could also be attributed to more efficient computation of massive number of parameters in LoTeNet resulting in generalised decision rules, as reported in Table 2 where we notice that the number of parameters in LoTeNet are higher than other methods (except MLP baseline). In Tables 1 and 2, we report the GPU memory requirement for each of the models. LoTeNet requires only a fraction of the GPU memory utilised by the corresponding baseline methods, including the Densenet  or Rotation Eq-CNN models (Veeling et al., 2018). One reason for this drastic reduction in GPU memory utilisation is because tensor networks do not maintain massive intermediate feature maps. This behaviour is similar to feedforward neural networks and unlike CNNs which use a large chunk of GPU memory mainly to store intermediate feature maps (Rhu et al., 2016). Further, as LoTeNet is based on contracting input data into smaller tensors in a hierarchical manner the memory consumption with successive contracted layers does not increase. This can be an important feature in medical imaging applications as larger images and larger batch sizes can be processed without compromising the quality of the learned decision rules.

Comparison with Neural Networks
Feed-forward neural networks: The class of tensor network models adapted for machine learning tasks (Stoudenmire and Schwab, 2016;Efthymiou et al., 2019) are conceptually most similar to feed-forward neural networks such as the MLP, as both classes of models operate on vector inputs (see Figure 8). The difference, however, arises in the type of decision boundaries optimised by each class of model as illustrated in Figure 1. MLP based feed-forward networks obtain non-linear decision boundaries using several layers of neurons and non-linear activation functions (such as sigmoid, ReLU etc.). Tensor network models, including the proposed model, optimise linear decision boundaries in exponentially high dimensional spaces without using any non-linear activation functions. Figure 8: Conceptual comparison of a feedforward neural network such as an MLP with the proposed LoTeNet model. Notice the functional correspondence between linear and MPS layers, and non-linear activation functions and the squeeze operations between the MLP and LoTeNet models, respectively. The input image is flattened to a vector to input to the MLP model, whereas for the LoTeNet model we squeeze small image regions before flattening as described in Section 3.1.
Convolutional neural networks: LoTeNet uses MPS blocks on image patches and aggregates these tensor representations in a hierarchical manner to learn the decision rule. While this could appear similar to the use of strided convolution kernels in CNNs, it is indeed closer to feed-forward neural networks. The use of an MPS block per patch perhaps can be interpreted as a form of strided convolution but only in how the local operation is performed on each patch. The primary difference is in the weight sharing; the weights of MPS blocks are not shared across the image. Each k × k image region has an MPS block acting on it. This is indicated in Figure 5, where we report the number of MPS blocks used per layer.

Limitations and Future Work
The computation time reported in Table 2 shows it to be higher for LoTeNet (45.3s) compared to the CNN (12.8s) or MLP (4.1s) baselines as efficient implementations of tensor contractions are not natively supported in frameworks such as PyTorch. There are ongoing efforts to further improve efficient implementations which can accelerate tensor network operations (Fishman et al., 2020;Novikov et al., 2020). Considering that LoTeNet usually converges between 10-20 epochs, the increase in computation time with the current implementation is not substantial as the model can be trained even on 3D data easily under an hour.
One possible drawback of having a single hyperparameter β controlling the MPS approximations is the lack of granular control of the model capacity. In the 2D OASIS experiments (Table 2), we further investigated the lower average balanced accuracy score by trying out different β values. This experiment revealed that LoTeNet was either under-fitting (β < 4) or over-fitting (4 ≤ β ≤ 8) to the training set. This could be attributed to the fact that the number of parameters grow quadratically with β and these discrete jumps make it harder to obtain models of optimal capacity for some tasks.
Another concern with tensor network based models is their tendency to be over parameterised as the number of parameters scale as {d · N · β 2 }. The linear dependency on the number of pixels, N , for volumetric data results in a model with massive number of parameters at the outset which is further aggravated by the quadratic dependency on β. While in our experiments, this has not had an adverse effect on the performance, we have observed the model is capable of over-fitting on small training data quite easily. This can be a concern when dealing with datasets with few samples.
One way of handling these problems would be to introduce weight-sharing between MPS blocks at each layer so that the number of parameters can be reduced and controlled better. A possible future work is to investigate the effect of sharing MPS blocks per layer and study if equivariance to translation can be induced in tensor network based models (Cohen and Welling, 2016).
Finally, the issue of extending LoTeNet type models beyond classification or regression tasks to more diverse medical imaging tasks such as segmentation and registration are open research questions waiting to be explored further.

Conclusion
We have presented a tensor network model for classification of 2D and 3D medical images, of higher spatial resolutions. Using LoTeNet we have shown that aggregating MPS approximations on small image regions in a hierarchical manner retains global structure of the image data. With experiments on three datasets we have shown such models can yield competitive performance on a variety of classification tasks. We have measured the GPU memory requirement of the proposed model and shown it to be a fraction of the memory utilisation of state-of-the-art CNN based models when trained to attain comparable performance. Finally, we have demonstrated that a single model architecture with β = 5 (tuned on one of the datasets) can perform well on other datasets also, which can be an appealing feature as it reduces dataset specific architecture search.