# Coordinate-adaptive integration of PDEs on tensor manifolds

Alec Dektor  Daniele Venturi<EMAIL_ADDRESS>Lawrence Berkeley National
Laboratory, Berkeley (CA) 94720, USA.  Department of Applied Mathematics,
University of California Santa Cruz, Santa Cruz (CA) 95064, USA.

###### Abstract

We introduce a new tensor integration method for time-dependent PDEs that
controls the tensor rank of the PDE solution via time-dependent diffeomorphic
coordinate transformations. Such coordinate transformations are generated by
minimizing the normal component of the PDE operator relative to the tensor
manifold that approximates the PDE solution via a convex functional. The
proposed method significantly improves upon and may be used in conjunction
with the coordinate-adaptive algorithm we recently proposed in [14], which is
based on non-convex relaxations of the rank minimization problem and
Riemannian optimization. Numerical applications demonstrating the
effectiveness of the proposed coordinate-adaptive tensor integration method
are presented and discussed for prototype Liouville and Fokker-Planck
equations.

††journal: ArXiv

##  1 Introduction

Developing efficient numerical methods to solve high-dimensional partial
differential equations (PDEs) is a central task in many areas of engineering,
physical sciences and mathematics. Such PDEs are often written in the form of
an abstract initial/boundary value problem

$\begin{cases}\displaystyle\frac{\partial u(\bm{x},t)}{\partial
t}=G_{\bm{x}}(u(\bm{x},t)),\qquad t\in[0,T],\\\
u(\bm{x},0)=u_{0}(\bm{x}),\end{cases}$ (1)

which governs the time evolution of a quantity of interest $u(\bm{x},t)$
(high-dimensional field) over a compact domain $\Omega\subseteq\mathbb{R}^{d}$
($d\gg 1$) and has temporal dynamics generated by the nonlinear operator
$G_{\bm{x}}$. The subscript “$\bm{x}$” in $G_{\bm{x}}$ indicates that the
operator can explicitly depend on the variables $\bm{x}\in\Omega$. For
instance,

$G_{\bm{x}}(u)=\bm{f}(\bm{x})\cdot\nabla u+\nabla\cdot(\sigma(\bm{x})\nabla
u)+R(u),$ (2)

where $R(u)$ is a nonlinear reaction term. The PDE (1) may involve tens,
hundreds, or thousands of independent variables, and arises naturally in a
variety of applications of kinetic theory, e.g., the Liouville equation, the
Fokker-Planck equation, the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) PDF
hierarchy [26, 3, 8, 9] or the Lundgren-Monin-Novikov (LMN) PDF hierarchy [24,
16, 20, 37]. These equations allow us to perform uncertainty quantification in
many different physical systems including non-neutral plasmas, turbulent
combustion, hypersonic flows, and stochastic particle dynamics.

Several general-purpose algorithms, e.g., based on tensor networks [22, 31,
11, 12, 10, 13, 14] and physics-informed machine learning [29, 30, 21], have
recently been proposed to integrate the PDE (1). Tensor networks can be seen
as factorizations of entangled objects such as multivariate functions or
operators into networks of simpler objects which are amenable to efficient
representation and computation. The vast majority of tensor algorithms
currently available to approximate functions, operators and PDEs on tensor
manifolds rely on canonical polyadic (CP) decomposition [4, 6, 7, 11], Tucker
tensors, or tensors corresponding to binary trees such as functional tensor
train (FTT) [5, 13, 27] and hierarchical Tucker tensors [15, 17, 32, 34]. A
compelling reason for using binary tensor trees is that they allow the tensor
expansion to be constructed via spectral theory for linear operators, in
particular the hierarchical Schmidt decomposition [12, 13, 18].

Regardless of the chosen tensor format, the efficiency of tensor algorithms
for high-dimensional PDEs depends heavily on the rank of the solution and the
PDE operator in the chosen format. Indeed, the computational cost of tensor
approximation algorithms scales linearly with the dimension of the system,
e.g., the number of independent variables $d$ in the PDE (1), and polynomially
with the rank. To address the rank-related unfavorable scaling, we recently
developed a new tensor rank reduction algorithm based on coordinate
transformations that can significantly increase the efficiency of high-
dimensional tensor approximation algorithms [14]. Given a multivariate
function or operator, the algorithm determines a coordinate transformation so
that the function, the operator, or the operator applied to the function in
the new coordinate system has smaller tensor rank. In [14] we considered
linear coordinate transformations, which yield a new class of functions that
we called tensor ridge functions. Tensor ridge functions can be written
analytically as

$v_{\bm{s}}(\bm{x},t)=u_{\bm{r}}(\bm{A}\bm{x},t),$ (3)

where $u_{\bm{r}}$ is a given FTT tensor, and $v_{\bm{s}}$ is a FTT tensor
with smaller rank, i.e., $|\bm{s}|\leq|\bm{r}|$. To compute the unknown matrix
$\bm{A}$ (which enables tensor rank reduction), we developed a Riemannian
gradient descent algorithm on some matrix manifold (e.g.,
$\bm{A}\in\text{SL}_{d}(\mathbb{R})$) that minimizes the non-convex cost
functional

$C(\bm{A})=S\left[u_{\bm{r}}(\bm{A}\bm{x},t)\right],$ (4)

where $S$ is the Schatten 1-norm. The numerical results we obtained in the
recent paper [14] suggest that rank reduction based on linear coordinate
transformations can significantly speed up tensor computations for PDEs while
retaining accuracy.

In this paper we propose a new rank-reduction method for (1) based on
coordinate transformations which does not rely on non-convex relaxations of
the rank minimization problem, e.g., the minimization of (4). The main idea is
to select an infinitesimal generator of the coordinate transformation at each
time which minimizes the component of the PDE operator (expressed in intrinsic
coordinates) that is normal to the tensor manifold where we approximate the
PDE solution. This generates a time-dependent nonlinear coordinate
transformation (referred to as a coordinate flow) directly from the PDE via a
convex optimization problem. The Euler-Lagrange equations arising from this
time-dependent variational principle (TDVP) allow us to determine an optimal
coordinate flow that controls/minimizes the tensor rank of the solution during
time integration.

This paper is organized as follows. In Section 2 we develop the theory for
coordinate-adaptive tensor integration of PDEs via rank-reducing coordinate
flows. To this end, we begin with general nonlinear coordinate flows
(diffeomorphisms) and formulate a new convex functional that links the PDE
generator to the geometry of the tensor manifold and the flow (Section 2.1).
We then restrict our attention to linear coordinate flows, i.e., ridge
tensors, and develop the corresponding coordinate-adaptive tensor integration
scheme. In Section 3 we demonstrate the effectiveness of the new tensor
integrators by applying them to prototype Liouville and Fokker-Planck
equations. Our findings are summarized in Section 4.

##  2 Coordinate-adaptive time integration of PDEs on tensor manifolds

We begin by representing the solution $u(\bm{x},t)$ to the PDE (1) with
another function

$v(\bm{y}(\bm{x},t),t)=u(\bm{x},t),$ (5)

