Dynamically orthogonal tensor methods for high-dimensional nonlinear PDEs

TL;DR

This paper introduces dynamically orthogonal tensor methods that hierarchically decompose high-dimensional nonlinear PDEs, enabling efficient solution computation without rank reduction, demonstrated through numerical examples.

Contribution

The paper presents a novel hierarchical tensor decomposition approach with dynamic orthogonality for solving high-dimensional nonlinear PDEs, avoiding traditional rank reduction techniques.

Findings

Effective high-dimensional PDE solutions demonstrated

Hierarchical tensor formats enable efficient computation

New algorithms for mode management introduced

Abstract

We develop new dynamically orthogonal tensor methods to approximate multivariate functions and the solution of high-dimensional time-dependent nonlinear partial differential equations (PDEs). The key idea relies on a hierarchical decomposition of the approximation space obtained by splitting the independent variables of the problem into disjoint subsets. This process, which can be conveniently be visualized in terms of binary trees, yields series expansions analogous to the classical Tensor-Train and Hierarchical Tucker tensor formats. By enforcing dynamic orthogonality conditions at each level of binary tree, we obtain coupled evolution equations for the modes spanning each subspace within the hierarchical decomposition. This allows us to effectively compute the solution to high-dimensional time-dependent nonlinear PDEs on tensor manifolds of constant rank, with no need for rank reduction methods. We also propose new algorithms for dynamic addition and removal of modes within each subspace. Numerical examples are presented and discussed for high-dimensional hyperbolic and parabolic PDEs in bounded domains.

Figures

Equations

$$u:\Omega \to \mathbb{R}$$

$$H^k(\Omega )=\left\{u\in L^2(\Omega ):D^{\mathbf{\alpha }}u\in L^2(\Omega )\text{for all}|\mathbf{\alpha }|\le k\right\},k=0,1,2,\dots$$

$$D^{\mathbf{\alpha }}u=\frac{{\partial }^{|\mathbf{\alpha }|}u}{\partial x_1^{{\alpha }_1}\dots \partial x_d^{{\alpha }_d}},|\mathbf{\alpha }|={\alpha }_1+\cdots +{\alpha }_d.$$

$$⟨f,g{⟩}_{H^k(\Omega )}=\sum _{|\mathbf{\alpha }|\le k}{\int }_{\Omega }{D}^{\mathbf{\alpha }}f(\mathbf{x})D^{\mathbf{\alpha }}g(\mathbf{x})d{x}_1\dots d{x}_d.$$

$$\mathcal{H}\cong {\mathcal{H}}_1\otimes {\mathcal{H}}_2.$$

$$\Omega ={\Omega }^{(1,\dots ,p)}\times {\Omega }^{(p+1,\dots ,d)}$$

$$\underset{\text{in}\Omega }{\underbrace{(x_1,x_2,\dots ,x_d)}}=(\underset{\text{in}{\Omega }^{(1,\dots ,p)}}{\underbrace{(x_1,\dots ,x_p)}},\underset{\text{in}{\Omega }^{(p+1,\dots ,d)}}{\underbrace{(x_{p+1},\dots ,x_d)}}).$$

$$H^k(\Omega )\cong H^k\left({\Omega }^{(1,\dots ,p)}\right)\otimes H^k\left({\Omega }^{(p+1,\dots ,d)}\right).$$

$$⟨f,g{⟩}_{H^k\left({\Omega }^{(1,\dots ,p)}\right)}=\sum _{{\alpha }_1+\cdots +{\alpha }_p\le k}{\int }_{{\Omega }^{(1,\dots ,p)}}\frac{{\partial }^{{\alpha }_1+\cdots +{\alpha }_p}f}{\partial x_1^{{\alpha }_1}\dots \partial x_p^{{\alpha }_p}}\frac{{\partial }^{{\alpha }_1+\cdots +{\alpha }_p}g}{\partial x_1^{{\alpha }_1}\dots \partial x_p^{{\alpha }_p}}d{x}_1\cdots d{x}_p,$$

$$⟨f,g{⟩}_{H^k\left({\Omega }^{(p+1,\dots ,d)}\right)}=\sum _{{\alpha }_{p+1}+\cdots +{\alpha }_d\le k}{\int }_{{\Omega }^{(p+1,\dots ,d)}}\frac{{\partial }^{{\alpha }_{p+1}+\cdots +{\alpha }_d}f}{\partial x_{p+1}^{{\alpha }_{p+1}}\dots \partial x_d^{{\alpha }_d}}\frac{{\partial }^{{\alpha }_{p+1}+\cdots +{\alpha }_d}g}{\partial x_{p+1}^{{\alpha }_{p+1}}\dots \partial x_d^{{\alpha }_d}}d{x}_{p+1}\cdots d{x}_d.$$

$$u(x_1,\dots ,x_d)=\sum _{i,j=1}^{\infty }{a}_{ij}{\phi }_i^{(1,\dots ,p)}(x_1,\dots ,x_p){\phi }_j^{(p+1,\dots ,d)}(x_{p+1},\dots ,x_d),$$

$$U_u:H^k\left({\Omega }^{(1,\dots ,p)}\right)\to H^k\left({\Omega }^{(p+1,\dots ,d)}\right),$$

$$\left(U_u\psi \right)(x_{p+1},\dots ,x_d)=⟨u,\psi {⟩}_{H^k\left({\Omega }^{(1,\dots ,p)}\right)}.$$

$$⟨U_u\psi ,\phi {⟩}_{H^k\left({\Omega }^{(p+1,\dots ,d)}\right)}=⟨\psi ,U_u^{†}\phi {⟩}_{H^k\left({\Omega }^{(1,\dots ,p)}\right)},$$

$$U_u^{†}:H^k\left({\Omega }^{(p+1,\dots ,d)}\right)\to H^k\left({\Omega }^{(1,\dots ,p)}\right),$$

$$\left(U_u^{†}\phi \right)(x_1,\dots ,x_p)=⟨u,\phi {⟩}_{H^k\left({\Omega }^{(p+1,\dots ,d)}\right)},$$

$$R_u=U_u{U}_u^{†},R_u:H^k({\Omega }^{(p+1,\dots ,d)})\to H^k({\Omega }^{(p+1,\dots ,d)}),$$

$$L_u=U_u^{†}{U}_u,L_u:H^k({\Omega }^{(1,\dots ,p)})\to H^k({\Omega }^{(1,\dots ,p)}).$$

$$\left(R_u\phi \right)(x_{p+1},\dots ,x_d)=⟨r_u,\phi {⟩}_{H^k\left({\Omega }^{(p+1,\dots ,d)}\right)}\phi \in H^k({\Omega }^{(p+1,\dots ,d)}),$$

$$r_u(x_{p+1},\dots ,x_d,x_{p+1}^{\prime },\dots ,x_d^{\prime })=⟨u,u{⟩}_{H^k\left({\Omega }^{(1,\dots ,p)}\right)}.$$

$$\left(L_u\psi \right)(x_1,\dots ,x_p)=⟨l_u,\psi {⟩}_{H^k\left({\Omega }^{(1,\dots ,p)}\right)}\psi \in H^k({\Omega }^{(1,\dots ,p)}),$$

$$l_u(x_1,\dots ,x_p,x_1^{\prime },\dots ,x_p^{\prime })=⟨u,u{⟩}_{H^k\left({\Omega }^{(p+1,\dots ,d)}\right)}.$$

$$u=\sum _{k=1}^{\infty }{\lambda }_k{\psi }_k^{(1,\dots ,p)}{\psi }_k^{(p+1,\dots ,d)},$$

$$\left[\begin{array}[]{cc}U_{u}&0\\ 0&U_{u}^{\dagger}\end{array}\right]\left[\begin{array}[]{c}\psi_{k}^{(1,\ldots,p)}\\ \psi_{k}^{(p+1,\ldots,d)}\end{array}\right]=\lambda_{k}\left[\begin{array}[]{cc}0&1\\ 1&0\end{array}\right]\left[\begin{array}[]{c}\psi_{k}^{(1,\ldots,p)}\\ \psi_{k}^{(p+1,\ldots,d)}\end{array}\right].$$

$$⟨{\psi }_i^{(1,\dots ,p)}{\psi }_j^{(1,\dots ,p)}{⟩}_{H^k\left({\Omega }^{(1,\dots ,p)}\right)}=⟨{\psi }_i^{(p+1,\dots ,d)}{\psi }_j^{(p+1,\dots ,d)}{⟩}_{H^k\left({\Omega }^{(p+1,\dots ,d)}\right)}={\delta }_{ij}.$$

