{'fraction_non_alphanumeric': 0.05282505499594767, 'fraction_numerical': 0.012562232256570568, 'mean_word_length': 4.433941491034917, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Graph-Regularized Tensor Regression: A Domain-Aware Framework for Interpretable Multi-Way Financial Modelling

Analytics of financial data is inherently a Big Data paradigm, as such data are collected over many assets, asset classes, countries, and time periods. This represents a challenge for modern machine learning models, as the number of model parameters needed to process such data grows exponentially with the data dimensions; an effect known as the Curse-of-Dimensionality. Recently, Tensor Decomposition (TD) techniques have shown promising results in reducing the computational costs associated with large-dimensional financial models while achieving comparable performance. However, tensor models are often unable to incorporate the underlying economic domain knowledge. To this end, we develop a novel Graph-Regularized Tensor Regression (GRTR) framework, whereby knowledge about cross-asset relations is incorporated into the model in the form of a graph Laplacian matrix. This is then used as a regularization tool to promote an economically meaningful structure within the model parameters. By virtue of tensor algebra, the proposed framework is shown to be fully interpretable, both coefficient-wise and dimension-wise. The GRTR model is validated in a multi-way financial forecasting setting and compared against competing models, and is shown to achieve improved performance at reduced computational costs. Detailed visualizations are provided to help the reader gain an intuitive understanding of the employed tensor operations.


INTRODUCTION
Financial data are naturally multi-dimensional and multivariate; indeed, a typical financial data setting involves various market data features which span over time and across several assets. It therefore occurs routinely that when considering long time horizons, many assets, various geographies, and asset classes, the dimensionality of each mode can quickly explode. This results in a large number of model parameters which renders the training prone to over-fitting and impedes efficient learning, a phenomenon known as the Curse of Dimensionality [1].
An emerging framework which is suited to deal with the Curse of Dimensionality associated with financial data are tensor-based models [2,3,4]. Tensors are multi-dimensional generalization of vectors (order-1 tensors) and matrices (order-2 tensors). In particular, Tensor Decomposition (TD) techniques, similar to their matrix decomposition counterparts, have the ability to logarithmically compress large multi-dimensional tensors into an efficient low-rank format by exploiting intra-and inter-dependencies across data dimensions (tensor modes) [5,6,7]. This both preserves the underlying data structure and interpretability, while bypassing the Curse-of-Dimensionality, as the associated computational costs reduce from an original exponential one to a linear one in the data dimensions.
In addition, financial data exhibit irregular relational structures, which cannot be handled via traditional machine learning techniques. For instance, pairs of stocks can exhibit varying degrees of correlations, depending on their industries, market capitalizations, and value factors, but to mention a few. Currencies also exhibit varying relational structures depending on the interest rate differential and other macroeconomic factors between pairs of countries. Such irregular relational structures, which are well understood and researched within the quantitative finance community, cannot be straightforwardly integrated into standard machine learning methods, resulting in sub-optimal models.
To this end, Graph Signal Processing techniques have recently gained the interest of the data analytics community, as they generalize traditional machine learning techniques to data on irregular domains [8,9,10,11], achieving state-of-the-art performance in tasks where the underlying data relations can be captured as a graph, including financial modelling [4,12]. In this way, assets can be represented as nodes of a graph, while the edges connecting different nodes may be weighted depending on their relational information.
Financial modelling would therefore greatly benefit from a joint consideration of tensor and graph techniques, to simultaneously overcome the Curse of Dimensionality inherent to financial data while benefiting from domain knowledge about the underlying relational structures. To this end, we introduce a novel Graph-Regularized Tensor Regression (GRTR) framework, which allows for the regression of large multi-dimensional data while efficiently incorporating prior domain knowledge into the model parameter estimation. In addition to improved performance and reduced space complexity costs, we show that such a model is fully interpretable, both in the coefficient-wise sense and in the mode-wise sense. Particular effort has been made to provide intuitive understanding of the employed techniques through detailed visualizations of the underlying operations using Tensor Network (TN) diagrams, to provide readers with an insight into the virtues of tensor algebra when dealing with multiple data dimensions.
The rest of the paper is organized as follows. Related work is first discussed in Section 2. The tensor and graph preliminaries necessary to follow this work are presented in Section 3. Next, we introduce the proposed Graph Regularised Tensor Regression (GRTR) framework, derive the corresponding learning algorithm, and analyse the associated space and time complexities in Section 4. The proposed model is validated and compared against traditional algorithms in a synthetic data based experiment in Section 5, as well as in a real data based financial forecasting experiment in Section 6. The interpretability of the trained model is then discussed in detail in Section 7. Finally, we conclude and summarize our findings in Section 8.

RELATED WORK
The proposed framework can be considered under the tensor regression paradigm [13,14] whereby, contrary to standard tensor regression, graph regularization makes it possible to incorporate domain knowledge within the financial modelling setting. This is different from current graph regularized tensor decomposition models [15,16,17], which have been used to only help preserve local relationships among tensor samples. In other words, whereas such models are concerned with tensor representations of multi-dimensional data, we are concerned with regression tasks. The recently introduced Graph Tensor Networks framework [18,19,4] also employs graphs within the tensor network paradigm, however, it is a blackbox model with graphs serving as an added filter, whereas in this work we employ graph regularization as a mathematical constraint. This allows us to maintain model interpretability throughout the data processing chain, while identifying and promoting meaningful relationships between the regression variables. The proposed GRTR framework is also different from Tensor Decomposition based kernel regression methods discussed in [20,21], as hand-crafted feature mappings are not required. Research on the use of tensors in finance includes the work in [3,22,23,2], which elaborate on the benefits of tensors as an ideal data structure to compactly represent the inherently multi-dimensional financial data. Tensor decomposition were therefore a natural choice to form the backbone of the proposed GRTR framework.

PRELIMINARIES
This section introduces the necessary background for this work. We refer the readers to [5] and [8] for an in-depth overview of tensor algebra and graph signal processing,

Tensors and Basic Tensor Operations
Definitions. A real-valued tensor is a multidimensional array, denoted by a calligraphic font, e.g., X ∈ R I1×···×I N , where N is the order of the tensor (i.e., number of dimensions), and I n (1 ≤ n ≤ N ) is the size (i.e., the dimensionality) of its n th mode. Matrices (denoted by bold capital letters, e.g., X ∈ R I1×I2 ) can be seen as second order tensors (N = 2), vectors are denoted by bold lower-case letters, e.g., x ∈ R I , and can be seen as order-1 tensors (N = 1), and scalars (denoted by lower-case letters, e.g., x ∈ R) are tensors of order N = 0. A specific entry of a tensor X ∈ R I1×···×I N is given by x i1,...,i N ∈ R. The following conventions for basic linear/multilinear operations are employed throughout the paper. Vectorization and matricization. Given an order-N tensor, X ∈ R I1×···×I N , its vectorization reshapes the multidimensional array into a large vector, x ∈ R I1I2···I N . The moden matricization unfolds the tensor into a matrix, X (n) ∈ R In×I1I2···In−1In+1···I N , as illustrated in Figure 2.
Outer product. The outer product of vectors a ∈ R I and b ∈ R J is given by Kronecker product. The Kronecker product of two matrices A ∈ R I×J and B ∈ R K×L is denoted by Khatri-Rao product. The Khatri-Rao product of two matri- Tensor contraction. The contraction of an N th -order tensor, A ∈ R I1×···×I N , and an M th -order tensor B ∈ R J1×···×J M , over the n th and m th modes respectively, where I n = J m , results in an (N + M − 2) th -order tensor with entries c i1,...,in−1,in+1,...,i N ,j1,...,jm−1,jm+1,...,j M = In in=1 a i1,...,in−1,in,in+1,...,i N b j1,...,jm−1,in,jm+1,...,j M . Tensor inner product. The inner product of two N thorder tensors, A ∈ R I1×···×I N and B ∈ R I1×···×I N , results in a scalar, c = A, B ∈ R, such that c = i1,...,i N a i1,...,i N b i1,...,i N . Rank-1 tensor. An N th -order tensor, A ∈ R I1×···×I N , is said to be of rank-R = 1, if it can be expressed as an outer product of N vectors, A = a (1) • · · · • a (N ) = N n=1 a (n) , where a (n) ∈ R In . Rank-1 tensor inner product. The inner product between two N th -order tensors, A ∈ R I1×···×I N and B ∈ R I1×···×I N , can be expressed as a sequence of tensor contractions, . Tensor-by-vector product. A tensor-by-vector product is a special case of tensor contraction between an N thorder tensor, A ∈ R I1×···×I N , and an 1 st -order tensor, b ∈ R In , which results in a (N − 1)-th order tensor, C ∈ R I1×···×In−1×In+1×···×I N , which is exactly 1 order smaller than A. This operation is mathematically equivalent to performing the following steps in a sequence: (i) Matricize the tensor, A (n) unf old ←−−−− A; (ii) Perform a linear combination (matrix-by-vector multiplication) in the n-th sub- . Therefore, the tensor-by-vector product can be interpreted as performing a linear combination in the n-th subspace. Tensor network diagram. For visualization purpose and to gain physical insight into the multi-dimensional nature of tensor operations, we adopt the Tensor Network (TN) notation, whereby a tensor is diagrammatically represented as a circular vertex with outgoing edges; the number of edges equals the order of the tensor (see [5] for more information on this graphical notation), as illustrated in Figure 1. If two vertices are connected by an edge, this connection represents a tensor contraction between two tensors over their shared common mode. Canonical Polyadic Decomposition. Canonical Polyadic Decomposition (CPD) of rank, R, decomposes an order-N tensor, W ∈ R I1×I2×···×I N , into a sum of outer products of Fig. 3: Visualization of the Canonical Polyadic Decomposition (CPD) for an order-3 tensor, W ∈ R I 1 ×I 2 ×I 3 , according to (1). The sizes of tensor modes, I1, I2, and I3, are labelled in gray.

Fig. 4:
Alternative visualization of the Canonical Polyadic Decomposition for an order-3 tensor, W ∈ R I 1 ×I 2 ×I 3 , according to (2) in: (Top) standard matrix notation and (Bottom) TN notation. The sizes of tensor modes, I1, I2, and I3, are labelled in gray.
N vectors, and is given by where u (n) r ∈ R In denotes the r-th mode-n factor vector, as illustrated diagrammatically in Figure 3.
Alternatively, for each mode-n, we can stack the vectors, u (n) r for r = 1, . . . , R, to form the factor matrix, . This allows us to express the CPD as a series of tensor contractions, given by as illustrated in Figure 4 in both standard matrix notation and TN notation.
The mode-n matricization of a tensor in its CPD format, W (n) ∈ R In×I1···In−1In+1···I N , can be expressed as a series of Khatri-Rao products of the factors matrices, as Similarly, the vectorization of a tensor in its CPD format, vec(W) ∈ R I1···I N , can be expressed as where 1 denotes a vector of unities.

Graphs and Signals on Graphs
. . , I and i 2 = 1, . . . , I. A graph can be fully described by the weighted adjacency matrix, Alternatively, a graph can be described in terms of its Laplacian matrix, [8].
A graph signal can be represented by a vector, u ∈ R I , which associates a scalar value to each of its I vertices. The smoothness of such a graph signal with respect to the underlying graph topology can then be quantified as s = u T Lu, Fig. 6: Tensor network diagram of the inner product, W, X , within an order-3 tensor regression. The weight tensor, W ∈ R I 1 ×I 2 ×···×I N , on the left is factorized by the CPD according to equation (2). The connected tensor modes represent tensor contractions along these modes.
which measures the change in signal value between neighbouring (i.e., connected) vertices [9]. Therefore, a graph signal is said to be smooth if it does not change abruptly across connected nodes, as illustrated in Figure 5. For a given set of R different graph signals, the total smoothness over all such signals can be computed as s = r u T r Lu r . Equivalently, the given set of R graph signals can be stacked as a set of column vectors into a matrix, U ∈ R I×R . The smoothness can then be computed as

Tensor Regression
For a set of tensor-structured inputs, X ∈ R I1×I2×···×I N , and a set of corresponding output labels, y ∈ R, the tensor regression [24] can be expressed as Analogously to standard regression, Tensor Regression (TR) aims to learn: (i) the regression constant, b ∈ R, and (ii) the weight tensor, W ∈ R I1×I2×···×I N , in the CPD format, such that the prediction mean-squared-error (MSE) in (6) is minimized for a given dataset. Remark 1. Instead of directly learning the large weight tensor, W, tensor regression learns its low-rank CPD factor matrices, U (n) (one for each of the N input dimensions), such that equation (2) holds. The tensor regression framework in (6) is illustrated in Figure 6 in the TN notation.
The learning of the CPD factor matrices is often achieved through the Alternating Least Squares (ALS) procedure [24], which iteratively updates each factor matrix by minimizing the MSE with respect to mode-n, n = 1, . . . , N , until convergence.
Remark 2. For a learned weight tensor, W, the tensor regression framework in (6) is mathematically equivalent to a classical linear regression, y = w T x + b, where w and x are respectively the vectorized versions of W and X , given by The main difference (and advantage) of the tensor regression framework over a standard multivariate regression approach is that it maintains the multi-dimensional structure of the original data, allowing us to exploit its multi-way, low-rank CPD structure during the training (i.e., parameter estimation) stage. This reduces both the computational complexity and the likelihood of over-fitting, while maintaining interpretability.

Graph Regularized Tensor Regression
In practice, when learning the factor matrices, U (n) ∈ R In×R , in addition to minimizing the MSE, the parameter estimation is often regularized by the L 2 -norm of the factor matrices. However, such a regularization scheme is ignorant to the intramodal relations present in the data. In other words, given some domain knowledge about the relational structure within mode-n (e.g., sector relations among I n stocks), a naive L 2regularization scheme cannot incorporate such domain knowledge into the model training process.
To this end, we propose an extension to the classical tensor regression paradigm to enable the incorporation of domain knowledge into the learning process. Given a graph Laplacian matrix, L (n) ∈ R In×In , which encodes domain knowledge about the relational structure within mode-n, we aid the learning of the mode-n factor matrix, U (n) ∈ R In×R , through a corresponding mode-n graph Laplacian regularization, as discussed in Section 3.2. In this way, domain-knowledge about data modes informs the learning, which both reduces the likelihood of over-fitting, and accelerates the learning process.
Remark 3. The column vectors of the mode-n CPD factor matrix, u (n) r ∈ R In for r = 1, . . . , R, can be viewed as graph signals. In this way, the estimation of the factor matrix can be interpreted as the estimation of R graph signals that are smooth with respect to the graph topology encoded by the mode-n Laplacian matrix, L (n) ∈ R In×In . Example 1. Consider an order-3 input tensor, X ∈ R T ×S×F , where F features for S stocks are indexed over T time-steps. For this data, we wish to learn a tensor regression model, such that the mode-2 (the stocks mode) factor matrix, U (2) ∈ R S×R , incorporates the sector knowledge of the S stocks, as illustrated in Figure 7. To this end, we construct a graph with S vertices representing S different stocks, where the edges connect a pair of stocks if they belong to the same sector. Given this graph, we can interpret the column vectors of the mode-2 factor matrix, u R , as a set of R graph signals. Therefore, the learning of the factor matrix, U (2) , can be regularised by ensuring that its column vectors, u r , r = 1, . . . , R, are smooth with respect to the underlying stock graph (i.e., the graph signal smoothness, s = tr U (2) T L (2) U (2) , is small).
To formally introduce the concept of graph regularization in the classical tensor regression framework, we define a moden graph Laplacian matrix, L (n) , for each of the N tensor modes, and assign a corresponding regularization constant, λ (n) . The mode-wise graph Laplacian regularization terms are then included into the MSE loss function, calculated over m = 1, . . . , M samples, to yield the following loss function to be minimized Regularization term (8) The constant, b ∈ R, that minimizes the loss function, L, is estimated based on the derivative of the loss function with respect to this constant using the chain-rule, resulting in The factor matrices that minimize the loss function, are estimated by the Alternating Least Squares (ALS) framework, which iteratively minimizes the loss function sequentially onemode-at-a-time. To this end, we can apply the properties discussed in Section 3.1 and matricize the tensor contraction with respect to the n-th mode as where U (−n) = 1 i=N,i =n U (i) . This mode-n matricized formulation allows for the computation of the partial derivative with respect to mode-n factor matrix, by leveraging on the well-known matrix derivatives rules Through the chain rule, the above results can be included into the calculation of the derivative of the loss function with R ], can be viewed as a set of R graph signals (its column vectors), which can be learned so as to be smooth with respect to the underlying stocks graph. This makes it possible to regularize the learning of the factor matrix and results, as desired, in similar weights for stocks from the same sector.
respect to the mode-n factor matrix as where (m) = y (m) − W (n) , X