defined on a time-dependent curvilinear coordinate system $\bm{y}(\bm{x},t)$
[2, 33, 35] with $\bm{y}(\bm{x},0)=\bm{x}$. It is assumed that
$\bm{y}(\bm{x},t)$ is a diffeomorphism111Recall that if $\bm{y}(\bm{x},t)$ is
a diffeomorphism then there exists a unique differentiable inverse flow
$\bm{x}(\bm{y},t)$ such that $\bm{x}(\bm{y}(\bm{x},t),t)=\bm{x}$ for all
$t\in[0,T]$ ($T<\infty$). which we will refer to as coordinate flow. To
simplify notation we may not explicitly write the dependence of $\bm{y}$ on
$\bm{x}$ (and vice-versa the dependence of $\bm{x}$ on $\bm{y}$), or the
dependence of $\bm{x}$, $\bm{y}$ on $t$. However, it is always assumed that
$\bm{x}$ and $\bm{y}$ are related via a (time-dependent) diffeomorphism. Next,
we define two distinct time derivatives of $v(\bm{y}(\bm{x},t),t)$. The first
represents a change in $v$ at fixed location $\bm{y}$, i.e., a derivative of
$v$ with respect to time with $\bm{y}$ constant, which we denote by $\partial
v(\bm{y}(\bm{x},t),t)/\partial t$. The second represents a change in $v$ along
the coordinate flow $\bm{y}$, i.e., a derivative of $v$ with respect to time
with $\bm{x}$ constant, which we denote by $Dv(\bm{y}(\bm{x},t),t)/Dt$. The
derivative $Dv(\bm{y}(\bm{x},t),t)/Dt$ is known as material derivative [2,
25], or convective derivative [35, 36]. Of course these time derivatives are
related via the equation

$\frac{Dv(\bm{y}(\bm{x},t),t)}{Dt}=\frac{\partial
v(\bm{y}(\bm{x},t),t)}{\partial t}+\frac{\partial\bm{y}(\bm{x},t)}{\partial
t}\cdot\nabla_{\bm{y}}v(\bm{y}(\bm{x},t),t).$ (6)

Since $u(\bm{x},t)$ is a solution to (1), $v$ is related to $u$ by (5) and
$\bm{y}(\bm{x},0)=\bm{x}$, it follows immediately that $v$ satisfies the PDE

$\begin{cases}\displaystyle\frac{Dv(\bm{y},t)}{Dt}=G_{\bm{y}}(v(\bm{y},t)),\qquad
t\in[0,T],\\\ v(\bm{y},0)=u_{0}(\bm{y}),\end{cases}$ (7)

where the operator $G_{\bm{y}}$ can be derived by writing $G_{\bm{x}}$ in
coordinates $\bm{y}$ using standard tools of differential geometry [2, 35,
36]. Combining (6) and (7) we obtain222The PDE (8) represents the evolution
equation for the solution to (7) written along the coordinate flow
$\bm{y}(\bm{x},t)$.

$\begin{cases}\displaystyle\frac{\partial v(\bm{y},t)}{\partial
t}=Q_{\bm{y}}(v(\bm{y},t),\dot{\bm{y}}),\qquad t\in[0,T],\\\
v(\bm{y},0)=u_{0}(\bm{y}),\end{cases}$ (8)

where

$Q_{\bm{y}}(v(\bm{y},t),\dot{\bm{y}})=G_{\bm{y}}(v(\bm{y},t))-\dot{\bm{y}}(\bm{x},t)\cdot\nabla_{\bm{y}}v(\bm{y},t)$
(9)

and $\dot{\bm{y}}(\bm{x},t)=\partial\bm{y}(\bm{x},t)/\partial t$. Of course
$\dot{\bm{y}}(\bm{x},t)$ can be expressed in coordinates $\bm{y}$ via the
inverse map $\bm{x}(\bm{y},t)$. Moreover, the PDE domain $\Omega$ is mapped
into

$\Omega_{y}(t)=\left\\{\bm{z}\in\mathbb{R}^{d}:\bm{z}=\bm{y}(\bm{x},t),\quad\bm{x}\in\Omega\right\\}.$
(10)

by the coordinate flow $\bm{y}(\bm{x},t)$.

###  2.1 Dynamic tensor approximation

Denote by $\mathcal{M}_{\bm{r}}(t)\subseteq\mathcal{H}(\Omega_{y}(t))$ a
tensor manifold embedded in the Hilbert space

$\mathcal{H}(\Omega_{y}(t))=L^{2}(\Omega_{y}(t))$ (11)

of functions defined on the time-dependent domain
$\Omega_{y}(t)\subseteq\mathbb{R}^{d}$. Such manifolds include hierarchical
tensor formats such as the hierarchical Tucker and the tensor train formats
[1, 32]. We approximate the solution $v(\bm{y},t)$ to the PDE (8) with an
element $v_{\bm{r}}(\bm{y},t)$ belonging to $\mathcal{M}_{\bm{r}}(t)$. We
allow the tensor manifold to depend on $t$ so that the solution rank can be
chosen adaptively during time integration to ensure an accurate approximation.
At each point $v_{\bm{r}}\in\mathcal{M}_{\bm{r}}$ the function space
$\mathcal{H}(\Omega_{y})$ can be partitioned as a direct sum of two vector
subspaces [12]

$\mathcal{H}(\Omega_{y})=T_{v_{\bm{r}}}\mathcal{M}_{\bm{r}}\oplus
N_{v_{\bm{r}}}\mathcal{M}_{\bm{r}},$ (12)

where $T_{v_{\bm{r}}}\mathcal{M}_{\bm{r}}$ (resp.
$N_{v_{\bm{r}}}\mathcal{M}_{\bm{r}}$) denotes the tangent space (resp. normal
space) relative to the manifold $\mathcal{M}_{\bm{r}}$ at the point
$v_{\bm{r}}$. Assuming the solution to (8) at time $t$ can be represented by
an element of the tensor manifold $\mathcal{M}_{\bm{r}}(t)$, we utilize the
direct sum (12) to decompose the right hand side of the PDE (8) into a
tangential component and a normal component relative to
$\mathcal{M}_{\bm{r}}(t)$

$\displaystyle
Q_{\bm{y}}(v_{\bm{r}},\dot{\bm{y}})=\mathcal{P}_{T}({v}_{\bm{r}})Q_{\bm{y}}(v_{\bm{r}},\dot{\bm{y}})+\mathcal{P}_{N}({v}_{\bm{r}})Q_{\bm{y}}(v_{\bm{r}},\dot{\bm{y}}),\qquad
v_{\bm{r}}\in\mathcal{M}_{\bm{r}}(t),$ (13)

where $\mathcal{P}_{T}({v}_{\bm{r}})$ and
$\mathcal{P}_{N}({v}_{\bm{r}})=\mathcal{I}-\mathcal{P}_{T}({v}_{\bm{r}})$
denote, respectively, the orthogonal projections onto the tangent and normal
space of $\mathcal{M}_{\bm{r}}(t)$ at the point $v_{\bm{r}}$ (see Figure 1).
Projecting the initial condition in (8) onto $\mathcal{M}_{\bm{r}}(0)$ with a
truncation operator $\mathfrak{T}_{\bm{r}}(\cdot)$ and using only the
tangential component of the PDE dynamics for all $t\in[0,T]$ we obtain the
projected PDE

$\begin{cases}\displaystyle\frac{\partial v_{\bm{r}}(\bm{y},t)}{\partial
t}=\mathcal{P}_{T}(v_{\bm{r}})Q_{\bm{y}}(v_{\bm{r}}(\bm{y},t),\dot{\bm{y}}),\qquad
t\in[0,T],\\\
v_{\bm{r}}(\bm{y},0)=\mathfrak{T}_{\bm{r}}(u_{0}(\bm{y})),\end{cases}$ (14)