$$L_u{\psi }_k^{(1,\dots ,p)}$$

$$R_u{\psi }_k^{(p+1,\dots ,d)}$$

$${\psi }_k^{(p+1,\dots ,d)}=\frac{1}{{\lambda }_k}{U}_u{\psi }_k^{(1,\dots ,p)},{\psi }_k^{(1,\dots ,p)}=\frac{1}{{\lambda }_k}{U}_u^{†}{\psi }_k^{(p+1,\dots ,d)}.$$

$$H^k(\Omega )=$$

$$=$$

$$\cdots$$

$$u$$

Full text

††journal: arXiv

1 Introduction

High-dimensional partial differential equations (PDEs) arise in many areas of engineering, physical sciences and mathematics. Classical examples are equations involving probability density functions (PDFs) such as the Fokker-Plank equation [55], the Liouville equation [64, 18], or the Boltzmann equation [12, 44]. Other types of high-dimensional PDEs can be obtained as finite-dimensional approximations of functional differential equations [63], such as the Hopf equation of turbulence [30, 31, 46], the Schwinger-Dyson equation of quantum mechanics [32], or the Martin-Siggia-Rose formulation of classical statistical dynamics [45, 33, 34, 50]. Computing the solution to high-dimensional PDEs is a challenging problem that requires approximating high-dimensional functions, i.e., the solution to the PDE, and then developing appropriate numerical schemes to compute such functions accurately. Classical numerical methods based on tensor product representations are not viable in high-dimensions, as the number of degrees of freedom grows exponentially fast with the dimension. To address this problem there have been substantial research efforts in recent years on high-dimensional numerical approximation theory. Techniques such as sparse collocation [10, 17, 7, 23, 47], high-dimensional model representations (HDMR) [42, 11, 6] and, more recently, deep neural networks [51, 52, 66] and tensor methods [37, 5, 9, 28, 16, 39] were proposed to mitigate the exponential growth of the degrees of freedom, the computational cost and memory requirements.

In this paper, we develop a new dynamically orthogonal tensor method to approximate multivariate functions and the solution of high-dimensional time-dependent nonlinear PDEs. The key idea relies on a hierarchical decomposition of the function space in terms of a sequence of nested subspaces of smaller dimension. Such decomposition is induced by by splitting the independent variables of the problem recursively into two disjoint subsets which can conveniently be visualized by binary trees. In particular, we study two classes of trees which are analogous to the Tensor-Train (TT) [48] and Hierarchical Tucker (HT) [26] tensor formats. By enforcing dynamic orthogonality (DO) [57] or bi-orthogonality (BO) [13] conditions at each level of the TT or the HT binary tree, we obtain coupled evolution equations for the modes spanning each subspace within the hierarchy. This allows us to represent the time evolution of high-dimensional functions and compute the solution of high-dimensional time-dependent nonlinear PDEs on a tensor manifold with constant rank. This formulation has several advantages over classical numerical tensor methods. In particular, the hard-to-compute nonlinear projection [43, 38] that maps the solution of high-dimensional PDEs onto a tensor manifold with constant rank [60] here is represented explicitly by the hierarchical DO/BO propagator111The hierarchical DO/BO propagator is nonlinear even for linear PDEs. Such nonlinearity implicitly represnts the projection onto a tensor manifold with constant rank., i.e., by a system of coupled one-dimensional nonlinear PDEs. In other words, there is no need to perform tensor rank reduction [24, 26, 40], rank-constrained temporal integration [43, 38], or Riemannian optimization [59], when solving high-dimensional PDEs with the hierarchical subspace decomposition method we propose222Classical numerical tensor methods for high-dimensional PDEs with explicit time stepping schemes require rank-reduction to project the solution back into the tensor manifold with prescribed rank (see [63] §5.5), the so-called retraction step [59]. This can be achieved, e.g., by a sequence of suitable matricizations followed by hierarchical singular value decomposition [24, 26, 40], or by optimization [59, 39, 9, 20, 56, 35]. Rank reduction can be computationally intensive, especially if performed at each time step. Tensor methods with implicit time stepping suffer from similar issues. In particular, the nonlinear system that yields the solution at the next time step needs to be solved on a tensor manifold with constant rank by using, e.g., Riemannian optimization algorithms [59, 60, 22]..

This paper is organized as follows. In Section 2 we introduce a recursive bi-orthogonal decomposition method for time-independent multivariate functions, develop error estimates and provide simple examples of application. In Section 3 we extend the recursive bi-orthogonal decomposition to time-dependent functions, and develop a hierarchy of nested time-dependent orthogonal projections generalizing the DO and BO conditions to tensor formats with multiple levels. We also prove that the approximations resulting from the BO and the DO conditions are equivalent, in the sense that they span the same function spaces. In Section 4, we apply the recursive subspace decomposition method to compute the solution of high-dimensional nonlinear PDEs. In Section 5 we provide numerical examples demonstrating the accuracy and computational effectiveness of the recursive subspace decomposition method we propose. Specifically we study high-dimensional hyperbolic and parabolic PDEs. The main findings are summarized in Section 6.

2 Recursive bi-orthogonal decomposition of time-independent multivariate functions

Let $\Omega$ be a subset of $\mathbb{R}^{d}$ ($d\geq 2$) that contains an open set333If $\Omega\subseteq\mathbb{R}^{d}$ contains an open set then $\text{dim}(\Omega)=d$., and let

[TABLE]

be a multivariate function which we assume to be an element of a separable Hilbert space $\mathcal{H}(\Omega)$. A possible choice of such Hilbert space is the Sobolev space

[TABLE]

where $\bm{\alpha}=(\alpha_{1},\dots,\alpha_{d})$ is a multi-index and

[TABLE]

Note that (2) includes the classical Lebesgue space $L^{2}(\Omega)=H^{0}(\Omega)$. We equip (2) with the standard inner product

[TABLE]

If needed, this inner product can be weighted by a non-negative separable density $\rho_{1}(x_{1})\cdots\rho_{d}(x_{d})$. Any separable Hilbert space is isomorphic to $L^{2}$, and it can be represented as a tensor product of two Hilbert spaces [53, p.51], i.e.,

[TABLE]

The spaces $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ may be specified by partitioning the spatial variables $\{x_{1},\ldots,x_{d}\}$ into two disjoint subsets. This is equivalent to represent the domain $\Omega$ as a Cartesian product of two sub-domains (whenever possible). For instance, consider the partition

[TABLE]

induced by the following splitting of the spatial variables

[TABLE]

In this setting, the Sobolev space (2) admits the following decomposition

[TABLE]

The inner products within each subspace $H^{k}\left(\Omega^{(1,\ldots,p)}\right)$ and $H^{k}\left(\Omega^{(p+1,\ldots,d)}\right)$ can be defined, respectively, as

[TABLE]

and

[TABLE]

A representation of the multivariate function (1) in the tensor product space (8) has the general form

[TABLE]

where $\varphi_{i}^{(1,\ldots,p)}$ and $\varphi_{j}^{(p+1,\ldots,d)}$ are orthonormal basis functions in $H^{k}(\Omega^{(1,\ldots,p)})$ and $H^{k}(\Omega^{(p+1,\ldots,d)})$, respectively444Orthonormality is relative to the inner products in $H^{k}(\Omega^{(1,\ldots,p)})$ and $H^{k}(\Omega^{(p+1,\ldots,d)})$.. The superscripts in (11) denote which spatial components the function depends on. This will be the case throughout this paper and for notational simplicity, the spatial arguments will be often omitted when there is no ambiguity.

With the isomorphism (8) and the inner products (9)-(10) set, it is straightforward to develop an operator framework which guarantees the existence of a diagonalized bi-orthogonal representation of the field $u(x_{1},\ldots,x_{d})$. To this end, following Aubry et. al. [2, 3, 1] and Venturi [62] (see also [61, 65]), we define the integral operator

[TABLE]

The formal adjoint of $U_{u}$, denoted by $U^{\dagger}$, is a linear operator defined by the requirement

[TABLE]

for all $\psi\in H^{k}\left(\Omega^{(1,\ldots,p)}\right)$, and all $\varphi\in H^{k}\left(\Omega^{(p+1,\ldots,d)}\right)$. By using integration by parts and discarding boundary conditions (formal adjoint operator) we obtain

[TABLE]

The subscript “$u$” in $U_{u}$ and $U^{\dagger}_{u}$ identifies the kernel of the integral operators. We will shortly define a hierarchy of such operators and it will be important to distinguish them by their kernels. Next, we introduce the following correlation operators

