# Tensor approximation of functional differential equations

Abram Rodgers<EMAIL_ADDRESS>Daniele Venturi<EMAIL_ADDRESS>∗Advanced Supercomputing Division, NASA Ames Research Center,
†Department of Applied Mathematics, University of California Santa Cruz

###### Abstract

Functional Differential Equations (FDEs) play a fundamental role in many areas
of mathematical physics, including fluid dynamics (Hopf characteristic
functional equation), quantum field theory (Schwinger-Dyson equation), and
statistical physics. Despite their significance, computing solutions to FDEs
remains a longstanding challenge in mathematical physics. In this paper we
address this challenge by introducing new approximation theory and high-
performance computational algorithms designed for solving FDEs on tensor
manifolds. Our approach involves approximating FDEs using high-dimensional
partial differential equations (PDEs), and then solving such high-dimensional
PDEs on a low-rank tensor manifold leveraging high-performance parallel tensor
algorithms. The effectiveness of the proposed approach is demonstrated through
its application to the Burgers-Hopf FDE, which governs the characteristic
functional of the stochastic solution to the Burgers equation evolving from a
random initial state.

††preprint: APS/123-QED

##  I Introduction

Functional differential equations (FDEs) are equations involving operators
(nonlinear functionals), and derivatives/integrals of such operators with
respect functions (functional derivatives) and other independent variables
such as space and time [1, 2]. FDEs arise naturally in many different branches
of physics. A classical example in fluid dynamics is the Hopf-Navier-Stokes
equation [3, 4, 5, 6]

$\displaystyle\frac{\partial\Phi([\bm{\theta}],t)}{\partial
t}=\sum_{k=1}^{3}\int_{\Omega}\theta_{k}(\bm{x})$
$\displaystyle\left(i\sum_{j=1}^{3}\frac{\partial}{\partial
x_{j}}\frac{\delta^{2}\Phi([\bm{\theta}],t)}{\delta\theta_{k}(\bm{x})\delta\theta_{j}(\bm{x})}+\right.$
$\displaystyle\nu\nabla^{2}\frac{\delta\Phi([\bm{\theta}],t)}{\delta\theta_{k}(\bm{x})}\Biggl{)}d\bm{x},$
(1)

which governs the temporal evolution of the characteristic functional

$\Phi([\bm{\theta}],t)=\mathbb{E}\left\\{\exp\left(i\int_{\Omega}\bm{u}(\bm{x},t)\cdot\bm{\theta}(\bm{x})d\bm{x}\right)\right\\}.$
(2)

Here, $\Omega\subseteq\mathbb{R}^{3}$ is the spatial domain,
$\bm{u}(\bm{x},t)$ represents a stochastic solution of the Navier-Stokes
equations

$\begin{cases}\displaystyle\frac{\partial\bm{u}}{\partial
t}+(\bm{u}\cdot\nabla)\bm{u}=-\nabla p+\nu\nabla^{2}\bm{u},\vspace{0.2cm}\\\
\nabla\cdot\bm{u}=0,\end{cases}$ (3)

corresponding to a random initial velocity distribution [7, 8],
$\mathbb{E}\\{\cdot\\}$ is the expectation over the probability measure of
such random velocity distribution, and $\bm{\theta}(\bm{x})$ is a divergence-
free test function. Equation (1) involves derivatives of the functional
$\Phi([\bm{\theta}],t)$ with respect to the functions $\theta_{i}(\bm{x})$,
i.e., functional derivatives $\delta/\delta\theta_{i}(\bm{x})$ [2], and
derivatives with respect to the independent variables $x_{j}$ and $t$.

As is well-known, the characteristic functional (2) encodes the full
statistical information of the random field $\bm{u}(\bm{x},t)$ that solves the
Navier-Stokes equations (3), including multi-point statistical moments,
cumulants, and multi-point joint probability density functions. For instance,
by taking functional derivatives of $\Phi$ we can immediately express the mean
and auto-correlation function of $\bm{u}(\bm{x},t)$ as

$\displaystyle\mathbb{E}\left\\{u_{k}(\bm{x},t)\right\\}$
$\displaystyle=\frac{1}{i}\left.\frac{\delta\Phi([\bm{\theta}],t)}{\delta\theta_{k}(\bm{x})}\right|_{\bm{\theta}=\bm{0}},$
$\displaystyle\mathbb{E}\left\\{u_{k}(\bm{x},t)u_{j}(\bm{y},t)\right\\}$
$\displaystyle=\frac{1}{i^{2}}\left.\frac{\delta^{2}\Phi([\bm{\theta}],t)}{\delta\theta_{k}(\bm{x})\delta\theta_{j}(\bm{y})}\right|_{\bm{\theta}=\bm{0}}.$
(4)

Similarly, by evaluating (2) at
$\theta_{i}(\bm{x})=a_{i}\delta(\bm{x}-\bm{y})$ yields the one-point
characteristic function of the solution $\bm{u}(\bm{y},t)$ at an arbitrary
spatial location $\bm{y}\in\Omega$ as

$\displaystyle\phi(\bm{a},\bm{y},t)$
$\displaystyle=\mathbb{E}\left\\{\exp\left(i\int_{\Omega}\bm{u}(\bm{x},t)\cdot\bm{a}\delta(\bm{x}-\bm{y})d\bm{x}\right)\right\\}$
$\displaystyle=\mathbb{E}\left\\{e^{i\bm{a}\cdot\bm{u}(\bm{y},t)}\right\\}.$
(5)

The inverse Fourier transform (in the sense of distributions) of
$\phi(\bm{a},\bm{y},t)$ with respect to the variable $\bm{a}$ is the
probability density function PDF of the solution $\bm{u}(\bm{y},t)$ at the
spatial location $\bm{y}\in\Omega$ and time $t$. Any other statistical
property of $\bm{u}(\bm{x},t)$ can be derived from the characteristic
functional (2). For this reason, the FDE (1) was deemed by Monin and Yaglom
([4, Ch. 10]) to be “the most compact formulation of the general turbulence
problem”, which is the problem of determining the statistical properties of
the velocity and the pressure fields of the Navier-Stokes equations given
statistical information on the initial state.

Another well-known classical example of functional differential equation is
the Schwinger-Dyson (SD) equation of quantum field theory [9, 10]. The SD
equation describes the dynamics of the generating functional of the Green
functions of a quantum field theory, allowing us to propagate quantum field
interactions in a perturbation setting (e.g., with Feynman diagrams), or in a
strong coupling regime. The usage of FDEs saw rapid growth in the 1970s,
driven by the realization that techniques originally developed for quantum
field theory by Dyson, Feynman, and Schwinger could be, at least formally,
extended to other branches of mathematical physics. A pivotal contribution to
this evolution was the groundbreaking work by Martin, Siggia, and Rose [11].
This work served as a landmark in the field, showcasing the possibility of
applying quantum field theoretic methods, such as functional integrals and
diagrammatic expansions [12, 13, 14, 15], to problems in classical statistical
physics.

More recently, FDEs have appeared in both mean-field games and mean-field
optimal control [16, 17]. Mean-field games are optimization problems involving
a large (potentially infinite) number of interacting players. In some cases,
it is possible to reformulate such optimization problems in terms of a
nonlinear Hamilton-Jacobi-Bellman FDE in a Wasserstein space of probability
measures. The standard form of such equation is (see [18, p. 1])

$\displaystyle\frac{\partial F([\rho],t)}{\partial t}$
$\displaystyle+\mathcal{W}\left([\rho],\left[\frac{\delta
F([\rho],t)}{\delta\rho(\bm{x})}\right]\right)=0,$ $\displaystyle F([\rho],0)$
$\displaystyle=F_{0}([\rho]),$ (6)

where $\rho(\bm{x})$ is a $d$-dimensional probability density function
supported on $R\subseteq\mathbb{R}^{d}$, $\delta F/\delta\rho(\bm{x})$ is the
first-order functional derivative of $F$ relative to $\rho(\bm{x})$, and
$\mathcal{W}$ is the Hamilton functional

$\displaystyle\mathcal{W}\left([\rho],\left[\frac{\delta
F([\rho],t)}{\delta\rho(\bm{x})}\right]\right)=\int_{R}$
$\displaystyle{\Psi}\left(\bm{x},\nabla\frac{\delta
F([\rho],t)}{\delta\rho(\bm{x})}\right)\rho(\bm{x})d\bm{x}$
$\displaystyle+\mathcal{J}([\rho]).$ (7)

Here, $\Psi$ is the Hamilton’s function and $\mathcal{J}([\rho])$ is an
interaction potential. Mean-field theory is also useful in optimal feedback
control of nonlinear stochastic dynamical systems and in deep learning. For
instance, E and collaborators [17] established the mathematical groundwork for
the population risk minimization problem in deep learning by framing it as a
mean-field optimal control problem. This yields yields a generalized version
of the Hamilton-Jacobi-Bellman equation in a Wasserstein space of probability
measures, which is a nonlinear FDE of the form (6) – e.g., equation (20) in
[17] and equation (1.1) in [19].

In this paper we develop approximation theory and high-performance
computational algorithms designed for solving FDEs on tensor manifolds. Our
approach involves first approximating FDEs in terms of high-dimensional
partial differential equations (PDEs), and then computing the solution to such
high-dimensional PDEs on a low-rank tensor manifold leveraging high-
performance (parallel) tensor algorithms. To this end, we will build upon our
recent work on tensors, in particular step-truncation tensor methods [20, 21],
and demonstrate convergence to functional approximations of FDEs.

This paper is organized as follows. In Section II we briefly review necessary
and sufficient conditions for existence and uniqueness of solutions to FDEs.
In Section III we discuss functional approximation of FDEs in terms of high-
dimensional PDEs. In Section IV we briefly review numerical tensor methods
approximate the solution of high-dimensional PDEs on tensor manifolds. In
Section IV we demonstrate the new functional approximation methods for FDEs to
the Burgers-Hopf equation. The main findings are summarized in Section VI. We
also include Appendix A in which we describe the high-performance (parallel)
tensor algorithms we developed to solve high-dimensional PDEs corresponding to
FDEs on tensor manifolds.

##  II Existence and uniqueness of solutions to FDEs

Let us briefly describe what the problem is here, and our intended course of
action. To this end, consider the Hopf-Navier-Stokes functional differential
equation (1). Back in 1972 Monin and Yaglom stated in [4, p. 773] that:

“When we tried to develop a complete statistical description of turbulence
with the aid of the Hopf equation for the characteristic functional we found
that no general mathematical formalism for solving linear equations in
functional derivatives was available. There are also no rigorous theorems on
the existence and uniqueness of the solution to such equations.”