with solution $v_{\bm{r}}(\bm{y},t)$ which remains on the tensor manifold
$\mathcal{M}_{\bm{r}}(t)$ for all $t\in[0,T]$. We have previously shown [12]
that the approximation error

$\epsilon(t)=\left\|v(\bm{y},t)-v_{\bm{r}}(\bm{y},t)\right\|_{\mathcal{H}(\Omega_{y}(t))}$
(15)

can be controlled by selecting the rank $\bm{r}(t)$ at each $t$ so that the
norm of the normal component in (13) remains bounded by some user-defined
threshold $\varepsilon$, i.e.,

$\left\|\mathcal{P}_{N}(v_{\bm{r}})Q_{\bm{y}}(v_{\bm{r}}(\bm{y},t),\dot{\bm{y}})\right\|_{\mathcal{H}(\Omega_{y}(t))}<\varepsilon.$
(16)

Therefore the normal component determines the increase in solution rank
required to maintain an accurate approximation during time integration. On the
other hand, rank decrease during time integration can be interpreted as the
tangential component of the PDE operator guiding the solution to a region of
the manifold $\mathcal{M}_{\bm{r}}(t)$ with higher curvature. In fact, recall
that the smallest singular value of $v_{\bm{r}}$ is inversely proportional to
the curvature of the manifold $\mathcal{M}_{\bm{r}}(t)$ at the point
$v_{\bm{r}}$ [23].

###  2.2 Variational principle for rank-reducing coordinate flows

Figure 1: Sketch of the tensor manifold
$\mathcal{M}_{\bm{r}}(t)\subset\mathcal{H}(\Omega_{y}(t))$. Shown are the
normal and tangent components of $Q_{\bm{y}}(v_{\bm{r}},\dot{\bm{y}})$ (PDE
operator in equation (9)) at $v_{\bm{r}}\in\mathcal{M}_{\bm{r}}(t)$. The
manifold $\mathcal{M}_{\bm{r}}(t)$ has multilinear rank that depends on the
coordinate flow $\bm{y}(\bm{x},t)$ (see [14]). The variational principle (17)
minimizes the component of $Q_{\bm{y}}$ normal to the tensor manifold
$\mathcal{M}_{\bm{r}}(t)$ at each time, i.e., the component of $Q_{\bm{y}}$
that is responsible for the rank increase, relative to arbitrary variations of
the coordinate flow generator. This mitigates the tensor rank increase when
integrating (14) using the rank-adaptive tensor methods we recently developed
in [12, 31].

Leveraging the geometric structure of the tensor manifold
$\mathcal{M}_{\bm{r}}(t)$ described in Section 2.1, we aim at generating a
coordinate flow $\bm{y}(\bm{x},t)$ that minimizes the rank of the PDE solution
$v_{\bm{r}}(\bm{y},t)$ to (14) while maintaining an accurate approximation to
(7). Since the approximation error (15) is controlled by the normal component
of the PDE operator we select an infinitesimal generator for the coordinate
flow $\dot{\bm{y}}$ that minimizes such normal component, i.e., that minimizes
the left hand side of (16). This idea results in a convex optimization problem
over the Lie algebra $\mathfrak{g}$ of infinitesimal coordinate flow
generators $\dot{\bm{y}}$ at time $t$

$\displaystyle\min_{\dot{\bm{y}}\in\mathfrak{g}}\left\|\mathcal{P}_{N}(v_{\bm{r}})Q_{\bm{y}}(v_{\bm{r}},\dot{\bm{y}})\right\|_{\mathcal{H}(\Omega_{y}(t))}^{2}=\min_{\dot{\bm{y}}\in\mathfrak{g}}\left\|\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}})-\dot{\bm{y}}(\bm{x},t)\cdot\nabla_{\bm{y}}v_{\bm{r}}\right]\right\|_{\mathcal{H}(\Omega_{y}(t))}^{2}.$
(17)

Note that this functional is not optimizing the projection of
$\dot{\bm{y}}\cdot\nabla v_{\bm{r}}$ onto the tangent space of the tensor
manifold at $v_{\bm{r}}$. In principle, it is possible to control such
projection, and eventually have the tangential component of $Q_{\bm{y}}$
pointing towards a region of the tensor manifold with higher curvature. This
results in smaller singular values of the tensor decomposition of the solution
as time increases333Recall that the curvature of a tensor manifold is
inversely proportional to the inverse of the singular values [23]., which may
induce further tensor rank reduction via thresholded tensor truncation.

###### Proposition 2.1

Let $\bm{f}$ be a solution of the convex optimization problem (17). Then
$\bm{f}$ satisfies the linear system of equations

$\mathcal{P}_{N}(v_{\bm{r}})\left[\bm{f}\cdot\nabla v_{\bm{r}}\right]\nabla
v_{\bm{r}}=\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}})\right]\nabla
v_{\bm{r}}.$ (18)

Proof: Due to the convexity of (17) any critical point of the cost functional

$\displaystyle
C(\dot{\bm{y}})=\left\|\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}})-\dot{\bm{y}}(\bm{x},t)\cdot\nabla_{\bm{y}}v_{\bm{r}}\right]\right\|_{\mathcal{H}(\Omega_{y}(t))}^{2}$
(19)

in the Lie algebra $\mathfrak{g}$ is necessarily a global solution to the
optimization problem (17). The first variation of $C$ with respect to
$\dot{y}_{i}$ is easily obtained as

$\displaystyle\delta_{\dot{y}_{i}}C\eta_{i}$
$\displaystyle=2\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}})-\dot{\bm{y}}\cdot\nabla
v_{\bm{r}}\right]\mathcal{P}_{N}(v_{\bm{r}})\left[-\frac{\partial
v_{\bm{r}}}{\partial y_{i}}\eta_{i}\right]d\bm{y}$ (20)
$\displaystyle=2\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}})-\dot{\bm{y}}\cdot\nabla
v_{\bm{r}}\right]\left(-\frac{\partial v_{\bm{r}}}{\partial
y_{i}}\eta_{i}\right)d\bm{y},$

where $\eta_{i}\in\mathfrak{g}$ is an arbitrary perturbation. To obtain the
second equality in (20) we used the fact that the orthogonal projection
$\mathcal{P}_{N}(v_{\bm{r}})$ is symmetric with respect to the inner product
in $L^{2}(\Omega_{y}(t))$ and idempotent. Setting (20) equal to zero for
arbitrary $\eta_{i}$ yields the Euler-Lagrange equations

$\mathcal{P}_{N}(v_{\bm{r}})\left[\dot{\bm{y}}\cdot\nabla
v_{\bm{r}}\right]\frac{\partial v_{\bm{r}}}{\partial
y_{i}}=\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}})\right]\frac{\partial
v_{\bm{r}}}{\partial y_{i}},\qquad i=1,2,\ldots,d.$ (21)

Writing the preceding system of equations in vector notation proves the
Proposition.

$\square$

Note that if $\nabla v_{\bm{r}}$ is non-zero then from the optimality
conditions we obtain $\mathcal{P}_{N}(v_{\bm{r}})\left[\bm{f}\cdot\nabla
v_{\bm{r}}\right]=\mathcal{P}_{N}(v_{\bm{r}})G_{\bm{y}}(v_{\bm{r}})$. In this
case the normal component of $v_{\bm{r}}(\bm{y},t)$, i.e., the right hand side
of (16), is zero along the optimal coordinate flow and therefore the rank of
$v_{\bm{r}}$ never increases during time integration.