Remark 4. Equation 12
considers the general case where each mode of the input tensor, X ∈ R I1×I2×···×I N , is associated with an underlying graph domain. However, if the n-th mode does not have an associated Laplacian matrix, then the corresponding regularization parameter, λ (n) , can simply be set to 0 to ignore the regularization term in the loss function.
Based on the partial derivatives in (12), a gradient-descent optimization framework can be employed to iteratively update the factor matrices until convergence. The learning algorithm is summarized in Algorithm 1.

Choice of Model Rank
A key hyper-parameter choice within the GRTR model is the model rank. A larger rank increases the expressive power of the model, but it also increases the space complexity of the model which can lead to over-fitting. In practice, the optimal choice of model rank can be determined via grid-search, whereby the rank with the smallest cross-validation error is used for the final model.

Choice of Regularization Parameter
The regularization parameter, λ (n) , is another important hyperparameter. As λ (n) increases, the regularization term begins to dominate the loss function over the error term in (8), which reduces the likelihood of in-sample over-fitting and results in better out-of-sample performance. In practice, the optimal choice of regularization parameter can also be determined empirically via grid search.

Example 2.
To illustrate the effect of the graph regularization parameter, consider a GRTR model of rank R = 1 , whereby the mode-n factor matrix, U (n) ∈ R In×R , is reduced to a column vector, u (n) ∈ R In . This factor vector is effectively a graph signal, which contains a scalar value for each of the I n vertices in the mode-n graph implied by the corresponding Laplacian matrix, L (n) ∈ R In×In . For large values of λ (n) , this leads to highly smooth vectors, u (n) ∈ R In , where the difference between entries, |u j |, is small if the i-th and j-th vertices are closely connected in the underlying mode-n graph. This effect is illustrated in Figure 8.