Since then, there were advancements in the theory of existence and uniqueness
of solutions to the FDE (1). In particular, it was found the existence of a
solution to (1) is strongly related to the existence of a Hopf-Leray weak
solution for the Cauchy problem of the Navier–Stokes equations, which is a
well-established result [22, 23, 24, 25] (see also [26, Theorem 3.1]).
Specifically, it can be proved that if there exists a Hopf-Leray solution to
the Navier-Stokes equations then there exists a statistical solution to the
Hopf-Navier-Stokes equation (1) (see [27, §5 in Ch. 4]). A similar results was
obtained by Foias in [8, p. 254] for initial velocity conditions described by
probability measures with finite second moments.

Regarding uniqueness of the solution to the Hopf-Navier-Stokes equation (1),
it turns out that this is again connected to the uniqueness of Hopf-Leray weak
solutions to the Navier-Stokes equation. Indeed, Vishik and Fursikov proved in
[27, Theorem 6.1] that if there exists a unique Hopf-Leray solution to the
Navier-Stokes equation then any space-time statistical solution of the Hopf-
Navier-Stokes equation is uniquely determined by the initial probability
measure (see also [8, pp. 323-344]). In two spatial dimensions it is well-
known that there exists a unique global Hopf-Leray solution, and therefore a
unique statistical solution. The proof was formalized by Gishlarkaev in [28]
for the Hopf-Navier-Stokes equation using techniques developed in [27]. The
uniqueness of the Hopf-Leray solution for the three-dimensional Navier-Stokes
equations is still an open problem, and a rather active area of research. For
example, it can be shown that the Hopf-Leray solution is unique within at
least a short-time integration period if we choose the initial condition in
non-conventional function spaces [29, 30, 24, 23, 25]. Also, the Hopf-Leray
solution is unique and global in time if we choose the initial data small
enough [31, 25]. Recently, Galdi proved in [32] that very weak solutions to
the Cauchy problem for the Navier-Stokes equations must be Hopf-Leray
solutions if their initial data are solenoidal with finite kinetic energy. A
smooth (mild or strong) solution of the Navier-Stokes equations is of course
also a Hopf-Leray (weak) solution. More precisely, it can be shown that a mild
solution to Navier-Stokes in $H^{s}$ agrees with a weak Hopf-Leray solution
almost everywhere in $[0,T^{*})$. This result, which belongs to the so-called
weak-strong uniqueness methods [33], implies that weak Hopf-Leray solutions
are unique as long as a mild (or a strong) solution exists. In particular,
Leray-Hopf solutions are unique whenever they are regular enough to be strong
solutions. Once the strong solution is integrated past $t\sim T^{*}$ (maximum
Cauchy development), the is no further guarantee of uniqueness [34, 25]. We
also recall that if the energy inequality that characterizes Hopf-Leray weak
solutions is dropped then it is possible to construct multiple weak solutions
of the Navier-Stokes equations with finite kinetic energy [35]. On the other
hand, weak solutions obeying the energy inequality are necessarily unique and
smooth if they lie in $L^{p}([0,T])\otimes L^{q}(\Omega)$ with $2/p+3/q=1$ and
$q>3$. These conditions are known as Prodi-Serrin-Ladyzhenskaya conditions,
and they have been recently generalized in various ways (see [23, 36, 37]).

##  III Approximation of FDEs by high-dimensional PDEs

Consider a well-posed initial/boundary value problem for an FDE, i.e., a
problem that admits a unique solution in some function space as discussed in
Section II. How do we compute the solution? This is a longstanding problem in
mathematical physics and a rather new research area of computational
mathematics.

The FDEs (1) and (6) can be seen as PDEs in an infinite number of independent
variables. Such infinite-dimensional PDEs may be approximated by PDEs in a
finite (though possibly extremely large) number of variables using the
functional analytic methods described in [2, 1, 19, 18]. For instance, we have
shown in [1, 2] that equation (1) can be approximated by the high-dimensional
linear PDE

$\frac{\partial\phi}{\partial
t}=\sum_{p=1}^{n}a_{p}\left(\nu\sum_{k=1}^{n}B_{pk}\frac{\partial\phi}{\partial
a_{k}}+i\sum_{j,k=1}^{n}A_{pjk}\frac{\partial^{2}\phi}{\partial a_{k}\partial
a_{j}}\right),$ (8)

where

$\displaystyle a_{p}=\int_{\Omega}\bm{\theta}\cdot\bm{\Gamma}_{p}d\bm{x},$
$\displaystyle\qquad
B_{pk}=\int_{\Omega}{\bm{\Gamma}}_{p}\cdot\nabla^{2}{\bm{\Gamma}}_{k}d\bm{x},$
$\displaystyle A_{pjk}$
$\displaystyle=\int_{\Omega}{\bm{\Gamma}}_{p}\cdot\left[\left({\bm{\Gamma}}_{k}\cdot\nabla\right){\bm{\Gamma}}_{j}\right]d\bm{x},$
(9)

$\\{\bm{\Gamma}_{1}(\bm{x}),\ldots,\bm{\Gamma}_{n}(\bm{x})\\}$ is a
divergence-free basis, and $\phi(a_{1},\ldots,a_{n},t)$ is the characteristic
function approximating the characteristic functional $\Phi([\bm{\theta}],t)$
in the following sense:

$\phi(a_{1},\ldots,a_{n},t)=\Phi([P_{n}\bm{\theta}],t),$ (10)

where $P_{n}$ is an orthogonal projection onto
$\\{\bm{\Gamma}_{1},\ldots,\bm{\Gamma}_{n}\\}$. The divergence-free basis
$\\{\bm{\Gamma}_{1},\ldots,\bm{\Gamma}_{n}\\}$ can be constructed in terms of
div-free wavelets [38, 39], radial basis functions [40], trigonometric
polynomials [41], or eigenvalue problems for the Stokes operator with
appropriate boundary conditions [22]. Approximating the FDE (1) using the PDE
(8) bears a resemblance, in a certain aspect, to what is known as the “method
of lines” for discretizing a nonlinear partial differential equation into a
high-dimensional system of nonlinear ordinary differential equations (ODEs).
However, in the context of FDEs, the “lines” are determined by a hyper-
dimensional PDE instead of a system of ODEs.

In [2, Theorem 7.1] we proved that if we choose the test function
$\bm{\theta}(\bm{x})$ in a Sobolev sphere $\mathcal{S}(\Omega)$ of
$H^{1}(\Omega)$-div functions (which is a compact subset of
$L^{2}(\Omega)$-div) then $\phi$ converges uniformly to the characteristic
functional $\Phi$ as the number of independent variables $n$ goes to infinity,
i.e.,

$\lim_{n\rightarrow\infty}\sup_{\begin{subarray}{c}\bm{\theta}\in\mathcal{S}(\Omega)\\\
t\in[0,T]\end{subarray}}\left|\phi(a_{1},\ldots,a_{n},t)-\Phi([\bm{\theta}],t)\right|=0.$
(11)

Moreover, convergence of the high-dimensional PDE to the FDE can be
exponential. The result (11) leverages the Trotter-Kato approximation theorem
[42, p. 210] for abstract evolution equations in Banach spaces [43].

Gangbo et al. [19] proved a similar result for the nonlinear functional HJB
equation (6). Specifically, by evaluating the functional equation (6) on the
Wasserstein space $\mathcal{P}_{2}\left(\mathbb{R}^{d}\right)$ (space of
probability densities in $\mathbb{R}^{d}$ with finite second-order moment)
with elements

$\rho_{n}(\bm{x})=\frac{1}{n}\sum_{k=1}^{n}\delta(\bm{x}-\bm{x}_{k})\qquad\bm{x},\bm{x}_{i}\in\mathbb{R}^{d}$
(12)

then we can write the FDE (6) as a hyper-dimensional nonlinear PDE

$\frac{\partial f}{\partial
t}+\frac{1}{n}\sum_{k=1}^{n}\Psi\left(\bm{x}_{i},\frac{\partial
f}{\partial\bm{x}_{k}}\right)+j(\bm{x}_{1},\ldots,\bm{x}_{n})=0,$ (13)

where $f(\bm{x}_{1},\ldots,\bm{x}_{n},t)=F([\rho_{n}],t)$ (see [19,
Eq.(1.2)]). Moreover, if we choose the empirical measure (12) in a bounded
subset ${B}$ of the Wasserstein space
$\mathcal{P}_{2}\left(\mathbb{R}^{d}\right)$ then
$f(\bm{x}_{1},\ldots,\bm{x}_{n},t)$ converges uniformly to the (unique)
solution of (6) as $n$ goes to infinity ([19, Theorem 1.2]), i.e.,

$\lim_{n\rightarrow\infty}\sup_{\begin{subarray}{c}\rho_{n}\in{B}\\\
t\in[0,T]\end{subarray}}\left|f(\bm{x}_{1},\ldots,\bm{x}_{n},t)-F([\rho_{n}],t)\right|=0.$
(14)

The boundedness assumption of $B$ is essential for convergence. Similarly, the
compactness assumption of the function space $S$ in (11) is essential for
convergence (see [2, Theorems 7.1 and 8.3]).

To solve high-dimensional PDEs such as (8) and (13) numerically, we need a
computational paradigm for linear PDEs that can handle tens, hundreds, or
thousands of independent variables efficiently. Such computational paradigm is
discussed in the next section.

##  IV Tensor approximation of high-dimensional PDEs

Figure 1: Examples of tensor networks. The vertices represent tensor modes
(functions) used in the decomposition. The edges connecting to vertices
represent summation over an index (contraction) between two modes. The free
edges represent input variables $x_{j}$ ($j=1,\ldots,n$).

We have seen in Section III that FDEs can be approximated by high-dimensional
PDEs of the form

$\displaystyle\frac{\partial f({\bm{x}},t)}{\partial t}={\cal
G}\left(f({\bm{x}},t),{\bm{x}}\right),\qquad f({\bm{x}},0)=f_{0}({\bm{x}}),$
(15)

where $f:\Omega\times[0,T]\to\mathbb{R}$ is a $n$-dimensional time-dependent
scalar field defined on the domain $\Omega\subseteq\mathbb{R}^{n}$, $T$ is the
period of integration, and $\cal G$ is a nonlinear operator which may depend
on the variables ${\bm{x}}=(x_{1},\ldots,x_{n})\in\Omega$, and may incorporate
boundary conditions. For simplicity, we assume that the domain $\Omega$ is a
Cartesian product of $n$ one-dimensional domains $\Omega_{i}$

$\Omega=\Omega_{1}\times\cdots\times\Omega_{n},$ (16)

and that $f$ is an element of a Hilbert space $H(\Omega;[0,T])$. In these
hypotheses, we can leverage the isomorphism $H(\Omega;[0,T])\simeq
H([0,T])\otimes H(\Omega_{1})\otimes\cdots\otimes H(\Omega_{n})$ and represent
the solution of (15) as