With the infinitesimal generator $\bm{f}(\bm{y},t)$ solving the convex
optimization problem (17) available (i.e., the generator satisfying the linear
system (18)), we can now compute a coordinate transformation
$\bm{y}(\bm{x},t)$ via the dynamical system

$\begin{cases}\displaystyle\frac{\partial\bm{y}}{\partial
t}=\bm{f}(\bm{y},t),\\\ \bm{y}(\bm{x},0)=\bm{x}.\end{cases}$ (22)

Such coordinate transformation allows us to solve the PDE (14) along a
coordinate flow that minimizes the projection of $Q_{\bm{y}}$ onto the normal
component of the tensor manifold in which the solution lives, henceforth
minimizing the rank increase of $v_{\bm{r}}$ during time integration. The
corresponding evolution equation along the rank reducing coordinate flow is
obtained by coupling the ODE (22) with (14) as

$\begin{cases}\displaystyle\frac{\partial v_{\bm{r}}(\bm{y},t)}{\partial
t}=\mathcal{P}_{T}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}}(\bm{y},t))-\bm{f}(\bm{y},t)\cdot\nabla_{\bm{y}}v_{\bm{r}}(\bm{y},t)\right],\\\
v_{\bm{r}}(\bm{y}(\bm{x},0),0)=\mathfrak{T}_{\bm{r}}(u_{0}(\bm{y})).\end{cases}$
(23)

The coupled system of equations (17), (22), and (23) allows us to devise a
time integration scheme for the PDE (1) that leverages diffeomorphic
coordinate flows to control the component of $Q_{\bm{y}}$ that is normal to
tensor manifold $\mathcal{M}_{\bm{r}}$, and therefore control the rank
increase of the PDE solution $v_{\bm{r}}$ in time.

To describe the scheme in more detail let us discretize the temporal domain
$[0,T]$ into $N+1$ evenly spaced time instants

$t_{k}=k\Delta t,\qquad\Delta t=\frac{T}{N},\qquad k=0,1,\ldots,N.$ (24)

To integrate the solution and the coordinate flow from time $t_{k}$ to time
$t_{k+1}$ we first solve the linear system of equations (18). This gives us an
optimal infinitesimal coordinate flow generator $\bm{f}(\bm{y},t_{k})$ at time
$t_{k}$. Using the evolution equations (22) and (23) we then integrate both
the coordinate flow $\bm{y}(\bm{x},t)$ and the PDE solution along the
coordinate flow $v_{\bm{r}}(\bm{y}(\bm{x},t),t)$ from time $t_{k}$ to time
$t_{k+1}$. Finally we update the PDE operator $G_{\bm{y}}$ to operate in the
new coordinate system.

The proposed coordinate flow time marching scheme comes with two computational
challenges. The first challenge is computing the optimal infinitesimal
generator $\bm{f}(\bm{y},t_{k})$, i.e., solving the linear system (18) at each
time step. If we represent $\bm{f}(\bm{y},t_{k})$ in the span of a finite-
dimensional basis then (18) yields a linear system for the coefficients of the
expansion. The size of such linear system depends on the dimension of the
chosen space. For instance, in the case of linear coordinate flows we obtain a
$d^{2}\times d^{2}$ system. The second challenge is representing the PDE
operator in coordinates $\bm{y}(\bm{x},t)$, i.e., constructing $G_{\bm{y}}$.
For a general nonlinear coordinate transformation $\bm{y}(\bm{x},t)$ the
operator $G_{\bm{y}}$ includes the metric tensor of the transformation which
can significantly complicate the form of $G_{\bm{y}}$ (see, e.g., [2, 25]).
Hereafter we narrow our attention to linear coordinate flows and derive the
corresponding time integration scheme which does not suffer from the
aforementioned computational challenges.

####  2.2.1 Linear coordinate flows

Any coordinate flow $\bm{y}(\bm{x},t)$ that is linear in $\bm{x}$ at each time
$t$ can be represented by a time-dependent invertible $d\times d$ matrix
$\bm{\Gamma}(t)$

$\bm{y}(\bm{x},t)=\bm{\Gamma}(t)\bm{x}.$ (25)

Here we consider arbitrary matrices $\bm{\Gamma}(t)$ belonging to the
collection of all444It is also possible to generate linear coordinate flows
using subgroups of $\operatorname*{\text{GL}}_{d}(\mathbb{R})$ such as the
subgroup of matrices with determinant equal to 1 (i.e. volume preserving
linear maps), or the group of rigid rotations in $d$ dimensions. $d\times d$
invertible matrices with real entries, denoted by
$\operatorname*{\text{GL}}_{d}(\mathbb{R})$. The corresponding collection of
infinitesimal coordinate flow generators is the vector space of all $d\times
d$ matrices with real entries. In this setting the time-dependent tensor
approximation of $u(\bm{x},t)$ in coordinates $\bm{y}$ takes the form

$v_{\bm{r}}(\bm{\Gamma}(t)\bm{x},t)\approx u(\bm{x},t).$ (26)

This is known as a (time-dependent) generalized tensor ridge function [14,
28]. The time-dependent optimization problem (17) that generates the optimal
linear coordinate flow can be now written as

$\min_{\dot{\bm{\Gamma}}\in M_{d\times
d}(\mathbb{R})}C\left(\dot{\bm{\Gamma}}\bm{y}\right),$ (27)

where $C(\cdot)$ is the cost function defined in (19).

###### Proposition 2.2

Let $\bm{\Sigma}(t)$ be a solution of the optimization problem (27) at time
$t$. Then $\bm{\Sigma}(t)$ satisfies the symmetric $d^{2}\times d^{2}$ linear
system

$\bm{A}{\rm vec}(\bm{\Sigma}(t))=\bm{b},$ (28)

where ${\rm vec}(\bm{\Sigma}(t))$ is a vectorization of the $d\times d$ matrix
$\bm{\Sigma}(t)$,

$\bm{A}=\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[{\rm
vec}(\bm{c}(\bm{y}))\right]{\rm
vec}(\bm{c}(\bm{y}))^{\text{T}}d\bm{y},\qquad\bm{b}=\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}}(\bm{y},t))\right]{\rm
vec}(\bm{c}(\bm{y}))d\bm{y},$ (29)

and $\bm{c}(\bm{y})$ is the $d\times d$ matrix with entries

$c_{ij}(\bm{y})=\frac{\partial v_{\bm{r}}(\bm{y},t)}{\partial
y_{i}}y_{j},\qquad i,j=1,2,\ldots,d.$ (30)

Proof: First note that the cost function
$C\left(\dot{\bm{\Gamma}}\bm{y}\right)$ is convex in $\dot{\bm{\Gamma}}$ and
the search space $M_{d\times d}(\mathbb{R})$ is linear. Hence, any critical
point of the cost function is necessarily a global minimum. To find such a
critical point we calculate the derivative of $C(\dot{\bm{\Gamma}}\bm{y})$
with respect to each entry of the matrix $\dot{\bm{\Gamma}}$ directly