Complexity Analysis
We now proceed to analyse the space complexity (i.e., the number of trainable parameters) and the time complexity (i.e., the number of basic computations) of the proposed algorithm.
Once the learning is complete, inference can be performed by reconstructing the weight tensor, W, according to (1), which requires a space complexity of O I N , where I n = I for all n. In this way, the time complexity of the inference becomes the same as in standard linear regression, O I N . This is because tensor regression can be re-written in a vectorized format, as shown in (7).
Alternatively, the trained factor matrices, U (1) , . . . , U (N ) , can be kept in a low-rank CPD format, without reconstructing the full weight tensor, W. This reduces the space-complexity from an exponential one, O I N , to a linear one, O (N IR), in the tensor order, N . However, in this way the timecomplexity of inference would increase to O RI N , as this requires evaluation of a series of tensor contractions

EXPERIMENTS ON SYNTHETIC DATA
To validate the ability of the proposed model to capture interand intra-modal relations in data, as opposed to classical models, we conducted an experiment on synthetic data involving large and noisy tensor data. For reproducibility, all experiment code can be found on GitHub 1 .

Experiment Setting
Consider a regression problem which is fully determined by an unknown CPD weight tensor, W ∈ R I1×···×I N . This system generates a set of M observations, X (m) , y (m) , such that y (m) = W, X (m) + η (m) , where η is additive random noise. Given this set of noisy observations, the goal of the experiment is to estimate the ground-truth weight tensor,W ∈ R I1×···×I N , using regression methods, such that , is as small as possible.