[TABLE]

and

[TABLE]

Note that $L_{u}$ and $R_{u}$ are self-adjont relative to (9) and (10), respectively. Moreover, if $U_{u}$ is compact (e.g., if we consider a decomposition in $H^{0}=L_{2}$), then $U_{u}^{\dagger}$ is compact, and therefore $L_{u}$ and $R_{u}$ are compact. Hence, by the Riesz-Schauder theorem, they have the same discrete spectra (see e.g. [36, p.185]). By a direct calculation, it can be show that

[TABLE]

where the correlation function $r_{u}$ is defined by

[TABLE]

Similarly,

[TABLE]

where the correlation function $l_{u}$ is defined by

[TABLE]

It is a classical result in the spectral theory of compact operators (see e.g., [2, 1]) that there exists a canonical decomposition of the field (1) of the form

[TABLE]

where the modes $\psi_{k}^{(1,\ldots,p)}$ and $\psi_{k}^{(p+1,\ldots,d)}$ satisfy the eigenvalue problem

[TABLE]

Moreover, it can be shown that $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$, and

[TABLE]

The series (21) is usually called bi-orthogonal (or Schmidt) decomposition of the multivariate function $u$, and it converges in norm. The modes $\psi_{k}^{(1,\ldots,p)}$ are eigenfunctions of the operator $L_{u}$ with corresponding eigenvalue $\lambda_{k}^{2}$, while the modes $\psi_{k}^{(p+1,\ldots,d)}$ are eigenfunctions of the operator $R_{u}$ with corresponding eigenvalues $\lambda_{k}^{2}$, i.e.,

[TABLE]

In practice, since $\psi_{k}^{(1,\ldots,p)}$ and $\psi_{k}^{(p+1,\ldots,d)}$ are determined up to two unitary transformations [3], to compute (21) we can solve one of the two eigenvalue problems in (24) (the one with smaller dimension), and then use one of the dispersion relations (22), i.e.,

[TABLE]

2.1 Hierarchical subspace decomposition and tensor formats

To obtain a series expansion of the multivariate function $u(x_{1},\ldots,x_{d})$ in terms of univariate functions, we apply the bi-orthogonal decomposition method discussed in the previous Section recursively. The way in which the variables are split in each step of the recursive procedure, i.e., the choice of $p$ in (8), can be conveniently visualized by binary trees. In Figure 1 we provide two simple examples of such binary trees corresponding to the Tensor-Train (TT) [48] and the Hierarchical Tucker (HT) [26] tensor formats (see also [28, 63], and the references therein).

2.1.1 Tensor Train (TT) format

The Tensor-Train format singles out one variable at a time, resulting in a binary tree with depth $(d-1)$ when decomposing $d$-variate functions $u(x_{1},\ldots,x_{d})$ (see Figure 1). This corresponds to the following hierarchical subspace decomposition of the Sobolev space (2)

[TABLE]

in which we diagonalize each tensor product representation using the bi-orthogonal decomposition method. This yields the following TT expansion of the multivariate function $u(x_{1},\ldots,x_{d})$

[TABLE]

i.e.,

[TABLE]

Each of the bi-orthogonal modes can be obtained by solving a sequence of one-dimensional eigenvalue problems followed by projections. Specifically, the eigenvalue problems are

[TABLE]

(see Eq. (19)) while the corresponding projections are defined as

[TABLE]

2.1.2 Hierarchical Tucker (HT) format

The Hierarchical Tucker format splits variables into disjoint subsets of equal size, whenever possible. In the case $d$-variate functions $u(x_{1},\ldots,x_{d})$, where $d=2^{n}$ for some natural number $n$, the tree is balanced. In general, the Hierarchical Tucker tree is more shallow than a Tensor Train tree for the same number of variables $d$. In fact, the depth of the HT tree for $d=2^{n}$ is $n=\log_{2}(d)$, while the corresponding TT tree has depth $2^{n}-1$. The HT format is based on the following hierarchical decomposition of the Sobolev space $H^{k}(\Omega)$

[TABLE]

As before, we diagonalize each tensor product representation as we proceed splitting variables down the tree. This yields the following sequence of bi-orthogonal decompositions

[TABLE]

and the expansion

[TABLE]

Similar to the TT format, a sequence of eigenfunction problems followed by projections are used to obtain the modes spanning the hierarchical subspaces. However, in the HT case the eigenfunction problems are higher dimensional and not tractable for large $d$.

Remark

Clearly, we may decompose a multivariate function $u(x_{1},\ldots,x_{d})$ by splitting variables in various ways at different levels of the decompositions. Any binary tree which has leaves containing one index leads to a series expansion in terms of functions of one spatial variable. Separable Hilbert spaces defined on a Cartesian product of one-dimensional domains always allow such reduction.

2.2 Error analysis

In this Section we develop an error analysis for the recursive biorthogonal decomposition we discussed in Section 2.1. To this end we first state a Lemma which will be useful in Section 2.4 for establishing a thresholding criterion to truncate the infinite sums in (30) and (38).

Lemma 2.1**.**

If $u\in H^{k}(\Omega)$ admits the bi-orthogonal expansion

[TABLE]

then $\displaystyle\sum_{i=1}^{\infty}\lambda_{i}^{2}=\|u\|^{2}_{H^{k}(\Omega)}$.

Proof.

This result follows immediately from the orthonormality of the modes $\psi_{i}^{(1,\ldots,p)}$ and $\psi_{i}^{(p+1,\ldots,d)}$ relative to the inner products (9)-(10).

∎

Next, we analyze the error in the $H^{k}(\Omega)$ norm between the Tensor Train series expansion (30) and the truncated expansion

[TABLE]

where $r_{1},r_{2},\ldots,r_{d-1}$ are truncation ranks. To simplify indexing and array bounds for truncated TT expansions such as (40), we will omit the array indices in the rank arrays and write, e.g., $\psi^{(j_{1},\ldots,j_{p})}_{i_{1}\cdots i_{k}},\ k=1,\ldots,r_{k}$ instead of $k=1,\ldots,r_{k}(i_{1},\ldots,i_{k-1})$, since the rank array indices are clear from the subscripts of the mode $\psi^{(j_{1},\ldots,j_{p})}_{i_{1}\cdots i_{k}}$. In this simplified notation, the truncated TT expansion (40) can be written as

[TABLE]

Proposition 2.1**.**

Let $u\in H^{k}(\Omega)$. The error incurred by truncating the infinite expanion (30) to the finite expansion (41) is given by

[TABLE]

Proof.

Let us rewrite (30) as

[TABLE]

and split each infinite sum into the superimposition of a finite sum and an infinite sum, i.e.,

[TABLE]

Expanding the products in (44) yields the following expression

[TABLE]

Using the orthogonality of each set of modes we obtain

[TABLE]

i.e.,

[TABLE]

∎

Remark

Proposition 2.1 can be generalized to tensor formats corresponding to arbitrary binary trees, e.g., the HT format (38). In some sense, the equality (42) represents the infinite-dimensional version of well-known finite-dimensional results which bound the overall squared approximation error of multilinear singular value decompositions in the 2-norm by the sum (over the whole tree) of squares of deleted singular values. These types of results were first proven by De Lathauwer et al. in [41], and later generalized by Grasedyck [24] (see also Schneider and Uschmajew [58]).

Remark

Recent error estimates by Griebel and Li [27] on the decay rate of singular values allow us to develop sharp upper bounds for (42) depending only on the multivariate rank, the smoothness of the function $u(x_{1},\ldots,x_{d})$ and other computable quantities. To obtain such estimates, it is sufficient to bound each eigenvalue $\lambda_{i_{1}}^{2}$, $\lambda_{i_{1}i_{2}}^{2}$, etc., with the corresponding sharp upper bound recently obtained in [27].

2.3 Computational aspects of TT and HT series expansions

To compute a recursive bi-orthogonal decomposition of the multivariate function $u(x_{1},\ldots x_{d})$ we first need to identify the self-adjoint operator (19) (or (17)) at each level of the binary tree (see Figure 1). In other words, we need to compute the kernel (20), (or (18)) and then solve the corresponding eigenfunction problem. In the HT tensor format, computing such kernel requires evaluating a multivariate integral of dimension $d/2$ (at the first level of the tree) and then solving an eigenfunction problem of dimension $d/2$, which can be extremely challenging when $d$ is large. In the TT tensor format, this problem can be mitigated substantially. In fact, at the first level of the TT tree we have that the kernel of $L_{u}$ is a $(d-1)$-dimensional integral which can be evaluated, e.g., by using Quasi-Monte Carlo or more general lattice cubature rules [21]. The corresponding eigenvalue problem (31) and the projection (32) are both one-dimensional. From a numerical viewpoint, this is extremely advantageous, as we can accurately solve one dimensional eigenvalue problems in a collocation or a Galerkin setting at a low computational cost, once the kernels of the operators in (31) are available.