$\displaystyle\frac{\partial
C(\dot{\bm{\Gamma}}\bm{y})}{\partial\dot{\Gamma}_{ij}}$
$\displaystyle=2\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{\Gamma}}(v_{\bm{r}})-(\dot{\bm{\Gamma}}\bm{y})\cdot\nabla
v_{\bm{r}}\right]\mathcal{P}_{N}(v_{\bm{r}})\left[-y_{j}\frac{\partial
v_{\bm{r}}}{\partial y_{i}}\right]d\bm{y}$ (31)
$\displaystyle=2\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{\Gamma}}(v_{\bm{r}})-(\dot{\bm{\Gamma}}\bm{y})\cdot\nabla
v_{\bm{r}}\right]\left(-y_{j}\frac{\partial v_{\bm{r}}}{\partial
y_{i}}\right)d\bm{y},$

for all $i,j=1,2,\ldots,d$, where in the second line we used the fact that the
orthogonal projection $\mathcal{P}_{N}(v_{\bm{r}})$ is symmetric with respect
to the inner product in $L^{2}(\Omega_{y}(t))$ and idempotent. The critical
points are then obtained by setting $\partial
C(\dot{\bm{\Gamma}})/\partial\dot{\Gamma}_{ij}=0$ for all $i,j=1,2,\ldots,d$,
resulting in a linear system of equations for the entries of the infinitesimal
linear coordinate flow generator $\dot{\bm{\Gamma}}$

$\sum_{k,p=1}^{d}\dot{\Gamma}_{kp}\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[y_{p}\frac{\partial
v_{\bm{r}}}{\partial y_{k}}\right]y_{j}\frac{\partial v_{\bm{r}}}{\partial
y_{i}}d\bm{y}=\int_{\Omega_{y}(t)}\mathcal{P}_{N}(v_{\bm{r}})\left[G_{\bm{y}}(v_{\bm{r}})\right]y_{j}\frac{\partial
v_{\bm{r}}}{\partial y_{i}}d\bm{y},$ (32)

for $i=1,2,\ldots,d$. We obtain (28) by writing the linear system (32) in
matrix notation.

$\square$

The linear system of equations (32) can also be obtained directly from the
Euler Lagrange equations (21) by substituting
$\dot{\bm{y}}=\dot{\bm{\Gamma}}\bm{y}$ and projecting onto $y_{j}$ for
$j=1,2,\ldots,d$.

A few remarks are in order regarding the linear system (28). The matrix
$\bm{A}$ is a $d^{2}\times d^{2}$ symmetric matrix which is determined by
$(d^{4}+d^{2})/2$ entries and the vector $\bm{b}$ has length $d^{2}$.
Computing each entry of the matrix $\bm{A}$ and the vector $\bm{b}$ requires
evaluating a $d$-dimensional integral. Therefore computing $\bm{A}$ and
$\bm{b}$ to set up the linear system at each time $t$ requires evaluating a
total of $(d^{4}+3d^{2})/2$ $d$-dimensional integrals. Since $v_{\bm{r}}$ is a
low-rank tensor these high-dimensional integrals can computed by applying one-
dimensional quadrature rules to the tensor modes. By orthogonalizing the low-
rank tensor expansion (e.g., tensor train) it is possible to reduce the number
one-dimensional integrals needed to compute $\bm{A}$. If the matrix $\bm{A}$
is singular then the solution to the linear system of equations is not unique.
In this case any solution will suffice for generating a coordinate flow.

For linear coordinate transformations it may be advantageous to consider the
related functional

$\displaystyle C_{2}(\dot{\bm{\Gamma}}\bm{y})$
$\displaystyle=\left\|\frac{\partial v_{\bm{r}}}{\partial
t}\right\|_{\mathcal{H}(\Omega_{y}(t))}^{2}$ (33)
$\displaystyle=\left\|Q_{\bm{y}}(v_{\bm{r}})\right\|_{\mathcal{H}(\Omega_{y}(t))}^{2}$
$\displaystyle=\left\|\mathcal{P}_{T}(v_{\bm{r}})Q_{\bm{y}}(v_{\bm{r}},\dot{\bm{y}})\right\|_{\mathcal{H}(\Omega_{y}(t))}^{2}+\left\|\mathcal{P}_{N}(v_{\bm{r}})Q_{\bm{y}}(v_{\bm{r}},\dot{\bm{y}})\right\|_{\mathcal{H}(\Omega_{y}(t))}^{2}.$

The linear coordinate flow minimizing (33) aims at “undoing” the action of the
operator $G_{\bm{y}}$ in coordinates $\bm{y}$. More precisely the PDE solution
$u(\bm{x},t)$ as seen in coordinate $\bm{y}(\bm{x},t)=\bm{\Gamma}(t)\bm{x}$,
i.e., $v_{r}(\bm{\Gamma}(t)\bm{x},t)\approx u(\bm{x},t)$, varies in time as
little as possible. In the third equality of (33) we used the fact that
$\mathcal{P}_{T}(v_{\bm{r}})$ and $\mathcal{P}_{N}(v_{\bm{r}})$ are orthogonal
projections. Following similar steps used to prove Proposition 2.2 it is
straightforward to show that the critical points $\bm{\Sigma}(t)$ of (33)
satisfy the linear system (28) with

$\bm{A}=\int_{\Omega_{y}(t)}{\rm vec}(\bm{c}(\bm{y})){\rm
vec}(\bm{c}(\bm{y}))^{\text{T}}d\bm{y},\qquad\bm{b}=\int_{\Omega_{y}(t)}G_{\bm{y}}(v_{\bm{r}}(\bm{y},t)){\rm
vec}(\bm{c}(\bm{y}))d\bm{y}.$ (34)

With the optimal coordinate flow generator $\bm{\Sigma}(t)$ available we can
compute the matrix $\bm{\Gamma}(t)$ appearing in (25) by solving the matrix
differential equation

$\begin{cases}\displaystyle\frac{d\bm{\Gamma}(t)}{dt}=\bm{\Sigma}(t)\bm{\Gamma}(t),\\\
\bm{\Gamma}(0)=\bm{I}_{d\times d}.\end{cases}$ (35)

Correspondingly, the PDE (23) can be written along the optimal rank-reducing
linear coordinate flow (25) as

$\begin{cases}\displaystyle\frac{\partial
v_{\bm{r}}(\bm{\Gamma}(t)\bm{x},t)}{\partial
t}=G_{\bm{y}}(v_{\bm{r}}(\bm{y},t),\bm{y})-\left(\bm{\Sigma}(t)\bm{y}(\bm{x},t)\right)\cdot\nabla_{\bm{y}}v_{\bm{r}}(\bm{y}(\bm{x},t),t),\\\
v_{\bm{r}}(\bm{x},0)=\mathfrak{T}_{\bm{r}}(u_{0}(\bm{x})).\end{cases}$ (36)

Using the linear system (28), the matrix differential equation (35), and the
PDE (36) we can specialize the time integration scheme described in Section
2.2 for nonlinear coordinate flows to linear coordinate flows. In Algorithm 1
we provide pseudo-code for PDE time stepping with linear coordinate flows. In
the Algorithm ${\rm TDVP}$ denotes a subroutine that solves the optimization
problem (27) (e.g., by solving the symmetric $d^{2}\times d^{2}$ linear system
of equations (28)) and ${\rm LMM}$ denotes any time stepping scheme.

Input:

$v_{0}$ $\rightarrow$ initial condition