Data
To carry out the proposed regression experiment, N = 4 ground-truth factor matrices, U (n) ∈ R 10×5 , were sampled randomly from a uniform distribution with rank R = 5 and dimensions I n = 10, n = 1, . . . , 4. The factor matrices were then contracted according to (2) to generate the groundtruth weight tensor of order-4, W ∈ R 10×10×10×10 . A set of M = 125 input tensors, X (m) ∈ R 10×10×10×10 , were then sampled randomly from a Gaussian distribution with zeromean and unit-variance, which were then contracted with the weight tensor to generate a corresponding set of regression labels, y (m) = W, X (m) + η (m) , where η (m) is additive noise sampled from a Gaussian distribution with zero-mean and variance, σ 2 η , such that σ η = 0.5σ y . The M = 125 observations were then split according to a standard 80% − 20% train-test split.
Remark 5. The input tensors, X (m) ∈ R 10×10×10×10 , are order-4 tensors with modal dimensions I n = 10, giving effectively 10 4 = 10, 000 features for every sample. However, only 100 samples were available for the training set, which is 2 orders of magnitude smaller than the number of features. The Curse-of-Dimensionality inherent to this experiment, coupled with a low signal-to-noise ratio of σy ση = 2, is ubiquitous in financial Big Data. This poses significant challenges to traditional models, both in terms of complexity and over-fitting.

