Tensor-Train Format Solution with Preconditioned Iterative Method for High Dimensional Time-Dependent Space-Fractional Diffusion Equations with Error Analysis

Abstract

In this paper, a first order implicit finite difference scheme with Krylov subspace linear system solver is employed to solving time-dependent space-fractional diffusion equations in high dimensions where the initial condition and source term are in tensor-train (TT) format with low TT-ranks. In the time-marching process, TT-format of the solution is maintained and the increment of TT-ranks due to addition is moderated by rounding. The error introduced by rounding is shown to be consistent with the first order finite difference scheme. On the other hand, the linear systems involved in the solution process are shown to possess Toeplitz-like structure so that the complexity and required memory for Krylov subspace solver can be optimized. Further reduction in complexity is made by utilizing a circulant preconditioner which accelerates the convergence rate of Krylov subspace method dramatically. Numerical examples for problems up to 20 dimensions are presented.

Appendices

Examples of FTT-format

The $x^{(k)}$ will denote the variable in the kth direction, $f^{[k]}$ will denote a function of variable $x^{(k)}$, $x_{i_{k}}^{(k)}$ denotes the $i_{k}^{th}$ grid point in the kth direction, and $f_{i_{k}}^{[k]}:=f^{[k]}(x_{i_{k}}^{(k)})$.

• 1. Product-separable function: $f(x^{(1)},\ldots ,x^{(d)})=f^{[1]}(x^{(1)})\cdots f^{[d]}(x^{(d)})$.

  FTT-format: $f=\begin{bmatrix}f^{[1]}(x^{(1)})\end{bmatrix}\begin{bmatrix}f^{[2]}(x^{(2)})\end{bmatrix}\cdots \begin{bmatrix}f^{[d-1]}(x^{(d-1)})\end{bmatrix}\begin{bmatrix}f^{[d]}(x^{(d)})\end{bmatrix}$.

  Tensor form: ${\mathscr {F}}(i_{1},\ldots ,i_{d})=f_{i_{1}}^{[1]}\cdots f_{i_{d}}^{[d]}$.

  TT-format: ${\mathscr {F}}(i_{1},\ldots ,i_{d})=\begin{bmatrix}f_{i_{1}}^{[1]}\end{bmatrix}\begin{bmatrix}f_{i_{2}}^{[2]}\end{bmatrix}\cdots \begin{bmatrix}f_{i_{d-1}}^{[d-1]}\end{bmatrix}\begin{bmatrix}f_{i_{d}}^{[d]}\end{bmatrix}$.

  TT-ranks: At most 1.

• 2. Sum-separable function Oseledets [45]: $f(x^{(1)},\ldots ,x^{(d)})=f^{[1]}(x^{(1)})+\cdots +f^{[d]}(x^{(d)})$.

  FTT-format: $f=\!\begin{bmatrix}f^{[1]}(x^{(1)})&1\end{bmatrix}\!\begin{bmatrix}1&0\\ f^{[2]}(x^{(2)})&1\end{bmatrix}\!\cdots \!\begin{bmatrix}1&0\\ f^{[d-1]}(x^{(d-1)})&1\end{bmatrix}\!\begin{bmatrix}1\\ f^{[d]}(x^{(d)})\end{bmatrix}$.

  Tensor form: ${\mathscr {F}}(i_{1},\ldots ,i_{d})=f_{i_{1}}^{[1]}+\cdots +f_{i_{d}}^{[d]}$.

  TT-format: ${\mathscr {F}}(i_{1},\ldots ,i_{d})=\begin{bmatrix}f_{i_{1}}^{[1]}&1\end{bmatrix}\begin{bmatrix}1&0\\ f_{i_{2}}^{[2]}&1\end{bmatrix}\cdots \begin{bmatrix}1&0\\ f_{i_{d-1}}^{[d-1]}&1\end{bmatrix}\begin{bmatrix}1\\ f_{i_{d}}^{[d]}\end{bmatrix}$.

  TT-ranks: At most 2.