2.4 Thresholding hierarchical bi-orthogonal series expansions

In this Section we develop a new thresholding criterion to truncate the series expansion (30) and (38) to finite rank. To this end, we first notice that the amplitude of each term in the series is represented by products of eigenvalues from each level of the binary tree, since all eigenmodes are normalized. With this in mind, it is clear that a reasonable criterion to truncate the bi-orthogonal series expansion (30) (or (38)) to finite rank is to ensure that each of these eigenvalue products remains above a specified threshold $\sigma$. Hereafter, we develop this criterion for the TT format (30). The same technique can be applied to any other tensor format. We begin by setting some threshold value $\sigma$ for which we enforce $\lambda_{i_{1}}\lambda_{i_{1}i_{2}}\cdots\lambda_{i_{1}\cdots i_{d-1}}\geq\sigma$. In the first level of the TT tree (Figure (1)), i.e., Eq. (26), we keep all modes with eigenvalues $\lambda_{i_{1}}\geq\sigma$, of which there will be a finite number $r_{1}$ because of property (23). Then we proceed to the second level of the TT tree and decompose $\psi^{(2,\ldots,d)}_{i_{1}}$ ($1\leq i_{1}\leq r_{1}$) as in (27). Here we set new thresholds $\sigma_{i_{1}}=\sigma/\lambda_{i_{1}}$ and keep all modes $\psi^{(2)}_{i_{1}i_{2}},\psi^{(3,\ldots,d)}_{i_{1}i_{2}}$ with eigenvalues $\lambda_{i_{1}i_{2}}\geq\sigma_{i_{1}}$. Proceeding recursively in this way down to the $j^{\text{th}}$ level of the TT tree we have the thresholds $\sigma_{i_{1}\cdots i_{j}}=\sigma_{i_{1}\cdots i_{j-1}}/\lambda_{i_{1}\cdots i_{j}}$. It is reasonable to disregard modes corresponding to eigenvalues smaller than $\sigma$ in the first bi-orthogonal decomposition since

[TABLE]

for all $j=2,\ldots,d-1$. Indeed, Lemma 2.1 implies that $\lambda_{i_{1}\cdots i_{j}}\leq 1$ for all $j=2,\ldots,d-1$ from which (46) immediately follows. Another desirable consequence of Lemma 2.1 is that $\sigma_{i_{1}\cdots i_{j-1}}\leq\sigma_{i_{1}\cdots i_{j}}$ for all $j=2,\ldots,d-1$. As a result, bi-orthogonal decompositions at different levels of the binary trees are truncated to a different number of modes. Let us now summarize the thresholding algorithm for Tensor Train formats. On the first level of the tree we decompose $u(x_{1},\ldots,x_{d})$ as in (26), for which we keep $r_{1}$ modes, identified by the criterion $\lambda_{i_{1}}\geq\sigma$. For each of the modes $\psi_{i_{1}}^{(2,\ldots,d)}$ ($i_{1}=1,\ldots,r_{1}$), we perform the decomposition (27) on the second level of the binary tree with mode-specific thresholds $\sigma_{i_{1}}=\sigma/\lambda_{i_{1}}$. Hence, the bi-orthogonal decomposition of $\psi_{i_{1}}^{(2,\ldots,d)}$, has $r_{2}(i_{1})$ modes, i.e., the ranks in the second level are described by the vector $r_{2}$. For each of the modes555The total number of modes $\psi^{(3,\ldots,d)}_{i_{1}i_{2}}$ is $r_{1}\displaystyle\sum_{i_{1}=1}^{r_{1}}r_{2}(i_{1})$. $\psi^{(3,\ldots,d)}_{i_{1}i_{2}}$, the decomposition at the third level of the tree is performed with thresholds $\sigma_{i_{1}i_{2}}=\sigma/\lambda_{i_{1}i_{2}}$. This yields the truncation rank $r_{3}(i_{1},i_{2})$ in the decomposition of $\psi_{i_{1}i_{2}}^{(3,\ldots,d)}$. Thus, the bi-orthogonal ranks for the third level of the TT tree are described by a matrix. In general, on the $j^{\text{th}}$ level of the tree, the rank is described by a tensor $r_{j}(i_{1},\ldots,i_{j-1})$ of dimension $j-1$.

2.4.1 An example: recursive bi-orthogonal decomposition of a 3D function

In this Section, we apply the recursive bi-orthogonal decomposition method to a simple a three-dimensional function666For two- and three-dimensional functions TT and HT tensor formats are equivalent. defined on the cube $\Omega=[-1,1]^{3}$. Specifically, we consider

[TABLE]

This function is shown in Figure 2 together with the binary tree representing the tensor format we use in the recursive bi-orthogonal decomposition.

We discretize (47) on a three-dimensional tensor product grid with $50$ Gauss-Legendre collocation points [29] in each variable (125000 points total). Regarding the function space in which we perform the decomposition, in this example we set $k=0$ in (2), i.e., we consider the classical $L^{2}(\Omega)=H^{0}(\Omega)$ function space. In this setting, the kernel of the integral operator $L_{u}$ in (31) (first level of the binary tree), reduces to

[TABLE]

This integral is computed with the Gauss-Legendre quadrature rule corresponding to the chosen grid points. The $x_{1}$-modes are solutions of the eigenvalue problem

[TABLE]

We discretize (49) using Gauss-Legendre spectral collocation with 50 points. This yields $50$ leading eigenvalues and corresponding eigenfunctions. Following the thresholding technique discussed in Section 2.4, we set $\sigma=10^{-5}$ to determine how many level-1 eigenvalues and eigenfunctions to keep. It turns out that only $9$ eigenvalues are larger than $\sigma$ which determines the first truncation rank as $r_{1}=9$. These $9$ eigenvalues $\{\lambda_{1},\ldots,\lambda_{9}\}$ constitute the level $1$ spectrum which is shown in Figure 4(a). The corresponding eigenfunctions are $\psi^{(1)}_{1},\ldots,\psi^{(1)}_{9}$. The modes $\psi^{(2,3)}_{1},\ldots,\psi^{(2,3)}_{9}$ can now be obtained through projection as in (32)

[TABLE]

which we compute with Gauss-Legendre quadrature. For each of the $9$ modes $\psi^{(2,3)}_{i_{1}}$ we follow the same procedure we used to compute $\psi^{(1)}_{i_{1}}$, i.e., we build the kernels

[TABLE]

and then solve the eigenvalue problems

[TABLE]

to obtain $50$ eigenvalues and eigenfunctions of for each ${i_{1}}=1,\ldots,9$. To decide how many eigenvalues and eigenfunctions to keep we use the thresholds $\sigma_{i_{1}}=\sigma/\lambda_{i_{1}}$ (see Figure 4(b)). This yields the following vector of level-$2$ ranks $r_{2}=[11,11,11,11,11,11,10,6,0]$. The level-$2$ spectra $\lambda_{i_{1}i_{2}}$ is shown in Figure 4(c) for the level-$1$ modes above the threshold $\sigma=10^{-5}$.

3 Recursive bi-orthogonal decomposition of time-dependent multivariate functions

Let us consider the time-dependent multivariate function

[TABLE]

At any fixed $t\in[0,T]$ we assume that $u(x_{1},\dots,x_{d},t)$ is an element of the Sobolev space (2). With the goal of solving high-dimensional nonlinear evolution equations (PDEs), in this Section we develop a recursive bi-orthogonal decomposition theory for (53). To decompose the function $u(x_{1},\ldots,x_{d},t)$ with recursive bi-orthogonal methods we have at least two different options. For instance, if $u$ is square integrable in $t$, then we may treat $t$ the same as a spatial variable $x_{j}$. In this way, we essentially include the time variable in the inner products (9)-(10), and decompose the time-dependent field (53) using the methods of Section 2.1.

An alternative approach relies on introducing time dependence in all modes appearing in the tensor expansion of $u(x_{1},\ldots,x_{d},t)$, and then deriving problem-dependent evolution equations for each mode. To this end, one can enforce, e.g., a dynamic orthogonality (DO) condition or a bi-orthogonality (BO) condition [14, 15, 57, 19, 4] at each level of the binary tree representing the tensor expansion. This generalizes the DO/BO framework originally proposed by Sapsis and Lermusiaux [57], and Cheng et. al. [14], which holds for binary trees with only one level, to binary trees with an arbitrary number of levels. The DO/BO method was originally proposed to solve initial/boundary value problems for nonlinear PDEs with parametric uncertainty modeled by the random vector $\bm{\xi}(\omega)$. The key idea was to introduce time-redundancy in a Karhunen-Loève-type expansion of the solution for the purpose of efficiently representing the time evolution of the stochastic modes and the corresponding space-time modes. As we mentioned above, the classical DO/BO expansion corresponds to a binary tree with only one level, where the random field $u(\bm{x},\bm{\xi},t)$ is decomposed as

[TABLE]

In this Section we generalize this idea, and apply it recursively to multivariate functions not necessarily dependent on random parameters, until we obtain an expansion in terms of one-dimensional functions (whenever possible). To illustrate the hierarchical DO/BO method, we consider the TT series expansion (30), and introduce redundant time-dependence in all eigenvalues and eigenmodes. This yields the following representation of (53)

[TABLE]

Other tensor expansions corresponding to different binary trees, e.g., the HT expansion sketched in Figure 1, can be generalized in a similar way. Hereafter, we derive the DO and BO evolution equations for the time-dependent modes in (55), and show that these two ways of propagating $u(x_{1},\ldots,x_{d},t)$ forward in time on a (smooth) low-dimensional manifold of constant rank [60] are equivalent in the sense that the finite dimensional function spaces containing the DO and BO components are essentially the same.

In the following Sections every function is time dependent, so $t$ is omitted from the function arguments. Superscripts indicate spatial dependencies as in Section 2, so spatial arguments of functions are also omitted when there is no ambiguity. Angled brackets $\langle\cdot,\cdot\rangle$ here denote the $L^{2}=H^{(0)}$ inner product over all spatial components for which the two arguments are defined.

3.1 DO-TT propagator

To derive the Dynamically-Orthogonal Tensor-Train (DO-TT) propagator, let us first consider the level-1 expansion

[TABLE]

where

[TABLE]

By differentiating (56) with respect to time we obtain

[TABLE]

Clearly, if $u(x_{1},\ldots,x_{d},t)$ is given then $\partial u/\partial t$ is known, and therefore the left hand side of (58) is fully determined. On the other hand, in the context of nonlinear evolution equations, $\partial u/\partial t$ is represented by the right hand side of the PDE. Applying $\langle\cdot,\psi_{k_{1}}^{(1)}\rangle$ to (58) and utlizing the DO conditions (see [57])

[TABLE]

yields

[TABLE]

The evolution equations for the modes $\psi_{k_{1}}^{(1)}$ can be obtained by applying $\langle\cdot,\Psi_{k_{1}}^{(2,\ldots d)}\rangle$ to (58) and (60), and the using simple algebra. This yields

[TABLE]

Equations (60) and (61) can be conveniently expressed in a matrix-vector form as

[TABLE]

where

[TABLE]

[TABLE]

Note that $\bm{M}^{(1)}(t)$ and $\bm{C}(t)$ depend on $\Psi^{(2,\ldots,d)}$ and $\partial u/\partial t$, while $\bm{N}^{(2,\ldots,d)}(t)$ depends on $\bm{\psi}^{(1)}$ and $\partial u/\partial t$. Therefore, given $\partial u/\partial t$ we have that the system (62) is closed. At this point, we move to the next level of the TT binary tree and derive DO evolution equations for the level-$2$ modes. By following the same steps as in the derivation of the level-$1$ system (62), we obtain

[TABLE]

where we defined

[TABLE]

Proceeding recursively, it is possible to obtain evolution equations for each mode in the TT binary tree sketched in Figure 1. Specifically, we have

[TABLE]

Equations (68)-(69) can be conveniently expressed in a matrix-vector form as

[TABLE]

where,

[TABLE]

[TABLE]

Remark

The “non-leaf” modes $\psi_{k_{1}\cdots k_{j}}^{(j+1,\ldots,d)}$ can be constructed at any time in terms of the “leaf” modes as

[TABLE]

To build the time-dependent multivariate function (55) it is sufficient to integrate only the evolution equations corresponding to the leaf modes.

3.2 An example: DO-TT decomposition of a time-dependent 3D function

Let us consider the time-dependent multivariate function

[TABLE]

where $\Omega=[-1,1]^{3}$ is the standard three-dimensional cube. We are interested in decomposing (74) with the DO-TT method. To this end, we consider the $L^{2}$ inner product and first perform a recursive bi-orthogonal decomposition of the initial state $u(x_{1},x_{2},x_{3},0)$ using the methods of Section 2.1. To this end, we consider $50$ Gauss-Legendre quadrature points and set the eigenvalue threshold to $\sigma=10^{-5}$. By following the same steps as in the example 2.4.1, this yields multivariate ranks $r_{1}=2$ and $r_{2}=[1\,1]$. This allows us to approximate the initial condition as

[TABLE]

The time derivative of $u$ is easily obtained as

[TABLE]

A substitution of (76) intro the multi-level DO evolution equations (70) yields

[TABLE]

for $j=1,2$. In a spectral collocation setting, all modes are represented by their values at the 50 Gauss-Legendre collocation points. Hence, the DO evolution equations reduce to a system of ODEs. All integrals in (77)-(79) are computed by using the one-dimensional Gauss-Legendre quadrature rule. In Figure 5 we plot the temporal evolution of all $6$ modes appearing in the system (77)-(79). The time-dependent $L^{2}(\Omega)$ error between the DO-TT expansion

[TABLE]

and the function (74) is plotted in Figure 6 versus time. We notice that (74) is a separable function with rank $r_{1}=2$, $r_{2}=[1,1]$ at each time $t$. This means that it can be represented exactly on a low-dimensional tensor manifold of constant rank at each time. This is the reason why the $L^{2}$ error in Figure 6 is of order $10^{-8}$ for all $t\geq 0$. In general, the multivariate function $u(x_{1},\ldots,x_{d},t)$ is not separable and may become rougher/wavier as time increases. In these cases, the multivariate rank $(r_{1},\ldots,r_{d-1})$ of the DO-TT expansion usually needs to be increased in time to accurately represent represent $u(x_{1},\ldots,x_{d},t)$. Methods for increasing rank will be addressed in Section 4.1 in the context of adaptive DO-TT approximation of the solution to high-dimensional PDEs.

3.3 BO-TT propagator

Alongside the multi-level DO propagator (70), we derive the Bi-Orthogonal Tensor-Train (BO-TT) propagator. The key idea is to replace the dynamic orthogonality condition between the modes at each level of the binary tree, with a bi-orthogonality condition (see [14]). To this end, let us consider the following sequence of bi-orthogonal decompositions

[TABLE]

where the modes $\varphi_{i_{1}\cdots i_{j-1}k}^{(j)}(t)$ and $\varphi_{i_{1}\cdots i_{j-1}k}^{(j+1,\ldots,d)}(t)$ satisfy the orthogonality conditions

[TABLE]

for all $t\geq 0$. In contrast with the multi-level DO modes, here the one-dimensional BO modes carry the eigenvalues. By extending the derivation given in [14, 15] to the case of multi-level binary trees, it is straightforward to obtain the following BO-TT evolution equations

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

for $j=2,\ldots,d-2$. The matrices $\bm{S}_{i_{1}\cdots i_{j}}$ are skew-symmetric, as it can be verified by differentiating (82) with respect to time.

3.4 Equivalence between DO-TT and BO-TT expansions

The DO-TT and the BO-TT series expansions we discussed in Section 3.1 and Section 3.3 are equivalent in the sense that they approximate a time-dependent multivariate function using components from the same finite-dimensional function space. In this Section we prove such equivalence. The approach we take is to prove the equivalence for one level of the TT binary tree, following similar steps taken in [19], and then proceed inductively to show the equivalence throughout the whole TT binary tree. Clearly, switching from DO to BO at any level of the binary tree affects all expansions in the child nodes. Let $\bm{\varphi}^{(j)}(t)=[\varphi^{(j)}_{1}(t),\ldots,\varphi^{(j)}_{r_{j}}(t)]$ and $\bm{\varphi}^{(j+1,\ldots,d)}(t)=[\varphi^{(j+1,\ldots,d)}_{1}(t),\ldots,\varphi^{(j+1,\ldots,d)}_{r_{j}}(t)]$ (row vectors) be the BO modes at the $j$th-level of the TT binary tree. Consider the transformation

[TABLE]

where $[\bm{\Lambda}_{j}(t)]_{ik}=\langle\varphi^{(j)}_{i},\varphi^{(j)}_{k}\rangle$, and $\bm{P}_{j}(t)$ satisfies the matrix differential equation777$\bm{P}_{j}(t)$ is a time-dependent $r_{j}\times r_{j}$ matrix with real coefficients.

[TABLE]

with

[TABLE]

Lemma 3.1**.**

The matrix $\bm{P}_{j}(t)$ defined by the initial value problem (90) is orthogonal for all $t\geq 0$.

Proof.

The matrix $\bm{F}_{j}(t)=-\bm{\Lambda}_{j}^{-\frac{1}{2}}(t)\bm{\Sigma}_{j}(t)\bm{\Lambda}_{j}^{-\frac{1}{2}}(t)$ is skew-symmetric for all $t\geq 0$, since $\bm{\Sigma}_{j}$ is skew-symmetric (see Eq. (91)). Therefore, we have

[TABLE]

This implies that $\bm{P}_{j}(t)^{T}\bm{P}_{j}(t)=\bm{P}_{j}(0)^{T}\bm{P}_{j}(0)=\bm{I}$, i.e., $\bm{P}_{j}(t)$ is orthogonal for all $t\geq 0$. We notice that the same conclusion holds if we replace $\bm{P}_{j}(0)$ with any orthogonal matrix, not just the identity.

∎

Hereafter we prove that the linear transformation (89)-(90) defines the mapping between the DO and BO modes at the $j^{\text{th}}$-level of the TT binary tree.

Theorem 3.1**.**

The linear transformation defined by (89)-(90) is invertible and it defines a new set of modes $\{\psi^{(j)}_{1}(t),\ldots,\psi^{(j)}_{r_{j}}(t)\}$ such that, for all $t\geq 0$,

• (i)

$\{\psi^{(j)}_{1},\ldots,\psi^{(j)}_{r_{j}}\}$* is orthonormal,*

• (ii)

$\displaystyle\sum_{k=1}^{r_{j}}\psi^{(j)}_{k}(t)\psi^{(j+1,\ldots,d)}_{k}(t)=\sum_{k=1}^{r_{j}}\varphi^{(j)}_{k}(t)\varphi^{(j+1,\ldots,d)}_{k}(t)$,

• (iii)

$\{\psi^{(j)}_{1},\ldots,\psi^{(j)}_{r_{j}}\}$* satisfies the DO condition $\displaystyle\langle\frac{\partial\psi^{(j)}_{i}}{\partial t}\psi^{(j)}_{k}\rangle=0\qquad\forall i,k=1,\ldots,r_{j}$.*

Proof.

The transformations defined in (89) are invertible by Lemma 3.1. To prove $(i)$, we notice that

[TABLE]

Multiply by $\bm{P}_{j}^{T}\bm{\Lambda}_{j}^{-\frac{1}{2}}$ and $\bm{\Lambda}_{j}^{-\frac{1}{2}}\bm{P}_{j}$ to the left and the right hand sides, respectively. This yields,

[TABLE]

which proves $(i)$. To prove $(ii)$, it is sufficient to apply the transformation (89). In fact,

[TABLE]

To prove $(iii)$, we first differentiate $\bm{\varphi}^{(j)}=\bm{\psi}^{(j)}\bm{P}_{j}^{T}\bm{\Lambda}_{j}^{\frac{1}{2}}$ with respect to time to obtain888This equality follows from (90) and the identity

$\bm{S}_{j}=\bm{\Sigma}_{j}+\frac{1}{2}\frac{\partial\bm{\Lambda}_{j}}{\partial t}.$

[TABLE]

At this point, we have

[TABLE]

i.e.,

[TABLE]

because $\bm{\Sigma}_{j}$ is the skew-symmetric part of $\bm{S}_{j}$. We know that $\bm{P}_{j}$ and $\bm{\Lambda}_{j}$ are nonzero, and therefore it must be the case that

[TABLE]

i.e., the modes $\psi^{(j)}_{k}$ are dynamically orthogonal.

∎

With the equivalence between the BO and DO modes established at the $j^{\text{th}}$ level of the TT binary tree, we can discuss the effects of this transformation on the remaining parts of the tree. An immediate consequence is that every mode belonging to levels above the $j^{\text{th}}$ one remains unchanged. However, all modes below the $j^{\text{th}}$-level need to be recomputed. Moreover, in the case of the BO-TT representation, one has to make sure that there are no eigenvalue crossings in any of the bi-orthogonal decompositions [19].

4 Dynamically orthogonal tensor methods for high-dimensional nonlinear PDEs

In this Section we develop dynamically orthogonal tensor methods to compute to solution of initial/boundary value problems involving nonlinear PDEs of the form

[TABLE]

where $G$ is a nonlinear operator, $B$ is a linear boundary operator, $\Omega$ is a bounded subset of $\mathbb{R}^{d}$ which can be represented as a Cartesian products of $d$ one-dimensional domains999The numerical results we present in this paper are for spatial domains $\Omega$ that can be represented as Cartesian products of one-dimensional domains. Obviously, this is not always the case. Handling high-dimensional functions and PDEs in complex geometries is not a trivial task. For instance, a four-dimensional sphere may be discretized by a set of four-dimensional cubes, i.e., tesseracts. Each tesseract consists of eight cubical cells, 24 faces, 32 edges and 16 vertices. Connecting such tesseracts in a finite-element fashion is not straightforward. Similarly, mapping high-dimensional complex domains into separable domains (whenever possible) is not straightforward., while $u_{0}(\bm{x})$ and $h(\bm{x},t)$ are, respectively, the initial condition and the boundary condition. To compute the solution of (93) we substitute any of the tensor series expansion we discussed in Section 3, e.g., (55), into (93) and derive a coupled system of nonlinear evolution equations for the one-dimensional modes $\psi^{(j)}_{i_{1}\cdots i_{j}}(t)$. The derivation of such evolution equations is identical to the derivation given in Section 3, with $\partial u/\partial t$ replaced by $G(u)$. Specifically, the DO-TT system (68)-(69) takes the form

[TABLE]

where

[TABLE]

If $G(u)$ is separable, i.e., if it can be written as

[TABLE]

where $G_{j}^{(i)}$ ($i=1,\ldots,r_{G}$) are nonlinear operators acting only on functions of $x_{j}$, then after replacing $u$ with the series expansion (55) all inner products in (94)-(95) are essentially one-dimensional integrals. A simple example of a separable nonlinear operator is $G(u)=-u\cdot\nabla u+\nabla^{2}u$, where $\nabla$ and $\nabla^{2}$ are, respectively, $d$-dimensional gradient and Laplace operators. Other examples of separable linear operators will be given in Section 5. Projecting (93) recursively onto the tensor-train modes, e.g., as in Section 3.1, allows us to compute the solution on a tensor manifold with constant rank. This is achieved by solving a coupled system of one-dimensional nonlinear PDEs, e.g., the system (94). However, the solution of (93) may not have an accurate representation on a tensor manifold with constant rank for all times. Hence, we may need to add or remove modes adaptively as time integration proceeds.

Remark

Classical numerical tensor methods to solve high-dimensional PDEs with explicit time stepping schemes require rank-reduction to project the solution back into a tensor manifold [60] with specified rank (see [63] §5.5). This can be achieved, e.g., by a sequence of singular value decompositions [24, 40], or by optimization [39, 59, 9, 20, 56, 35]. Rank reduction can be computationally intensive, especially if performed at each time step. Numerical tensor methods with implicit time stepping suffer from similar issues. In particular, the nonlinear system that yields the solution at the next time step needs to be solved on a tensor manifold with constant rank, e.g., by using Riemannian optimization algorithms [59, 60, 22]. The dynamically orthogonal tensor method we propose operates in a different way. In particular, the hard-to-compute nonlinear projection [43, 38] that maps the solution of high-dimensional PDEs onto a tensor manifold with constant rank [60] here is represented explicitly by the hierarchical DO/BO propagator, i.e., by a system of coupled one-dimensional nonlinear PDEs. In other words, there is no need to perform tensor rank reduction [24, 26, 40], rank-constrained temporal integration [43, 38], or Riemannian optimization [59], when solving high-dimensional PDEs with the dynamically orthogonal tensor method we propose.

4.1 Adaptive addition and removal of modes

The solution to the PDE (93) may not be accurately approximated by elements of a tensor manifold with fixed rank at all times. Therefore, it may be desirable to increase or decrease the tensor rank of the solution as time integration proceeds. Removing modes is straightforward since one can simply truncate the BO-TT or DO-TT decomposition we discussed in Section 3.3 and Section 3.1, respectively, at any level of the binary tree to a decomposition of smaller rank. This obviously affects all child nodes and corresponding modes in the binary tree.

Adding modes is a more subtle task. If the hierarchical rank of the tensor representing the solution is too large at any level of the binary tree, then the matrices at the left hand side of DO-TT and BO-TT propagators (Eqs. (70) and (84)), i.e., $\bm{C}_{i_{1}\cdots i_{j}}(t)$ and $\bm{\Lambda}_{i_{1}\cdots i_{j}}(t)$, may become nearly singular. On the other hand, if the hierarchical rank is too small, the tensor approximation may not be accurate. The energy by which each mode contributes to the series expansion of the solution can be tracked by the eigenvalues of $\bm{\Lambda}_{i_{1}\cdots i_{j}}$ in the BO-TT setting, and by the eigenvalues of $\bm{C}_{i_{1}\cdots i_{j}}(t)$ in the DO-TT setting. Once a new mode is added with zero energy, one or more of the matrices $\bm{C}_{i_{1}\cdots i_{j}}$ or $\bm{\Lambda}_{i_{1}\cdots i_{j}}$ become singular. To overcome this issue, and be able to continue integrating the system Babaee et al. [4] proposed an energy threshold criterion combined with matrix pseudo-inverse. This technique is effective, although it does slightly pollute the solution at the time of adding the mode. The process of activating a new mode from a state with zero energy, can be rigorously addressed by using the theory of fast-slow systems [8]. In fact, both DO-TT and BO-TT propagators (Eqs. (70) and (84)) become very stiff systems when adding a new mode, which requires appropriate temporal-integrators. Hereafter we propose two new algorithms to add modes in the DO-TT/BO-TT propagators:

Algorithm 1

In this algorithm, we add a pair of modes (left and right) satisfying the orthogonality conditions of DO or BO (whichever condition is currently being enforced) with zero energy. Evolve the modes which do not require inverting the now singular matrix (right hand modes in the DO setting and left hand modes in the BO setting) for a short amount of time while keeping the other modes constant. At this point, the energy of the new mode comes up to some value $\lambda_{\epsilon}$ which is significantly smaller than the energy of the more developed modes. Now when the matrix is inverted to continue mode propagation of the modes which were fixed, we obtain a slow-fast system [8]. The evolution of these modes remains slow-fast until the energy amongst all modes become more balanced.

Algorithm 2

In this algorithm, we switch from DO-TT or BO-TT propagators to an explicit time stepping scheme involving numerical tensors [63, 37, 25] for a number of time steps. This naturally increases the rank of the solution tensor (see [63], §5.5) based on the structure of the PDE. At this point, we perform an orthogonalization of such tensor [24, 40], and restart the DO-TT/BO-TT propagators using such new orthogonal set of expansion modes. In other words, given a DO-TT (or BO-TT) representation of the solution of (93) at time $t$, we convert such representation into a numerical tensor in any format [39, 5, 24]. Then we perform a number of time steps of (93) with an explicit temporal integration scheme using the numerical tensor representation, without performing any rank reduction. This naturally sends the solution tensor into a tensor manifold with larger multivariate rank. At this point we orthogonalize the solution tensor to obtain a new set of bi-orthogonal modes, which can be further truncated using the thresholding technique explained in Section 2.4, and then restart the DO-TT (or BO-TT) propagator with the new initial condition.

Remark

Algorithm 1 relies on enlarging the finite dimensional function space which the approximate solution lives in. Algorithm 2 is different in that the new bi-orthogonal modes obtained from the numerical tensor need not lie in the same finite dimensional function space as the previous set of modes. In other words, we are re-representing the solution in a different finite dimensional function space (tensor manifold with larger multivariate rank). Moreover, we are not explicitly adding modes with low energy, but we are letting the PDE itself increase the solution rank using numerical tensor techniques. Algorithm 2 may or may not result in a slow-fast system. Numerical examples demonstrating the effectiveness of Algorithm 2 will be provided in Section 5.1.

5 Numerical examples

In this Section, we demonstrate the accuracy and computational efficiency of the dynamically orthogonal tensor method we propose in this paper to compute the solution of high-dimensional PDEs. Specifically, we study the DO-TT representation, and apply it to hyperbolic and parabolic PDEs in periodic hypercubes with dimension ranging from 2 to 50. The reason for the choice of such PDEs is that they admit analytical solutions, which we will use to rigorously assess the accuracy and convergence rate of the proposed methods. We also demonstrate the adpative algorithm we developed in Section 4.1 to dynamically enrich the time-evolving basis at each level of the TT binary tree.

5.1 Hyperbolic PDEs

In this Section we study the DO-TT propagator of several linear hyperbolic PDEs in periodic domains.

5.1.1 Two-dimensional hyperbolic PDE

Let us begin with the two-dimensional initial/boundary value problem

[TABLE]

in the spatial domain $\Omega=[0,2\pi]^{2}$ with periodic boundary conditions.

As is well known, the PDE (97) can be reduced to the trivial ODE $du/dt=0$ along the flow generated by the dynamical system (see, e.g., [54])

[TABLE]

In Figure 7 we plot the phase portrait we obtained by solving (98) numerically for different initial conditions $(x_{01},x_{02})\in\Omega$. With the flow $\{x_{1}(t,x_{01},x_{02}),x_{1}(t,x_{01},x_{02})\}$ available, we can write the analytical solution to (97) as

[TABLE]

where $\{x_{01}(x_{1},x_{2},t),x_{02}(x_{1},x_{2},t)\}$ denotes the inverse flow. The semi-analytical solution (99) is plotted in the first row of Figure 8 at different times. Next, we compute the solution to (97) using the DO-TT method we proposed in Section 4. To this end, we first perform a bi-orthogonal decomposition of the initial condition $u(x_{1},x_{2},0)$ in (97) with threshold set to $\sigma=10^{-13}$ (see Section 2.4). This yields $17$ modes $(\psi^{(1)}_{i_{1}},\psi^{(2)}_{i_{1}})_{i_{1}=1}^{17}$, each of which is collocated on an evenly-spaced grid with $257$ nodes in $[0,2\pi]$. The DO-TT system associated with the PDE (97) is

[TABLE]

In a Fourier spectral collocation setting [29], the partial derivatives $\partial/\partial x_{1}$, $\partial/\partial x_{2}$ and the inner products can be easily represented by one-dimensional spectral differentiation matrices, and one-dimensional Fourier quadrature rules [29]. This allows us to transforms the PDE system (100) into a system of nonlinear ODEs with $2\times r_{1}\times 256$ equations. It is worthwhile emphasizing that the DO-TT system (100) is nonlinear even though the PDE (97) is linear. The nonlinearity is related to the fact that we implicitly project the solution back on a tensor manifold with constant rank $r_{1}$ at each time. As we pointed out in Section 4, projecting on a tensor manifold with constant rank is a nonlinear operation. The system (100) is solved numerically by inverting the matrix $C_{jk}(t)=\langle\psi^{(2)}_{j}\psi^{(2)}_{k}\rangle$ (assuming it is non-singular) at each time step, and an explicit RK4 scheme with time step $\Delta t=10^{-3}$. We ran one simulation with constant rank $r_{1}=17$, two adaptive simulations using the pseudo-inverse (PI) technique proposed in [4] for addition of modes, and one adaptive simulation using our Algorithm 2 in Section 4.1. In Figure 8 we compare the time snapshots of the constant rank DO-TT solution ($r_{1}=17$) with the semi-analytical solution (99).

In Figure 12 we plot the time-dependent $L^{2}(\Omega)$ errors between various DO-TT simulations and the semi-analytical solution (99). It is seen that the simulation with constant rank has an error slope that increases substantially around $t=0.5$. This suggests that $17$ DO-TT modes are no longer sufficient to represent the solution (99) for $t>0.5$. This can also be also seen from the fact that the DO-TT spectrum tends to flatten out in time (Figure 9), suggesting that each of the $17$ modes is picking up more and more energy. To overcome this problem, and therefore control the error growth in time, we implemented the adaptive algorithm for mode addition/removal proposed in [4], and compared it with our Algorithm 2 in Section 4.1 (Algorithm 2). As it is seen in Figure 11, such algorithm can indeed control the temporal growth of the DO-TT spectrum, hence the overall error (see Figure 12). Each time modes are added using our Algorithm 2 (Section 4.1), there is a re-orthogonalization process which can yield a discontinuity in the temporal evolution of each mode (see Figure 11). In practice, the DO-TT system is re-started from a new initial condition after such re-orthogonalization. This does is not affect the solution, nor creates any temporal discontinuity or error jump (see Figure 12).