Incorporation of Domain Knowledge through Graphs
The proposed GRTR framework allows for the incorporation of domain knowledge about the underlying ground-truth weight tensor. This is achieved by defining a mode-n graph which captures the corresponding intra-modal structure. For the proposed experiment, the mode-n graph adjacency matrices, A (n) ∈ R 10×10 , n = 1, . . . , 4, were computed such that where u (n) i ∈ R 5 and u (n) j ∈ R 5 denote respectively the ith and j-th row vector of the mode-n factor matrix, U (n) ∈ R 10×5 . In this way, the entries of the mode-n graph adjacency matrix encode the similarity between row vectors of the ground-truth mode-n factor matrix. The mode-n Laplacian matrices, L (n) ∈ R 10×10 , were then computed as discussed in Section 3.2, which were used to regularize the learning process and to guide the estimated weights closer to the ground-truth.

Results
Simulation results are illustrated in Table 1. Linear models, such as LR and L2LR, needed to estimate 10,000 parameters from 100 noisy samples, which led to over-fitting in the training set due to the noisy and large-dimensional nature of data. This resulted in negative explained variance scores in the test set and large mean-squared-errors in the estimation of the ground-truth weight tensor. Classical tensor models, such as TR and L2TR, were able to better capture the inherent multiway relations in data by virtue of their tensor formulation, which reduced the space complexity to only 200 parameters. However, they still struggled due to large noise and small number of observations. The proposed GRTR, aided by the domain knowledge encoded in the mode-n graphs, achieved the best overall scores by a large margin, while incurring the same small complexity cost of other tensor models.