$\Delta t$ $\rightarrow$ temporal step size

$N_{t}$ $\rightarrow$ total number of time steps

Output:

$\bm{\Gamma}$ $\rightarrow$ linear coordinate transformation at final time
$t_{f}$

$v(\bm{x},t_{f})=u(\bm{\Gamma}\bm{x},t_{f})$ $\rightarrow$ Solution at time
$t_{f}$ on transformed coordinate system

Runtime:

  * $\bm{\Gamma}=\bm{I}_{d\times d}$

  * for $k=0$ to $N_{t}$

    * 1

$\dot{\bm{\Gamma}}_{k}={\rm TDVP}(v_{k},G_{\bm{y}})$

* 2

$v_{k+1}={\rm LMM}(v_{k},G_{\bm{y}},\Delta t)$

* 3

$\bm{\Gamma}_{k+1}={\rm LMM}(\bm{\Gamma}_{k},\dot{\bm{\Gamma}}_{k},\Delta t)$

end

Algorithm 1 PDE integrator with rank-reducing linear coordinate flow.

##  3 Numerical examples

We demonstrate the proposed coordinate-adaptive time integration scheme based
on linear coordinate flows (Algorithm 1) on several prototype PDEs projected
on functional tensor train (FTT) manifolds $\mathcal{M}_{\bm{r}}$ [13].

###  3.1 Liouville equation

Consider the $d$-dimensional Liouville equation

$\begin{cases}\displaystyle\frac{\partial u(\bm{x},t)}{\partial
t}=\bm{f}(\bm{x})\cdot\nabla u(\bm{x},t),\vspace{0.2cm}\\\
u(\bm{x},0)=u_{0}(\bm{x}),\end{cases}$ (37)

where $\bm{f}(\bm{x})$ is a divergence-free vector field. We set the initial
condition to be a product of independent Gaussians

$u_{0}(\bm{x})=\frac{1}{m}\exp\left(-\sum_{j=1}^{d}\frac{1}{\beta_{j}}\left(x_{j}+t_{j}\right)^{2}\right),$
(38)

where $\beta_{j}>0$, $t_{j}\in\mathbb{R}$ and $m$ is the normalization factor

$m=\left\|\exp\left(-\sum_{j=1}^{d}\frac{1}{\beta_{j}}\left(x_{j}+t_{j}\right)^{2}\right)\right\|_{L^{1}\left(\mathbb{R}^{d}\right)}.$

####  3.1.1 Two-dimensional simulation results

We first consider the Liouville on a two-dimensional flat torus
$\Omega\subset\mathbb{R}^{2}$ with velocity vector
$\bm{f}(\bm{x})=(f_{1}(\bm{x}),f_{2}(\bm{x}))$ consisting of

$f_{1}(\bm{x})=\frac{\partial\psi}{\partial x_{2}},\qquad
f_{2}(\bm{x})=-\frac{\partial\psi}{\partial x_{1}}$ (39)

generated via the two-dimensional stream function [38]

$\psi(x_{1},x_{2})=\Theta(x_{1})\Theta(x_{2})$ (40)

with

$\Theta(x)=\frac{\cos(\alpha x/L)}{\cos(\alpha/2)}-\frac{\cosh(\alpha
x/L)}{\cosh(\alpha/2)},$ (41)

$L=30$ and $\alpha=4.73$. This yields a two-dimensional measure-preserving
(divergence-free) velocity vector $\bm{f}(\bm{x})$. In the initial condition
(38) we set $\bm{\beta}=\begin{bmatrix}1/4&2\end{bmatrix}$ and
$\bm{t}=\begin{bmatrix}3&3\end{bmatrix}$ resulting in a rank-$1$ initial
condition.

(a) (b) (c)

Figure 2: Two-dimensional advection PDE. (a) Initial condition. (b) Solution
in Cartesian coordinates at time $t=35$. (c) Solution in time-dependent
adaptive coordinates at time $t=35$.

Note that this is the same problem we recently studied in [14] using
coordinate-adaptive tensor integration based on non-convex relaxations of the
rank minimization problem. We discretize the PDE with the Fourier pseudo-
spectral method [19] using $200$ points in each spatial variable $x_{i}$ and
integrate the PDE from $t=0$ to $t=35$ using the two-steps Adams-Bashforth
scheme (AB2) with time step size $\Delta t=10^{-3}$. To demonstrate rank
reduction with coordinate flows, we run three distinct simulations. In the
first simulation we integrate the PDE on a full two-dimensional tensor product
grid in Cartesian coordinates for all time resulting in our benchmark
solution. In the second simulation we construct a low-rank representation of
the solution in fixed Cartesian coordinates with the rank-adaptive step-
truncation time integration scheme [31] built upon the AB2 integrator using
relative tensor truncation tolerance $\delta=10^{-8}$. The third simulation
demonstrates the low-rank tensor format in adaptive coordinates generated by
Algorithm 1 (using cost function 33) with a step-truncation time integration
scheme built upon the AB2 integrator using relative tensor truncation
tolerance $\delta=10^{-8}$.

In Figure 2 we plot the low-rank solution at time $t=0$ (a), $t=35$ computed
in fixed Cartesian coordinates (b) and $t=35$ computed in adaptive coordinates
(c). We notice that the time-dependent coordinate system in the adaptive
simulation captures the affine components (i.e., translation, rotation, and
stretching) of the dynamics generated by the PDE operator. In Figure 3 (d) we
plot the rank versus time of the solution computed in fixed Cartesian
coordinates and the solution computed in adaptive time-dependent coordinates.
We also plot the norm of the normal components, i.e., the left hand side of
(16), versus time in Cartesian coordinates and in adaptive coordinates. Since
the linear coordinate flow is minimizing the normal component of the PDE
operator we observe that the red dashed line increases at a slower rate than
the blue dashed line. Once the norm of the normal component reaches a
threshold depending on the relative tensor truncation tolerance $\delta$ set
on the solution and time step size $\Delta t$, the rank-adaptive criterion
discussed [12] is triggered which results in an increase of the solution rank.
At the subsequent time step the normal component is reduced due to the
addition of a tensor mode.

(a) (b)

Figure 3: Two-dimensional advection PDE. (a) Rank versus time and (b) norm of
the normal component of PDE dynamics versus time for the low-rank solution in
Cartesian coordinates (blue solid line) and in adaptive coordiantes (red
dashed line).

At time $t=35$ we check the $L^{\infty}$ error of the low-rank tensor solution
computed in Cartesian coordinates and the low-rank solution computed using
adaptive coordinates (for the adaptive solution we map it back to Cartesian
using a two-dimensional interpolant) and find that both errors are bounded by
$1.5\times 10^{-5}$.

####  3.1.2 Three-dimensional simulation results

Next we consider the Liouville equation (37) on a three-dimensional flat torus
$\Omega\subset\mathbb{R}^{3}$ with velocity vector

$\bm{f}(\bm{x})=\begin{bmatrix}x_{2}\\\ -x_{1}\\\ x_{2}\\\ \end{bmatrix}.$
(42)

This allows us to test our coordinate-adaptive algorithm for a problem in
which we know the optimal ridge matrix. We set the parameters of the initial
condition (38) as $\bm{\beta}=\begin{bmatrix}4&1/4&4\end{bmatrix}$ and
$\bm{t}=\begin{bmatrix}1&-1&1\end{bmatrix}$ resulting in the FTT rank
$\bm{r}(0)=\begin{bmatrix}1&1&1&1\end{bmatrix}$. Due to the choice of
coefficients (42), the analytical solution to the PDE (37) can be written as a
ridge function in terms of the PDE initial condition