$f(\bm{x},t)\approx\sum_{i_{1}=1}^{m_{1}}\cdots\sum_{i_{n}=1}^{m_{n}}f_{i_{1}\ldots
i_{n}}(t)h_{i_{1}}(x_{1})\cdots h_{i_{1}}(x_{1}),$ (17)

where $h_{i_{j}}(x_{j})$ are one-dimensional orthonormal basis functions of
$H(\Omega_{i})$. Substituting (17) into (15) and projecting onto an
appropriate finite-dimensional subspace of $H(\Omega)$ yields the semi-
discrete form

$\frac{d{\bm{f}}}{dt}={\bm{G}}({\bm{f}}),\qquad{\bm{f}}(0)={\bm{f}}_{0}$ (18)

where ${\bm{f}}:[0,T]\rightarrow{\mathbb{R}}^{m_{1}\times
m_{2}\times\dots\times m_{n}}$ is a multivariate array with coefficients
$f_{i_{1}\ldots i_{d}}(t)$, and $\bm{G}$ is the finite-dimensional
representation of the nonlinear operator $\cal G$. The number of degrees of
freedom associated with the solution to the Cauchy problem (18) is
$N_{\text{dof}}=m_{1}m_{2}\cdots m_{n}$ at each time $t\geq 0$, which can be
extremely large even for moderately small dimension $n$. For instance, the
solution of the Boltzmann-BGK equation on a six-dimensional ($n=6$) flat torus
[44, 45] with $m_{i}=128$ basis functions in each position and momentum
variable yields $N_{\text{dof}}=128^{6}=4398046511104$ degrees of freedom at
each time $t$. This requires approximately $35.18$ Terabytes per temporal
snapshot if we store the solution tensor $\bm{f}$ in a double precision IEEE
754 floating point format. Several general-purpose algorithms have been
developed to mitigate such an exponential growth of degrees of freedom, the
computational cost, and the memory requirements. These algorithms include,
e.g., sparse collocation methods [46, 47, 48] and techniques based on deep
neural networks [49, 50].

In a parallel research effort that has its roots in quantum field theory and
quantum entanglement, researchers have recently developed a new generation of
algorithms based on tensor networks and low-rank tensor techniques to compute
the solution of high-dimensional PDEs [51, 52, 53, 54, 55, 56]. A tensor
network is a factorization of an entangled object such as a multivariate
function or an operator into simpler objects, e.g., low-dimensional functions
or separable operators, which are amenable to efficient representation and
computation. A tensor network can be visualized in terms of graphs (see Figure
1).

Figure 2:  Sketch of implicit and explicit step-truncation integration
methods. Given a tensor $\bm{f}_{k}$ with multilinear rank $\bm{s}$ on the
tensor manifold ${\cal M}_{\bm{s}}$, we first perform an explicit time-step
with any conventional time-stepping scheme. The explicit step-truncation
integrator then projects $\bm{\Psi}_{\Delta t}(\bm{G},\bm{f}_{k})$ to a tensor
manifold with rank $\bm{s}$. The multilinear rank $\bm{s}$ is chosen
adaptively based on desired accuracy and stability constraints [20]. On the
other hand, the implicit step-truncation method takes $\bm{\Psi}_{\Delta
t}(\bm{G},\bm{f}_{k})$ as input and generates a sequence of fixed-point
iterates ${\bm{f}}^{[j]}$ shown as dots connected with blue lines. The last
iterate is then projected onto a low rank tensor manifold, illustrated here
also as a red line landing on ${\cal M}_{\bm{s}}$. This operation is
equivalent to the compression step in the HT/TT-GMRES algorithm described in
[57].

The vast majority of tensor algorithms currently available to approximate
functions, operators and PDEs on tensor spaces is based on canonical polyadic
(CP) decompositions [56, 53, 44, 58], Tucker tensors, or tensors corresponding
to binary trees such as tensor train (TT) [59, 60, 61] and hierarchical Tucker
(HT) tensors [62, 63, 64, 65]. A compelling reason for using binary tensor
trees is that they allow to construct the tensor expansion by leveraging the
spectral theory for linear operators, in particular the hierarchical Schmidt
decomposition [66, 59, 67, 68, 69].

###  IV.1 Step-truncation tensor methods

A new class of algorithms to integrate (18) on a low-rank tensor manifold was
recently proposed in [20, 21]. These algorithms are known as step-truncation
methods and they are based on integrating the solution ${\bm{f}}(t)$ of the
ODE (18) off the tensor manifold for a short time using any conventional
explicit time-stepping scheme, and then mapping it back onto the manifold
using a tensor truncation operation. To briefly describe these methods, let us
discretize the ODE (18) in time with a one-step method on an evenly-spaced
temporal grid as

$\displaystyle{\bm{f}}_{k+1}={\bm{\Psi}}_{\Delta
t}({\bm{G}},{\bm{f}}_{k}),\qquad{\bm{f}}_{0}={\bm{f}}(0),$ (19)

where ${\bm{f}}_{k}$ denotes an approximation of ${\bm{f}}(k\Delta t)$ for
$k=0,1,\ldots$, and ${\bm{\Psi}}_{\Delta t}$ is an increment function. To
obtain a step-truncation integrator, we simply apply a truncation operator
$\mathfrak{T}_{\bm{s}}(\cdot)$, i.e., a nonlinear projection onto a tensor
manifold ${\cal M}_{\bm{s}}$ of HT (or TT) tensors with multilinear rank
$\bm{s}$ [70], to the scheme (19) (see Figure 2). This yields

$\displaystyle{\bm{f}}_{k+1}={\mathfrak{T}}_{\bm{s}}\left[{\bm{\Psi}}_{\Delta
t}({\bm{G}},{\bm{f}}_{k})\right].$ (20)

The need for tensor rank-reduction when iterating (19) can be easily
understood by noting that tensor operations such as the application of an
operator to a tensor and the addition between two tensors naturally increase
tensor rank [71].

We now present the rank-adaptive step-truncation algorithms used here to
integrate FDEs on tensor manifolds. The first one is a truncated version of
the Euler forward scheme, i.e.,

${\bm{f}}_{k+1}={\mathfrak{T}}_{{\bm{r}}}[{\bm{f}}_{k}+\Delta
t{\mathfrak{T}}_{{\bm{s}}}[{\bm{G}}({\bm{f}}_{k})]].$ (21)

We have shown that this scheme is convergent to first-order in $\Delta t$
provided the rank $\bm{r}$ is properly updated at each time step (see [20] for
details).

Figure 3: Sketch of the logical flow to transform a random nonlinear initial
value problem (IVP) for the Burgers equation into the (linear) Burgers-Hopf
FDE, which is linear. The Burgers-Hopf FDE is subsequently approximated using
the methods outlined in Section III. This yields an $N$-dimensional
characteristic function equation of the form (8). By taking the inverse
Fourier transform of such characteristic function equation and integrating it
with respect to the phase variable we obtain the $N$-dimensional CDF equation
shown above. The CDF equation is then solved using the proposed step-
truncation tensor methods.

The second step-truncation scheme is derived from the explicit midpoint rule,
and it is defined as

${\bm{f}}_{k+1}={\mathfrak{T}}_{{\bm{\alpha}}}\left[{\bm{f}}_{k}+\Delta
t{\mathfrak{T}}_{{\bm{\beta}}}\left[{\bm{G}}\left({\bm{f}}_{k}+\frac{\Delta
t}{2}{\mathfrak{T}}_{{\bm{\gamma}}}({\bm{G}}({\bm{f}}_{k}))\right)\right]\right],$
(22)

where $\bm{\alpha}$, $\bm{\beta}$, and $\bm{\gamma}$ are once again selected
adaptively as time integration advances.

Both (21) and (22) are explicit step-truncation methods. These schemes are
very simple to implement and have proven successful in integrating initial
value problems for a variety of PDEs [20]. However, the combination of
dimensionality, non-linearity and stiffness may introduce time-step
restrictions which could make explicit time integration on tensor manifolds
computationally infeasible. To overcome this problem, we recently developed a
new class of implicit rank-adaptive step-truncation algorithms for temporal
integration PDEs on tensor manifolds [21]. These algorithms are based on an
inexact Newton’s method with tensor train GMRES iterations [57] to solve the
algebraic equation arising from time discretization on the tensor manifold.

To describe these methods, let us consider the standard Euler backward scheme

${\bm{f}}_{k+1}={\bm{f}}_{k}+\Delta t{\bm{G}}({\bm{f}}_{k+1}),$ (23)

and the associated root-finding problem

$\bm{H}_{k}(\bm{f}_{k+1})={\bm{f}}_{k+1}-{\bm{f}}_{k}-\Delta
t{\bm{G}}({\bm{f}}_{k+1})={\bm{0}}.$ (24)

If $\bm{G}$ is linear, e.g., for PDEs of the form (8), this reduces to a
linear inversion problem

$\left({\bm{I}}-\Delta t{\bm{G}}\right){\bm{f}}_{k+1}={\bm{f}}_{k}.$ (25)

To solve the linear system (25) in a tensor format, we apply the TT-GMRES
method proposed by Dolgov in [57]. For the implicit midpoint method, the logic
is identical. We apply the above steps to the linear system

$\left({\bm{I}}-\frac{1}{2}\Delta
t{\bm{G}}\right){\bm{f}}_{k+1}=\left({\bm{I}}+\frac{1}{2}\Delta
t{\bm{G}}\right)\bm{f}_{k}.$ (26)

For the fully nonlinear case, one may couple the TT-GMRES iteration to an
inexact Newton method and retain all convergence properties (see [21, Appendix
A] for further details). In [21] we have recently show that under mild
conditions on the tensor truncation error, the implicit step-truncation Euler
and implicit step-truncation midpoint converge, respectively, with order one
and order two in $\Delta t$.

##  V Burgers-Hopf FDE

Local Truncation Error

Figure 4:  Local truncation errors (calculated as $O(\Delta t^{p+1})$ for and
order $p$ method) of the proposed step-truncation tensor methods. These errors
are computed by comparing one time step of step-truncation method with its
Richardson extrapolation. It is seen that explicit/implicit Euler and
explicit/implicit midpoint step-truncation methods have accuracy of order one
and order two, respectively.

$t=0.0$ $t=0.04$ $t=0.08$ $t=0.12$

MC CDF

ST CDF

MC PDF

Figure 5:  Two joint CDF of $u(0,t)$ and $u(\pi,t)$. The CDF is computed by
generating numerical solutions to (31) for $N=20$ using the proposed step-
truncation tensor methods, and then marginalizing the solution in the
remaining 18 variables. We also show a Monte-Carlo estimate of the joint CDF
obtained by sampling $5\times 10^{6}$ solutions to (27).

Consider Burger’s equation on the unit circle $\Omega=[0,2\pi]$

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
x}=\gamma\frac{\partial^{2}u}{\partial x^{2}}+R(u),$ (27)

with a reaction term $R(u)=\alpha\sin(\pi u)$ and random ($2\pi$-periodic)
initial condition $u(x,0)$. Let

$\Phi([\theta(x)],t)=\mathbb{E}\left\\{\exp\left(i\int_{0}^{2\pi}u(x,t)\theta(x)dx\right)\right\\}$
(28)

be the characteristic functional of the stochastic solution to (27). It is
straightforward to show that $\Phi([\theta(x)],t)$ satisfies the following
Burgers-Hopf FDE

$\displaystyle\frac{\partial\Phi([\theta],t)}{\partial t}=$
$\displaystyle\int_{0}^{2\pi}\Biggl{[}\frac{i}{2}\frac{\partial}{\partial
x}\frac{\delta^{2}\Phi([\theta],t)}{\delta\theta(x)^{2}}+\gamma\frac{\partial^{2}}{\partial
x^{2}}\frac{\delta\Phi([\theta],t)}{\delta\theta(x)}$
$\displaystyle+\alpha\sum_{k=1}^{\infty}\frac{(-\pi)^{k}}{(2k+1)!i^{2k}}\frac{\delta^{2+1}\Phi([\theta],t)}{\delta\theta(x)^{2k+1}}\Biggl{]}\theta(x)dx,$
(29)

which admits a unique solution due to the uniqueness of the solution to the
weak form of the Burgers equation (see Section II).

By approximating the FDE (29) in terms of a high-dimensional PDEs of the form
of (8) and taking the inverse Fourier transform (in the sense of
distributions) yields the Liouville-type linear hyperbolic conservation law

$\displaystyle\frac{\partial p}{\partial
t}=-\sum_{i=1}^{N}\frac{\partial}{\partial u_{j}}\Biggl{[}$
$\displaystyle\sum_{j=1}^{N}\Biggl{(}-u_{i}D_{ij}^{(1)}u_{j}+$
$\displaystyle\quad\delta_{ij}R(u_{j})+\gamma
D_{ij}^{(2)}u_{j}\Biggl{)}p\Biggl{]},$ (30)

where $\delta_{ij}$ is the Kronecker delta and $p(u_{1},u_{2},\dots u_{N},t)$
is the joint probability density function of the solution to (27) at all
evenly-spaced grid points $x_{j}=2\pi(j-1)/N$ ($j=1,\ldots,N$), i.e.,
$u_{j}=u(x_{j},t)$.

In (30), $D_{ij}^{(1)}$ and $D_{ij}^{(2)}$ are differentiation matrices. In
our numerical experiments, we use the second-order finite difference matrix on
a periodic domain. From here, it is straightforward to write down a step-
truncation tensor method to solve the high-dimensional PDE (30).

However, this may not yield an accurate scheme due to the properties of the
solution. In fact, it is well-known that PDF equations of the form (30) can
develop shocks [72], or can have solutions converging to Dirac delta functions
if the time asymptotic solution is deterministic. In these cases, it is
convenient to first transform the kinetic equation into a cumulative
distribution function (CDF) equation via an integral transform [73], and then
solve such CDF equation. This allows us to transform PDF shocks into CDF
ramps, and PDF Dirac deltas into CDF shocks. In other words, the CDF equations
produces solutions with features that are well-known to the numerical analysis
community (shocks, ramps and solutions with range bounded in $[0,1]$) and
therefore can be handled with numerical schemes for hyperbolic conservation
laws [74].

Figure 6:  Left: Highest rank of the Tensor Train cores for tensor solutions
of the CDF to equation (31) computed with the explicit step-truncation
midpoint method for dimension $N=\\{8,12,16,20\\}$. Right: Convergence of the
solution to the CDF equation (31) the Burgers-Hopf FDE as we increase the
dimension $N$ from $8$ to $20$. The error plot represents represents the
$L^{2}$ error between the joint CDFs shown in Figure 5.

By integrating (30) with respect to the phase variables
$\\{u_{1},\ldots,u_{N}\\}$ and using we obtain the following CDF equation

$\displaystyle\frac{\partial F}{\partial t}$
$\displaystyle=\sum_{i,j=1}^{N}(u_{i}D_{ij}^{(1)}u_{j}-\gamma
D_{ij}^{(2)}u_{j}-\delta_{ij}R(u_{j}))\frac{\partial F}{\partial u_{i}}$
$\displaystyle\quad-\sum_{\begin{subarray}{c}{i,j=1}\\\ {i\neq
j}\end{subarray}}^{N}(D_{ij}^{(1)}u_{i}-\gamma
D_{ij}^{(2)})\int_{-\infty}^{u_{j}}\frac{\partial F}{\partial u_{i}}du_{j}$
$\displaystyle\quad-\alpha\sum_{j=1}^{N}\sin(\pi u_{j})\frac{\partial
F}{\partial u_{j}},$ (31)

where

$F(u_{1},\ldots,u_{n},t)=\int_{-\infty}^{u_{1}}\cdots\int_{-\infty}^{u_{n}}p(v_{1},\ldots,v_{n},t)dv_{1}\cdots
dv_{n}.$ (32)

With the CDF $F(u_{1},\ldots,u_{N},t)$ available it is straightforward to
compute the PDF

$p(u_{1},\ldots,u_{N},t)=\frac{\partial^{n}F(u_{1},\ldots,u_{N},t)}{\partial
u_{1}\cdots\partial u_{N}}.$ (33)

To apply the proposed step-truncation tensor methods to equation (31), we
discretize the phase space in terms of a $N$-dimensional hypercube centered at
the origin which covers the support of $p(u_{1},u_{2},\dots u_{N},t)$. For the
specific application at hand such hypercube is $\Sigma=[-1.3,1.3]^{N}$. On the
boundary of $\Sigma$, it is straightforward to show that the derivative of $F$
normal to the hypercube boundary is zero. Therefore we may apply outflow
boundary conditions, i.e., set $\partial F/\partial u_{j}=0$ at the boundary
of $\Sigma$. We discretize each variable $u_{i}$ in $\Sigma$ on an evenly-
spaced grid with $n_{i}=64$ points. Also, we approximate the integrals in (31)
using the trapezoidal rule, and the partial derivatives with centered second-
order finite differences. For the initial condition, we set the CDF to be a
product of $N$ zero-mean Gaussian CDFs with standard deviation $\sigma=0.25$,
i.e.,

$F(u_{1},\ldots,u_{N},0)=\prod_{k=1}^{N}g(u_{k}),$ (34)

where

$g(u)=\frac{1}{2}\left[1+\text{erf}\left(\frac{u}{\sigma\sqrt{2}}\right)\right]$
(35)

and $\text{erf}(x)$ is the standard error function.

###  V.1 Numerical results

To solve the high-dimensional CDF equation (31) approximating the Burgers-Hopf
FDE (29) we developed a high-performance (parallel) tensor code implementing
the step truncation tensor methods discussed in Section IV.1. The code can be
found at the GitHub repository: https://github.com/akrodger/paratt.

We first perform a numerical study to verify temporal accuracy of order 1 and
order 2 of the proposed step-truncation tensor methods applied to the high-
dimensional CDF equation (31) (see Figure 3). To stabilize explicit step-
truncation methods, we add a numerical diffusion term proportional to the
round-off error in $\Delta t$ to the right hand side of (31), in a similar
manner to the Lax-Wendroff method [75]. To guarantee that we are accurately
capturing second and first order errors in time, we must ensure that the semi-
discrete PDE we are solving for our reference solution matches that of our
step-truncation integrators. To this end, we compute the local truncation
errors of each scheme by computing its difference with a Richardson
extrapolation (see [76, II.4]) of it to 1 order higher, i.e., we take two
steps of size $\Delta t/2$ and combine it with a step of size $\Delta t$ to
create a more accurate estimate for the purposes of comparing it to a single
step of size $\Delta t$. To this end, we calculate the local truncation errors
by computing the difference between the prediction of each scheme and a
Richardson extrapolation of one order higher (as outlined in [76, II.4]). In
other words, we take two steps of size $\Delta t/2$ and integrate them with a
step of size $\Delta t$ to generate a more precise estimation for comparison
with a single step of size $\Delta t$. The results are shown in Figure 4. It
is seen that the explicit/implicit Euler and explicit/implicit midpoint step-
truncation integrators have accuracy of order one and order two, respectively.
To quantify the errors, we constructed a benchmark CDF via Monte-Carlo
simulation. In practice, we sampled $5\times 10^{6}$ initial conditions
$u(x,0)$ from the Gaussian distribution described above, computed the solution
$u(x,t)$ of (27) using the Fourier pseudospectral method [77], and then used a
kernel density estimator (KDE) [78] to obtain the one-point and the two-points
CDFs. This allows us to study convergence of the solutions to the CDF equation
(31) as we increase the dimension $N$ of the PDE. This convergence study is
crucial for establishing convergence of the solution of the CDF equation (31)
to the solution of the Burgers-Hopf FDE. For this analysis, we selected the
explicit step-truncation midpoint tensor method as the temporal integrator,
due to its favorable performance and ease of implementation.

$t=0.0$ $t=0.04$

$t=0.08$ $t=0.12$

Figure 7: CDF of $u(0,t)$ computed by solving equation (31) with the proposed
step-truncation tensor method. It is seen that as we increase the dimension
$N$ (number of independent variables in the CDF equation (31)) the tensor
solution converges to the Monte-Carlo benchmark (black continuous line).

In Figure 5 plot the joint CDF of $u(0,t)$ and $u(\pi,t)$ at
$t=\\{0,0.04,0.08,0.012\\}$ computed with Monte Carlo, and the proposed step-
truncation midpoint tensor method applied to a 20-dimensional ($N=20$) CDF
equation (31). Clearly there is visual agreement between the results form the
tensor simulation and the results form Monte-Carlo simulation. In Figure 6 we
demonstrate convergence of the solution to the CDF equation (31) the Burgers-
Hopf FDE as we increase the dimension $N$ from $8$ to $20$. The error plot in
Figure 6 represents the $L^{2}$ error between the joint CDFs shown in Figure
5. Runtime increases considerably with dimension though, especially for
moderately high solution ranks. For example once rank surpasses $100$, the
$N=20$ case requires approximately 45 minutes per time step on an Intel
i9-7980xe workstation. The storage of the CDF tensor solution in 20 dimensions
requires approximately 690 megabytes per temporal snapshot for tensor ranks of
the order of $150$. This is a significant compression of a multidimensional
double precision floating point array of size $64^{20}$, which would normally
require approximately $10^{22}$ petabytes. That is to say that we capture
$O(10^{-2})$ accuracy of the data while achieving a data compression ratio of
$1.44\times 10^{26}\%$ and data space savings of approximately $100\%$. Of
course, 690 megabytes is not negligible on practical computing systems. It is
essentially nothing compared to the decompressed data though. In Figure 7 we
compare the CDF of $u(0,t)$ computed by solving the CDF equation (31) with
$N=\\{8,12,16,20\\}$ and then marginalizing, with a benchmark CDF obtained by
Monte-Carlo simulation of (27). It is seen that as we increase the dimension
$N$ the marginal of the CDF tensor solution approaches the Monte-Carlo
benchmark. While the proposed numerical tensor schemes are accurate and
provably convergent, they may not preserve important properties of the CDF
solution, such as positivity, monotonicity, and range in $[0,1]$. The
development of a structure-preserving tensor integration schemes is indeed a
rather unexplored research area. One possibility is to enforce structure at
the level of the tensor truncation operation, i.e., perform rounding
conditional to some required property. For instance, if we are interested in
enforcing non-negativity in a step-truncation tensor scheme, then we can
simply replace the truncation operator based on recursive QR decomposition
(see Appendix A), with the distributed non-negative tensor truncation
discussed in [79]. Another approach to enforce structure can be built directly
at the level of the tensor equations. This approach was recently proposed by
Einkemmer et al. in [80, 81, 82] in the context of Vlasov and Vlasov-Poisson
kinetic equations. The resulting dynamical low-rank algorithm conserves mass,
momentum, and energy.

##  VI Summary

The numerical simulation of functional differential equations (FDEs) remains
an open problem in many respects. In this paper, we developed approximation
theory and new high-performance computational algorithms to solve FDEs on
tensor manifolds. Our approach involves initially approximating the given FDE
with a high-dimensional partial differential equation (PDE) using the
functional approximation methods described in [2, 1, 19]. Subsequently, we
compute the numerical solution of such PDE using parallel rank-adaptive step-
truncation methods [20, 21]. This is a novel and effective strategy for
computing numerical solutions to FDEs, which paves the way for further
research and practical applications across various branches of physics and
engineering. We demonstrated the effectiveness of the proposed new approach
through its application to the Burgers-Hopf FDE, which governs the
characteristic functional of the stochastic solution to the Burgers equation
evolving from a random initial state. Our numerical results and convergence
studies demonstrate that, despite their computational complexity, tensor
methods allow us to compute accurate solutions to the Burgers-Hopf FDE.

###### Acknowledgements.

This research was supported by the U.S. Army Research Office grant
W911NF1810309. Dr. Rodgers would like to acknowledge the Transformational
Tools and Technologies $(\text{T}^{3})$ Project for partially funding his
Ph.D. research through the NASA Pathways internship program.

##  Appendix A Tensor algorithms

In this Appendix, we describe the high-performance (parallel) tensor
algorithms we developed to solve high-dimensional PDEs on tensor manifolds, in
particular the CDF equation discussed in Section V. To this end, we begin with
a brief review of the algebraic representation of the tensor train format
[83].

###  A.1 Tensor train decomposition

A tensor in ${\bm{f}}\in{\mathbb{R}}^{n_{1}\times n_{2}\times\cdots\times
n_{d}}$ is in a tensor train (TT) format if there is an array of positive
integers (called the TT rank)
${\bm{r}}=(r_{0},r_{1},r_{2},\dots,r_{d-1},r_{d})$ with $r_{0}=1=r_{d}$ and a
list of order 3 tensors (called the TT cores)
${\bm{C}}=({\bm{C}}_{1},{\bm{C}}_{2},\dots,{\bm{C}}_{d})$ where
${\bm{C}}_{k}\in{\mathbb{R}}^{r_{k-1}\times n_{k}\times r_{k}}$ such that the
entries of ${\bm{f}}$ may be written as the iterated matrix product

$\displaystyle{\bm{f}}[i_{1},i_{2},\dots,i_{d}]={\bm{C}}_{1}[1,$
$\displaystyle i_{1},:]{\bm{C}}_{2}[:,i_{2},:]\cdots$
$\displaystyle{\bm{C}}_{d-1}[:,i_{d-1},:]{\bm{C}}_{d}[:,i_{d},1],$ (36)

where the colon operator indicates multidimensional array slicing, akin to its
usage in MATLAB. Though rather involved at first glance, the above expression
may be derived by writing down a multivariate function series expansion by the
method of separation of variables, rearranging the expression into a finite
sequence of infinite matrix products, truncating the series so that the matrix
products are finite, then discretizing in space so that the multivariate
function is sampled on a tensor product grid. From this perspective, it
becomes apparent that the tensor train cores represent a two- dimensional
array of functions of a single variable and the ranks are the number of
functions present in a series expansion approximation. A more concrete
exploration of this perspective is given in [84, 85].

Due the the iterated matrix product definition of the format, it becomes
apparent that a number of arithmetic operations may be represented in the
format by producing a new TT tensor with different ranks. In particular if
${\bm{f}},{\bm{g}}$ have TT core lists ${\bm{C}},{\bm{D}}$ then by a simple
inductive argument, we have

$\displaystyle{\bm{f}}$
$\displaystyle[i_{1},i_{2},\dots,i_{d}]+{\bm{g}}[i_{1},i_{2},\dots,i_{d}]=$
$\displaystyle\quad\begin{bmatrix}{\bm{C}}_{1}[1,i_{1},:]\ \big{|}\
{\bm{D}}_{1}[1,i_{1},:]\end{bmatrix}\begin{bmatrix}{\bm{C}}_{2}[:,i_{2},:]&\big{|}&{\bm{0}}\\\
{\bm{0}}&\big{|}&{\bm{D}}_{2}[:,i_{2},:]\end{bmatrix}$
$\displaystyle\cdots\begin{bmatrix}{\bm{C}}_{d-1}[:,i_{d-1},:]&\big{|}&{\bm{0}}\\\
{\bm{0}}&\big{|}&{\bm{D}}_{d-1}[:,i_{d-1},:]\end{bmatrix}\begin{bmatrix}{\bm{C}}_{1}[:,i_{d},1]\\\
{\bm{D}}_{1}[:,i_{d},1]\end{bmatrix}.$ (37)

Data: A tensor $\bm{f}$ in TT format with cores
$({\bm{C}}_{1},{\bm{C}}_{2},\dots,{\bm{C}}_{d})$.

Result: A tensor ${\bm{g}}={\bm{f}}$ in TT format with left orthogonal cores.

Reserve memory for each core ${\bm{D}}_{c}$ with sizes of ${\bm{C}}_{c}$.

${\bm{D}}_{1}={\bm{C}}_{1}$

for  _$c=1,2,\dots,d-1$_ do

$[{\bm{Q}}_{c},{\bm{R}}_{c}]=QR({\mathtt{V}}({\bm{D}}_{c}))$

${\bm{H}}_{c+1}={\bm{R}}_{c}{\mathtt{H}}({\bm{C}}_{c+1})$

${\bm{D}}_{c}[:]={\bm{Q}}_{c}[:]$

${\bm{D}}_{c+1}[:]={\bm{H}}_{c+1}[:]$

end for

Set the cores of $\bm{g}$ as $({\bm{D}}_{1},{\bm{D}}_{2},\dots,{\bm{D}}_{d})$

Algorithm 1 Left Orthogonalization

Thus the sum ${\bm{w}}={\bm{f}}+{\bm{g}}$ is a TT tensor with cores obtained
by concatenating those of ${\bm{f}}$ and ${\bm{g}}$. Scalar multiplication is
straightforward, just scale any of the cores of the tensor.

$\displaystyle\alpha{\bm{f}}[i_{1},i_{2},\dots,i_{d}]=(\alpha{\bm{C}}_{1}$
$\displaystyle[1,i_{1},:]){\bm{C}}_{2}[:,i_{2},:]\cdots$
$\displaystyle{\bm{C}}_{d-1}[:,i_{d-1},:]{\bm{C}}_{d}[:,i_{d},1].$ (38)

In order to perform more sophisticated operations such as tensor truncations,
we must frequently reshape the tensor cores into lower dimensional arrays. Two
particularly useful reshapings are the horizontal flattening

$\displaystyle{\mathtt{H}}:{\mathbb{R}}^{a\times b\times c}$
$\displaystyle\longrightarrow{\mathbb{R}}^{a\times(bc)}$
$\displaystyle{\bm{C}}[i,j,k]$ $\displaystyle\longmapsto{\bm{H}}[i,j+bk],$

and the vertical flattening,

$\displaystyle{\mathtt{V}}:{\mathbb{R}}^{a\times b\times c}$
$\displaystyle\longrightarrow{\mathbb{R}}^{(ab)\times c}$
$\displaystyle{\bm{C}}[i,j,k]$ $\displaystyle\longmapsto{\bm{V}}[i+aj,k].$

The above maps are two dimensional array analogs of the vectorization of a
tensor, in which we have the coordinates are listed as
${\bm{C}}[i,j,k]={\bm{v}}[i+aj+abk]$. Note that when the entries of a tensor
are stored contiguously in computer memory, no memory movement or copying
needs to be done to interpret a tensor as its flattening or vectorization. To
undo the flattening in an algorithm, we denote copying all the entries as
${\bm{C}}[:]={\bm{H}}[:]$ or ${\bm{C}}[:]={\bm{V}}[:]$.

Data: A tensor $\bm{f}$ in TT format with cores
$({\bm{C}}_{1},{\bm{C}}_{2},\dots,{\bm{C}}_{d})$ and a desired accuracy
$\varepsilon$.

Result: A tensor ${\bm{g}}$ so that $\|{\bm{g}}-{\bm{f}}\|\leq\varepsilon$ in
TT format with right orthogonal cores.

Reserve memory for each core ${\bm{D}}_{c}$ with sizes of ${\bm{C}}_{c}$.

Set $({\hat{\bm{C}}}_{1},{\hat{\bm{C}}}_{2},\dots,{\hat{\bm{C}}}_{d})$ as
orthogonal by Alg. 1.

${\hat{\varepsilon}}=\varepsilon/\sqrt{d-1}$

${\bm{D}}_{1}={\hat{\bm{C}}}_{1}$

for  _$c=d,d-1,\dots,2$_ do

$[{\bm{U}}_{c},{\bm{\Sigma}}_{c},{\bm{V}}_{c}^{\top}]=SVD({\mathtt{H}}({\bm{D}}_{c}),{\hat{\varepsilon}})$
(Truncated SVD)

${\bm{H}}_{c-1}={\mathtt{V}}({\hat{\bm{C}}}_{c-1}){\bm{U}}_{c}{\bm{\Sigma}}_{c}$

${\bm{D}}_{c}[:]={\bm{U}}_{c}[:]$

${\bm{D}}_{c-1}[:]={\bm{H}}_{c-1}[:]$

end for

Set the cores of $\bm{g}$ as $({\bm{D}}_{1},{\bm{D}}_{2},\dots,{\bm{D}}_{d})$

