Coordinate-adaptive integration of PDEs on tensor manifolds

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 JCP (2023) Vol. 491, 112378, 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.


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 which governs the time evolution of a quantity of interest u(x, t) (high-dimensional field) over a compact domain Ω ⊆ R d (d ≫ 1) and has temporal dynamics generated by the nonlinear operator G x .The subscript "x" in G x indicates that the operator can explicitly depend on the variables x ∈ Ω.For instance, 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 physicsinformed 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 where u r is a given FTT tensor, and v s is a FTT tensor with smaller rank, i.e., |s| ≤ |r|.To compute the unknown matrix A (which enables tensor rank reduction), we developed a Riemannian gradient descent algorithm on some matrix manifold (e.g., A ∈ SL d (R)) that minimizes the non-convex cost functional 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 timedependent 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.

Coordinate-adaptive time integration of PDEs on tensor manifolds
We begin by representing the solution u(x, t) to the PDE (1) with another function defined on a time-dependent curvilinear coordinate system y(x, t) [2,33,35] with y(x, 0) = x.It is assumed that y(x, t) is a diffeomorphism1 which we will refer to as coordinate flow.To simplify notation we may not explicitly write the dependence of y on x (and vice-versa the dependence of x on y), or the dependence of x, y on t.However, it is always assumed that x and y are related via a (time-dependent) diffeomorphism.Next, we define two distinct time derivatives of v(y(x, t), t).The first represents a change in v at fixed location y, i.e., a derivative of v with respect to time with y constant, which we denote by ∂v(y(x, t), t)/∂t.The second represents a change in v along the coordinate flow y, i.e., a derivative of v with respect to time with x constant, which we denote by Dv(y(x, t), t)/Dt.
The derivative Dv(y(x, t), t)/Dt is known as material derivative [2,25], or convective derivative [35,36] Since u(x, t) is a solution to (1), v is related to u by (5) and y(x, 0) = x, it follows immediately that v satisfies the where the operator G y can be derived by writing G x in coordinates y using standard tools of differential geometry [2,35,36].Combining ( 6) and ( 7) we obtain2 where and ẏ(x, t) = ∂y(x, t)/∂t.Of course ẏ(x, t) can be expressed in coordinates y via the inverse map x(y, t).Moreover, the PDE domain Ω is mapped into by the coordinate flow y(x, t).

Dynamic tensor approximation
Denote by M r (t) ⊆ H(Ω y (t)) a tensor manifold embedded in the Hilbert space of functions defined on the time-dependent domain Ω y (t) ⊆ 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(y, t) to the PDE (8) with an element v r (y, t) belonging to M 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 r ∈ M r the function space H(Ω y ) can be partitioned as a direct sum of two vector subspaces [12] H where T vr M r (resp.N vr M r ) denotes the tangent space (resp.normal space) relative to the manifold M r at the point v r .Assuming the solution to (8) at time t can be represented by an element of the tensor manifold M 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 M r (t) where P T (v r ) and P N (v r ) = I −P T (v r ) denote, respectively, the orthogonal projections onto the tangent and normal space of M r (t) at the point v r (see Figure 1).Projecting the initial condition in (8) onto M r (0) with a truncation operator T r (•) and using only the tangential component of the PDE dynamics for all t ∈ [0, T ] we obtain the projected Figure 1: Sketch of the tensor manifold Mr(t) ⊂ H(Ωy(t)).Shown are the normal and tangent components of Qy(vr, ẏ) (PDE operator in equation ( 9)) at vr ∈ Mr(t).The manifold Mr(t) has multilinear rank that depends on the coordinate flow y(x, t) (see [14]).The variational principle ( 17) minimizes the component of Qy normal to the tensor manifold Mr(t) at each time, i.e., the component of Qy 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].
with solution v r (y, t) which remains on the tensor manifold M r (t) for all t ∈ [0, T ].We have previously shown [12] that the approximation error ϵ(t) = ∥v(y, t) − v r (y, t)∥ H(Ωy(t)) (15) can be controlled by selecting the rank r(t) at each t so that the norm of the normal component in ( 13) remains bounded by some user-defined threshold ε, i.e., 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 M r (t) with higher curvature.In fact, recall that the smallest singular value of v r is inversely proportional to the curvature of the manifold M r (t) at the point v r [23].