EXPERIMENTS ON REAL DATA
We next validate the proposed model in a real-world multi-way financial forecasting experiment.

Data
Stocks within the S&P500 index between 2010 and 2020 were considered. Some stocks were removed due to data quality issues, leaving 450 stocks in total for the experiment. For each of these stocks, 6 daily market data features were obtained: (i) adjusted closing price, (ii) closing price, (iii) highest price, (iv) lowest price, (v) opening price, and (vi) trading volume. The raw market data were then processed to ensure stationary features by computing the log-returns, r t = log(p t ) − log(p t−1 ), where p t is the value of the market data at the t-th time-step. Using a rolling window of 5 time-steps, a set of tensor-structured inputs was generated at each point in time, X (t) ∈ R T ×S×F , where F = 6 features for S = 450 stocks were extracted over T = 5 past time-steps. For each of the order-3 input tensors at a particular time-step, the corresponding output label was set to be the return of the S&P500 index at the next time-step. The data was then partitioned using a standard 50%-30%-20% train-validation-test split, and standardized to have zero mean and unit variance within the training set.

Domain Knowledge through Graphs
The performance of the proposed GRTR model will vary significantly depending on the design of the graph Laplacian matrix. Designing the optimal graph for equity returns forecasting is beyond the scope of the paper and is left to the domain experts. Instead, a sector graph was used in this experiment for illustrative purposes.
Using the sector information from mode-2 (i.e., the stocks mode), a graph with S vertices representing S different stocks was constructed, where the edges connect a pair of stocks if they belong to the same sector. We then computed the corresponding Laplacian matrix for the mode-2, L (2) ∈ R S×S , as discussed in Section 3.2. The proposed graph Laplacian matrix encodes knowledge about the sector information of the mode-2 stocks, which is expected to improve the learning process, as it guides the learning of the mode-2 factor matrix to be consistent with the correlation structure exhibited by the stock market.

Models
As in Section 5, the proposed GRTR was evaluated against: (i) Linear Regression (LR), (ii) Linear Regression with L2 regularization (L2LR), (iii) standard Tensor Regression (TR), and (iv) Tensor Regression with L2 regularization (L2TR). For the LR and L2LR models, the ground-truth weight tensor was vectorized as w = vec (W) ∈ R 13,500 .  Table 2: Experimental results for the considered models. The first and second rows show respectively the TRain set ACCuracy (TR-ACC) and TEst set ACCuracy (TE-ACC) achieved by the considered models. The final row shows the required number of parameters (i.e., space complexity). The proposed GRTR model achieved the best results while requiring only 461 parameters, which is closely followed by the TR and the L2TR models. Linear models once again over-fitted the training set due to the noisy and large-dimensional nature of the problem.

Results
The results for the proposed experiment are reported in Table 2, which validates the proposed GRTR framework, as it obtained the best scores both in terms of computational complexity costs and out-of-sample accuracy scores.

INTERPRETABILITY
By virtue of multi-linear algebra, the GRTR model offers enhanced interpretability, both in terms of individual regression coefficients and different regression dimensions, as elaborated below.

Coefficient-Wise Interpretability
Similar to the classical linear regression paradigm, the proposed GRTR model allows us to recover the exact coefficient responsible for modelling a given input variable. For instance, the input value, x i1,i2,...,i N within X ∈ R I1×I2×···×I N , is directly linked to the coefficient, w i1,i2,...,i N of the weight tensor, W ∈ R I1×I2×···×I N , which can be reconstructed from its factor matrices according to (2).
Example 3. Consider the order-3 input tensor from the experiment in Section 6, X ∈ R T ×S×F . The effect of the input variable, x 1,10,2 (i.e., the f = 2-nd feature of the s = 10-th stock at the t = 1-st past time-step), on the output, is directly linked to the corresponding entry in the weight tensor, w 1,10,2 .
In the proposed real data experiment, this value was equal to -0.00017.

Mode-Wise Interpretability
By construction, the large weight tensor, W, is decomposed by mode-wise factor matrices, U (n) , which effectively breaksdown the large multi-dimensional regression problem into a set of smaller, physically meaningful, mode-wise sub-problems.  Fig. 9: Illustration of a rank R = 1 order N = 3 tensor regression operation. The data tensor, X ∈ R T ×S×F , contracts with the factor vectors, u (n) 1 , n = 1, 2, 3, to generate the overall scalar output, y = X × 1 1 u 1 . Each mode-n contraction represents a linear combination in the n-th subspace of X ∈ R T ×S×F , which reduces the order of the tensor until a scalar is obtained. In the above example, the original order-3 tensor is reduced to a matrix through a tensor-by-vector contraction (defined in Section 3), which is then reduced to a vector through a matrix-by-vector multiplication, which is then reduced to a scalar through a vector-by-vector inner product.
This allows us to interpret the obtained model coefficients from a mode-wise point of view.
To break down the tensor regression into a series of interpretable mode-wise operations, the tensor regression equation can be rewritten as a sum of tensor contractions by applying the rank-1 tensor inner product property discussed in Section 3.1, as Note that the computation of the tensor-by-vector contraction, Z = X × 1 n u (n) r , can be interpreted as performing a simple linear combination in the n-th subspace of X ∈ R I1×I2×···×I N , as discussed in Section 3.1. Therefore, the sequence of contractions in equation (15), r , can be interpreted as performing a sequence of linear combinations over every n-th subspace of X ∈ R I1×I2×···×I N .
Example 4. Consider the order-3 data tensor from the experiment in Section 6, X ∈ R T ×S×F . For R = 1, equation (15) simplifies to one single contraction, y = X × 1 1 u + b, which can be visualized as in Figure 9. Therefore, tensor regression performs three linear combinations, one in each mode. The time weights, u (1) r ∈ R T , the stocks weights, u (2) r ∈ R S , and the feature weights, u Mode-3 (Features) Fig. 10: Estimated rank R = 1 tensor-regression factor vectors from the proposed experiment. A factor vector is estimated for the time mode, u (1) ∈ R T , the stock mode, u (2) ∈ R S , and feature mode, u (3) ∈ R F . Note that the factor vector values for the time-mode can be interpreted as a weighted moving average operation in the time-subspace, as larger weighting are placed towards more recent data (t = 1) compared to older data (t = 5).
features to the final regression output. For illustration, Figure  10 shows the estimated weights of the GRTR model with rank R = 1. In particular, the time-domain weights can be interpreted as a weighted moving average operation over the time-dimension, with more weighting placed towards recent data (time-step t = 1) compared to older data (time-step t = 5).

CONCLUSION
We have introduced a tensor regression framework with graph-based regularization, which makes it possible to efficiently incorporate domain knowledge into large-scale, multi-dimensional regression through Laplacian graph regularization. The so introduced model has been shown to be physically meaningful and to provide improved results against standard vector and tensor regression models, both in terms of out-of-sample performance and computational costs. By virtue of tensor algebra, the proposed model is also shown to enable enhanced interpretability and control over the learning process. The proposed framework is general and opens avenues for exploring numerous possibilities for formulating economically meaningful graphs from cross-asset variables; these can also can guide learning models towards higher performance and meaningful predictions.