Algorithm 2 Truncation of a TT tensor

The general procedure of truncating the ranks of a TT tensor comes in two
phases. The first phase is to collect all the norm of the tensor into a single
core through a sequence of matrix factorizations, leaving all other cores to
have Frobenius norm 1. Then a sequence of SVDs are applied, multiplying a non-
unitary factor from a tensor train core to its neighboring cores. The first
phase is called orthogonalization of a TT tensor and we now present an
algorithm for it.

There are two variants of orthogonalization we will discuss. One is Left-to-
Right orthogonalization and the other Right-to-Left orthogonalization. They
are the same algorithm mirrored across the middle core of the tensor train and
each step in one is a transpose of its counterpart in the other. A TT tensor
is left orthogonal if for $c=1,2,\dots,d-1$, we have
${\mathtt{V}}({\bm{C}}_{c})^{\top}{\mathtt{V}}({\bm{C}}_{c})={\bm{I}}$.
Similarly, a TT tensor is right orthogonal if for $c=2,3,\dots,d$, we have
${\mathtt{H}}({\bm{C}}_{c}){\mathtt{H}}({\bm{C}}_{c})^{\top}={\bm{I}}$. One
may transform the cores of a TT tensor to be left orthogonal while the full
tensor is unchanged by applying iterated matrix factorizations. All we must do
is sequentially flatten, orthogonalize by QR decompoosition, then unflatten
each core, until the final one is reached. This process is summarized in
Algorithm 1. For the right orthogonal case, the procedure is largely
identical, though we use LQ factorization, the lower-triangular
orthogonalizing decomposition, instead.

We now describe one algorithm for truncation of a tensor train given an
orthogonalized TT tensor. The process is essentially the same as
orthogonalization, we apply a sequence of orthogonal matrix factorizations to
each core, multiplying one of the terms into that core’s left or right
neighbor. We then replace the current core with the entries of an orthogonal
factor. The difference in this algorithm is that we apply an SVD, rather than
QR or LQ factorization. The term “truncation” comes from truncating the
singular value decomposition series expansion of a matrix up to a desired
tolerance. This truncation is described as the rounding algorithm in [86]. We
present the right orthogonal variant of it here as Algorithm 2. There is no
clear general rule for preferring one left or right variants truncation
algorithm. They have identical computational costs and similar formulations.
It is possible that differing runtime will occur based on row-major versus
column major layout of the matrices, though this is highly dependent on the
specific memory mapping and workspaces used in the SVD, QR, and LQ calls.

###  A.2 Parallel algorithms for tensor train arithmetic

Among the arithmetic operations performed on the tensor train during the
evolution of the solution to a differential equation, the most costly is the
truncation. This is due to the the fact that it requires data access and
editing of every core twice. Due to the sequential nature of the loops in
Algorithms 1 and 2 we see that parallelizing the inner loop would be most
effective. To this end, we follow the approach of [87]. Our extension of the
algorithms presented in [87] can be found at the GitHub repository:
https://github.com/akrodger/paratt. Since the rank of the tensor train is
expected to change frequently during program execution, we must split the
memory stored on a total of $P$ compute nodes in a manner independent of the
tensor rank.

To this end, we introduce a ${P\times d}$ memory partition matrix ${\bm{M}}$
with positive integer values so that for each core
${\bm{C}}_{k}\in\mathbb{R}^{r_{k-1}\times n_{k}\times r_{k}}$, we have
$\sum_{p=1}^{P}{\bm{M}}[p,k]=n_{k}$. This matrix describes a block-tensor
storage layout for the tensor train core list. More precisely, for each core
${\bm{C}}_{k}$, we define an array of core blocks
$({\bm{C}}_{k}^{1},{\bm{C}}_{k}^{2},\dots,{\bm{C}}_{k}^{P})$ with array sizes
defined via ${\bm{C}}_{k}^{p}\in\mathbb{R}^{r_{k-1}\times{\bm{M}}[p,k]\times
r_{k}}$ so that

${\bm{C}}_{k}[i,:,j]=\begin{bmatrix}{\bm{C}}_{k}^{1}[i,:,j]\\\
{\bm{C}}_{k}^{2}[i,:,j]\\\ \vdots\\\ {\bm{C}}_{k}^{P}[i,:,j]\end{bmatrix}.$
(39)

Each distributed memory compute node with index $p$ stores the core list
$({\bm{C}}_{1}^{p},$ ${\bm{C}}_{2}^{p},$ $\dots,$ ${\bm{C}}_{d}^{p})$. This
memory layout is designed so that

$\displaystyle{\mathtt{V}}({\bm{C}}_{k})$
$\displaystyle=\begin{bmatrix}{\mathtt{V}}({\bm{C}}_{k}^{1})\\\
{\mathtt{V}}({\bm{C}}_{k}^{2})\\\ \vdots\\\
{\mathtt{V}}({\bm{C}}_{k}^{P})\end{bmatrix}\quad\text{and}$
$\displaystyle{\mathtt{H}}({\bm{C}}_{k})$
$\displaystyle=\begin{bmatrix}{\mathtt{H}}({\bm{C}}_{k}^{1})&|&{\mathtt{H}}({\bm{C}}_{k}^{2})&|&\cdots&|&{\mathtt{H}}({\bm{C}}_{k}^{P})\end{bmatrix}.$

We may therefore compute a flattening of the cores without any memory
movements or cross-node communications by reinterpreting the array storage
offsets in column major layout. Additionally, sums and scalar multiplications
of tensors may also be computed in parallel without communications.

This parallel layout also allows for straightforward applications of finite
difference stencils. Consider a finite difference stencil which requires $s$
many points where the number of 1 dimensional ghost cells required is
$g=(s-1)/2$. Let ${\bm{D}}_{FD}$ be the linear operator for this stencil. To
apply a partial derivative in variable $k$ to a TT tensor, we replace the
entries of core ${\bm{C}}_{k}$ with the entries of ${\bm{D}}_{FD}$ applied to
equation (39) for each $i,j$. To perform the parallel version of this stencil,
we first perform a nonblocking communication to deliver the ghost cells of
compute node $p$ to their neighboring nodes $p\pm 1$. For nodes $p=1$ or
$p=P$, we instead rely on discrete boundary condition formulas, communicating
again between nodes $p=1$ and $p=P$ if periodic boundary conditions are
required. After boundary condition communications, every processor has an
appropriate local version of the finite difference stencil, which is identical
save for the definition of the ghost cells.

Numerical integration is similar, though instead of sharing ghost cells of
compute node $p$ with compute node $p\pm 1$, we need only pass the data with
compute node $p+1$, due to the one-sided stencil of a cumulative summation
formula.

###  A.3 Parallel tall-skinny QR factorization

In order to truncate a TT tensor with the described block memory layout, we
must perform a distributed memory QR factorization on a matrix of the form

${\bm{A}}=\begin{bmatrix}{\bm{A}}_{1}\\\ {\bm{A}}_{2}\\\ \vdots\\\
{\bm{A}}_{P}\end{bmatrix},$

or equivalently, perform an LQ factorization on the horizontally concatenated
transpose of $\bm{A}$. For mathematical simplicity, we only present the QR
variant, though our GitHub code available at
https://github.com/akrodger/paratt has both variants.

$[{\bm{R}}_{12};{\bm{R}}_{34}]={\bm{Q}}_{1234}{\bm{R}}_{1234}$$[{\bm{R}}_{1};{\bm{R}}_{2}]={\bm{Q}}_{12}{\bm{R}}_{12}$${\bm{A}}_{1}={\bm{Q}}_{1}{\bm{R}}_{1}$${\bm{A}}_{2}={\bm{Q}}_{2}{\bm{R}}_{2}$$[{\bm{R}}_{3};{\bm{R}}_{4}]={\bm{Q}}_{34}{\bm{R}}_{34}$${\bm{A}}_{3}={\bm{Q}}_{3}{\bm{R}}_{3}$${\bm{A}}_{4}={\bm{Q}}_{4}{\bm{R}}_{4}$
Figure 8:  Parallel TSQR algorithm with 4 compute notes. A binary tree is
formed which outlines the communication pattern for calculating the ${\bm{R}}$
factor. At each level, the node with larger ID sends its ${\bm{R}}$ factor to
the tree sibling it connects to. The sibling then does a QR factorization with
a concatenated pair of child ${\bm{R}}$ factors. This process repeats until
the final ${\bm{R}}$ factor is found at the root of the tree stored on
processor $p=1$. We then broadcast this matrix to all other processors. To get
the ${\bm{Q}}$ factor, we store the orthogonal ${\bm{Q}}_{i}$ for each tree
node $i$ and multiply the parent’s ${\bm{Q}}_{i}$ on the right side of the
child’s ${\bm{Q}}$ factor, traversing the tree from the root to the leaves.
This results in a distributed memory orthogonal matrix
${\bm{Q}}=[{\bm{Q}}_{1};{\bm{Q}}_{2};{\bm{Q}}_{3};{\bm{Q}}_{4}].$

One such factorization amenable to this layout is the Tall-Skinny QR (TSQR)
factorization [88]. This algorithm uses a binary tree structure to define a
communication pattern for decomposing the QR factorization of $\bm{A}$ into a
collection of smaller QR factorizations, avoiding communications if possible.
By factoring the blocks into their own QR factorizations, we see that the Q
factor may be presented as a sequence of products of block diagonal matrices.
This process is then applied to a recursive binary tree to track the nesting
of the orthogonal Q factor. The final result is an orthogonal matrix
${\bm{Q}}=[{\bm{Q}}_{1};{\bm{Q}}_{2};\dots;{\bm{Q}}_{P}]$ distributed in the
same block layout as $\bm{A}$ and an upper triangular matrix $\bm{R}$ copied
across all compute nodes. See Figure 8 for a visual representation of the
binary tree storage structure for the case $P=4$. The transposed version of
this algorithm is called the Wide-Fat LQ (WFLQ) factorization and has the same
tree data storage structure, but every matrix factorization is transposed.

###  A.4 Parallel orthogonalization and truncation

We now introduce a parallelization of Algorithms 1 and 2. It can be seen that
parallelizing the orthogonalization algorithms is as simple as replacing the
QR factorizations with their TSQR variants. However, the truncation is
requires a bit of manipulation. In order to parallelize this, we first note
the following relationship of the SVD and QR factorizations,

$\displaystyle{\bm{Q}}{\bm{R}}$
$\displaystyle={\bm{A}}={\bm{U}}{\bm{\Sigma}}{\bm{V}}^{\top}$
$\displaystyle{\bm{R}}$
$\displaystyle={\bm{Q}}^{\top}{\bm{U}}{\bm{\Sigma}}{\bm{V}}^{\top}.$