• 3. Sine of a sum Oseledets [45]: $f(x^{(1)},\ldots ,x^{(d)})=\sin (x^{(1)}+\cdots +x^{(d)})$.

  FTT-format: $f= \begin{bmatrix}\sin x^{(1)}&\cos x^{(1)}\end{bmatrix}\begin{bmatrix}\cos x^{(2)}&-\,\sin x^{(2)}\\ \sin x^{(2)}&\cos x^{(2)}\end{bmatrix}\cdots \begin{bmatrix}\cos x^{(d-1)}&-\,\sin x^{(d-1)}\\ \sin x^{(d-1)}&\cos x^{(d-1)}\end{bmatrix}\begin{bmatrix}\cos x^{(d)}\\ \sin x^{(d)}\end{bmatrix}$.

  Tensor form: ${\mathscr {F}}(i_{1},\ldots ,i_{d})=\sin (x_{i_{1}}^{(1)}+\cdots +x_{i_{d}}^{(d)})$.

  TT-format: ${\mathscr {F}}(i_{1},\ldots ,i_{d})= \begin{bmatrix}\sin x_{i_{1}}^{(1)}&\cos x_{i_{1}}^{(1)}\end{bmatrix}\begin{bmatrix}\cos x_{i_{2}}^{(2)}&-\,\sin x_{i_{2}}^{(2)}\\ \sin x_{i_{2}}^{(2)}&\cos x_{i_{2}}^{(2)}\end{bmatrix}\cdots \begin{bmatrix}\cos x_{i_{d-1}}^{(d-1)}&-\,\sin x_{i_{d-1}}^{(d-1)}\\ \sin x_{i_{d-1}}^{(d-1)}&\cos x_{i_{d-1}}^{(d-1)}\end{bmatrix}\begin{bmatrix}\cos x_{i_{d}}^{(d)}\\ \sin x_{i_{d}}^{(d)}\end{bmatrix}$.

  TT-ranks: At most 2.

• 4. “Laplace-like” function Oseledets [44]:

  $f(x^{(1)},\ldots ,x^{(d)})=g^{[1]}(x^{(1)})f^{[2]}(x^{(2)})\cdots f^{[d]}(x^{(d)})+f^{[1]}(x^{(1)})g^{[2]}(x^{(2)})\cdots f^{[d]}(x^{(d)}) +\cdots +f^{[1]}(x^{(1)})f^{[2]}(x^{(2)})\cdots g^{[d]}(x^{(d)}))$.

  FTT-format: $f=\begin{bmatrix} g^{[1]}(x^{(1)})&f^{[1]}(x^{(1)}) \end{bmatrix} \begin{bmatrix} f^{[2]}(x^{(2)})&0\\ g^{[2]}(x^{(2)})&f^{[2]}(x^{(2)}) \end{bmatrix}$

  $\qquad \qquad \qquad \qquad \quad \cdots \begin{bmatrix}f^{[d-1]}(x^{(d-1)})&0\\ g^{[d-1]}(x^{(d-1)})&f^{[d-1]}(x^{(d-1)})\end{bmatrix} \begin{bmatrix}f^{[d]}(x^{(d)})\\ g^{[d]}(x^{(d)})\end{bmatrix}$.

  Tensor form:

  ${\mathscr {F}}(i_{1},\ldots ,i_{d})=g^{[1]}(x_{i_{1}}^{(1)})f^{[2]}(x_{i_{2}}^{(2)})\cdots f^{[d]}(x_{i_{d}}^{(d)})+f^{[1]}(x_{i_{1}}^{(1)})g^{[2]}(x_{i_{2}}^{(2)})\cdots f^{[d]}(x_{i_{d}}^{(d)})$

  $\qquad \qquad \qquad \qquad \quad +\cdots +f^{[1]}(x_{i_{1}}^{(1)})f^{[2]}(x_{i_{2}}^{(2)})\cdots g^{[d]}(x_{i_{d}}^{(d)}))$.

  TT-format:

  ${\mathscr {F}}(i_{1},\ldots ,i_{d})=\begin{bmatrix}g^{[1]}(x_{i_{1}}^{(1)})&f^{[1]}(x_{i_{1}}^{(1)})\end{bmatrix}\begin{bmatrix}f^{[2]}(x_{i_{2}}^{(2)})&0\\ g^{[2]}(x_{i_{2}}^{(2)})&f^{[2]}(x_{i_{2}}^{(2)})\end{bmatrix}$

  $\qquad \qquad \qquad \qquad \quad \cdots \begin{bmatrix}f^{[d-1]}(x_{i_{d-1}}^{(d-1)})&0\\ g^{[d-1]}(x_{i_{d-1}}^{(d-1)})&f^{[d-1]}(x_{i_{d-1}}^{(d-1)})\end{bmatrix} \begin{bmatrix}f^{[d]}(x_{i_{d}}^{(d)})\\ g^{[d]}(x_{i_{d}}^{(d)})\end{bmatrix}$.

  TT-ranks: At most 2.

For general functions, a detailed study on continuous analogue of TT-decomposition can be found in Gorodetsky et al. [19].