$u(\bm{x},t)=u_{0}\left(e^{t\bm{B}}\bm{x}\right),$ (43)

where

$\bm{B}=\begin{bmatrix}0&1&0\\\ -1&0&0\\\ 0&1&0\end{bmatrix}.$ (44)

Thus (43) is a tensor ridge solution to the 3D Lioville equation (37) with the
same rank as the initial condition, i.e., in this case there exists a tensor
ridge solution at each time with FTT rank equal to
$\begin{bmatrix}1&1&1&1\end{bmatrix}$.

To demonstrate the performance of the rank-reducing coordinate flow algorithm,
we ran two low-rank simulations of the PDE (37) up to $t=1$ and compared the
results with the analytic solution (43). The first simulation is computed
using the rank-adaptive step-truncation methods discussed in [12, 31] with
relative truncation accuracy $\delta=10^{-5}$, Cartesian coordinates, and time
step size $\Delta t=10^{-3}$. The second simulation demonstrates the
performance of the rank-reducing coordinate flow (Algorithm 1) and is computed
with the same step-truncation scheme as the simulation in Cartesian
coordinates. For this example the coordinate flow is constructed by minimizing
the functional (33).

In Figure 4(a) we plot the $1$-norm of the FTT solution rank versus time. We
observe that the FTT solution rank in Cartesian coordinates grows rapidly
while the FTT-ridge solution generated by the proposed Algorithm 1 recovers
the optimal time-dependent low-rank tensor ridge solution, which remains
constant rank for all time at $\left\|\bm{r}(t)\right\|_{1}=4$. We also map
the FTT ridge solution back to Cartesian coordinates using a three-dimensional
trigonometric interpolant every $100$ time steps and compare both low-rank FTT
solutions to the benchmark solution. In Figure 4(b) we plot the $L^{\infty}$
error versus time. The solution computed in adaptive coordinates does not
change in time and therefore is about one order of magnitude more accurate
than the solution computed in Cartesian coordinates. In this case the proposed
algorithm outperforms our previous algorithm proposed in [14] which does not
recover the optimal linear coordinate transformation due to the non-convexity
of the cost function used to generate the time-dependent coordinate system.

(a) (b)

Figure 4: Three-dimensional advection equation (37). (a) Rank versus time and
(b) $L^{\infty}$ error of the low-rank solutions relative to the benchmark
solution.

For this example the PDE operator has separation rank (CP-rank) $3$ in
Cartesian coordinates and CP-rank $9$ in a general linearly transformed
coordinate system.

###  3.2 Fokker-Planck equation

Finally we consider the Fokker-Planck equation

$\begin{cases}\displaystyle\frac{\partial u(\bm{x},t)}{\partial
t}=-\nabla\cdot\left(\bm{f}(\bm{x})u(\bm{x},t)\right)+\sum_{i,j=1}^{d}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}\left(D_{ij}(\bm{x})u(\bm{x},t)\right),\vspace{0.2cm}\\\
u(\bm{x},0)=u_{0}(\bm{x}),\end{cases}$ (45)

in three spatial dimensions ($d=3$) on a flat torus
$\Omega\subseteq\mathbb{R}^{3}$. We set the drift vector $\bm{f}$ as in (42),
and diffusion tensor

$\bm{D}=\sigma\begin{bmatrix}\exp\left(-x_{2}^{2}\right)&0&0\\\
0&\exp\left(-x_{3}^{2}\right)&0\\\
0&0&\exp\left(-x_{1}^{2}\right)\end{bmatrix}$ (46)

with $\sigma=1/4$. We use the same initial condition as in Section 3.1.2.

Once again, to demonstrate the performance of the rank-reducing coordinate
flow algorithm, we ran three simulations of the PDE (45) up to $t=1$. The
first simulation is computed with a three-dimensional Fourier pseudo-spectral
method on a tensor product grid with 200 points per dimension and AB2 time
stepping with $\Delta t=10^{-4}$. This yields an accurate solution benchmark.
The second simulation is computed using the rank-adaptive step-truncation
methods discussed in [12, 31] using relative truncation accuracy
$\delta=10^{-5}$, fixed Cartesian coordinates, and time step size $\Delta
t=10^{-3}$. The third simulation demonstrates the performance of the rank-
reducing coordinate flow (Algorithm 1) and is computed with the same step-
truncation scheme as the simulation in Cartesian coordinates. For this example
the coordinate flow is constructed by minimizing the functional (33).

$t=0.1$ $t=0.5$ $t=1.0$

Cartesian coordinates

adaptive coordinates

Figure 5: Three-dimensional Fokker-Planck equation (45). Time snapshots of the
solution at $t=0.1,0.5,1.0$ in Cartesian coordinates (top row) and adaptive
coordinates (bottom row). The adaptive coordinate system generated by
Algorithm 1 effectively captures the affine effects (i.e., rotation,
stretching, translation) of the transformation generated by the Fokker-Planck
equation (45).

In Figure 5 we plot three time snapshots of the low-rank solutions to (45) in
fixed Cartesian coordinates (top row) and rank-reducing adaptive coordinates
(bottom row). We observe that the rank-reducing adaptive coordinate
transformation effectively captures the affine effects of the transformation
generated by the PDE (45) even in the presence of diffusion. In Figure 6(a) we
plot the $1$-norm of the FTT solution rank versus time. We observe that the
FTT solution rank in Cartesian coordinates grows rapidly while the FTT-ridge
solution generated by Algorithm 1 grows significantly slower. We also map the
FTT ridge solution back to Cartesian coordinates using a three-dimensional
trigonometric interpolant every $100$ time steps and compare both low-rank FTT
solutions to the benchmark solution. In Figure 6(b) we plot the $L^{\infty}$
error of Cartesian and coordinate-adaptive low-rank solutions relative to the
benchmark solution versus time. We observe that the process of solving the PDE
in transformed coordinates and mapping the transformed solution back to
Cartesian coordinates incurs negligible error.

(a) (b)

Figure 6: Three-dimensional Fokker-Planck equation (45). (a) Rank versus time
and (b) $L^{\infty}$ error of the low-rank solutions relative to the benchmark
solution.

##  4 Conclusions

We proposed a new tensor integration method for time-dependent PDEs that
controls the rank of the PDE solution in time by using diffeomorphic
coordinate transformations. Such coordinate transformations are generated by
minimizing the normal component of the PDE operator (written in intrinsic
coordinates) relative to the tensor manifold that approximates the PDE
solution. This minimization principle is defined by a convex functional which
can be computed efficiently and has optimality guarantees. The proposed method
significantly improves upon and may be used in conjunction with the
coordinate-adaptive algorithms we recently proposed in [14], which are based
on non-convex relaxations of the rank minimization problem and Riemannian
optimization. We demonstrated the proposed coordinate-adaptive time-
integration algorithm for linear coordinate transformations on prototype
Liouville and Fokker-Planck equations in two and three dimensions. Our
numerical results clearly demonstrate that linear coordinate flows can capture
the affine component (i.e., rotation, translation, and stretching) of the
transformation generated by PDE operator very effectively. In general, one
cannot expect linear coordinate flow (or even nonlinear coordinate flows) to
fully control the rank of the solution generated by an arbitrary nonlinear
PDE. Yet, the proposed method allows us to solve certain classes of PDEs at a
computational cost that is significantly lower than standard temporal
integration on tensor manifolds in Cartesian coordinates. In general overall
computational cost of the proposed method is not only determined by the tensor
rank of the solution but also the rank of the PDE operator $G_{\bm{y}}$, as
discussed in [14]. Further research is warranted to determine if the
efficiency of the proposed coordinate-adaptive methodology can be improved by
including conditions on the tangent projection of $\dot{\bm{y}}\cdot\nabla
v_{\bm{r}}$ in (17), or by simultaneously controlling the PDE solution and
operator rank during time integration.