Data: A tensor $\bm{f}$ in distributed memory TT format with cores
$([{\bm{C}}_{1}^{p}],[{\bm{C}}_{2}^{p}],\dots,[{\bm{C}}_{d}^{p}])$.

Result: A tensor ${\bm{g}}={\bm{f}}$ in distributed memory TT format with left
orthogonal cores.

Reserve memory for each block ${\bm{D}}_{c}^{p}$ with sizes of
${\bm{C}}_{c}^{p}$.

${\bm{D}}_{1}^{p}={\bm{C}}_{1}^{p}$

for  _$c=1,2,\dots,d-1$_ do

$[{\bm{Q}}_{c}^{p},{\bm{R}}_{c}]=TSQR({\mathtt{V}}({\bm{D}}_{c}^{p}))$

${\bm{H}}_{c+1}^{p}={\bm{R}}_{c}{\mathtt{H}}({\bm{C}}_{c+1}^{p})$

${\bm{D}}_{c}^{p}[:]={\bm{Q}}_{c}^{p}[:]$

${\bm{D}}_{c+1}^{p}[:]={\bm{H}}_{c+1}^{p}[:]$

end for

Set the cores of $\bm{g}$ as
$([{\bm{D}}_{1}^{p}],[{\bm{D}}_{2}^{p}],\dots,[{\bm{D}}_{d}^{p}])$

Algorithm 3 Parallel Left Orthogonalization

Thus ${\bm{R}}$ has the same singular values as ${\bm{A}}$, and so truncating
an SVD of ${\bm{R}}$ produces the same approximation error as truncating
${\bm{A}}$ directly. Since the TSQR algorithm is designed to produce a copy of
${\bm{R}}$ on each compute node, we apply a shared memory SVD redundantly on
this copied ${\bm{R}}$ to determine the truncated tensor. This process is
formalized in Algorithms 3 and 4. These algorithms are designed so that the
only communications required are those that are involved in the computation of
the QR or LQ factorizations. Each algorithm is mathematically equivalent to
its serial counterpart. The core difference is the use of a superscript
$p=1,2,\dots,P$, to denote which matrices take on different values within the
different compute nodes. Matrices which are numerically identical across all
nodes lack this superscript.

Data: A tensor $\bm{f}$ in distributed memory TT format with cores
$([{\bm{C}}_{1}^{p}],[{\bm{C}}_{2}^{p}],\dots,[{\bm{C}}_{d}^{p}])$ and a
desired accuracy $\varepsilon$.

Result: A tensor ${\bm{g}}$ so that $\|{\bm{g}}-{\bm{f}}\|\leq\varepsilon$ in
TT format with right orthogonal cores.

Reserve memory for each block ${\bm{D}}_{c}^{p}$ with sizes of
${\bm{C}}_{c}^{p}$.

Set
$([{\hat{\bm{C}}}_{1}^{p}],[{\hat{\bm{C}}}_{2}^{p}],\dots,[{\hat{\bm{C}}}_{d}^{p}])$
as orthogonal by Alg. 3.

${\hat{\varepsilon}}=\varepsilon/\sqrt{d-1}$

${\bm{D}}_{d}^{p}={\hat{\bm{C}}}_{d}^{p}$

for  _$c=d,d-1,\dots,2$_ do

$[{\bm{L}}_{c},{\bm{Q}}_{c}^{p}]=WFLQ({\mathtt{H}}({\bm{D}}_{c}^{p}))$

$[{\bm{U}}_{c},{\bm{\Sigma}}_{c},{\bm{V}}_{c}^{\top}]=SVD({\bm{L}}_{c},{\hat{\varepsilon}})$
(Truncated SVD)

${\bm{H}}_{c-1}^{p}={\mathtt{V}}({\hat{\bm{C}}}_{c-1}^{p}){\bm{U}}_{c}{\bm{\Sigma}}_{c}$

${\bm{\hat{Q}}}_{c}^{p}={\bm{V}}_{c}^{\top}{\bm{Q}}_{c}^{p}$

${\bm{D}}_{c}^{p}[:]={\bm{\hat{Q}}}_{c}^{p}[:]$

${\bm{D}}_{c-1}^{p}[:]={\bm{H}}_{c-1}^{p}[:]$

end for

Set the cores of $\bm{g}$ as
$([{\bm{D}}_{1}^{p}],[{\bm{D}}_{2}^{p}],\dots,[{\bm{D}}_{d}^{p}])$

Algorithm 4 Parallel TT Truncation

## References

  * Venturi [2018] D. Venturi, The numerical approximation of nonlinear functionals and functional differential equations, Physics Reports 732, 1 (2018).
  * Venturi and Dektor [2021] D. Venturi and A. Dektor, Spectral methods for nonlinear functionals and functional differential equations, Res. Math. Sci. 8, 1 (2021).
  * Hopf [1952] E. Hopf, Statistical hydromechanics and functional calculus, J. Rat. Mech. Anal. 1, 87 (1952).
  * Monin and Yaglom [2007] A. S. Monin and A. M. Yaglom,  _Statistical Fluid Mechanics, Volume II: Mechanics of Turbulence_ (Dover, 2007).
  * Ohkitani [2020] K. Ohkitani, Study of the Hopf functional equation for turbulence: Duhamel principle and dynamical scaling, Phys. Rev. E 101, 013104 (2020).
  * Rosen [1971] G. Rosen, Functional calculus theory for incompressible fluid turbulence, J. Math. Phys. 12, 812 (1971).
  * Vakhania  _et al._ [1987] N. Vakhania, V. Tarieladze, and S. Chobanyan,  _Probability distributions on Banach spaces_ , 1st ed. (Springer, 1987).
  * Foias [1973] C. Foias, Statistical study of Navier-Stokes equations, part I, Rend. Sem. Mat. Univ. Padova 48, 219 (1973).
  * Peskin and Schroede [2018] M. E. Peskin and D. V. Schroede,  _An introduction to quantum field theory_ (CRC Press, 2018).
  * Zinn-Justin [2002] J. Zinn-Justin,  _Quantum field theory and critical phenomena_ , 4th ed. (Oxford Univ. Press, 2002).
  * Martin  _et al._ [1973] P. C. Martin, E. D. Siggia, and H. A. Rose, Statistical dynamics of classical systems, Phys. Rev. A 8, 423 (1973).
  * Phythian [1977] R. Phythian, The functional formalism of classical statistical dynamics, J. Phys A: Math. Gen. 10, 777 (1977).
  * Jensen [1981] R. V. Jensen, Functional integral approach to classical statistical dynamics, J. Stat. Phys. 25, 183 (1981).
  * Phythian [1975] R. Phythian, The operator formalism of classical statistical dynamics, J. Phys A: Math. Gen. 8, 1423 (1975).
  * Jouvet and Phythian [1979] B. Jouvet and R. Phythian, Quantum aspects of classical and statistical fields, Phys. Rev. A 19, 1350 (1979).
  * Ruthotto  _et al._ [2020] L. Ruthotto, S. Osher, W. Li, L. Nurbekyan, and S. W. Fung, A machine learning framework for solving high-dimensional mean field game and mean field control problems, PNAS 117, 9183 (2020).
  * E  _et al._ [2019] W. E, J. Han, and Q. Li, A mean-field optimal control formulation of deep learning, Res. Math. Sci. 6, 10 (2019).
  * Chow  _et al._ [2019] Y. T. Chow, W. Li, S. Osher, and W. Yin, Algorithm for Hamilton–Jacobi equations in density space via a generalized Hopf formula, J. Sci. Comp. 80, 1195 (2019).
  * Gangbo  _et al._ [2021] W. Gangbo, S. Mayorga, and A. Swiech, Finite dimensional approximations of Hamilton-Jacobi-Bellman equations in spaces of probability measures, SIMA. J. Math. Anal. 53, 1320 (2021).
  * Rodgers  _et al._ [2022] A. Rodgers, A. Dektor, and D. Venturi, Adaptive integration of nonlinear evolution equations on tensor manifolds, J. Sci. Comput. 92, 1 (2022).
  * Rodgers and Venturi [2023] A. Rodgers and D. Venturi, Implicit integration of nonlinear evolution equations on tensor manifolds, J. Sci. Comput 97, 1 (2023).
  * Berselli and Spirito [2021] L. Berselli and S. Spirito, On the existence of Leray-Hopf weak solutions to the Navier-Stokes equations, Fluids 6, 6010042 (2021).
  * Prodi [1959] G. Prodi, Un teorema di unicità per le equazioni di Navier–Stokes, Annali di Matematica 48, 173 (1959).
  * Leray [1934] J. Leray, Sur le mouvement d’un liquide visqueux emplissant l’espace, Acta Math. 63, 193 (1934).
  * Fabes  _et al._ [1972] E. Fabes, B. Jones, and N. Rivière, The initial value problem for the navier-stokes equations with data in $L_{p}$, Arch. Rational Mech. Anal. 45, 222 (1972).
  * Temam [1984] R. Temam,  _Navier-Stokes equations: theory and numerical analysis_ (AMS Chelsea Publishing, 1984).
  * Vishik and Fursikov [1988] M. J. Vishik and A. V. Fursikov,  _Mathematical problems of statistical hydromechanics_ , 2nd ed. (Kluwer Academic Publishers, 1988).
  * Gishlarkaev [2010] V. I. Gishlarkaev, Uniqueness of a solution to the Cauchy problem for the Hopf equation in the two-dimensional case, Journal of Mathematical Sciences 169, 64 (2010).
  * Barker [2018] T. Barker, Uniqueness results for weak Leray-Hopf solutions of the Navier-Stokes system with initial values in critical spaces, J. Math. Fluid Mech. 20, 133 (2018).
  * Barker [2020] T. Barker, About local continuity with respect to $L_{2}$ initial data for energy solutions of the Navier-Stokes equations, Mathematische Annalen , https://doi.org/10.1007/s00208 (2020).
  * Kato [1984] T. Kato, Strong $L^{p}$-solutions of the Navier-Stokes equation in $\mathbb{R}^{m}$, with applications to weak solutions, Math. Z. 187, 471 (1984).
  * Galdi [2019] G. P. Galdi, On the relation between very weak and Leray-Hopf solutions to Navier–Stokes equations, Proc. Amer. Math. Soc. 147, 5349 (2019).
  * Dubois [2003] S. Dubois, Uniqueness for some Leray–Hopf solutions to the Navier–Stokes, Journal of Differential Equations 189, 99 (2003).
  * Fujita and Kato [1964] H. Fujita and T. Kato, On the Navier–Stokes initial value problem. I, Arch. Rational Mech. Anal. 16, 269 (1964).
  * Buckmaster and Vicol [2019] T. Buckmaster and V. Vicol, Nonuniqueness of weak solutions to the Navier–Stokes equation, Annals of Mathematics 189, 101 (2019).
  * Kang and Protas [2022] D. Kang and B. Protas, Searching for singularities in Navier–Stokes flows based on the Ladyzhenskaya–Prodi–Serrin conditions, J Nonlinear Sci 32, 1 (2022).
  * Serrin [1963] J. Serrin, The initial value problem for the Navier–Stokes equations, Proc. Sympos., Madison, Wis. , 69 (1963).
  * Deriaz and Perrier [2006] E. Deriaz and V. Perrier, Divergence-free and curl-free wavelets in two dimensions and three dimensions: application to turbulent flows, Journal of Turbulence 7, 1 (2006).
  * Deriaz and Perrier [2008] E. Deriaz and V. Perrier, Direct numerical simulation of turbulence using divergence-free wavelets, Multiscale Model. Simul. 7, 1101 (2008).
  * Fuselier and Wright [2017] E. J. Fuselier and G. B. Wright, A radial basis function method for computing Helmholtz–Hodge decompositions, IMA J. Numer. Anal. 37, 774 (2017).
  * Sacchi-Landriani and Vandeven [1989] G. Sacchi-Landriani and H. Vandeven, Polynomial approximation of divergence-free functions, Mathematics of Computation 185, 103 (1989).
  * Engel and Nagel [1999] K.-J. Engel and R. Nagel,  _One-parameter semigroups for linear evolution equations_ , Vol. 194 (Springer, 1999).
  * Guidetti  _et al._ [2004] D. Guidetti, B. Karasozen, and S. Piskarev, Approximation of abstract differential equations, Journal of Mathematical Sciences 122, 3013 (2004).
  * Boelens  _et al._ [2020] A. M. P. Boelens, D. Venturi, and D. M. Tartakovsky, Tensor methods for the Boltzmann-BGK equation, J. Comput. Phys. 421, 109744 (2020).
  * di Marco and Pareschi [2014] G. di Marco and L. Pareschi, Numerical methods for kinetic equations, Acta Numerica 23, 369 (2014).
  * Bungartz and Griebel [2004] H. J. Bungartz and M. Griebel, Sparse grids, Acta Numerica 13, 147 (2004).
  * Barthelmann  _et al._ [2000] V. Barthelmann, E. Novak, and K. Ritter, High dimensional polynomial interpolation on sparse grids, Advances in Computational Mechanics 12, 273 (2000).
  * Narayan and Jakeman [2014] A. Narayan and J. Jakeman, Adaptive Leja sparse grid constructions for stochastic collocation and high-dimensional approximation, SIAM J. Sci. Comput. 36, A2952 (2014).
  * Raissi and Karniadakis [2018] M. Raissi and G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, J. Comput. Phys. 357, 125 (2018).
  * Raissi  _et al._ [2019] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378, 606 (2019).
  * Khoromskij [2015] B. N. Khoromskij, Tensor numerical methods for multidimensional PDEs: Theoretical analysis and initial applications, in  _CEMRACS 2013–modelling and simulation of complex systems: stochastic and deterministic approaches_ , ESAIM Proc. Surveys, Vol. 48 (EDP Sci., Les Ulis, 2015) pp. 1–28.
  * Bachmayr  _et al._ [2016] M. Bachmayr, R. Schneider, and A. Uschmajew, Tensor networks and hierarchical tensors for the solution of high-dimensional partial differential equations, Foundations of Computational Mathematics 16, 1423 (2016).
  * Boelens  _et al._ [2018] A. M. P. Boelens, D. Venturi, and D. M. Tartakovsky, Parallel tensor methods for high-dimensional linear PDEs, J. Comput. Phys. 375, 519 (2018).
  * Dektor and Venturi [2023] A. Dektor and D. Venturi, Tensor rank reduction via coordinate flows, J. Comp. Phys. 491, 112378 (2023).
  * Dektor and Venturi [2024] A. Dektor and D. Venturi, Coordinate-adaptive integration of PDEs on tensor manifolds, Comm. Appl. Math. Comput. , https://doi.org/10.1007/s42967 (2024).
  * Cho  _et al._ [2017] H. Cho, D. Venturi, and G. E. Karniadakis, Numerical methods for high-dimensional kinetic equations, in  _Uncertainty quantification for kinetic and hyperbolic equations_ , edited by S. Jin and L. Pareschi (Springer, 2017) pp. 93–125.
  * Dolgov [2013] S. V. Dolgov, TT-GMRES: solution to a linear system in the structured tensor format, Russian Journal of Numerical Analysis and Mathematical Modelling 28, 149 (2013).
  * Beylkin and Mohlenkamp [2002] G. Beylkin and M. J. Mohlenkamp, Numerical operator calculus in higher dimensions, PNAS 99, 10246 (2002).
  * Dektor and Venturi [2021a] A. Dektor and D. Venturi, Dynamic tensor approximation of high-dimensional nonlinear PDEs, J. Comput. Phys. 437, 110295 (2021a).
  * Bigoni  _et al._ [2016] D. Bigoni, A. P. Engsig-Karup, and Y. M. Marzouk, Spectral tensor-train decomposition, SIAM J. Sci. Comput. 38, A2405 (2016).
  * Oseledets [2011a] I. V. Oseledets, Tensor-train decomposition, SIAM J. Sci. Comput. 33, 2295– (2011a).
  * Schneider and Uschmajew [2014] R. Schneider and A. Uschmajew, Approximation rates for the hierarchical tensor format in periodic Sobolev spaces, J. Complexity 30, 56 (2014).
  * Grasedyck and Löbbert [2018] L. Grasedyck and C. Löbbert, Distributed hierarchical SVD in the hierarchical Tucker format, Numer. Linear Algebra Appl. 25, e2174 (2018).
  * Etter [2016] A. Etter, Parallel ALS algorithm for solving linear systems in the hierarchical Tucker representation, SIAM J. Sci. Comput. 38, A2585–A2609 (2016).
  * Uschmajew and Vandereycken [2013a] A. Uschmajew and B. Vandereycken, The geometry of algorithms using hierarchical tensors, Linear Algebra Appl. 439, 133 (2013a).
  * Dektor  _et al._ [2021] A. Dektor, A. Rodgers, and D. Venturi, Rank-adaptive tensor methods for high-dimensional nonlinear PDEs, J. Sci. Comput. 88, 1 (2021).
  * Griebel and Li [2019] M. Griebel and G. Li, On the decay rate of the singular values of bivariate functions, SIAM J. Numer. Anal. 56, 974 (2019).
  * Rodgers and Venturi [2020] A. Rodgers and D. Venturi, Stability analysis of hierarchical tensor methods for time-dependent PDEs, J. Comput. Phys. 409, 109341 (2020).
  * Kato [1995] T. Kato,  _Perturbation theory for linear operators_ , Classics in Mathematics (Springer-Verlag, Berlin, 1995) pp. xxii+619, reprint of the 1980 edition.
  * Uschmajew and Vandereycken [2013b] A. Uschmajew and B. Vandereycken, The geometry of algorithms using hierarchical tensors, Linear Algebra Appl. 439, 133 (2013b).
  * Kressner and Tobler [2014] D. Kressner and C. Tobler, Algorithm 941: htucker – a Matlab toolbox for tensors in hierarchical Tucker format, ACM Transactions on Mathematical Software 40, 1 (2014).
  * Cho  _et al._ [2014] H. Cho, D. Venturi, and G. E. Karniadakis, Statistical analysis and simulation of random shocks in Burgers equation, Proc. R. Soc. A 2171, 1 (2014).
  * Venturi  _et al._ [2013] D. Venturi, D. M. Tartakovsky, A. M. Tartakovsky, and G. E. Karniadakis, Exact PDF equations and closure approximations for advective-reactive transport, J. Comput. Phys. 243, 323 (2013).
  * Hesthaven [2018] J. S. Hesthaven,  _Numerical Methods for Conservation Laws: From Analysis to Algorithm_ , 1st ed. (SIAM, 2018).
  * Lax and Wendroff [1958] P. Lax and B. Wendroff,  _Systems of conservation laws_ , Tech. Rep. (Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 1958).
  * Hairer  _et al._ [1993] E. Hairer, G. Wanner, and S. P. Nørsett,  _Solving Ordinary Differential Equations I, Nonstiff Problems_ , second revised edition. ed., Springer Series in Computational Mathematics,, Vol. 8 (Springer Berlin Heidelberg, Berlin, Heidelberg, 1993).
  * Hesthaven  _et al._ [2007] J. S. Hesthaven, S. Gottlieb, and D. Gottlieb,  _Spectral methods for time-dependent problems_ (Cambridge Univ. Press, 2007).
  * Botev  _et al._ [2010] Z. Botev, J. Grotowski, and D. Kroese, Kernel density estimation via diffusion, Annals of Statistics , 2916 (2010).
  * Bhattarai  _et al._ [2020] M. Bhattarai, G. Chennupati, E. Skau, R. Vangara, H. Djidjev, and B. Alexandrov, Distributed non-negative tensor train decomposition, arXiv:2008.01340 , 1 (2020).
  * Einkemmer and Joseph [2021] L. Einkemmer and I. Joseph, A mass, momentum, and energy conservative dynamical low-rank scheme for the Vlasov equation, J. Comput. Phys. 443, 110495 (2021).
  * Einkemmer and Lubich [2018] L. Einkemmer and C. Lubich, A low-rank projector-splitting integrator for the Vlasov-Poisson equation, SIAM J. Sci. Comput. 40, B1330 (2018).
  * Einkemmer and Lubich [2019] L. Einkemmer and C. Lubich, A quasi-conservative dynamical low-rank algorithm for the Vlasov equation, SIAM J. Sci. Comput. 41, B1061 (2019).
  * Oseledets and Tyrtyshnikov [2010] I. Oseledets and E. Tyrtyshnikov, TT-cross approximation for multidimensional arrays, Linear Algebra and its Applications 432, 70 (2010).
  * Dektor and Venturi [2020] A. Dektor and D. Venturi, Dynamically orthogonal tensor methods for high-dimensional nonlinear PDEs, J. Comput. Phys. 404, 109125 (2020).
  * Dektor and Venturi [2021b] A. Dektor and D. Venturi, Dynamic tensor approximation of high-dimensional nonlinear PDEs, J. Comput. Phys. 437, 110295 (2021b).
  * Oseledets [2011b] I. V. Oseledets, Tensor-train decomposition, SIAM Journal on Scientific Computing 33, 2295 (2011b), https://doi.org/10.1137/090752286 .
  * Daas  _et al._ [2022] H. A. Daas, G. Ballard, and P. Benner, Parallel algorithms for tensor train arithmetic, SIAM J. Sci. Comput. 44, C25 (2022).
  * Demmel  _et al._ [2012] J. Demmel, L. Grigori, M. Hoemmen, and J. Langou, Communication-optimal parallel and sequential QR and LU factorizations, SIAM Journal on Scientific Computing 34, A206 (2012).