Examples of TT-Format Operations

For tensors ${\mathscr {A}}$ and ${\mathscr {B}}$ of the same size and with TT-formats

$$\begin{aligned} {\mathscr {A}}(i_{1},\ldots ,i_{d})=A_{1}(i_{1})\cdots A_{d}(i_{d}),\quad {\mathscr {B}}(i_{1},\ldots ,i_{d})=B_{1}(i_{1})\cdots B_{d}(i_{d}), \end{aligned}$$

• 1. Their sum ${\mathscr {C}}={\mathscr {A}}+{\mathscr {B}}$ has TT-format

  ${\mathscr {C}}(i_{1},\ldots ,i_{d})=A_{1}(i_{1})\cdots A_{d}(i_{d})+B_{1}(i_{1})\cdots B_{d}(i_{d})$

  $\quad =\begin{bmatrix}A_{1}(i_{1})&B_{1}(i_{1}) \end{bmatrix}\!\begin{bmatrix}A_{2}(i_{2})&0\\ 0&B_{2}(i_{2})\end{bmatrix}\!\cdots \!\begin{bmatrix}A_{d-1}(i_{d-1})&0\\ 0&B_{d-1}(i_{d-1})\end{bmatrix}\!\begin{bmatrix}B_{d}(i_{d})\\ A_{d}(i_{d})\end{bmatrix}$,

  with TT-ranks being added.

• 2. Their Hadamard product (entrywise product) has TT-format

  ${\mathscr {C}}(i_{1},\ldots ,i_{d})=A_{1}(i_{1})\cdots A_{d}(i_{d})B_{1}(i_{1})\cdots B_{d}(i_{d})=\begin{bmatrix}A_{1}(i_{1})\otimes B_{1}(i_{1})\end{bmatrix}\cdots \begin{bmatrix}A_{d}(i_{d})\otimes B_{d}(i_{d})\end{bmatrix}$,

  with TT-ranks being multiplied.

• 3. Their inner product (sum of entrywise products) is obtained by

  $\langle {\mathscr {A}},{\mathscr {B}}\rangle =\sum \limits _{i_{1},\ldots ,i_{d}}A_{1}(i_{1})\cdots A_{d}(i_{d})B_{1}(i_{1})\cdots B_{d}(i_{d})=\left( \sum \limits _{i_{1}}A_{1}(i_{1})\otimes B_{1}(i_{1})\right) \cdots \left( \sum \limits _{i_{d}}A_{d}(i_{d})\otimes B_{d}(i_{d})\right)$.

  with TT-ranks being multiplied. In practice, if the TT-cores are of $r\times N\times r$, the operation cost can be made ${\mathscr {O}}(dNr^{3})$ by the tensor structure of $(A\otimes B){\mathbf {x}}$ Oseledets [44].

Matrix-Core Multiplication

Suppose G is a TT-core of $r_{0}\times N\times r_{1}$, then the fibers of G are vectors as denoted in the following.

$$\begin{aligned} \text {mode-1 fibers: }G(\cdot ,\beta ,\gamma )&:=[G(1,\beta ,\gamma ),G(2,\beta ,\gamma ),\cdots ,G(r_{0}, \beta ,\gamma )]^{\intercal },\\&1\le \beta \le N, ~ 1\le \gamma \le r_{1}; \\ \text {mode-2 fibers: }G(\alpha ,\cdot ,\gamma )&:=[G(\alpha ,1,\gamma ),G(\alpha ,2,\gamma ), \cdots ,G(\alpha ,N,\gamma )]^{\intercal },\\&1\le \alpha \le r_{0}, ~ 1\le \gamma \le r_{1}; \\ \text {mode-3 fibers: }G(\alpha ,\beta ,\cdot )&:=[G(\alpha ,\beta ,1),G(\alpha ,\beta ,2),\cdots , G(\alpha ,\beta ,r_{1})]^{\intercal },\\&1\le \alpha \le r_{0}, ~ 1\le \beta \le N. \end{aligned}$$

We may also put the fibers into the matrices