5.1.2 Four-dimensional hyperbolic PDE

Let us consider the following four-dimensional linear hyperbolic PDE

[TABLE]

in the spatial domain $\Omega=[0,2\pi]^{4}$, with periodic boundary conditions. In equation (101) $c_{ij}$ are real numbers and $f_{j}(x_{j})$ are real-valued functions. The system of DO-TT evolution equations (94) can be explicitly written for the linear PDE (101). Such system is rather complicated, and therefore not presented here. For numerical demonstration, we set the coefficient matrix $c_{ij}$ as

[TABLE]

and consider the following functions

[TABLE]

This yields non-trivial dynamics ($c_{ij}$ is not diagonal). We solve the DO-TT system (94) numerically using a Fourier spectral collocation method with $20$ Fourier points in each variable $x_{i}$, and RK4 time integration with $\Delta t=10^{-3}$. To this end, we first decompose the four-dimensional initial condition in (101) with the bi-orthogonal method we discussed in Section 2.1. Specifically, we consider the space $H^{0}(\Omega)=L_{2}(\Omega)$ and set the eigenvalue threshold to $\sigma=10^{-10}$. This yields the following hierarchical ranks

[TABLE]

which we keep constant throughout the simulation. In other words, here we do not perform adaptive addition/removal of DO-TT modes as we did in the previous PDE example (Eqs. (97) and (100)). In Figure 13 we plot the time evolution of a few representative DO-TT modes.

It is seen that as time evolves the spatial frequency of such modes increases which suggests that the hyperbolic dynamics activate higher spatial frequencies in the spectral representation of the PDE solution. This is shown in Figure 14, where we plot one section of the solution to the initial value problem (101) we obtained with the method of characteristics (benchmark solution), and the solution we obtained with the DO-TT propagator. It is seen that even with the very few hierarchical ranks summarized in (102) we were able to resolve the 4D solution to a reasonable accuracy (see Figure 15(a)).

5.2 Fifty-dimensional hyperbolic PDE

Let us consider the following fifty-dimensional hyperbolic problem

[TABLE]

subject to periodic boundary conditions in the hyper-cube $\Omega=[0,2\pi]^{50}$. It is easy to show that the solution to (103) is rank-one for all $t\geq 0$. Correspondingly, we look for a rank-one DO-TT solution of the form

[TABLE]

By substituting (104) into (103) and imposing DO orthogonality conditions at each level of the binary tree yields the DO-TT propagator

[TABLE]

Note that this system is nonlinear. Specifically, the first 49 equations are uncoupled, while the 50th one is coupled with the entire system. We set the coefficients $f_{j}(x_{j})=j$, i.e., constant, and the components of the initial condition as

[TABLE]

which satisfy the bi-orthogonality condition. The analytical solution to (103) is easily obtained as

[TABLE]

As before, we solved the DO-TT system (105) using a Fourier spectral collocation method with $60$ Fourier points in each variable $x_{j}$, and RK4 time integration with $\Delta t=10^{-3}$. The temporal evolution of a few representative DO-TT modes is shown in Figure 16.

The $L^{2}(\Omega)$ error between the analytical solution (107) and the DO-TT approximation (104) can be represented in terms of one-dimensional integrals as follows

[TABLE]

This error is plotted in Figure 15(b) versus time. It is seen that the DO-TT propagator in this case is numerically exact.

5.3 Parabolic PDEs

In this Section we study DO-TT solution of constant-coefficients parabolic PDEs in periodic hypercubes. In contrast with the hyperbolic problems discussed in the previous Section, the solution in this case does not increase in rank throughout propagation. Rather, the solution smooths out over time and the rank tends to decrease. This implies that in order for the covariance matrices in the DO-TT system to remain invertible during propagation, a threshold must be set to dynamically remove modes, in particular when their energy becomes very small.

5.3.1 Four-dimensional parabolic PDE

Let us consider the four-dimensional initial value problem

[TABLE]

in the spatial domain $\Omega=[0,2\pi]^{4}$, with periodic boundary conditions. Note that here we set the same initial condition as in the hyperbolic PDE (101). In this case, however, the solution decays to zero because of diffusion. Hence, we expect that the DO-TT solution rank decays in time, instead of increasing. Applying the Fourier transform $\mathcal{F}[\cdot]$ to (109), yields the linear ODE

[TABLE]

which can be solved analytically or numerically to obtain a benchmark solution to (109). With such benchmark solution available, we can study the accuracy of the DO-TT propagator. To this end, we first compute the bi-orthogonal decomposition of the initial condition in (109) as we have done in Section 5.1.2, with threshold set to $\sigma=10^{-10}$. This yields the the hierarchical ranks (102), as before. In Figure 17 we plot $L^{2}(\Omega)$ error of the DO-TT approximation of the solution to (109) relative to the benchmark solution obtained via Fourier transform. It is seen that, contrary to hyperbolic problems, diffusion promotes a low-rank structure of the solution.

This means that the energy of the DO-TT modes tends to decrease in time, which implies that the covariance matrix at the left hand side of the DO-TT propagator can become numerically singular in time. To avoid such singularities we introduce a threshold $\epsilon=10^{-10}$ on the eigenvalues of the level-1 TT binary tree, i.e., $\lambda_{i_{1}}$, to determine whether we should decrease the leading rank $r_{1}$. Specifically, if $\lambda_{r_{1}}<\epsilon$ then $r_{1}=r_{1}-1$. For this example, thresholding only the first level is sufficient to avoid singular covariance matrices at all levels of the TT binary tree. In Figure 17 we see that each time we remove one mode we have that the relative error jumps. Note however, that it immediately decays again with a rate that is inversely proportional to the number of the DO-TT modes kept in the series expansion.

5.4 Fifty-dimensional parabolic PDE

The last problem we consider is a fifty-dimensional constant-coefficient diffusion problem with rank-one initial condition

[TABLE]

in the spatial domain $\Omega=[0,2\pi]^{5}0$, with periodic boundary conditions. The analytical solution can be obtained by using the method of separation of variables [49] as

[TABLE]

Hence, it is a rank-one solution. Correspondingly, we look for a rank-one DO-TT solution of the form

[TABLE]

where the modes $\psi^{(j)}(t)$ satisfy the DO-TT system of evolution equations

[TABLE]

We set the initial condition in (111) as

[TABLE]

which satisfy the orthogonality conditions. As before, we solved the system (114) using Fourier spectral collocation with 60 nodes in each variable $x_{j}$ and RK4 time integration with $\Delta t=10^{-3}$. In Figure 18 we plot the temporal evolution of a few representative DO-TT modes. As expected, they all decay to zero.

Finally in Figure 19 we plot the time-dependent $L^{2}(\Omega)$ error between the DO-TT solution (113) and the analytical solution (112). Such error can be expressed in terms of one-dimensional integrals as follows

[TABLE]

6 Summary

In this paper, we presented a new dynamically orthogonal tensor decomposition method to approximate multivariate functions and the solution to high-dimensional time-dependent nonlinear PDEs. The key idea relies on a hierarchical decomposition of the approximation space obtained by splitting the independent variables of the problem into disjoint subsets. This process, which can be conveniently be visualized in terms of binary trees, yields series expansions analogous to the classical tensor-train and hierarchical Tucker tensor formats. By enforcing dynamic orthogonality conditions at each level of binary tree, we obtained evolution equations for the orthogonal modes spanning each of the nested subspaces in the hierarchical decomposition. This allowed us to represent the temporal dynamics of high-dimensional functions, and compute the solution to high-dimensional time-dependent PDEs on a tensor manifold with constant rank. We also proposed a new algorithm for dynamic addition and removal of modes within each subspace, and demonstrated its effectiveness in numerical applications to hyperbolic and parabolic PDEs. The mathematical techniques and algorithms we presented in this paper can be readily applied to more general high-dimensional nonlinear systems and PDEs, such as the finite-dimensional approximation of nonlinear functionals and functional differential equations [63].

Acknowledgements This research was supported by the U.S. Army Research Office grant W911NF1810309.