## Acknowledgements

This research was supported by the U.S. Air Force Office of Scientific
Research grant FA9550-20-1-0174, and by the U.S. Army Research Office grant
W911NF1810309.

## References

  * [1] A. Nouy A. Falco, W. Hackbusch.  Geometric structures in tensor representations (final release).  arXiv:1505.03027v2, pages 1–50, 2015.
  * [2] R. Aris.  Vectors, tensors and the basic equations of fluid mechanics.  Dover publications, 1989.
  * [3] S.x D. Baalrud and J. Daligault.  Mean force kinetic theory: A convergent kinetic theory for weakly and strongly coupled plasmas.  Phys. Plasmas, 26:082106, 2019.
  * [4] G. Beylkin and M. J. Mohlenkamp.  Numerical operator calculus in higher dimensions.  PNAS, 99(16):10246–10251, 2002.
  * [5] D. Bigoni, A. P. Engsig-Karup, and Y. M. Marzouk.  Spectral tensor-train decomposition.  SIAM J. Sci. Comput., 38(4):A2405–A2439, 2016.
  * [6] A. M. P. Boelens, D. Venturi, and D. M. Tartakovsky.  Parallel tensor methods for high-dimensional linear PDEs.  J. Comput. Phys., 375:519–539, 2018.
  * [7] A. M. P. Boelens, D. Venturi, and D. M. Tartakovsky.  Tensor methods for the Boltzmann-BGK equation.  J. Comput. Phys., 421:109744, 2020.
  * [8] C. Brennan and D. Venturi.  Data-driven closures for stochastic dynamical systems.  J. Comp. Phys., 372:281–298, 2018.
  * [9] C. Cercignani.  The Boltzmann equation and its applications.  Springer, 1988.
  * [10] H. Cho, D. Venturi, and G. E. Karniadakis.  Numerical methods for high-dimensional probability density function equation.  J. Comput. Phys., 315:817–837, 2016.
  * [11] H. Cho, D. Venturi, and G. E. Karniadakis.  Numerical methods for high-dimensional kinetic equations.  In S. Jin and L. Pareschi, editors, Uncertainty quantification for kinetic and hyperbolic equations, pages 93–125. Springer, 2017.
  * [12] A. Dektor, A. Rodgers, and D. Venturi.  Rank-adaptive tensor methods for high-dimensional nonlinear PDEs.  J. Sci. Comput., 88(36):1–27, 2021.
  * [13] A. Dektor and D. Venturi.  Dynamic tensor approximation of high-dimensional nonlinear PDEs.  J. Comput. Phys., 437:110295, 2021.
  * [14] A. Dektor and D. Venturi.  Tensor rank reduction via coordinate flows.  Journal of Computational Physics, 491:112378, 2023.
  * [15] A. Etter.  Parallel ALS algorithm for solving linear systems in the hierarchical Tucker representation.  SIAM J. Sci. Comput., 38(4):A2585–A2609, 2016.
  * [16] R. Friedric, A. Daitche, O. Kamps, J. Lülff, M. Voßkuhle, and M. Wilczek.  The Lundgren-Monin-Novikov hierarchy: kinetic equations for turbulence.  Comptes Rendus Physique, 13(9–10):929–953, 2012.
  * [17] L. Grasedyck and C. Löbbert.  Distributed hierarchical SVD in the hierarchical Tucker format.  Numer. Linear Algebra Appl., 25(6):e2174, 2018.
  * [18] M. Griebel and G. Li.  On the decay rate of the singular values of bivariate functions.  SIAM J. Numer. Anal., 56(2):974–993, 2019.
  * [19] J. S. Hesthaven, S. Gottlieb, and D. Gottlieb.  Spectral methods for time-dependent problems, volume 21 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, 2007.
  * [20] I. Hosokawa.  Monin-Lundgren hierarchy versus the Hopf equation in the statistical theory of turbulence.  Phys. Rev. E, 73:067301(1–4), 2006.
  * [21] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3:422–440, 2021.
  * [22] 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, volume 48 of ESAIM Proc. Surveys, pages 1–28. EDP Sci., Les Ulis, 2015.
  * [23] E. Kieri, C. Lubich, and H. Walach.  Discretized dynamical low-rank approximation in the presence of small singular values.  SIAM J. Numer. Anal., 54(2):1020–1038, 2016.
  * [24] T. S. Lundgren.  Distribution functions in the statistical theory of turbulence.  Phys. Fluids, 10(5):969–975, 1967.
  * [25] H. Luo and T. R. Bewley.  On the contravariant form of the Navier-Stokes equations in time-dependent curvilinear coordinate systems.  J. Comput. Phys., 199:355–375, 2004.
  * [26] D. Montgomery.  A BBGKY framework for fluid turbulence.  Physics of Fluids, 19:802–810, 1976.
  * [27] I. V. Oseledets.  Tensor-train decomposition.  SIAM J. Sci. Comput., 33(5):2295––2317, 2011.
  * [28] A. Pinkus.  Ridge Functions.  Cambridge University Press, 2015.
  * [29] M. Raissi and G. E. Karniadakis.  Hidden physics models: Machine learning of nonlinear partial differential equations.  J. Comput. Phys., 357:125–141, 2018.
  * [30] 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–707, 2019.
  * [31] A. Rodgers, A. Dektor, and D. Venturi.  Adaptive integration of nonlinear evolution equations on tensor manifolds.  J. Sci. Comput., 92(39):1–31, 2022.
  * [32] R. Schneider and A. Uschmajew.  Approximation rates for the hierarchical tensor format in periodic Sobolev spaces.  J. Complexity, 30(2):56–71, 2014.
  * [33] J. Serrin.  The initial value problem for the Navier–Stokes equations.  Proc. Sympos., Madison, Wis., pages 69–98, 1963.
  * [34] A. Uschmajew and B. Vandereycken.  The geometry of algorithms using hierarchical tensors.  Linear Algebra Appl., 439(1):133–166, 2013.
  * [35] D. Venturi.  Convective derivatives and Reynolds transport in curvilinear time-dependent coordinate systems.  J. Phys. A: Math. Gen., 42:125203(1–16), 2009.
  * [36] D. Venturi.  Conjugate flow action functionals.  J. Math. Phys., 54:113502(1–19), 2013.
  * [37] D. Venturi.  The numerical approximation of nonlinear functionals and functional differential equations.  Physics Reports, 732:1–102, 2018.
  * [38] D. Venturi, M. Choi, and G.E. Karniadakis.  Supercritical quasi-conduction states in stochastic rayleigh–bénard convection.  International Journal of Heat and Mass Transfer, 55(13):3732–3743, 2012.