Accelerated Learning of High Dimensional Functions with a Tensor-Featured Training Network

# Accelerated Learning of High Dimensional Functions with a Tensor-Featured Training Network
Source: [https://arxiv.org/html/2608.10351](https://arxiv.org/html/2608.10351)
Yuehaw Khoo∗Haizhao Yang††The University of Maryland, College Park, College Park MD 20742 USA∗The University of Chicago, Chicago, IL 60637 USA

###### Abstract

In this work we present a method to accelerate the optimization of learning high dimensional functions using deep neural network \(DNN\)\. This optimization procedure introduces contextual features into the first layer of a DNN\. The parameters of DNN are optimized via standard gradient descent while keeping the input\-feature basis fixed\. After optimization of the DNN parameters, the feature layer is provided a chance to update and change before DNN optimization resumes\. The feature layer has two types of functions: those that can be evaluated quickly in a matrix\-free way on the domain \(i\.e\. rank\-1 features\) and more complex features that must first be decomposed using tensor network \(TN\) decomposition strategies \(tensor features\)\. In particular, we study the effect of adding features which distill pretrained DNN into TNs using a discretize and decompose strategy\. To efficiently decompose high\-dimensional functions constructed from discretized DNN, we leverage a randomized tensor decomposition strategy\. Using randomization, we are able to reduce the storage cost of decomposing high dimensional functions by at least 8 orders of magnitude\. Using this approach, we are able to efficiently train models between 5 and 40 dimensions\.

## IIntroduction

Methods to interpolate continuous functions and solve partial differential equations \(PDEs\) computationally are an interesting and on\-going challenge in scientific computing\. Canonical methods to represent continuous functions have require users to construct nontrivial grid\-based meshes or introduce complicated numerical techniques\. Not only are these strategies complicated from an implementation standpoint, but they can also require large amounts of computational resources\. Grid\-based strategies, in particular, suffer from the curse of dimensionality which states that the cost of storing and accessing elements of a multidimensional array grows exponentially with the number of dimensions in the array\. As a means to break the curse of dimensionality and solve large and high\-dimensional problems scientific machine learning \(ML\) has recently emerged\.\[[12](https://arxiv.org/html/2608.10351#bib.bib1),[18](https://arxiv.org/html/2608.10351#bib.bib2)\]

Scientific ML is a powerful set of tools which may be used to compactly represent high\-dimensional and unknown functions\. These tools work by learning the underlying structure of functions using a non\-linear network of tuneable parameters\. Because the scientific ML models do not attempt to memorize values of a function at specific points on a mesh, these models are commonly called mesh\-free methods\. By construction, these models are extremely expressive and the number of optimizable components in the models is not strongly correlated to the dimensionality of the function\. Though scientific ML models are exceptionally good at representing a diverse range of functions, the optimization of a ML model can be an extremely inefficient process, requiring hundreds of compute hours even on today’s high\-throughput GPU processing units\. This is especially true for deep neural network \(DNN\) optimizations, like physics informed NNs \(PINNs\), as their optimization is non\-convex\.\[[8](https://arxiv.org/html/2608.10351#bib.bib3),[23](https://arxiv.org/html/2608.10351#bib.bib4)\]In general, training ML models become more difficult and slowly converging with increasing function dimensionality\. There have been many ideas to accelerate the optimization process, some examples include sampling techniques to target high\-residual sub\-domains in the training process\[[5](https://arxiv.org/html/2608.10351#bib.bib6)\]and methods to construct random features using the kernel method\.\[[15](https://arxiv.org/html/2608.10351#bib.bib7)\]

In this work, we are interested in using ideas from tensor\-networks \(TN\) methods to improve the efficiency of training DNNs\. TNs are becoming increasingly popular in scientific computing as methods to efficiently represent high\-dimensional functions and data\.\[[20](https://arxiv.org/html/2608.10351#bib.bib5)\]Effectively, tensor network methods seek to factorize high\-dimensional functions into a sums of outer products of low\-dimensional component functions\. Using this factorization technique, tensor network methods are able to break the curse of dimensionality\. Unlike the mesh\-free ML models, TN methods do require a collocation grid\. As a result, while the exact storage complexity of the method depends on the topology of the TN, in general there is a strong correlation between cost of storing a TN, the dimension of the function being decomposed and number of collocation points per dimension\. How well a TN can represent a function explicitly depends on the type of function and the topology of the TN\. Though these properties make TN in general more difficult to apply to problems than ML models, TN methods do shine in their ability to be efficiently optimized\. This is because most tensor\-network algorithms rely on linear algebra techniques such as the singular value decomposition whereas machine learning models rely on non\-convex gradient descent techniques\. Therefore, in this work we seek to determine methods to apply the beneficial properties of TN methods, i\.e\. reliable optimization and interpretability, to improve and accelerate the optimization of DNN\.

Our idea is similar to a warm\-start initialization technique and works in the following way\. We generate a collection of functions and apply these functions to input data to a DNN\. We add the output of these functions to the input layer of the DNN and optimize parameters of the DNN subject to the input data and these new features\. A most straight\-forward feature would be a previously optimized DNN\. We demonstrate in this work, why this idea fails and instead use a TN decomposition strategy to improve on this ides\. Furthermore, we identify two distinct subclasses of basis functions, those that can be applied quickly in a matrix\-free way and those that can be approximated via a TN decomposition\. We, therefore, introduce a 2\-step optimization process for DNNs where, given a set of input features, DNN parameters are optimized recursively by canonical gradient descent techniques for a fixed number of iterations\. After parameter optimization, the basis functions are evaluated and updated to improve convergence\. After the basis functions are modified, a new DNN is initialized and parameter training is resumed\. This work draws parallels with random feature training as the DNN now could be considered the combination of it’s input layer and neuron layers\. Like random feature training, we fix the parameters associated with input level features during DNN training\. One major difference is that the interaction of features is facilitated implicitly by the DNN parameter layers\. This gives the feature layer of the DNN an outer\-product structure similar to a rank\-1 TN decomposition\. Therefore, we call this techniqueTensor\-Featuredtraining\.

As pointed out previously, we will introduce DNNs which are discretized and decomposed on their domain\. To do this decomposition we will utilize the canonical polyadic decomposition \(CPD\)\.\[[7](https://arxiv.org/html/2608.10351#bib.bib8),[3](https://arxiv.org/html/2608.10351#bib.bib10),[6](https://arxiv.org/html/2608.10351#bib.bib9)\]We limit this work to a single tensor decomposition strategy to simplify the analysis of the procedure\. In practice any flavor or combination of tensor decompositions could be enabled\. Furthermore, one can contextualize the tensor decompositions, like the CPD and the tensor train decomposition, as method for computing a high dimensional Gaussian process using a kernel matrix which is constructed as a tensor product of one\-dimensional kernel matrices\.\[[25](https://arxiv.org/html/2608.10351#bib.bib15),[27](https://arxiv.org/html/2608.10351#bib.bib16)\]As a result, we hope to demonstrate the following: first that the introduction of tensor features acts as a regularizer in the optimization of DNN and; second that the tensor decomposition acts as a smoothing process for the DNN via non\-local features introduced to the tensor kernel matrix\. Because tensor decompositions can be extremely expensive in both computer time and memory, we will leverage advancements in randomized linear algebra to efficiently decompose high dimensional functions\.

The remainder of the text is laid out in the following manner\. In[SectionII](https://arxiv.org/html/2608.10351#S2), we will discuss the theoretical background of the work including the randomized canonical polyadic decomposition and the relationship between tensor decompositions and the gaussian process for high\-dimensional functions\. In[SectionIII](https://arxiv.org/html/2608.10351#S3)we describe the outlines of our experiments to learn functions which are solutions to PDEs in various degrees of dimensionality spanning from55to4040\. In[SectionIV](https://arxiv.org/html/2608.10351#S4)we provide results for our optimization experiments, compare our result to canonical optimization strategies and discuss methods to improve the optimization process for DNN to learn unknown functions\. Finally in[SectionV](https://arxiv.org/html/2608.10351#S5)we summarize the findings and discuss future directions and improvements for this work\.

## IITheoretical Background

### II\.1The Canonical Polyadic Decomposition

The CPD is a higher\-order extension of the singular value decomposition and represents a function ofddvariables \(expressed as an order\-ddtensor\) as a sum of outer products of single variables, i\.e\. given𝒯∈ℝI1,I2,…,Id\\mathcal\{T\}\\in\\mathbb\{R\}^\{I\_\{1\},I\_\{2\},\\dots,I\_\{d\}\}

ta,b,…,n​=CPD​∑i=1RCPλi​𝐚i∘𝐛i∘⋯∘𝐧i\\displaystyle t\_\{a,b,\\dots,n\}\\overset\{\\mathrm\{CPD\}\}\{=\}\\sum\_\{i=1\}^\{R\_\{\\mathrm\{CP\}\}\}\\lambda\_\{i\}\{\\bf a\}\_\{i\}\\circ\{\\bf b\}\_\{i\}\\circ\\dots\\circ\{\\bf n\}\_\{i\}\(1\)where∘\\circdefines the vector outer product,𝐚i∈ℝI1,𝐛i∈ℝI1,…,𝐧i∈ℝId\{\\bf a\}\_\{i\}\\in\\mathbb\{R\}^\{I\_\{1\}\},\{\\bf b\}\_\{i\}\\in\\mathbb\{R\}^\{I\_\{1\}\},\\dots,\{\\bf n\}\_\{i\}\\in\\mathbb\{R\}^\{I\_\{d\}\},λi\\lambda\_\{i\}is a scaling constant that allows the vectors\[ai,…​ni\]\[a\_\{i\},\\dots n\_\{i\}\]be unit normalized, andRCPR\_\{\\mathrm\{CP\}\}is the CP rank, i\.e\. the smallest set of vectors such that[Eq\.1](https://arxiv.org/html/2608.10351#S2.E1)satisfies the equality\. All the vectors which span the same mode of the tensor can be group into a factor matrix, i\.e\.𝐀=\[𝐚1,…,𝐚RCP\]∈ℝIa×RCP\{\\bf A\}=\[\{\\bf a\}\_\{1\},\\dots,\{\\bf a\}\_\{R\_\{\\mathrm\{CP\}\}\}\]\\in\\mathbb\{R\}^\{I\_\{a\}\\times R\_\{\\mathrm\{CP\}\}\}\. Therefore, the CPD reduces the storage complexity of add\-dimensional function fromNdN^\{d\}whereN≈I1≈I2,…N\\approx I\_\{1\}\\approx I\_\{2\},\.\.\.toN​D​RCPNDR\_\{\\mathrm\{CP\}\}\. Although finding the exact CP rank is formally anNP\-hard problem, it is straightforward and easy to discover an approximate CPD of𝒯\\mathcal\{T\}given a desired rankRR\. By construction, the CPD defines a set of 1\-dimensional basis functions such that a function ofddinputs as a separable set of products of functions of single inputs\.

The most standard way to optimize the rank\-RRCPD is via the alternating least squares \(ALS\) optimization\.\[[13](https://arxiv.org/html/2608.10351#bib.bib11),[2](https://arxiv.org/html/2608.10351#bib.bib12)\]The ALS method splits the optimization of the CPD intodd\-linear least square problems\. Unfortunately, the CPD\-ALS optimization is non\-convex\. However, based on our experience and results, we find that the non\-convex nature of the CPD\-ALS has minimal effect on discovering accurate decomposition’s of tensor for most cases in physics and chemistry\.\[[16](https://arxiv.org/html/2608.10351#bib.bib17)\]As we will demonstrate in[SectionIV](https://arxiv.org/html/2608.10351#S4), this is also true for this work\. A more serious issue with the CPD\-ALS is that the optimization suffers from the curse of dimensionality, i\.e\. the computational cost of the method scales exponentially with the dimension of the function being decomposed\.

In general, the task of the CPD\-ALS is to iteratively solve a collection of extremely overdetermined least squares problems\. Because these problems are over\-determined, a direct result of the matricization of higher order tensors along a single mode of the tensor, the CPD\-ALS a perfect candidate for randomized linear algebra\. There are two main ways to apply randomization to the CPD\-ALS, via sampled least squares\[[1](https://arxiv.org/html/2608.10351#bib.bib14),[14](https://arxiv.org/html/2608.10351#bib.bib13)\]or via sketched and sampled least squares\.\[[26](https://arxiv.org/html/2608.10351#bib.bib18),[21](https://arxiv.org/html/2608.10351#bib.bib19)\]In this work, we will only consider the sampled least squares method because of its low\-cost and structure preserving properties\. In the literature, it has been demonstrated that the number of samples required for the CPD randomized ALS \(CPD\-RALS\) is proportional to the rank of the CPD and has no direct dependence on the dimensionality of the function being decomposed\.\[[1](https://arxiv.org/html/2608.10351#bib.bib14),[4](https://arxiv.org/html/2608.10351#bib.bib23)\]Therefore, in this work, we will assume that the CP rank of every tensor is practically small compared to the number of interpolation points along a single dimension and does not grow \(significantly\) with increasing tensor dimension\. The RALS procedure therefore reduces the computational complexity of the CPD from exponential in function dimension to polynomial and independent of function dimension\. Using the CPD\-RALS, the most time intensive component of the algorithm is the cost of sampling the target function which grows asd​S​NdSNwhereddis the dimension of the function,SSis the number of samples per mode andNNis the number of interpolation points along a single dimension\. To minimize the overall cost of this step, we will assume that single set of samples gathered from a uniform distribution is sufficient for optimizing the CPD\-RALS\.

In this work we will not explicitly derive the sampled CPD\-RALS method\. We utilize a modified version of the canonical version developed by Battaglino et\. al\. where samples are determined via the leverage score distribution of each factor matrix\. Our method deviates from the canonical method in the following way\. We found that the magnitude of elements in the sampled Khatri\-Rao product \(KRP\) associated CPD\-RALS decay’s rapidly with dimensionality of the tensor\. To accommodate this issue, we introduce an intermediate normalization where rows of the KRP are normalized to 1 such that

A​x\\displaystyle Ax=B\\displaystyle=B\(2\)\(l​A^\)​\(x​l−1\)\\displaystyle\(l\\hat\{A\}\)\(xl^\{\-1\}\)=B\\displaystyle=BA^​x¯\\displaystyle\\hat\{A\}\\bar\{x\}=B\.\\displaystyle=B\.The intermediate normalization tensorllcan, without loss of generality, be factored into the column\-wise normalization ofx¯\\bar\{x\}and, therefore, absorbed intoλ\\lambdain[Eq\.1](https://arxiv.org/html/2608.10351#S2.E1)\. We have found that this process, in general, stabilizes the LS solver in the CPD\-RALS procedure\.

### II\.2The CPD\-ALS as Gaussian Processes

Recently, it has been shown that there is a connection between popular tensor decomposition strategies, kernel machines\[[22](https://arxiv.org/html/2608.10351#bib.bib25)\]and Gaussian process\(GP\)\.\[[24](https://arxiv.org/html/2608.10351#bib.bib24)\]The GP is a mechanism for modeling functions and allows the user to incorporate prior knowledge of the function \(i\.e\. random inputs on the domain and associated noisy observations\)\. GP are fully described by two components, a meanμ∈ℝ\\mu\\in\\mathbb\{R\}which is usually chosen to be zero and a covariance or kernel matrixk∈ℝΩi×Ωik\\in\\mathbb\{R\}^\{\\Omega\_\{i\}\\times\\Omega\_\{i\}\}whereΩi\\Omega\_\{i\}is the extent of the domain in theiith dimension\. Given a set of points and their associated observations, one can determine random functions on the domain by sampling the multivariate gaussian associated with this mean and kernel matrix\. The kernel matrix is commonly approximated using a finite number of basis functions\[[17](https://arxiv.org/html/2608.10351#bib.bib26),[19](https://arxiv.org/html/2608.10351#bib.bib27)\]

k​\(x,x′\)=φ​\(x\)T​Λ​φ​\(x′\)\\displaystyle k\(x,x^\{\\prime\}\)=\\varphi\(x\)^\{T\}\\Lambda\\varphi\(x^\{\\prime\}\)\(3\)whereφ​\(x\)∈ℝΩi×M\\varphi\(x\)\\in\\mathbb\{R\}^\{\\Omega\_\{i\}\\times M\}are basis functions andΛ∈ℝM×M\\Lambda\\in\\mathbb\{R\}^\{M\\times M\}are basis function weights\. One can recognize that each CPD factor matrix are explicitly constructed as a collection ofRR, 1\-dimensional basis functions for a on set of observed points on a predefined collocation grid\.

The CDP\-ALS iteratively constructs each 1\-dimensional basis matrix \(CPD factor matrix\) in the following way\. For convenience we limit this explanation to a three dimensional function𝒯∈ℝIa×Ib×Ic\\mathcal\{T\}\\in\\mathbb\{R\}^\{I\_\{a\}\\times I\_\{b\}\\times I\_\{c\}\}but this idea can be extended to arbitrary dimension without loss of generality\. First, the loss function is defined for theIaI\_\{a\}th mode

f​\(A\)=‖Ta−A∗​\[B⊗C\]T‖22\\displaystyle f\(A\)=\\\|T\_\{a\}\-A^\{\*\}\[B\\otimes C\]^\{T\}\\\|^\{2\}\_\{2\}\(4\)whereTaT\_\{a\}is the so called matriciation of the third order tensor𝒯\\mathcal\{T\}that permutes and reshapes the tensor to beTa∈ℝIa×Ib​IcT\_\{a\}\\in\\mathbb\{R\}^\{I\_\{a\}\\times I\_\{b\}I\_\{c\}\}, a short and long matrix\. Also⊗\\otimesdefines the Khatri\-Rao tensor product \(KRP\) which is the column\-wise Kronecker product of two tensors\. The KRP defines the following productWr=br∘crW\_\{r\}=b\_\{r\}\\circ c\_\{r\}where is a matrixWr∈ℝIb×IcW\_\{r\}\\in\\mathbb\{R\}^\{I\_\{b\}\\times I\_\{c\}\}for each value ofr∈ℝRr\\in\\mathbb\{R\}^\{R\}\. By definitionBBandCCare also collections of 1D basis which were either chosen at random or drawn randomly a different GP\. Next, one solves forA∗A^\{\*\}by taking the gradient of the loss function and solving the following least squares problem

Ta​\(B⊗C\)=A∗​\(B⊗C\)T​\(B⊗C\)\.\\displaystyle T\_\{a\}\(B\\otimes C\)=A^\{\*\}\(B\\otimes C\)^\{T\}\(B\\otimes C\)\.\(5\)This equation consists of two parts:Ta​\(B⊗C\)∈ℝIa×RT\_\{a\}\(B\\otimes C\)\\in\\mathbb\{R\}^\{I\_\{a\}\\times R\}which down\-folds the high dimensional observations of the function into a single dimensional set of mean values for each basis basis function by taking the inner product of the higher order tensor and all other degrees of freedom\. The second part is the construction of the squared KRP,\(B⊗C\)T​\(B⊗C\)\(B\\otimes C\)^\{T\}\(B\\otimes C\)\. Using Tensor Algebra, this expression can be rewritten as\(BT​B\)∗\(CT​C\)\(B^\{T\}B\)\*\(C^\{T\}C\)where∗\*defines the Hadamard or element\-wise product of two matrices\. Written in this way one can see thatBT​B∈ℝR×RB^\{T\}B\\in\\mathbb\{R\}^\{R\\times R\}andCT​C∈ℝR×RC^\{T\}C\\in\\mathbb\{R\}^\{R\\times R\}are both approximations to 1\-dimensional kernel matrices\. Therefore the squared KRP can be viewed as a product kernel\[[19](https://arxiv.org/html/2608.10351#bib.bib27)\]of many 1\-dimensional kernel matrices\. As a result, the above solution to the CPD\-ALS can be viewed as a method that leverages GP to iteratively draw and refine 1\-dimensional basis functions for high dimensional functions\. By leveraging the basis function approximation and kernel product structure, the CPD\-ALS is able to draw GP basis functions with information across the entire domain of the function with a relatively low cost\. Unfortunately, one can see that the construction of the left hand side of[Eq\.5](https://arxiv.org/html/2608.10351#S2.E5)suffers from the curse of dimensionality\.

In an effort to minimize the cost of the CPD\-ALS, the CPD\-RALS modifies the above process\. The gradient of the CPD\-RALS loss function for theIaI\_\{a\}th mode can be written as

\(Ta​Sa\)​\[SaT​\(B⊗C\)\]=A∗​\[SaT​\(B⊗C\)\]T​\[SaT​\(B⊗C\)\]\\displaystyle\(T\_\{a\}S\_\{a\}\)\[S^\{T\}\_\{a\}\(B\\otimes C\)\]=A^\{\*\}\[S^\{T\}\_\{a\}\(B\\otimes C\)\]^\{T\}\[S^\{T\}\_\{a\}\(B\\otimes C\)\]\(6\)whereSa∈ℝIb​Ic×IsS\_\{a\}\\in\\mathbb\{R\}^\{I\_\{b\}I\_\{c\}\\times I\_\{s\}\}is a sparse, column\-wise selection matrix chosen specifically for theIaI\_\{a\}th subproblem\. By constructing the randomized solver in this way we may significantly reduce the cost of inner\-product down\-folding process fromIa​Ib​Ic​RI\_\{a\}I\_\{b\}I\_\{c\}RtoIa​Is​RI\_\{a\}I\_\{s\}R, provided thatIs∝RI\_\{s\}\\propto RandR≪Ib​IcR\\ll I\_\{b\}I\_\{c\}\.\[[1](https://arxiv.org/html/2608.10351#bib.bib14),[4](https://arxiv.org/html/2608.10351#bib.bib23)\]However, this sampling process changes the form of our kernel matrix The sampling matrix first down\-folds the high\-dimensional space orthogonal to theIaI\_\{a\}into a relatively small set of observations\. In this work, we choose these observations by first constructing a leverage score distribution of the kernel basis function\. By taking advantage of the KRP structure of the kernel basis functions, we may randomly sample the leverage score distribution of each basis matrix independently\. The row\-wise leverage score distribution of each factor matrix measures how far away observations on each quadrature point are from all the other observations and thus the leverage scores are correlated with the degree of misprediction of a least squares model\. Therefore, by importance sampling the KRP using the leverage score distribution, we are attempting to determine a compact basis of observations which contribute most strongly to the basis\-approximated kernel matrix\. After choosing this compact set of basis points, we now have a single set of kernel basis functions\. In this process you can recognize that our normalization process in[Eq\.2](https://arxiv.org/html/2608.10351#S2.E2)is effectively a modification of our basis function weights,Λ\\Lambda\.

## IIIComputational Experiments

As this work serves as a preliminary study into this improved optimization process, we make a number of restrictions to simplify and limit the analysis\. First we only consider functions of dimension,d∈\[5,40\]d\\in\[5,40\]\. Unless otherwise stated, we will use the following parameters by default in this manuscript\. In all cases we consider a DNN to have 3 layers and with 100 features per layer\. We apply the following activation functions to each layer\[relu,relu,identity\]\[\\mathrm\{relu\},\\mathrm\{relu\},\\mathrm\{identity\}\]of the DNN\. To optimize the parameters of the DNN we utilize an ADAM gradient descent algorithm\. We stop the parameter update of each DNN after 3000 epoch\. We consider one function that is time independent and one function which contain a time dependence\. For our test and training sets, we use a spatial domain of \(\-1,1\) and a time domain of \(0,1\) and determine points on this domain using a regular grid of 10,000 points per mode\. The training set contains 1000 random points and the validation set contains 5000 random points\. For this study a single batch of training points is drawn for all DNN optimizations and this set is not updated\. This decision comes from an analysis of the structure of the randomized CPD\-ALS optimization where we find no significant outliers in leverage scores across any single spatial dimension for all tested functions\. This decision was validated by decomposing each known functions via the CPD\-RALS procedure using a single uniformly drawn random sample\. Resampling may become necessary in future works when considering more difficult learning settings such as increased dimensionality\.

For the CPD, unless stated otherwise, we decompose DNN’s on a grid of 500 points per mode and to a rank of3030\. For the CPD\-RALS procedure, we fix the number of samples to30003000per least squares subproblem\. This reduces the size of target tensor in each least squares subproblem from500×500d−1500\\times 500^\{d\-1\}to500×3000500\\times 3000\. In the smallest tested case, this is a 20 million times reduction in storage cost of the CPD optimization process\. The CPD was optimized using 100 ALS iterations; however, early stopping conditions could be successfully applied to reduce the cost of the decomposition\. When necessary, to evaluate points which are not explicitly on the quadrature grid of the CPD, a simple nearest neighbor interpolation strategy is applied\.

The experiments are primarily run on a Mac M1 Max equipped with 8 high\-throughput CPU cores each with a 3\.2GHz processing speed and 64 GB of unified RAM\. DNN were constructed and optimized using the Julia\-based Flux\.jl\[[9](https://arxiv.org/html/2608.10351#bib.bib20),[10](https://arxiv.org/html/2608.10351#bib.bib21)\]library\. CPD were constructed and optimized using the Julia\-based ITensorCPD\.jl library\[[11](https://arxiv.org/html/2608.10351#bib.bib22)\]

## IVResults and Discussion

### IV\.1Iterative Tensor\-Featured Training via Tensor Decomposed DNN\.

First we investigate the effect of using a pretrained DNN and as a tensor feature\. For this section we attempt to learn a function which is the solution to the the nonlinear elliptic equation inside a bounded domain,

u​\(x\)=sin⁡\(π2​\(1−\|x\|2\)2\.5\),\\displaystyle u\(x\)=\\sin\(\\frac\{\\pi\}\{2\}\(1\-\|x\|\_\{2\}\)^\{2\.5\}\),\(7\)using the DNN structure described previously\. We consider the problem in three dimensionsd=5,10,d=5,10,and1515\. In this experiment, we first use canonical optimization strategies to train a DNN over a 3000 training epoch\. Throughout this work we will referred to this canonical training asConventional Training\. After this training we initialize two new DNN each with 1 additional feature added to their input layers\. For one DNN, this new feature is the output of the first DNN evaluated at the input training coordinates\. We call thisNN Featuredtraining\. In the other DNN, we provide the approximated output of the first DNN via CPD interpolated at the input training coordinates\. As explained early, for off\-grid points, we use a nearest neighbor interpolation\. We denote this second DNN as aTensor\-Featuredtrained model\.

![Refer to caption](https://arxiv.org/html/2608.10351v1/x1.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x2.png)\(b\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x3.png)\(c\)

Figure 1:Comparison of the convergence of different training procedures to learn[Eq\.7](https://arxiv.org/html/2608.10351#S4.E7)in \(a\)55, \(b\)1010, and \(c\)1515dimensions\.In[Fig\.1](https://arxiv.org/html/2608.10351#S4.F1)we compare accuracy of the three training strategies \(conventional, NN featured and Tensor\-Featured\) on the set of validation points over 3000 gradient updates\. For reference, we also show the error of the CPD interpolation of the conventionally trained NN on the same set of validation points\. One can see, across the board, decomposing the conventionally optimized DNN into a TN can significantly improve the accuracy of learned function\. As pointed out previously, we believe this is because the process of optimizing the tensor decomposition using non\-local information which effectively smooths the DNN function over the domain\.

With respect to Tensor\-Featured training, in[Fig\.1](https://arxiv.org/html/2608.10351#S4.F1)one can also see that using the DNN directly as a feature does not practically provide any improvement over the conventional training\. On the other hand with the CPD interpolation as a feature, we find an improvement in accuracy over the conventional training\. Unfortunately, although there is an improvement in accuracy with the Tensor\-Featured training in the validation set, the improvement is not dramatic and, in all cases, the CPD approximation alone outperforms the Tensor\-Featured DNN\. We believe this phenomena is not a failure in the Tensor\-Featured learning but, instead, associated with the accuracy of the conventionally trained model relative to the ground truth\.

To clarify this statement, consider the heatmaps in[Fig\.2](https://arxiv.org/html/2608.10351#S4.F2)which illustrates true and predicted values of[Eq\.7](https://arxiv.org/html/2608.10351#S4.E7)on a 2000 point interpolation grids with all but 2 modes of the tensor fixed to the 500th interpolation point\.

![Refer to caption](https://arxiv.org/html/2608.10351v1/x4.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x5.png)\(b\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x6.png)\(c\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x7.png)\(d\)

Figure 2:Heatmaps of \(a\) the true function \([Eq\.7](https://arxiv.org/html/2608.10351#S4.E7)\), \(b\) the function learned by the conventionally trained DNN, \(c\) the CPD approximation of the conventionally trained DNN and \(d\) the Tensor\-Featured trained DNN\. Heatmaps are evaluated on a 2000 point interpolation grid in all modes and all but 2 modes are fixed to the 500th point on their interpolation grid\. The DNN training is fixed to 1,000 samples in the training set\.One can see that, though[Fig\.1](https://arxiv.org/html/2608.10351#S4.F1)implies that the CPD approximation is a better representation of the true function, none of these models look significantly more accurate in this representation\. With this example, we demonstrate that because the reference function for the CPD \(the conventionally trained network\) is not representation of the ground truth, its addition to the feature list in the Tensor\-Featured training doesn’t significantly improve the DNN training\. Therefore, we hypothesize that, if the conventional training is more similar of the ground truth, then the Tensor\-Featured training procedure will be significantly accelerated\.

![Refer to caption](https://arxiv.org/html/2608.10351v1/x8.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x9.png)\(b\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x10.png)\(c\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x11.png)\(d\)

Figure 3:\(a\) The convergence of different training procedures to learn[Eq\.7](https://arxiv.org/html/2608.10351#S4.E7)withd=15d=15\. Heatmaps of \(b\) the function learned by the conventionally trained DNN, \(c\) the CPD approximation of the conventionally trained DNN and \(d\) the Tensor\-Featured trained DNN\. Heatmaps are evaluated on a 2000 point interpolation grid in all modes and all but 2 modes are fixed to the 500th point on their interpolation grid\. The DNN training is fixed to 10,000 samples in the training set\.In[Fig\.3](https://arxiv.org/html/2608.10351#S4.F3), we improve the accuracy of the conventionally trained DNN by increasing the training sample batch size from 1,000 to 10,000\. With a qualitatively more accuracy reference, the CPD approximation generates a function that captures more character of the true function\. Finally, equipped with this CPD, the Tensor\-Featured training is greatly improved and the mean squared error in the validation set falls dramatically\. In general, these results demonstrate that introduced features provide context for DNN training\. This 2\-step optimization process of training DNN with a fixed feature and decomposing partially trained DNN into CPD can be repeated iteratively in an attempt refine the Tensor\-Featured DNN\.

For example, we train the 15 dimensional example from above in a 3 step optimization procedure\. First a DNN is trained conventionally using 10,000 samples and 3000 training epoch\. After the training session, the DNN is discretized and decomposed into a CPD\. This CPD is passed as a feature into a new DNN which trained for 3000 iterations using the same 10,000 training points with the value of the CPD interpolation at each training point added to the feature list\. After this second training session, the resulting DNN is again decomposed into a CPD and this CPD is passed to the training of a third DNN\. In[Fig\.4](https://arxiv.org/html/2608.10351#S4.F4)we show the same 2D heatmaps for the CPD of the Tensor Featured trained DNN \(i\.e\.[Fig\.3\(d\)](https://arxiv.org/html/2608.10351#S4.F3.sf4)\) and the DNN constructed from the second round of Tensor\-Featured training\. While it is hard to tell if this network is an improvement on the first iteration, the data becomes clear when considering the relative L2 error and absolute mean squared error on the 2D slices in[Table1](https://arxiv.org/html/2608.10351#S4.T1)\. From these results one can see that this multi\-step training process, iteratively reduces the error in the trained DNNs\.

Table 1:Relative L2 and absolute mean squared error in DNN predicted functions on a 2D slice discretized by 2000 points in each dimension\.![Refer to caption](https://arxiv.org/html/2608.10351v1/x12.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x13.png)\(b\)

Figure 4:Heatmaps of \(a\) the CPD approximation of the first Tensor\-Featured trained DNN and \(b\) the second iteration Tensor\-Featured trained DNN\. Heatmaps are evaluated on a 2000 point interpolation grid in all modes and all but 2 modes are fixed to the 500th point on their interpolation grid\. The DNN training is fixed to 10,000 samples in the training set\.Because the tensor features provide context to the DNN optimization, we require that the initial training be relatively close \(in some measure\) to the global extrema\. If the original network is not within epsilon of the solution, the Tensor\-Featured training doesn’t help dramatically, if at all\. In the example above, we were able to improve the canonical training solution by increasing the number of batched training point samples\. In general, it may not be possible to improve the accuracy of conventional training so easily, for example there may be a limitation on the total number of available samples\. We are therefore interested in developing a strategy that introduces tensor features into the DNN training process from the start\. So far, we have only considered decomposing partially trained DNN however practically, we are able to introduce any one or many features into the Tensor\-Featured training process\. In the next section we will explore the idea of adding tensor features beyond TN decomposed DNN\.

### IV\.2Leveraging Tensor\-Featured Context to Accelerate Training

One of the biggest challenges of the Tensor\-Featured training procedure as presented in the previous section is need for features to provide meaningful context to a DNN’s training\. If a target function is unknown and the convergence of conventional DNN training is slow then the Tensor\-Featured training as presented thus far will not, in general, improve DNN optimization\. Therefore, to improve optimization, we consider introducing alternative and additional contextual features\.

In this section we define two types of features\. The first type of features are elementary functions which can be evaluated quickly for points in the domain of the unknown function\. Practically these features could be discretized and decomposed using TN methods but, because the functions can be evaluated quickly, we choose to simply access elements of the functions in a matrix\-free manner\. The first class of features we will callRank\-1 features\. The second type of feature are more complex functions which may be either slow to evaluate for points in the functions domain or otherwise limit the Tensor\-Featured optimization\. For this second class of features, we will leverage tensor decomposition strategies to address these issues\. An example of this type of feature are the partially trained DNN from the previous section\. We call this second class of featuresTN features\.

In this section we will attempt to learn a more difficult function, the solution to the initial boundary value problem of the hyperbolic \(wave\) equation

u​\(t,x\)=\(exp⁡\(t2\)−1\)​sin⁡\(π2​\(1−\|x\|\)2\.5\)\\displaystyle u\(t,x\)=\(\\exp\(t^\{2\}\)\-1\)\\sin\(\\frac\{\\pi\}\{2\}\(1\-\|x\|\)^\{2\.5\}\)\(8\)withd=15d=15\. In this study, we only consider the following rank\-1 features

- •y=∑ixiy=\\sum\_\{i\}x\_\{i\}
- •y=‖x‖2y=\\\|x\\\|\_\{2\}
- •y=sin⁡\(‖x‖2\)y=\\sin\(\\\|x\\\|\_\{2\}\)
- •y=cos⁡\(‖x‖2\)y=\\cos\(\\\|x\\\|\_\{2\}\)
- •y=∑isin⁡\(xi\)y=\\sum\_\{i\}\\sin\(x\_\{i\}\)
- •y=∑icos\(xi\)\)y=\\sum\_\{i\}\\cos\(x\_\{i\}\)\)
- •y=exp⁡\(‖x‖2\)y=\\exp\(\\\|x\\\|\_\{2\}\)

which mapℝn→ℝ\\mathbb\{R\}^\{n\}\\rightarrow\\mathbb\{R\}\. Our decision to use these rank\-1 features is relatively arbitrary\. Any function is acceptable as long as it mapℝk→ℝm\\mathbb\{R\}^\{k\}\\rightarrow\\mathbb\{R\}^\{m\}wherek≤nk\\leq nandℝm\\mathbb\{R\}^\{m\}is the dimensionality of function being learned\. In follow\-up studies we will investigate methods to automatically choose a good set of Rank\-1 features\. Furthermore, it should be noted that we will only consider the CPD of a previously trained DNN as a TN feature though in future work other TN features will be investigated\.

![Refer to caption](https://arxiv.org/html/2608.10351v1/x14.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x15.png)\(b\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x16.png)\(c\)

Figure 5:Comparison of the convergence of different training procedures to learn[Eq\.8](https://arxiv.org/html/2608.10351#S4.E8)\. \(a\) Comparison of trainings with at most a single feature added to the DNN\. \(b\) Comparison of trainings with single or multiple features added to the DNN\.In[Fig\.5](https://arxiv.org/html/2608.10351#S4.F5), we consider the effect of using different features in the Tensor\-Featured training\. In[Fig\.5\(a\)](https://arxiv.org/html/2608.10351#S4.F5.sf1)we add just a single to the training and find that the feature does influence the training efficiency\. One can see that adding a feature that does not align with the problem does not significantly improve the DNN optimization performance\. But adding a feature that does align with the problem improves the training\.

In[Fig\.5\(b\)](https://arxiv.org/html/2608.10351#S4.F5.sf2), we consider the effect of adding multiple features in parallel to the Tensor\-Featured training\. In this figure, we incrementally introduce the features from the rank\-1 feature list above in the order of the list\. One can see that introducing multiple features seems to systematically improve the accuracy over the conventional training and the best single feature training in[Fig\.5\(a\)](https://arxiv.org/html/2608.10351#S4.F5.sf1)\. In our tests we found that most of the functions in the rank\-1 feature list do not significantly improve the DNN training but together they outperform any single feature\. However, we note that this process can fail when a feature is added that puts the DNN close to a local extrema\. For example in[Fig\.5\(c\)](https://arxiv.org/html/2608.10351#S4.F5.sf3), we add the conventionally trained DNN into feature list and find that the optimization quickly converge back to the Conventional Training solution\. Though one can see that replacing the conventionally trained DNN with its CPD approximation can resolve this issue\.

![Refer to caption](https://arxiv.org/html/2608.10351v1/x17.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x18.png)\(b\)

Figure 6:Comparing the true values of[Eq\.8](https://arxiv.org/html/2608.10351#S4.E8)on the validation set to those predicted by trained DNN sorted from most negative to most positive\. \(a\) Fixes the training batch size to 1,000 points and \(b\) fixes the training batch size to 10,000 points\.Using Tensor\-Featured training with rank\-1 features, we can now improve the efficiency of training a DNN with no prior knowledge of the function\. Furthermore, the rank\-1 features do not depend on a pretrained DNN, like the CPD approximation did in the previous section\. In[Fig\.6](https://arxiv.org/html/2608.10351#S4.F6), we compare the accuracy of the conventional training to the Tensor\-Featured training with only rank\-1 features by plotting the true and predicted function values of points in the validation set\. Note that the points are ordered via the value of the true function from most negative to most positive\. One can see that the Tensor\-Featured training using only rank\-1 features more accurately predicts nearly all points in the validation set, compared to conventional training\. Furthermore, the Tensor\-Featured DNN training converges more quickly when the training batch size is increased\.

After training the DNN with rank\-1 features, we can discretize and decompose the trained DNN as a TN\. We compare the accuracy of the Tensor\-Featured trained DNN to those interpolated by the CPD of this DNN in[Fig\.7\(a\)](https://arxiv.org/html/2608.10351#S4.F7.sf1)\. From this figure, one can see that the decomposition process has little impact on the accuracy of predicting values in the validation set\. Furthermore, by decomposing the DNN we, effectively, cannibalize all of the features from the preliminary DNN optimization into a single function\. We demonstrate this in[Fig\.7\(b\)](https://arxiv.org/html/2608.10351#S4.F7.sf2)where we attempt to improve on the first Tensor\-Featured training\. In the curve marked ‘NN Feature‘ we use exactly the DNN trained via the Rank\-1 Tensor Featured Training into the optimization of a new DNN\. To do this optimization, we must evaluate the original DNN on the batch of training and validation points\. This evaluation may be relatively slow, especially if there were many input features in this DNN or the DNN is very deep\. The curve marked ‘TN Feature‘ uses only CPD interpolation of the Rank\-1 Featured Training\. Similar to the previous results, the CPD approximation of the DNN is more accurate than the DNN\. For the curve marked ‘Rank\-1 and TN Features‘ we combine the CPD with the original Rank\-1 features and find that the new Tensor\-Featured Training can become more accurate than the original Tensor\-Featured Training\. Because the original features were well captured by the CPD approximation, we are also free to modify the feature list in the DNN retraining\. For example for the curve marked ‘Modified Rank\-1 and TN Features‘ we remove features 1, 2, 5 and 6 from the input and find a result that outperforms the CPD approximation alone\. In these optimization curves we often see oscillation in the errors and will investigate this phenomenon in future studies\. This two step process of training DNN and updating features can be repeated iteratively until some convergence is reached\.

Finally in[Fig\.8](https://arxiv.org/html/2608.10351#S4.F8), we the accuracy of training a DNN to learn[Eq\.8](https://arxiv.org/html/2608.10351#S4.E8)in 40 dimensions with either conventional training or Tensor\-Featured training using only rank\-1 features\. One can see that the Tensor\-Featured training more accurately learns the objective function in just 3000 epoch\. Though comparing this figure to[Fig\.6\(b\)](https://arxiv.org/html/2608.10351#S4.F6.sf2), one can see that increasing the dimensionality of the problem slows the convergence of both training procedures\. In[Table2](https://arxiv.org/html/2608.10351#S4.T2), we show the cost comparing conventional training to Tensor\-Featured with only rank\-1 features and see that the Tensor\-Featured training is only marginally more expensive in both optimization time and storage complexity\.

Table 2:Costs associated with training different networks to represent[Eq\.8](https://arxiv.org/html/2608.10351#S4.E8)in 40 dimensions![Refer to caption](https://arxiv.org/html/2608.10351v1/x19.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10351v1/x20.png)\(b\)

Figure 7:\(a\) Comparing the true values of[Eq\.8](https://arxiv.org/html/2608.10351#S4.E8)on the validation set to those predicted by a Tensor\-Featured trained DNN and the CPD of the Tensor\-Featured trained DNN sorted from most negative to most positive\. \(b\) Comparison of the convergence of different training procedures to learn[Eq\.8](https://arxiv.org/html/2608.10351#S4.E8)\. All trainings utilize a training batch size of 10,000 points\.![Refer to caption](https://arxiv.org/html/2608.10351v1/x21.png)Figure 8:\(a\) Comparing the true values of[Eq\.8](https://arxiv.org/html/2608.10351#S4.E8)whered=40d=40on the validation set to those predicted by a Tensor\-Featured trained DNN sorted from most negative to most positive\. Trainings utilize a training batch size of 10,000 points\.

## VConclusion and Outlook

In this work we introduce an strategy to add context\-oriented features into the input layer of a DNN and leverage tensor decomposition strategy to develop a 2\-step optimization strategy for DNN training\. We are motivated by the idea of introducing pretrained DNN as input features similar to a warm\-start initialization strategies and by the outer product structure of tensor decompositions\. We find that DNN optimized with a pretrained DNN feature quickly converge to the same local minima as their pretrained feature\. Therefore, to improve this optimization strategy, we discretize the DNN over the domain of the problem on a quadrature grid and decompose the resulting higher\-order tensor using efficient randomized tensor decomposition strategies\. We restrict this work by only considering the CPD and to optimize the CPD we utilize a leveraged\-score based column selected randomized ALS scheme\. This CPD\-RALS strategy reducing the overall memory requirements of the decomposition by between 8 and 34 orders of magnitude\. We recognize that tensor decomposition algorithms like the ALS have a relationship to multidimensional gaussian process strategies\. This global optimization strategy appears to have a smoothing effect on functions that are encoded in trained DNN\. We introduce this TN approximation as a feature in the DNN optimization as a means to efficiently kick optimizations out of local minima\.

Our initial testing shows mixed results when the only feature introduced is the TN decomposition of a DNN\. Primarily, we found that, if the trained DNN is relatively close to a global extrema, Tensor\-Featured training can significantly improve the optimization of the subsequent DNN\. However, if the trained DNN is not close to a global extrema then the tensor feature does not improve the training, i\.e\. the context provided by the tensor feature is irrelevant\. We recognize that it may be difficult to train a DNN to be sufficiently close to a global extrema, which limits the utility of the Tensor\-Featured training as described above\. As a results, we introduce a more robust strategy to improve optimization\.

For this robust strategy, we recognize it is acceptable to introduce any number of functions that acts on any portion of the domain as features to the DNN optimization\. We therefore introduce an improved method where many functions, which we believe may provide context, are added as features to the DNN optimization\. The category of features can be split into functions which are easy to apply to points on the domain, which we call rank\-1 features, and functions which must first be decomposed as a TN to be applied efficiently, i\.e\. partially optimized DNN\. Using a combination of rank\-1 and TN features, we are able to significantly accelerate the convergence of model error in the validation set compared to conventional DNN training procedures\.

In follow\-up work, we will consider methods to efficiently choose rank\-1 features, determine the location and number of quadrature grid points for TN decomposition strategies and pick which TN decomposition strategy is necessary for particular problems\. We are also interested in studying how the product structure of these features influences the optimization of parameters in a DNN\. For example, when does adding a feature improve optimization, how do different features influence optimization of DNN parameters, and is there an interplay between activation functions and rank\-1 features? Finally, in this work we found some difficulties when optimizing extremely high dimensional functions using the randomized CPD\-ALS\. We are, therefore, interested in determining methods to improve the CPD optimization such as leveraging the symmetry and sparsity of high\-dimensional PDEs\.

## References

- \[1\]C\. Battaglino, G\. Ballard, and T\. G\. Kolda\(2018\)A practical randomized cp tensor decomposition\.SIAM Journal on Matrix Analysis and Applications39\(2\),pp\. 876–901\.Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p3.4),[§II\.2](https://arxiv.org/html/2608.10351#S2.SS2.p3.9)\.
- \[2\]G\. Beylkin and M\. J\. Mohlenkamp\(2002\-07\)Numerical operator calculus in higher dimensions\.Proc\. Natl\. Acad\. Sci\.99\(16\),pp\. 10246–10251\.External Links:[Document](https://dx.doi.org/10.1162/089976699300016863),[Link](http://www.pnas.org/cgi/doi/10.1073/pnas.112329799)Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p2.2)\.
- \[3\]J\. D\. Carroll and J\. Chang\(1970\)Analysis of individual differences in multidimensional scaling via an n\-way generalization of “eckart\-young” decomposition\.Psychometrika35\(3\),pp\. 283–319\.Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p5.1)\.
- \[4\]I\. Fakih, L\. Grigori, and K\. Pierce\(2026\-07\)Accelerating the Canonical Polyadic Alternating Least Squares Optimization via a Randomized Interpolative Decomposition\.Arxiv\.Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p3.4),[§II\.2](https://arxiv.org/html/2608.10351#S2.SS2.p3.9)\.
- \[5\]Y\. Gu, H\. Yang, and C\. Zhou\(2021\-09\)SelectNet: Self\-paced learning for high\-dimensional partial differential equations\.Journal of Computational Physics441,pp\. 110444\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2021.110444),ISSN 00219991Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p2.1)\.
- \[6\]R\. A\. Harshmanet al\.\(1970\)Foundations of the parafac procedure: models and conditions for an “explanatory” multi\-modal factor analysis\.UCLA working papers in phonetics16\(1\),pp\. 84\.Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p5.1)\.
- \[7\]F\. L\. Hitchcock\(1927\)The expression of a tensor or a polyadic as a sum of products\.Journal of Mathematics and Physics6\(1\-4\),pp\. 164–189\.Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p5.1)\.
- \[8\]Z\. Hu, K\. Shukla, G\. E\. Karniadakis, and K\. Kawaguchi\(2024\-08\)Tackling the curse of dimensionality with physics\-informed neural networks\.Neural Networks176,pp\. 106369\.External Links:[Document](https://dx.doi.org/10.1016/j.neunet.2024.106369),ISSN 08936080,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0893608024002934)Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p2.1)\.
- \[9\]M\. Innes, E\. Saba, K\. Fischer, D\. Gandhi, M\. C\. Rudilosso, N\. M\. Joy, T\. Karmali, A\. Pal, and V\. Shah\(2018\)Fashionable modelling with flux\.CoRRabs/1811\.01457\.External Links:[Link](https://arxiv.org/abs/1811.01457),1811\.01457Cited by:[§III](https://arxiv.org/html/2608.10351#S3.p3.1)\.
- \[10\]M\. Innes\(2018\)Flux: elegant machine learning with julia\.Journal of Open Source Software\.External Links:[Document](https://dx.doi.org/10.21105/joss.00602)Cited by:[§III](https://arxiv.org/html/2608.10351#S3.p3.1)\.
- \[11\]\(2026\)ITensorCPD library\.Note:Accessed: 2026\-01\-01External Links:[Link](https://github.com/kmp5VT/ITensorCPD.jl)Cited by:[§III](https://arxiv.org/html/2608.10351#S3.p3.1)\.
- \[12\]M\. I\. Jordan and T\. M\. Mitchell\(2015\-07\)Machine learning: trends, perspectives, and prospects\.Science349,pp\. 255–260\.External Links:[Document](https://dx.doi.org/10.1126/science.aaa8415),ISSN 0036\-8075Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p1.1)\.
- \[13\]P\. M\. Kroonenberg and J\. de Leeuw\(1980\-03\)Principal component analysis of three\-mode data by means of alternating least squares algorithms\.Psychometrika45\(1\),pp\. 69–97\.External Links:[Document](https://dx.doi.org/10.1121/1.1976063),[Link](http://link.springer.com/10.1007/BF02293599)Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p2.2)\.
- \[14\]B\. W\. Larsen and T\. G\. Kolda\(2022\)Practical leverage\-based sampling for low\-rank tensor decomposition\.SIAM Journal on Matrix Analysis and Applications43\(3\),pp\. 1488–1517\.Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p3.4)\.
- \[15\]C\. Liao\(2026\-01\)Solving partial differential equations with random feature models\.Communications in Nonlinear Science and Numerical Simulation152,pp\. 109343\.External Links:[Document](https://dx.doi.org/10.1016/j.cnsns.2025.109343),ISSN 10075704,[Link](https://linkinghub.elsevier.com/retrieve/pii/S100757042500752X)Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p2.1)\.
- \[16\]K\. Pierce, V\. Rishi, and E\. F\. Valeev\(2021\-03\)Robust Approximation of Tensor Networks: Application to Grid\-Free Tensor Factorization of the Coulomb Interaction\.J\. Chem\. Theory Comput\.17\(4\),pp\. 2217–2230\.External Links:ISSN 1549\-9618, 1549\-9626,[Document](https://dx.doi.org/10.1021/acs.jctc.0c01310),[Link](https://pubs.acs.org/doi/10.1021/acs.jctc.0c01310)Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p2.2)\.
- \[17\]J\. Quinonero\-Candela and C\. Rasmussen\(2005\-12\)A unifying view of sparse approximate gaussian process regression\.Journal of Machine Learning Research, v\.6, 1935\-1959 \(2005\)6,pp\.\.Cited by:[§II\.2](https://arxiv.org/html/2608.10351#S2.SS2.p1.4)\.
- \[18\]M\. Raissi, P\. Perdikaris, and G\.E\. Karniadakis\(2019\-02\)Physics\-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\.Journal of Computational Physics378,pp\. 686–707\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2018.10.045),ISSN 00219991,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0021999118307125)Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p1.1)\.
- \[19\]C\. Edward\. Rasmussen and C\. K\. I\.\. Williams\(2006\)Gaussian processes for machine learning\.MIT Press\.External Links:ISBN 026218253XCited by:[§II\.2](https://arxiv.org/html/2608.10351#S2.SS2.p1.4),[§II\.2](https://arxiv.org/html/2608.10351#S2.SS2.p2.18)\.
- \[20\]L\. Richter, L\. Sallandt, and N\. Nüsken\(2021\-07\)Solving high\-dimensional parabolic pdes using the tensor train format\.Arxiv\.Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p3.1)\.
- \[21\]V\. Rokhlin and M\. Tygert\(2008\-09\)A fast randomized algorithm for overdetermined linear least\-squares regression\.Proceedings of the National Academy of Sciences105,pp\. 13212–13217\.External Links:[Document](https://dx.doi.org/10.1073/pnas.0804869105),ISSN 0027\-8424,[Link](https://pnas.org/doi/full/10.1073/pnas.0804869105)Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p3.4)\.
- \[22\]E\. Stoudenmire and D\. Schwab\(2016\)Supervised learning with tensor networks\.InAdvances in Neural Information Processing Systems,D\. Lee, M\. Sugiyama, U\. Luxburg, I\. Guyon, and R\. Garnett \(Eds\.\),Vol\.29,pp\.\.External Links:[Link](https://proceedings.neurips.cc/paper_files/paper/2016/file/5314b9674c86e3f9d1ba25ef9bb32895-Paper.pdf)Cited by:[§II\.2](https://arxiv.org/html/2608.10351#S2.SS2.p1.4)\.
- \[23\]J\. F\. Urbán, P\. Stefanou, and J\. A\. Pons\(2025\-02\)Unveiling the optimization process of physics informed neural networks: how accurate and competitive can pinns be?\.Journal of Computational Physics523,pp\. 113656\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2024.113656),ISSN 00219991,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0021999124009045)Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p2.1)\.
- \[24\]F\. Wesel and K\. Batselier\(2024\-03\)Tensor network\-constrained kernel machines as gaussian processes\.Arxiv\.External Links:[Link](http://arxiv.org/abs/2403.19500)Cited by:[§II\.2](https://arxiv.org/html/2608.10351#S2.SS2.p1.4)\.
- \[25\]F\. Wesel and K\. Batselier\(2025\-03–05 May\)Tensor network\-constrained kernel machines as gaussian processes\.InProceedings of The 28th International Conference on Artificial Intelligence and Statistics,Y\. Li, S\. Mandt, S\. Agrawal, and E\. Khan \(Eds\.\),Proceedings of Machine Learning Research, Vol\.258,pp\. 2161–2169\.External Links:[Link](https://proceedings.mlr.press/v258/wesel25a.html)Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p5.1)\.
- \[26\]D\. P\. Woodruff\(2014\-10\)Sketching as a tool for numerical linear algebra\.Foundations and Trends® in Theoretical Computer Science10,pp\. 1–157\.External Links:[Document](https://dx.doi.org/10.1561/0400000060),ISSN 1551\-305X,[Link](https://www.emerald.com/fttcs/article/10/1-2/1/1332206/Sketching-as-a-Tool-for-Numerical-Linear-Algebra)Cited by:[§II\.1](https://arxiv.org/html/2608.10351#S2.SS1.p3.4)\.
- \[27\]Q\. Yuan, Z\. Xu, Y\. Chen, Y\. Xu, H\. Owhadi, and S\. Zhe\(2026\-03\)Tensor gaussian processes: efficient solvers for nonlinear pdes\.Arxiv\.External Links:[Link](http://arxiv.org/abs/2510.13772)Cited by:[§I](https://arxiv.org/html/2608.10351#S1.p5.1)\.