Variational principle for rank-reducing coordinate flows
Leveraging the geometric structure of the tensor manifold M r (t) described in Section 2.1, we aim at generating a coordinate flow y(x, t) that minimizes the rank of the PDE solution v r (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 ẏ 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 g of infinitesimal coordinate flow generators ẏ at time t Note that this functional is not optimizing the projection of ẏ • ∇v r onto the tangent space of the tensor manifold at v r .In principle, it is possible to control such projection, and eventually have the tangential component of Q 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 increases 3 , which may induce further tensor rank reduction via thresholded tensor truncation.
Proposition 2.1.Let f be a solution of the convex optimization problem (17).Then f satisfies the linear system of equations Proof: Due to the convexity of ( 17) any critical point of the cost functional in the Lie algebra g is necessarily a global solution to the optimization problem (17).The first variation of C with respect to ẏi is easily obtained as Ωy(t) Ωy(t) where η i ∈ g is an arbitrary perturbation.To obtain the second equality in (20) we used the fact that the orthogonal projection P N (v r ) is symmetric with respect to the inner product in L 2 (Ω y (t)) and idempotent.Setting (20) equal to zero for arbitrary η i yields the Euler-Lagrange equations Writing the preceding system of equations in vector notation proves the Proposition.□ Note that if ∇v r is non-zero then from the optimality conditions we obtain In this case the normal component of v r (y, t), i.e., the right hand side of ( 16), is zero along the optimal coordinate flow and therefore the rank of v r never increases during time integration.With the infinitesimal generator f (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 y(x, t) via the dynamical system Such coordinate transformation allows us to solve the PDE ( 14) along a coordinate flow that minimizes the projection of Q y onto the normal component of the tensor manifold in which the solution lives, henceforth minimizing the rank increase of v r during time integration.The corresponding evolution equation along the rank reducing coordinate flow is obtained by coupling the ODE ( 22) with ( 14) as 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 y that is normal to tensor manifold M r , and therefore control the rank increase of the PDE solution v 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 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 f (y, t k ) at time t k .Using the evolution equations ( 22) and ( 23) we then integrate both the coordinate flow y(x, t) and the PDE solution along the coordinate flow v r (y(x, t), t) from time t k to time t k+1 .Finally we update the PDE operator G 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 f (y, t k ), i.e., solving the linear system (18) at each time step.If we represent f (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 × d 2 system.The second challenge is representing the PDE operator in coordinates y(x, t), i.e., constructing G y .For a general nonlinear coordinate transformation y(x, t) the operator G y includes the metric tensor of the transformation which can significantly complicate the form of G 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.

Linear coordinate flows
Any coordinate flow y(x, t) that is linear in x at each time t can be represented by a time-dependent invertible Here we consider arbitrary matrices Γ(t) belonging to the collection of all4 d × d invertible matrices with real entries, denoted by GL d (R).The corresponding collection of infinitesimal coordinate flow generators is the vector space of all d × d matrices with real entries.In this setting the time-dependent tensor approximation of u(x, t) in coordinates 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 where C(•) is the cost function defined in (19).
Proposition 2.2.Let Σ(t) be a solution of the optimization problem (27) at time t.Then Σ(t) satisfies the symmetric where vec(Σ(t)) is a vectorization of the d × d matrix Σ(t), and c(y) is the d × d matrix with entries Proof: First note that the cost function C Γy is convex in Γ and the search space M d×d (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( Γy) with respect to each entry of the matrix Γ directly Ωy(t) for all i, j = 1, 2, . . ., d, where in the second line we used the fact that the orthogonal projection P N (v r ) is symmetric with respect to the inner product in L 2 (Ω y (t)) and idempotent.The critical points are then obtained by setting ∂C( Γ)/∂ Γij = 0 for all i, j = 1, 2, . . ., d, resulting in a linear system of equations for the entries of the infinitesimal linear coordinate flow generator Γkp Ωy(t) P N (v r ) y p ∂v r ∂y k y j ∂v r ∂y i dy = Ωy(t) for i = 1, 2, . . ., d.We obtain (28) by writing the linear system (32) in matrix notation.□ The linear system of equations ( 32) can also be obtained directly from the Euler Lagrange equations ( 21) by substituting ẏ = Γy and projecting onto y j for j = 1, 2, . . ., d.
A few remarks are in order regarding the linear system (28).The matrix A is a d 2 × d 2 symmetric matrix which is determined by (d 4 + d 2 )/2 entries and the vector b has length d 2 .Computing each entry of the matrix A and the vector b requires evaluating a d-dimensional integral.Therefore computing A and 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 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 A. If the matrix 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 The linear coordinate flow minimizing (33) aims at "undoing" the action of the operator G y in coordinates y.More precisely the PDE solution u(x, t) as seen in coordinate y(x, t) = Γ(t)x, i.e., v r (Γ(t)x, t) ≈ u(x, t), varies in time as little as possible.In the third equality of (33) we used the fact that P T (v r ) and P N (v r ) are orthogonal projections.Following similar steps used to prove Proposition 2.2 it is straightforward to show that the critical points Σ(t) of ( 33) satisfy the linear system (28) with With the optimal coordinate flow generator Σ(t) available we can compute the matrix Γ(t) appearing in (25) by solving the matrix differential equation Correspondingly, the PDE ( 23) can be written along the optimal rank-reducing linear coordinate flow (25) as 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 TDVP denotes a subroutine that solves the optimization problem ( 27) (e.g., by solving the symmetric d 2 ×d 2 linear system of equations ( 28)) and LMM denotes any time stepping scheme.

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 M r [13].

Liouville equation
Consider the d-dimensional Liouville equation where f (x) is a divergence-free vector field.We set the initial condition to be a product of independent Gaussians where .

Two-dimensional simulation results
We first consider the Liouville on a two-dimensional flat torus Ω ⊂ R 2 with velocity vector generated via the two-dimensional stream function [38]  with L = 30 and α = 4.73.This yields a two-dimensional measure-preserving (divergence-free) velocity vector f (x).
In the initial condition (38) we set β = 1/4 2 and t = 3 3 resulting in a rank-1 initial condition.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 ∆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 δ = 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 δ = 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 δ set on the solution and time step size ∆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.At time t = 35 we check the L ∞ 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 × 10 −5 .

Three-dimensional simulation results
Next we consider the Liouville equation (37) on a three-dimensional flat torus Ω ⊂ R 3 with velocity vector 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 β = 4 1/4 4 and t = 1 −1 1 resulting in the FTT rank r(0) = 1 1 1 1 .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 where 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 1 1 1 1 .
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 δ = 10 −5 , Cartesian coordinates, and time step size ∆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 Cartesian coordinates.For this example the coordinate flow is constructed by minimizing the functional (33).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 ∥r(t)∥ 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 ∞ 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.
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.

Fokker-Planck equation
Finally we consider the Fokker-Planck equation in three spatial dimensions (d = 3) on a flat torus Ω ⊆ R 3 .We set the drift vector f as in (42), and diffusion tensor with σ = 1/4.We use the same initial condition as in Section 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 ∆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 δ = 10 −5 , fixed Cartesian coordinates, and time step size ∆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).
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 ∞ 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.

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 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 ẏ • ∇v r in (17), or by simultaneously controlling the PDE solution and operator rank during time integration.

Figure 3 :
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).

Figure 4 :
Figure 4: Three-dimensional advection equation (37).(a) Rank versus time and (b) L ∞ error of the low-rank solutions relative to the benchmark solution.
f ) = u(Γx, t f ) → Solution at time t f on transformed coordinate system 1 Runtime: