Tensor train-Karhunen-Lo\`eve expansion for continuous-indexed random fields using higher-order cumulant functions

TL;DR

This paper introduces a novel tensor train-based framework for representing and computing with multidimensional non-Gaussian random fields, improving efficiency and adaptivity over traditional methods.

Contribution

It develops a new theoretical and algorithmic approach using tensor train decomposition of cumulant functions, eliminating the need for basis or quadrature point selection.

Findings

Accurate tensor train decomposition of cumulant functions achieved

Algorithm outperforms traditional methods in efficiency and adaptivity

Effective for non-Gaussian, non-homogeneous random fields

Abstract

The goals of this work are two-fold: firstly, to propose a new theoretical framework for representing random fields on a large class of multidimensional geometrical domain in the tensor train format; secondly, to develop a new algorithm framework for accurately computing the modes and the second and third-order cumulant tensors within moderate time. The core of the new theoretical framework is the tensor train decomposition of cumulant functions. This decomposition is accurately computed with a novel rank-revealing algorithm. Compared with existing Galerkin-type and collocation-type methods, the proposed computational procedure totally removes the need of selecting the basis functions or collocation points and the quadrature points, which not only greatly enhances adaptivity, but also avoids solving large-scale eigenvalue problems. Moreover, by computing with third-order cumulant functions, the new theoretical and algorithm frameworks show great potential for representing general non-Gaussian non-homogeneous random fields. Three numerical examples, including a three-dimensional random field discretization problem, illustrate the efficiency and accuracy of the proposed algorithm framework.

Tables

Table 1. Table 1: Time cost (seconds) of the quantities in Example 1

	${\tilde{C}}_{k,\mathrm{TT}}$	${\epsilon }_{gk}$	modes	${\epsilon }_{\mathrm{gf}k}$	total
$k=2$	107	1	440	1	549
$k=3$	4120	2	1467	1	5590

Table 2. Table 2: Control points and weights defining the geometry in Example 2

Subscripts $𝒊$		Control points $B_{𝒊}$		weights $w_{𝒊}$
(1,1)	(1,2)	(-0.5,-0.5,0)	(0.5,-0.5,1)	1	1
(2,1)	(2,2)	(-0.5,0.5,1)	(0.5,0.5,0)	1	1

Table 3. Table 3: Time cost (seconds) of the quantities in Example 2

	${\tilde{C}}_{k,\mathrm{TT}}$	${\epsilon }_{gk}$	modes	${\epsilon }_{\mathrm{gf}k}$	total
$k=2$	31	1	293	1	326
$k=3$	575	3	3962	1	4540

Table 4. Table 4: Control points and weights defining the geometry in Example 3

Subscripts $𝒊$		Control points $B_{𝒊}$		weights $w_{𝒊}$
(1,1,1)	(2,1,1)	$(R_i,0,0)$	$(R_o,0,0)$	1	1
(1,2,1)	(2,2,1)	$(R_i,R_i,0)$	$(R_o,R_o,0)$	$\frac{1}{\sqrt{2}}$	$\frac{1}{\sqrt{2}}$
(1,3,1)	(2,3,1)	$(0,R_i,0)$	$(0,R_o,0)$	1	1
(1,1,2)	(2,1,2)	$(R_i,0,R_i)$	$(R_o,0,R_o)$	$\frac{1}{\sqrt{2}}$	$\frac{1}{\sqrt{2}}$
(1,2,2)	(2,2,2)	$(R_i,R_i,R_i)$	$(R_o,R_o,R_o)$	$\frac{1}{2}$	$\frac{1}{2}$
(1,3,2)	(2,3,2)	$(0,R_i,R_i)$	$(0,R_o,R_o)$	$\frac{1}{\sqrt{2}}$	$\frac{1}{\sqrt{2}}$
(1,1,3)	(2,1,3)	$(0,0,R_i)$	$(0,0,R_o)$	1	1
(1,2,3)	(2,2,3)	$(0,0,R_i)$	$(0,0,R_o)$	$\frac{1}{\sqrt{2}}$	$\frac{1}{\sqrt{2}}$
(1,3,3)	(2,3,3)	$(0,0,R_i)$	$(0,0,R_o)$	1	1

Table 5. Table 5: Time cost (seconds) of the quantities in Example 3

	${\tilde{C}}_{k,\mathrm{TT}}$	${\epsilon }_{gk}$	modes	${\epsilon }_{\mathrm{gf}k}$	total
$k=2$	202	1	2014	1	2218
$k=3$	2821	11	31017	1	33850

Equations

$$\mathbf{x}(\xi )$$

$$\mathbf{x}(\xi ,\eta )$$

$$\mathbf{x}(\xi ,\eta ,\zeta )$$

$$R_i^p(\xi )$$

$$R_{i,j}^{p,q}(\xi ,\eta )$$

$$R_{i,j,k}^{p,q,r}(\xi ,\eta ,\zeta )$$

$$N_i^0$$

$$N_i^p(\xi )$$

$$\mathbf{x}(\mathbf{\xi })=\sum _{\mathbf{I}}{R}_{\mathbf{I}}^{\mathbf{p}}(\mathbf{\xi }){\mathbf{B}}_{\mathbf{I}}$$

$${\tilde{C}}_f(\xi ,{\xi }^{\prime })={\int }_0^1\tilde{C}(\xi ,\eta ;{\xi }^{\prime },\eta )\mathrm{d}\eta$$

$$\alpha (\xi ,\eta ;\theta )\approx \mathbf{f}(\xi ){\mathbf{\alpha }}_{-\xi }(\eta ,\theta )$$

$${\mathbf{\alpha }}_{-\xi }(\eta ,\theta )={\int }_0^1\mathbf{f}(\xi {)}^{\mathrm{T}}\alpha (\xi ,\eta ;\theta )\mathrm{d}\xi$$

$${\tilde{C}}_{-\xi }(\eta ,{\eta }^{\prime }{)}_{n_1\times n_1}={\int }_0^1{\int }_0^1\mathbf{f}(\xi {)}^{\mathrm{T}}\tilde{C}(\xi ,\eta ;{\xi }^{\prime },{\eta }^{\prime })\mathbf{f}({\xi }^{\prime })\mathrm{d}\xi \mathrm{d}{\xi }^{\prime }$$

$${\int }_0^1{\tilde{C}}_{-\xi }(\eta ,{\eta }^{\prime }){\mathbf{g}}_{:,j}({\eta }^{\prime })\mathrm{d}{\eta }^{\prime }={\mu }_j{\mathbf{g}}_{:,j}(\eta )$$

$$\alpha (\xi ,\eta ;\theta )\approx \mathbf{f}(\xi )\mathbf{g}(\eta )\mathbf{\gamma }(\theta )$$

$$\mathbf{\gamma }(\theta )={\int }_0^1{\int }_0^1\mathbf{g}(\eta {)}^{\mathrm{T}}\mathbf{f}(\xi {)}^{\mathrm{T}}\alpha (\xi ,\eta ;\theta )\mathrm{d}\xi \mathrm{d}\eta$$

$$⟨f_i,f_j⟩$$

$$⟨{\mathbf{g}}_{:,i},{\mathbf{g}}_{:,j}⟩$$

$$⟨{\gamma }_i,{\gamma }_j⟩$$

$$\alpha (\xi ,\eta ,\zeta ;\theta )\approx \mathbf{f}(\xi {)}_{1\times n_1}\mathbf{g}(\eta {)}_{n_1\times n_2}\mathbf{h}(\zeta {)}_{n_2\times n_3}\mathbf{\gamma }(\theta {)}_{n_3\times 1}$$

$${\tilde{C}}_f(\eta ,{\eta }^{\prime })={\int }_0^1{\int }_0^1\tilde{C}(\xi ,\eta ,\zeta ;{\xi }^{\prime },\eta ,\zeta )\mathrm{d}\eta \mathrm{d}\zeta$$

$${\tilde{C}}_{-\xi }(\eta ,{\eta }^{\prime }{)}_{n_1\times n_1}={\int }_0^1{\int }_0^1\mathbf{f}(\xi {)}^{\mathrm{T}}{\int }_0^1\tilde{C}(\xi ,\eta ,\zeta ;{\xi }^{\prime },{\eta }^{\prime },\zeta )\mathrm{d}\zeta \mathbf{f}({\xi }^{\prime })\mathrm{d}\xi \mathrm{d}{\xi }^{\prime }$$

$$\begin{array}{rl}{\tilde{C}}_{-\xi \eta }(\zeta ,{\zeta }^{\prime }{)}_{n_2\times n_2}=& {\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1\mathbf{g}(\eta {)}^{\mathrm{T}}\mathbf{f}(\xi {)}^{\mathrm{T}}\tilde{C}(\xi ,\eta ,\zeta ;{\xi }^{\prime },{\eta }^{\prime },{\zeta }^{\prime })\\ & \mathbf{f}({\xi }^{\prime })\mathbf{g}({\eta }^{\prime }))\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}{\xi }^{\prime }\mathrm{d}{\eta }^{\prime }\end{array}$$

$$\mathbf{\gamma }(\theta )={\int }_0^1{\int }_0^1{\int }_0^1\mathbf{h}(\zeta {)}^{\mathrm{T}}\mathbf{g}(\eta {)}^{\mathrm{T}}\mathbf{f}(\xi {)}^{\mathrm{T}}\alpha (\xi ,\eta ,\zeta ;\theta )\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta$$

$$⟨f_i,f_j⟩$$

$$⟨{\mathbf{g}}_{:,i},{\mathbf{g}}_{:,j}⟩$$

$$⟨{\mathbf{h}}_{:,i},{\mathbf{h}}_{:,j}⟩$$

$$⟨{\gamma }_i,{\gamma }_j⟩$$

$$F_k(\xi ,\eta )=\mathbf{f}(\xi ){\hat{\mathbf{g}}}_{:,k}(\eta )$$

$$\begin{array}{rl}& {\int }_0^1{\int }_0^1\tilde{C}(\xi ,\eta ;{\xi }^{\prime },{\eta }^{\prime })\mathbf{f}({\xi }^{\prime }){\hat{\mathbf{g}}}_{:,k}({\eta }^{\prime })\mathrm{d}{\xi }^{\prime }\mathrm{d}{\eta }^{\prime }={\lambda }_k\mathbf{f}(\xi ){\hat{\mathbf{g}}}_{:,k}(\eta )\\ \Rightarrow & {\int }_0^1\mathbf{f}(\xi {)}^{\mathrm{T}}{\int }_0^1{\int }_0^1\tilde{C}(\xi ,\eta ;{\xi }^{\prime },{\eta }^{\prime })\mathbf{f}({\xi }^{\prime }){\hat{\mathbf{g}}}_{:,k}({\eta }^{\prime })\mathrm{d}{\xi }^{\prime }\mathrm{d}{\eta }^{\prime }\mathrm{d}\xi ={\lambda }_k{\int }_0^1\mathbf{f}(\xi {)}^{\mathrm{T}}\mathbf{f}(\xi )\mathrm{d}\xi {\hat{\mathbf{g}}}_{:,k}(\eta )\\ \Rightarrow & {\int }_0^1{\tilde{C}}_{-\xi }(\eta ,{\eta }^{\prime }){\hat{\mathbf{g}}}_{:,k}({\eta }^{\prime })\mathrm{d}{\eta }^{\prime }={\lambda }_k{\hat{\mathbf{g}}}_{:,k}(\eta )\\ \Rightarrow & {\hat{\mathbf{g}}}_{:,k}(\eta )\mathrm{is}\mathrm{a}\mathrm{column}\mathrm{of}\mathbf{g}(\eta ),{\lambda }_k\in \mathrm{M}\\ \Rightarrow & E_1\subseteq E_2,\Lambda \in M.\end{array}$$

$$\begin{array}{rl}& {\epsilon }_k(\xi ,\eta )={\int }_0^1{\int }_0^1\tilde{C}(\xi ,\eta ;{\xi }^{\prime },{\eta }^{\prime })\mathbf{f}({\xi }^{\prime }){\mathbf{g}}_{:,k}({\eta }^{\prime })\mathrm{d}{\xi }^{\prime }\mathrm{d}{\eta }^{\prime }-{\mu }_k\mathbf{f}(\xi ){\mathbf{g}}_{:,k}(\eta )\\ \Rightarrow & {\int }_0^1\mathbf{f}(\xi {)}^{\mathrm{T}}{\epsilon }_k(\xi ,\eta )\mathrm{d}\xi ={\int }_0^1{\tilde{C}}_{-\xi }(\eta ,{\eta }^{\prime }){\mathbf{g}}_{:,k}({\eta }^{\prime })\mathrm{d}{\eta }^{\prime }-{\mu }_k{\mathbf{g}}_{:,k}(\eta )=0\\ \Rightarrow & {\epsilon }_k(\xi ,\eta )\equiv 0\\ \Rightarrow & {\int }_0^1{\int }_0^1\tilde{C}(\xi ,\eta ;{\xi }^{\prime },{\eta }^{\prime })\mathbf{f}({\xi }^{\prime }){\mathbf{g}}_{:,k}({\eta }^{\prime })\mathrm{d}{\xi }^{\prime }\mathrm{d}{\eta }^{\prime }={\mu }_k\mathbf{f}(\xi ){\mathbf{g}}_{:,k}(\eta )\\ \Rightarrow & E_2\subseteq E_1,\mathrm{M}\in \Lambda .\end{array}$$

$$\begin{array}{rl}{f}_{\mathbf{\Xi }}(\mathbf{u};\mathbf{\kappa })=& f_G(\mathbf{u})[1+\frac{1}{3!}\sum _{i,j,k}{\kappa }^{i,j,k}{h}_{ijk}(\mathbf{u})+\frac{1}{4!}{\kappa }^{i,j,k,l}{h}_{ijkl}(\mathbf{u})+\frac{1}{5!}{\kappa }^{i,j,k,l,m}{h}_{ijklm}(\mathbf{u})+\\ & \frac{1}{6!}\left({\kappa }^{i,j,k,l,m,n}+10{\kappa }^{i,j,k}{\kappa }^{l,m,n}\right)h_{ijklmn}(\mathbf{u})+\cdots ]\end{array}$$

$$\begin{array}{rl}{\breve{C}}_3(i,j,k)=& {\int }_{[0,1{]}^6}{\tilde{C}}_3(\xi ,\eta ;{\xi }^{\prime },{\eta }^{\prime };{\xi }^{\prime \prime },{\eta }^{\prime \prime })(\mathbf{g}(\eta {)}_{:,i}^{\mathrm{T}}\mathbf{f}(\xi {)}^{\mathrm{T}})(\mathbf{g}({\eta }^{\prime }{)}_{:,j}^{\mathrm{T}}\mathbf{f}({\xi }^{\prime }{)}^{\mathrm{T}})\\ & (\mathbf{g}({\eta }^{\prime \prime }{)}_{:,k}^{\mathrm{T}}\mathbf{f}({\xi }^{\prime \prime }{)}^{\mathrm{T}})\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}{\xi }^{\prime }\mathrm{d}{\eta }^{\prime }\mathrm{d}{\xi }^{\prime \prime }\mathrm{d}{\eta }^{\prime \prime }\end{array}$$

$$\begin{array}{rl}G(u_1,\dots ,u_a)\approx & G(u_1,{\mathcal{I}}^{>1})[G({\mathcal{I}}^{⩽1},{\mathcal{I}}^{>1}){]}^{-1}G({\mathcal{I}}^{⩽1},u_2,{\mathcal{I}}^{>2})\\ & [G({\mathcal{I}}^{⩽2},{\mathcal{I}}^{>2}){]}^{-1}\cdots G({\mathcal{I}}^{⩽a-1},u_a)\end{array}$$

$$\{{\mathcal{I}}^{⩽k\ast },{\mathcal{I}}^{>k\ast }\}=\arg \underset{\begin{array}{c}{\mathcal{I}}^{⩽k},{\mathcal{I}}^{>k}\end{array}}{\max }|\det G({\mathcal{I}}^{⩽k},{\mathcal{I}}^{>k})|$$

Full text

††journal: a journal

1 Introduction

Uncertainty quantification in engineering and applied sciences often requires random field descriptions of spatial variability of uncertain media and loads, etc. Essentially, random field representation is a problem of data compression and reconstruction. Information of a random field must be compressed in a sufficiently small number of deterministic functions and random variables which can reconstruct the random field with acceptable accuracy.

A frequently used approach is the Karhunen-Loève (K-L) expansion [1, 2] (also named as Principal Component Analysis (PCA) in the statistical community) which was originally proposed for representing stochastic processes. A continuous-time stochastic process is represented with countable modes and latent factors. To model spatial variability of a physical quantity in a multidimensional space, the most common approach is directly extending K-L expansion by treating the spatial coordinate as a macro parameter (named as trivial PCA throughout this paper). The modes are eigenfunctions of an integral operator induced by the covariance function and the latent factors are uncorrelated in the sense of second-order moments. A comprehensive review of the methods (including Galerkin-type and collocation-type ones) for computing the modes was made in [3]. Other procedures such as spectral element projection [4], multilevel finite element method [5], two-dimensional Haar wavelet Galerkin method [6] and high-order polynomial-based Ritz-Galerkin method [7] were also proposed in recent years. To tackle the difficulties of the previous methods in treating complex geometries, modes were projected onto a subspace spanned by isogeometric basis functions in recent papers [8, 9, 10]. Meanwhile, several recent works tried to avoid the computationally expensive eigenvalue problem in trivial PCA based on the idea of variable separation. Ghosh et al. [11] redefined the covariance function of a spatial-temporal random field, and computed both the temporal and spatial modes with CP and Tucker decomposition, respectively. Zentner et al. [12] and Guo et al. [13] independently proposed a hierarchical orthogonal decomposition method to split the temporal and spatial modes. Similar idea was also shown in [14], however, this work did not really achieve variable separation since the modes still contain multiple coordinates. The stepwise covariance matrix decomposition method [15] derives a Tucker-type representation and is only suitable for rank-1 covariance functions.

Latent factors also play an important role for describing the probabilistic structure of the random field. For non-Gaussian random fields which widely exist in practice, the latent factors are generally non-Gaussian and have higher-order dependencies, making the representation of their probabilistic structure a nontrivial task. To achieve this task, several methods are dedicated to match the prescribed single-point marginal distribution function and covariance function. Phoon et al. [16] and Dai et al. [17] proposed similar iterative procedures to update the marginal distribution of each latent factor. Kim et al. [18] iteratively updated the covariance function of the underlying Gaussian process based on the translation process theory. Other works tried to estimate the probabilistic structure of the latent factors with data-driven procedures [19, 20, 21, 22]. Since this paper is focused on representing a random field with prescribed statistics, we will not give further comments on these methods.

For the computation of modes, although the complex physical domains can be handled by the aforementioned isogeometric analysis (IGA)-based methods, however, the underlying architecture is still the trivial PCA in physical domain, making these methods suffer from the common drawbacks of the trivial PCA-based ones: (1) high computational and memory demand; (2) natural structure and correlation are broken, leading to loss of potentially more compact and useful representations [23]. Variable separation is indeed a promising way for reducing the computational scale of the modes and preserving the natural structure. However, existing methods can only represent random fields defined on domains which have Cartesian product decomposition format (i.e. cases where the random fields belong to structured data), while in many cases the random fields belong to unstructured data. Moreover, these methods generate a large number of latent factors which will bring heavy computational burden for subsequent stochastic analysis. In addition, all the existing methods in both categories requires manually predefining tensor product basis (or collocation points) to discretize the modes and selecting quadrature points to compute the stiffness matrices. Lack of automation in these selections limits generality and adaptivity [24].

It is well known that a general non-Gaussian random field should be described with multiple statistics rather than using the covariance function only. The most widely used description which combines the covariance function and single-point marginal distribution function cannot reflect the higher-order and nonlinear correlation structure of the random field. It was pointed out in [25] that joint distributions of multiple points can be non-Gaussian even if each marginal distribution is Gaussian. It is also well known that the cumulant functions with order three or higher are always zeros, making the cumulant functions be important measures of non-Gaussianity. Therefore, it will be better to develop cumulant descriptions of the latent factors based on prescribed cumulant functions.

This paper aims at overcoming the challenges in both aspects mentioned above by proposing novel theoretical and algorithm frameworks. Motivated by the aforementioned IGA-based works, we also use isogeometric transformation to deal with complex domains. However, rather than following the trivial PCA in the physical domain, we propose a new architecture by representing the random field on the parametric domain in tensor train(TT) format. The core of the architecture is a newly developed rank-revealing algorithm for separating variables of the cumulant functions. The need for predefining tensor product basis or collocation points is totally removed, which greatly enhances adaptivity. Moreover, higher-order cumulant tensors of the latent factors can be conveniently computed in a unified framework, which is beneficial for further dimension reduction.

The rest of this paper is organized in five sections. The proposed theoretical framework and the computational procedures are detailed in Section 2 and 3, respectively. In Section 4, we compare the proposed framework with three related methods. Particularly, the relationship with independent component analysis (ICA) will be discussed since this approach also uses higher-order cumulants for dimension reduction. In Section 5, three examples with increasing dimensions are employed to validate the performance of the proposed framework.

2 Theoretical framework of the proposed method

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space where $\Omega$ is a sample space, $\mathcal{F}$ is a $\sigma$-field on $\Omega$ and $\mathbb{P}$ is a probability measure. Let $(E,\mathcal{E})$ is the measurable space on the admissible set $E$ where each $f\in E$ is defined on a bounded domain $D\subseteq\mathbb{R}^{d},d=1,2,3$. A random field is defined as a measurable mapping $\omega:(\Omega,\mathcal{F},\mathbb{P})\to(E,\mathcal{E})$. To explicitly represent this abstract mapping, an auxiliary measurable space $(\mathbb{R}^{n},\mathcal{B}^{n}),n\in\mathbb{N}^{+}$ is needed and the corresponding two mappings $\Theta:(\Omega,\mathcal{F},\mathbb{P})\to(\mathbb{R}^{n},\mathcal{B}^{n})$ and $H:(\mathbb{R}^{n},\mathcal{B}^{n})\to(E,\mathcal{E})$ are to be found such that $\omega=H\circ\Theta$. The task is to capture the major part of information in $\omega$ with $n$ as small as possible.

2.1 Space transformation

The key to overcome the first difficulty is to represent the physical domain with a structured parametric domain. For a wide variety of curve-type, surface-type and solid-type domains, exact coordinate transformations can be constructed by using NURBS-based isogeometric transformation formulated as Eq.(1)

[TABLE]

where $R_{i}^{p},R_{i,j}^{p,q}$ and $R_{i,j,k}^{p,q,r}$ are NURBS basis functions defined as Eq.(2),

[TABLE]

$B_{i},B_{i,j}$ and $B_{i,j,k}$ are control points and $w_{i},w_{i,j}$ and $w_{i,j,k}$ are weights. Each $N_{i}^{p}$ is a $p$-degree ($(p+1)$-order) B-spline basis function. Given a knot vector $\Xi=(\xi_{1},\dotsc,\xi_{n+p+1})$ (a non-decreasing set of coordinates in the parametric space [0,1]), $N_{i}^{p}$ are defined recursively by the Cox-de Boor recursion formula in Eq.(3)

[TABLE]

where 0/0=0. Thus, given three ingredients: control points, weights and knot vectors, the isogeometric transformation can be constructed and formulated in a uniform format as Eq.(4).

[TABLE]

where $\bm{x}\in D\subseteq\mathbb{R}^{d},\ \bm{\xi}\in[0,1]^{m}$ and $m\leqslant d$.

2.2 Generalized Karhunen-Loève expansion in the parametric space

After representing the physical coordinates with parametric ones, the original random field $\omega(\bm{x},\theta)$ is transformed to the parametric space. Denote $\alpha(\bm{\xi},\theta)=\omega(\bm{x}(\bm{\xi}),\theta)$ and the corresponding covariance function $\tilde{C}(\bm{\xi},\bm{\xi}^{\prime})=C(\bm{x}(\bm{\xi}),\bm{x}(\bm{\xi}^{\prime}))$. We seek to represent $\alpha(\bm{\xi},\theta)$ with a generalized K-L expansion by improving the hierarchical SVD method in [12]. Taking the case of $d$=2 as an example, the first step is to split the $\xi$-modes $\bm{f}(\xi)$ = $(f_{i}(\xi))_{1\times n_{1}}$ from $\alpha(\xi,\eta;\theta)$ by computing the eigenfunctions of the kernel in Eq.(5).

[TABLE]

Thus, a rank-$n_{1}$ decomposition is derived as Eq.(6)

[TABLE]

where

[TABLE]

and the corresponding covariance function is a matrix-valued function in Eq.(8).

[TABLE]

Then, rather than decomposing each component of $\bm{\alpha}_{-\xi}(\eta,\theta)$ separately as in [12], $\bm{\alpha}_{-\xi}(\eta,\theta)$ is treated as whole object by regarding the index as a new coordinate. Next, the $\eta$-modes of $\alpha(\xi,\eta;\theta)$ (also the $\eta$-modes of $\bm{\alpha}_{-\xi}(\eta,\theta)$), denoting $\bm{g}(\eta)_{n_{1}\times n_{2}}$, are split by computing the eigenpairs of $\tilde{C}_{-\xi}(\eta,\eta^{\prime})$, see Eq.(9).

[TABLE]

Finally, we get the tensor train decomposition of $\alpha(\xi,\eta;\theta)$ as Eq.(10)

[TABLE]

where

[TABLE]

The component functions in Eq.(10) satisfy orthogonal conditions in Eq.(12)

[TABLE]

All the results can be directly extended to the case of $d$=3 as Eqs.(13) to (18).

[TABLE]

where each column of $\bm{f}(\xi)$ is an eigenfunction of the covariance kernel $\tilde{C}_{f}(\xi,\xi^{\prime})$ in Eq.(14),

[TABLE]

each column of $\bm{g}(\eta)$ is an eigenfunction of $\tilde{C}_{-\xi}(\eta,\eta^{\prime})$ in Eq.(15)

[TABLE]

and each column of $\bm{h}(\zeta)$ is an eigenfunction of $\tilde{C}_{-\xi\eta}(\zeta,\zeta^{\prime})$ in Eq.(16).

[TABLE]

[TABLE]

[TABLE]

where $\nu_{i}$ is the $i$th eigenvalue of the covariance kernel $\tilde{C}_{-\xi\eta}(\zeta,\zeta^{\prime})$ defined in Eq.(16).

Taking $m$=2 as an example, we have the following theorem:

Theorem 1**.**

The generalized K-L expansion in Eq.(10) is equivalent to the trivial PCA in the parametric space when $n_{1},\ n_{2}\to\infty$.

Proof.

Denote $E_{1}$ = $\{F_{k}^{(1)}(\xi,\eta)\}_{k=1}^{\infty}$ and $\Lambda$ = $\{\lambda_{k}\}_{k=1}^{\infty}$ as the sets of 2D modes and eigenvalues in the trivial PCA and $E_{2}$ = $\{F_{k}^{(2)}(\xi,\eta)=\bm{f}(\xi)\bm{g}_{:,k}(\eta)\}_{k=1}^{\infty}$ and $\rm M$ = $\{\mu_{k}\}_{k=1}^{\infty}$ as the set of 2D modes in the generalized K-L expansion. We only need to prove that $E_{1}$ = $E_{2}$ and $\Lambda$ = $\rm M$.

First, we prove $E_{1}\subseteq E_{2}$ and $\Lambda\in\rm M$. For each $k=1,2,\dotsc$, by projecting $F_{k}(\xi,\eta)$ on $\bm{f}(\xi)$,

[TABLE]

we derive

[TABLE]

Then, we prove $E_{2}\subseteq E_{1}$ and $\rm M\in\Lambda$. For each $k=1,2,\dotsc$, we derive

[TABLE]

Finally, we get $E_{1}$ = $E_{2}$ and $\Lambda$ = $\rm M$.

∎

Thus, Eq.(10) is mean-square convergent when $n_{1},\ n_{2}\to\infty$. in It is easy to verify that this theorem is valid for other values of $m$.

2.3 Representing higher-order cumulants of latent factors

It is well known that the probabilistic structure of a random field is uniquely defined by its family of finite-dimensional marginal distributions. For arbitrary $n$ points $\bm{\Xi}=(\bm{\xi}_{1},\dotsc,\bm{\xi}_{n})$ in parametric space and arbitrary $\bm{u}=(u_{1},\dotsc,u_{n})\in\mathbb{R}^{n}$, according to the A-type Gram-Charlier series, the corresponding marginal probability density function (PDF) can be represented as Eq.(19)

[TABLE]

where $f_{G}$ is the PDF of the Gaussian distribution with mean zero and covariance function $\tilde{C}(\bm{\xi},\bm{\xi}^{\prime})$, $h$ with subscripts are Hermite polynomials and $\bm{\kappa}$ is the family of cumulant functions consists of elements $\kappa^{i,j,k}$ = $\tilde{C}_{3}(\bm{\xi}_{i},\bm{\xi}_{j},\bm{\xi}_{k})$ (the third-order cumulant function) and so on. Meanwhile, $\bm{\kappa}$ is defined by the family of marginal PDFs by definition. Thus, $\bm{\kappa}$ is equivalent to the family of marginal PDFs by definition in terms of describing the probabilistic structure of a random field. In practical problems, only finite orders of cumulant functions are available. To absorb the information of high-order ($\geqslant$3) cumulant functions into the representation of $\alpha(\bm{\xi};\theta)$, taking the case of $d$=2 as an example, since the reconstructed random field has the form in Eq.(10), the relationship between the third-order cumulant functions and the third-order cumulants of latent factors is expressed as Eq.(20).

[TABLE]

This relationship can be directly extended to the cases of other values of $d$ and cumulant orders. Unfortunately, $\breve{C}_{3}(i,j,k)$ is generally not orthogonal decomposable [26], hence, we have to resort to HOSVD to further reduce the dimensionality of the latent factors. Finally, by combing the coordinate transformation in Eq.(4) with Eqs.(10) (or (13)) and (20), we get a parametric representation of the original random field $\omega(\bm{x};\theta)$ in the sense of given cumulants.

3 Computational procedure

In this section, we seek an algorithm framework which is general enough to overcome the difficulties in both aspects mentioned in Section 1.

3.1 Space transformation

Given three ingredients: control points, weights and knot vectors, to avoid redundant computations, geometry of the physical domain defined in Eq.(4) is represented with the fast isogeometric transformation algorithm [27].

3.2 Generalized Karhunen-Loève expansion in the parametric space

After representing the physical domain $D$ with a structured parametric domain $[0,1]^{m}$, the core of the new algorithm framework is representing the modes of $\alpha(\bm{\xi};\theta)$. This task is accomplished with the following two steps.

3.2.1 Tensor train decomposition of cumulant functions

According to the theoretical framework in the previous section, it will be very beneficial for fast computations if the cumulant functions have separable forms. Borrowing the idea in [28], we propose a tensor train decomposition algorithm to obtain a low-rank separable approximation to a $K$th ($K\geqslant 2$) cumulant function $\tilde{C}_{K}(\bm{\xi}_{1},\dotsc,\bm{\xi}_{K})$.

For the sake of simplicity, $\tilde{C}_{K}(\bm{\xi}_{1},\dotsc,\bm{\xi}_{K})$ is redefined as an auxiliary function $G:[0,1]^{a}\to\mathbb{R}\ (a=mK)$ where $\bm{u}=(u_{1},\dotsc,u_{a})\in[0,1]^{a}$ is a permutation of the original coordinates. The idea is reconstruct $G$ only by using some of its fibers. More precisely, we seek a separable representation of $G$ as Eq.(21)

[TABLE]

where $\mathcal{I}^{\leqslant k}$ and $\mathcal{I}^{>k}$ are the $k$th pair of interpolation sets (both have cardinality $r_{k}$ which is the $k$th rank). The optimal choice of the $k$th pair of interpolation sets is the solution to Eq.(22).

[TABLE]

However, the search for the maximum-volume submatrix is an NP (non-deterministic polynomial)-hard problem. Hence, we propose a heuristic algorithm to find a quasi-optimal choice of $\{\mathcal{I}^{\leqslant k},\mathcal{I}^{>k}\}_{k=1}^{a-1}$ and derive an adaptive tensor train decomposition of a cumulant function, see Algorithm 1.

After initialization, the interpolation sets for each two neighboring dimensions are progressively enriched by the left-to-right and right-to-left sweeps until the maximum approximation error is smaller than the prescribed threshold. The interpolation sets enriched in this way are two-side nested, i.e. $\forall k\in\{1,\dotsc,a-1\}$, $\mathcal{I}^{\leqslant k+1}\subseteq\mathcal{I}^{\leqslant k}\times[0,1]$, $\mathcal{I}^{>k}\subseteq[0,1]\times\mathcal{I}^{>k+1}$. Then, each fiber in Eq.(21) is adaptively represented in the chebfun format [29]. Next, fibers of each two neighboring dimensions form a chebmatrix as in lines 46, 48, 50.

For a special case where $K=2$ and $m=1$, Algorithm 1 degenerates into Algorithm 3.

3.2.2 Computation of modes in each direction using reconstructed covariance functions

Computation of modes in each direction is performed with $K=2$. When $m=1$, decomposition of the $\tilde{C}(\xi,\xi^{\prime})$ is trivial by using Algorithm 3, and we obtain a separable representation of $\tilde{C}(\xi,\xi^{\prime})$ as Eq.(27).

[TABLE]

Before computing the modes, accuracy of the TT decomposition should be checked by computing the empirical global relative error in Eq.(28)

[TABLE]

where $\{(\xi_{i},\xi^{\prime}_{i})\}_{i=1}^{N}$ is the set of test samples. If $\varepsilon_{\mathrm{g}}$ is larger than a prescribed tolerance $tol_{\mathrm{g}}$, then go to the next step; else, adjust values of $m_{1}$ or $tol$ in Algorithm 3 and re-compute the TT decomposition. Next,taking the multiplicative form in Eq.(27) into account, by following the idea of [30], we propose a specific algorithm for computing the $\xi$-modes $\bm{f}(\xi)_{1\times n_{1}}$ and the eigenvalues $\{\lambda_{i}\}_{i=1}^{n_{1}}$ see Algorithm 4.

For the QR decomposition in Algorithm 4, the first choice is the Householder triangularization [31](the default choice in the subsequent numerical experiments) due to its good numerical stability. Unfortunately, this method is too time-consuming for large-scale chebmatrices since the computational time of both the plus and inner-product operations are non-negligible. Hence, a second choice is the Cholesky decomposition of $\bm{M}^{\mathrm{T}}\bm{M}$. This method produces almost the same results as the first one with much higher efficiency, and still have good numerical stability for sufficiently large eigenvalues. When $m=2$, first, the auxiliary function $G$ is defined as: $G(u_{1},\dotsc,u_{4})$ = $G(\xi,\eta;\eta^{\prime},\xi^{\prime})$ = $C(\bm{x}(\xi,\eta);\bm{x}_{\mathrm{r}}(\eta^{\prime},\xi^{\prime}))$ where $\bm{x}_{\mathrm{r}}(\eta^{\prime},\xi^{\prime})$ = $\bm{x}(\xi^{\prime},\eta^{\prime})$, i.e. $G(\xi,\eta;\eta^{\prime},\xi^{\prime})$ = $\tilde{C}(\xi,\eta;\xi^{\prime},\eta^{\prime})$. Then, by applying Algorithm 1 on $G$, the tensor train decomposition of $\tilde{C}(\xi,\eta;\xi^{\prime},\eta^{\prime})$ is obtained as Eq.(29).

[TABLE]

After passing the global random test, covariance kernel in Eq.(5) is obtained as Eq.(30).

[TABLE]

The $\xi$-modes $\bm{f}(\xi)_{1\times n_{1}}$ are computed by using Algorithm 4. Next, the covariance kernel in Eq.(8) is obtained as Eq.(31).

[TABLE]

To compute the $\eta$-modes $\bm{g}(\eta)$ in Eq.(10), we join the supports of $\tilde{C}_{-\xi}(\eta,\eta^{\prime})$ in both dimensions, which is equivalent to joining each column of $\bm{M}_{3}(\eta)$ and each row of $\bm{M}_{4}(\eta^{\prime})$, respectively, deriving Eq.(32).

[TABLE]

Then, by applying Algorithm 4 on $\bm{M}_{3\mathrm{J}}(\eta_{\mathrm{J}})$ and $\bm{M}_{4\mathrm{J}}(\eta_{\mathrm{J}}^{\prime})$, we obtain a joint form of $\bm{g}(\eta)$, denoting $\bm{g}_{\mathrm{J}}(\eta)_{1\times n_{2}}$, and the eigenvalues $\{\mu_{i}\}_{i=1}^{n_{2}}$. Next, $\bm{g}(\eta)_{n_{1}\times n_{2}}$ is obtained by disjoining the support of each column of $\bm{g}_{\mathrm{J}}(\eta)$ with $n_{1}$ equal intervals.

The above procedure can be straightforwardly extended to the case of $m=3$. By letting $G(\xi,\eta,\zeta;\eta^{\prime},\xi^{\prime},\zeta^{\prime})$ = $\tilde{C}(\xi,\eta,\zeta;\xi^{\prime},\eta^{\prime},\zeta^{\prime})$, first, we get Eq.(33).

[TABLE]

Then, the covariance kernel in Eq.(14) is derived as Eq.(34),

[TABLE]

and the $\xi$-modes $\bm{f}(\xi)_{1\times n_{1}}$ are computed by applying Algorithm 4 on the $\bm{M}_{1}(\xi)$ and $\bm{M}_{2}(\xi^{\prime})$ above. Next, the covariance kernel in Eq.(15) is derived as Eq.(35)

[TABLE]

By joining supports, applying Algorithm 4 and disjoining supports sequentially, the $\eta$-modes $\bm{g}(\eta)$ are obtained. Next, the covariance kernel in Eq.(16) is derived as Eq.(36),

[TABLE]

By joining supports, applying Algorithm 4 and disjoining supports sequentially, the $\zeta$-modes $\bm{h}(\zeta)$ and the eigenvalues $\{\nu_{i}\}_{i=1}^{n_{3}}$ are obtained.

3.3 Representing higher-order cumulants of latent factors

Having obtained the modes in all directions and the covariance matrix of the latent factors, the final task is to efficiently represent higher-order cumulant tensors of the latent factors. Taking $K$ = 3 as an example, by recalling Eqs.(11) and (17), a single latent factor is expressed as Eq.(37) or uniformly written as Eq.(38).

[TABLE]

[TABLE]

By using Eq.(38), the third cumulant tensor of the latent factors is derived as Eq.(39).

[TABLE]

By letting the auxiliary function $G\equiv\tilde{C}_{3}$, a tensor train decomposition of $\tilde{C}_{3}$ (denoted as $\tilde{C}_{3,\mathrm{TT}}$) is obtained with Algorithm 1. Since this approximation is multiplicatively separable, taking $m$=3 as an example, for each index $I\ (I=i,j,k)$, the corresponding part in the decomposition has the form $\hat{\bm{G}}_{3I-2}(\xi)\hat{\bm{G}}_{3I-1}(\eta)\hat{\bm{G}}_{3I}(\zeta)$. Meanwhile, $F_{I}(\bm{\xi})$ =$(\bm{h}(\eta))^{\mathrm{T}}_{i,:}\bm{g}(\eta)^{\mathrm{T}}\bm{f}(\xi)^{\mathrm{T}}$. Thus, after passing the global random test, Eq.(39) can be transformed to Eq.(40)

[TABLE]

where $\times^{1}$ is the mode(3,1) contracted product and each $\bm{A}_{i}$ is a third-order tensor with each slice $\bm{A}_{i(:,j,:)}$ expressed in Eq.(41).

[TABLE]

After obtaining $\bm{\breve{C}}_{3}$, it is possible to further reduce the dimension of $\bm{\gamma}(\theta)$ by applying an orthogonal transformation represented by HOSVD bases of $\breve{C}_{3}(i,j,k)$. Since cumulant tensors are super-symmetric in theory, the modes are identical in all dimensions and need to be computed in only a single dimension. To achieve this goal, first, we derive Eq.(42)

[TABLE]

where subscripts $(i)$ is the $i$-mode unfolding of a tensor. It is clear in Eq.(42) that the summation with respect to $j,k$ is decoupled. Then the pursued orthogonal transformation matrix $\bm{U}_{[3]}$=$\bm{U}_{(:,1:n)}$ is obtained by spectral decomposition of matrix $\bm{\breve{C}}_{3(1)}\bm{\breve{C}}_{3(1)}^{\mathrm{T}}$. The number of columns $n$ in $\bm{U}$ is selected as the minimum value satisfying Eq.(43)

[TABLE]

where $\bm{\lambda}_{2}$ and $\bm{\lambda}_{3}$ are diagonal matrices containing eigenvalues of $\tilde{C}_{-\xi\eta}(\zeta,\zeta^{\prime})$ and mode-1 eigenvalues of $\bm{\breve{C}}_{3(1)}$, respectively. Next, based on the properties of TT decomposition [32], the cumulant tensor of the transformed latent factors $\bm{\gamma}_{[3]}(\theta)$ = $\bm{U}_{[3]}^{\mathrm{T}}\bm{\gamma}(\theta)$ is obtained by Eq.(44)

[TABLE]

where $\bm{\acute{A}}_{j(:,:,i)}$ = $\bm{A}_{j(:,:,i)}\bm{U}$. The final expression of the random field is obtained as Eq.(45)

[TABLE]

Accuracy of this expression for matching the prescribed second and third-order cumulant function is checked by computing the empirical global relative error in Eq.(46)

[TABLE]

where $\tilde{C}_{\mathrm{TTf},k}$ is the $k$-th order cumulant function of the final TT expression and is computed by Eq.(47)

[TABLE]

If $\min(\varepsilon_{\mathrm{gf}k})$ is smaller than a prescribed threshold, then finish; else, decrease $tol_{3}$ and re-compute each $\varepsilon_{\mathrm{gf}k}$. The above procedure is readily extended to other values of $m$.

3.4 Summary

For the sake of clarity, the whole algorithm framework is summarized as follows:

Do space transformation; 2. 2.

Compute TT decomposition of $\tilde{C}$ and perform global random test; 3. 3.

Compute modes in each direction; 4. 4.

Compute TT decomposition of $\tilde{C}_{3}$ and perform global random test; 5. 5.

Compute $\bm{\breve{C}}_{3}$ and its HOSVD; 6. 6.

Compress modes and latent factors, and perform the final global random test.

When multiple higher-order cumulant functions are simultaneously given, we need to evaluate their relative importance and find an orthogonal transformation which can best explain a functional of these cumulant functions. This is a non-trivial task. The simple criterion given in [25] lacks objectivity. Hence, this problem still needs further research.

4 Comparing with related works

4.1 Tensor train-PCA

A recent work worth mentioning is the tensor train-PCA [33]. Given $N$ samples of tensor data $\mathcal{X}_{i}\in\mathbb{R}^{I_{1}\times\cdots I_{n}}$, $i=1,\dotsc,N$, the goal is to find $\mathcal{U}_{j}\in\mathbb{R}^{r_{j-1}\times I_{j}\times r_{j}}$, $j=1,\dotsc,n$, such that the distance of the points to the space spanned by $\{\mathcal{U}_{j}\}_{j=1}^{n}$ is minimized. The algorithm in the original article requires concatenation of all the samples, hence is in essence a Full-to-TT algorithm [34], and is a data-driven realization of a part of our theoretical framework in Section 2.2 as well. The modes are computed by applying SVD on the unfolded data matrices, which is computationally infeasible for large values of $N$ and $I_{j}$, $j=1,\dotsc,n$. Different from the tensor train-PCA, our goal is to find modes of a continuously indexed random field in each direction with prescribed cumulant functions, as in Section 2.2. Estimating cumulant functions from raw data is out of the scope of this work. In addition, manipulations of large-scale arrays never occur in our work.

4.2 Independent component analysis

According to [35, 36], with the key assumption that the latent factors are mutually independent, independent component analysis (ICA) consists of the two sequential steps: (1) standardization (also called pre-whitening) of the original random field using the covariance function and (2) rotation of the standardized random field to the latent factors named independent components by using higher-order (usually the fourth) cumulant functions. Thus, the proposed theoretical framework has a similar form as that of ICA. However, rather than pursuing the independent latent factors, the proposed framework just aims at reducing the dimensionality of the latent factors and represent their higher-order non-Gaussian dependence.

4.3 Isogeometric analysis-based K-L expansion

IGA-based K-L expansion were proposed in recent papers to represent random fields on complex domains. However, the underlying architecture of these methods is still trivial PCA in physical domain. Modes are projected onto a subspace spanned by IGA basis functions, and the coordinates are computed with Galerkin or collocation method. Hence, these methods still suffer from the drawbacks of trivial PCA. Different from these methods, the underlying architecture of this work is the tensor train decomposition in parametric domain. Highlights of the proposed framework over the IGA-based methods is summarized as follows:

Since we only need a space transformation and the geometry is exactly preserved at the lowest level of refinement, enrichment of NURBS spaces is never needed. 2. 2.

The original unstructured random field is transformed into a structured one, making it feasible for using the powerful tensor train decomposition to compress the random field. The tensor train representation has the merits of both the trivial PCA (low stochastic dimensionality) and other tensor-based representations (moderate physical dimensionality). 3. 3.

It is totally unnecessary to predefine any basis functions or discrete point set to represent the modes, which avoids subjectivity. All modes are automatically represented in chebfun format in the sense of machine precision. 4. 4.

By adaptive variable separation of the covariance function, all integrals are transformed into uni-dimensional ones, which greatly reduces the computational scale of each mode. In addition, the tedious processes of forming and assembling the stiffness matrices are never needed. All the eigenvalues are obtained by solving only a few (1-3) small-scale ($O(10^{0-2})$) matrix SVD problems. 5. 5.

Since the random field is represented in the parametric domain rather than the physical counterpart, Jacobians are never needed in the integrals. 6. 6.

By adaptive variable separation of the random field and the higher-order cumulant functions, cumulants of latent factors can be conveniently computed to capture the higher-order non-Gaussian correlations of the random field, which has never been discussed in existing literature of K-L expansion.

Remark: Since the random field is represented in different spaces, results of the proposed algorithm framework are not directly comparable with that of existing IGA-based methods.

5 Case studies

In this section, we examine the performance of the proposed algorithms by using three examples with increasing parametric dimensionalities. The underlying isogeometric transformation is achieved by using the toolbox in [37]. All the numerical experiments were performed using MATLAB (Version 2019a) on a notepad (core i5 CPU and 12GB RAM).

5.1 Example 1

In the first example, to clearly demonstrate performance of the proposed method, we consider a simple random field with $D=[0,1]\subseteq\mathbb{R}$ and covariance function defined in Eq.(48)

[TABLE]

where $\lambda_{k}=4/(\pi^{2}(2k-1)^{2})$ and $f_{k}(x)=\sqrt{2}\sin(\xi/\sqrt{\lambda_{k}})$. Following the steps in Section 3, first, the isogeometric transformation is constructed with the following ingredients: (1) knot vector $\Xi$ = (0,0,1,1); (2) control points $x_{1}$ =0, $x_{2}$ =1; (3) weights $\bm{W}$ = (1,1). It is easy to derive that the transformation is $x=\xi$. Then, eigenpairs of the compound integral operator with kernel $C_{2}(x(\xi),x(\xi^{\prime}))$ should be computed. Since the parametric space coincide with the physical space, the eigenpairs are exactly $\lambda_{k}$ and $f_{k}(\xi)$. Thus, TT-rank of $\tilde{C}$ is exactly 80. Comparisons are made from different aspects as below.

5.1.1 Visualization of the reconstructed compound covariance function and spatial distribution of the errors

First, an experiment is made by letting $N$=400 and $tol=1.0000\times 10^{-6}$. Results are illustrated in Fig.1.

Fig.1(a) shows the reconstructed covariance function. Spatial distribution of errors in Fig.1(b) indicates that this reconstruction is almost perfect since errors are consistently of $O(10^{-12})$. In addition, the interpolation set is automatically selected as the points shown in Fig.1(c), which indicates the abilities of the proposed algorithm for detecting the regularity of the target function and adapting to its features.

5.1.2 Investigations on global relative errors and ranks

To investigate the effect of $m_{1}$ in line 9 of Algorithm 3 on variability of the global relative errors and the rank, we did four groups of experiment with $m_{1}$ = 50, 100, 200 and 400, respectively. Each group consisted of 100 independent tests. Each test was performed by letting $N$ = 1000 in Eq.(28) and $tol$ = $10^{-6}$ in Algorithm 3. Results are illustrated in Figs.2.

Fig.2(a) shows the number of outliers decreases from 14 to 5 as $m_{1}$ increases from 50 to 400. Despite the outliers, both the variability (represented by difference between the upper and lower adjacent value) and the median of $\varepsilon_{\mathrm{g}}$ decreases with $m_{1}$ as well. Moreover, the decreasing rates decreases with $m_{1}$. These three observations indicate that both the average accuracy and the robustness of the proposed algorithm increases with $m_{1}$, while the increasing rates tend to be slower. When $m_{1}=400$, very accurate (median=$O(10^{-12})$) and robust (upper adjacent value=$O(10^{-12})$ and rank $\equiv$ 80) can be obtained. Meanwhile, Fig.2(b) shows that the proposed algorithm almost always reveals the exact rank (except only a few outliers) with all values of $m_{1}$.

Next, two additional groups of experiment were made to investigate the effects of $tol$ on the variability of the global relative errors and the rank by setting $tol$ = $10^{-4}$ and $10^{-5}$ , respectively. $m_{1}$ is kept as 100 in both groups. Other settings are the same as that of the previous four groups. Results are illustrated in Fig.3.

Fig.3(a) shows that variability of $\varepsilon_{\mathrm{g}}$ significantly decreases as $tol$ decreases from $10^{-4}$ to $10^{-5}$, then almost remains constant as $tol$ decreased to $10^{-6}$. The median almost remains constant in all cases. Similar features are displayed in the figure of rank. Thus, to obtain sufficient robust results, $tol$ is suggested to be of $O(10^{-5})$ or smaller in this example.

Next, three independent tests were proceeded by setting $m_{1}$ = 400 and $tol$ = $10^{-6}$ to illustrate the accuracy and robustness of Algorithm 3 for computing the modes. Results are shown in Fig.4.

Fig.4 shows that for all orders of mode (from 1 to 80), the proposed algorithm consistently provides almost perfect reconstructions (the differences can not be recognized by naked eyes) in all the three tests, and so do the eigenvalues. These observations verify the high levels of accuracy and robustness of the proposed method from another point of view.

5.1.3 Dimension reduction of the latent factors

Before proceeding dimension reduction, we added an experiment with $K=2$, $m_{1}=400$ and $tol=10^{-6}$, and computed $\varepsilon_{\mathrm{g}2}$ = $3.8889\times 10^{-13}\ (N=1000)$. Then, rank-80 eigenpairs were computed with $\varepsilon_{\mathrm{gf}2}$ = $1.2715\times 10^{-2}\ (N=1000)$. Next, we assume the third-order cumulant function $C_{3}(x,y,z)$ as Eq.(49)

[TABLE]

and computed the TT decomposition $\tilde{C}_{3,\mathrm{TT}}$ with $\varepsilon_{\mathrm{g}3}$ = $1.1213\times 10^{-12}\ (N=1000)$ and ranks $\bm{r}$=$(80,80)$. Finally, by setting $tol_{3}=9.9990\times 10^{-1}$, the effective dimensionality of the latent factors is greatly reduced from 80 to 5 with $\varepsilon_{\mathrm{gf}3}$ = $1.2555\times 10^{-2}\ (N=1000)$. The time cost (seconds) of each quantity is listed in Table 1.

The results above indicate that the goal of matching the second and third-order cumulant functions within moderate time has been achieved.

5.2 Example 2

In this example, let the physical domain $D$ be a bilinear surface in Fig.5

where the knot vectors are $\Xi_{1}$ = (0,0,1,1) and $\Xi_{2}$ = (0,0,1,1), and the control points and weights are listed in Table 2.

The covariance function in physical domain is defined as Eq.(50)

[TABLE]

where $\sigma^{2}=1$, $b=1$ and $L=1$.

First, we decomposed the covariance function in the parametric space by using Algorithm 1 with $tol=1.0000\times 10^{-6}$ and each $m_{k}=800$ in Algorithm 2, and computed the global relative errors $\varepsilon_{\mathrm{g}2}=2.1065\times 10^{-7}(N=1000)$ and ranks $\bm{r}_{2}=(8,37,8)$. Then, we computed the eigenpairs in each direction with $\varepsilon_{\mathrm{gf}2}$ = $1.0847\times 10^{-2}\ (N=1000)$. To make comparisons, we also proceeded trivial PCA in the parametric space by using the FEM with rectangular bilinear elements and 40$\times$40 mesh. Main results are illustrated in Fig.6.

From Fig.6 we can see that the eigenvalues computed with the proposed method are highly consistent with the referenced values (the differences cannot be recognized by naked eyes). The dominant modes (whose corresponding eigenvalues are significantly larger than zero, the first four modes in this example) computed with the proposed method also show good agreement with the reference ones. Notice that the ’boundary effect’ arise in the FEM solutions which themselves are approximate ones, accuracy of the TT modes cannot be evaluated merely by the discrepancies from the FEM solutions. Moreover, the TT solution is more smooth near the boundary since the ’boundary effect’ hardly occurs. Thus, the TT solutions seem to be more accurate than the ones of FEM. On the whole, the above results are consistent with Theorem 1.

Next, by assuming the third-order cumulant function as Eq.(51)

[TABLE]

where $\sigma^{3}=1$, $b=1$ and $L=1$, we computed the TT decomposition of the counterpart in the parametric domain (denoted as $\tilde{C}_{3}(\xi,\eta;\xi^{\prime},\eta^{\prime};\xi^{\prime\prime},\eta^{\prime\prime})$) by setting $tol=1.0000\times 10^{-5}$ and each $m_{k}=800$ in Algorithm 2, and computed the global relative errors $\varepsilon_{\mathrm{g}3}=6.9247\times 10^{-6}(N=10000)$ and ranks $\bm{r}_{3}=(8,41,129,38,8)$. After HOSVD of $\bm{\breve{C}}_{3}$ ($tol_{3}=9.9990\times 10^{-1}$), the number of latent factors are reduced from 37 to 11 $\varepsilon_{\mathrm{gf}3}=3.9148\times 10^{-3},N=1000$. The time cost (seconds) of each quantity is listed in Table 3.

Computing the FEM solution above costs 5523s, which indicates that efficiency of the proposed algorithm is an order of magnitude higher than that of FEM. Thus, the goal of matching the second and third-order cumulant functions within moderate time has been achieved.

5.3 Example 3

In this example, let the physical domain $D$ be a hemispherical shell in Fig.7

where the inner radius $R_{\mathrm{i}}=1$ and the outer radius $R_{\mathrm{o}}=1.2$. The knot vectors are $\Xi_{1}$ = (0,0,1,1), $\Xi_{2}$ = (0,0,0,1,1,1) and $\Xi_{3}$ = (0,0,0,1,1,1), respectively. The control points and weights are listed in Table 4.

The covariance function in physical domain is defined as Eq.(50) where $\sigma^{2}=1$, $b=1$ and $L=1$.

First, we decomposed the covariance function in the parametric space by using Algorithm 1 with $tol=1.0000\times 10^{-5}$ and each $m_{k}=800$ in Algorithm 2, and computed the global relative errors $\varepsilon_{\mathrm{g}2}=3.3731\times 10^{-6}(N=1000)$ and ranks $\bm{r}_{2}=(4,21,60,21,4)$. Then, we computed the eigenpairs in each direction with $\varepsilon_{\mathrm{gf}2}$ = $1.3170\times 10^{-2}\ (N=1000)$. Next, by setting $tol=10^{-4}$ and each $m_{k}=800$ in Algorithm 2, we computed the TT decomposition of the third-order cumulant function and obtained $\varepsilon_{\mathrm{g}3}=8.8170\times 10^{-4}(N=1000)$ and ranks $\bm{r}_{3}=(3,20,60,110,164,56,35,8)$. After HOSVD of $\bm{\breve{C}}_{3}$ ($tol_{3}=9.9990\times 10^{-1}$), the number of latent factors are reduced from 60 to 17 $\varepsilon_{\mathrm{gf}3}=6.5547\times 10^{-3},N=1000$. The time cost (seconds) of each quantity is listed in Table 5.

The results above indicate that the goal of matching the second and third-order cumulant functions within moderate time has been achieved.

6 Conclusions

In this work, we developed novel theoretical and algorithm frameworks for extending the unidimensional Karhunen-Loève expansion to represent multidimensional random fields using higher-order cumulant functions. The algorithm framework can reveal the features of a random field in different directions in a highly automatic way due to the marriage of a rank-revealing tensor train decomposition and the Chebfun paradigm of continuous computation. Differences from some existing methods were also discussed. Numerical experiments indicate that the proposed algorithm framework is able to overcome the current challenges of both computing the modes and representing the cumulants of the latent factors. The efficiency is moderately dependent on the parametric dimensionality, and is generally more than an order of magnitude higher than that of FEM.

Finally, it is worth mentioning that the practical efficiency of an algorithm is also strongly affected by the computer languages. Computation of the third-order cumulant tensor of the latent factors contributes a major part of the time cost in each example. This phenomenon is mainly caused by the low efficiency of MATLAB on which the current Chebfun toolbox heavily depends. Hence, reconstructing the toolbox with a compiled language will also be beneficial for high-performance computing.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grant No.11572106) the National Key Research and Development Program of China (National Key Project No.2017YFC0703506), which are gratefully acknowledged by the authors.