$$\begin{aligned} G[1]&:=\left[ G(\cdot ,1,1)\big |\cdots \big |G(\cdot ,N,1)\big |\cdots \big | G(\cdot ,1,r_{1})\big |\cdots \big | G(\cdot ,N,r_{1})\right] \in M_{r_{0} \times Nr_{1}}; \\ G[2]&:=\left[ G(1,\cdot ,1)\big |\cdots \big | G(r_{0},\cdot ,1)\big |\cdots \big | G(1,\cdot ,r_{1})\big |\cdots \big | G(r_{0},\cdot ,r_{1})\right] \in M_{N\times r_{0}r_{1}}; \\ G[3]&:=\left[ G(1,1,\cdot )\big |\cdots \big | G(r_{0},1,\cdot )\big |\cdots \big | G(1,N,\cdot )\big |\cdots \big | G(r_{0},N,\cdot )\right] \in M_{r_{1}\times r_{0}N}. \end{aligned}$$

Then, given matrices A of $r'_{0}\times r_{0}$, B of $N'\times N$, C of $r'_{1}\times r_{1}$, we have the matrix-core multiplications

$$\begin{aligned} (G\times _{1}A)[1]&=AG[1]\in M_{r'_{0}\times Nr_{1}},\\ (G\times _{2}B)[2]&=BG[2]\in M_{N'\times r_{0}r_{1}},\\ (G\times _{3}C)[3]&=CG[3]\in M_{r'_{1}\times r_{0}N}, \end{aligned}$$

such that $G\times _{1}A$ is a TT-core of $r'_{0}\times N\times r_{1}$ and has all $Nr_{1}$ columns of AG[1] as mode-1 fibers; $G\times _{2}B$ is a TT-core of $r_{0}\times N'\times r_{1}$ has all $r_{0}r_{1}$ columns of BG[2] as mode-2 fibers; $G\times _{3}B$ is a TT-core of $r_{0}\times N\times r'_{1}$ and has all $r_{0}N$ columns of CG[3] as mode-3 fibers.

Rounding Algorithm

We first denote some notations.

• 1. For a TT-core G which contains matrices $G(1),G(2),\ldots ,G(N)$, denote the matrices

  $$\begin{aligned} \text {ro}(G)=[G(1)|G(2)|\cdots |G(N)]; \quad \text {co}(G)=[G(1)^{\intercal }|G(2)^{\intercal }|\cdots |G(N)^{\intercal }]^{\intercal }. \end{aligned}$$

• 2. For a matrix $A=[A_{1}|A_{2}|\cdots |A_{N}]$, where $A_{k}$ have the same size,

  and a matrix $B=[B_{1}^{\intercal }|B_{2}^{\intercal }|\cdots |B_{N}^{\intercal }]^{\intercal }$ where $B_{k}$ have the same size,

  denote

  $$\begin{aligned} \text {iro}(A,N):=G ~\text { and } ~\text {ico}(B,N):=H \end{aligned}$$

  as the TT-cores containing matrices $G(k)=A_{k}$ and $H(k)=B_{k}$ for $1\le k\le N.$

• 3. For a matrix A, denote

  $$\begin{aligned} \text {lq}(A)=\{L,Q\} \end{aligned}$$

  as a set of two matrices if $A=LQ$ where L is lower triangular and Q has orthonormal rows. Indeed, with QR decomposition $A^{\intercal }={\tilde{Q}}R$, then $L=R^{\intercal }$ and $Q={\tilde{Q}}^{\intercal }$.

• 4. For a matrix A with singular value decomposition (SVD)

  $$\begin{aligned} A=U\varSigma V^{\intercal }=\sigma _{1}{\mathbf {u}}_{1}{\mathbf {v}}_{1}^{\intercal }+\sigma _{2}{\mathbf {u}}_{2}{\mathbf {v}}_{2}^{\intercal } +\cdots +\sigma _{r}{\mathbf {u}}_{r}{\mathbf {v}}_{r}^{\intercal }, \end{aligned}$$

  given a $\delta >0$, denote the $\delta$-truncated SVD as

  $$\begin{aligned} \sigma _{1}{\mathbf {u}}_{1}{\mathbf {v}}_{1}^{\intercal }+\sigma _{2}{\mathbf {u}}_{2} {\mathbf {v}}_{2}^{\intercal }+\cdots +\sigma _{r'}{\mathbf {u}}_{r'} {\mathbf {v}}_{r'}^{\intercal }:=U'\varSigma 'V'^{\intercal }, \end{aligned}$$

  where

  $$\begin{aligned} r'=\min \left\{ k:\sqrt{\sigma _{k+1}^{2}+\cdots +\sigma _{r}^{2}}\le \delta \right\} . \end{aligned}$$

  Then, denote

  $$\begin{aligned} \text {svd}_{\delta }(A)=\{U',\varSigma 'V'^{\intercal }\} \end{aligned}$$

  as a set of two matrices.

With these notations, the rounding algorithm is presented as follows.