Tensor train-Karhunen-Loève expansion for continuous-indexed random fields using higher-order cumulant functions

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 Karhunen1947 ; Loeve1948

(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

Betz2014 . Other procedures such as spectral element projection Oliveira2014 , multilevel finite element method Xie2015 , two-dimensional Haar wavelet Galerkin method Azevedo2018 and high-order polynomial-based Ritz-Galerkin method Zhang2019 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 Rahman2018 ; Li2018 ; Jahanbin2019 . 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. Ghosh2014 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. Zentner2016 and Guo et al. Guo2016 independently proposed a hierarchical orthogonal decomposition method to split the temporal and spatial modes. Similar idea was also shown in Zheng2017 , however, this work did not really achieve variable separation since the modes still contain multiple coordinates. The stepwise covariance matrix decomposition method Li2019 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. Phoon2005 and Dai et al. Dai2019 proposed similar iterative procedures to update the marginal distribution of each latent factor. Kim et al. Kim2015a

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

Ganapathysubramanian2008 ; Ma2011 ; Poirion2014 ; Olivier2016 . 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 Lu2013 . 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 Gorodetsky2019 .

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 Morton2009

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 $(Ω, F, P)$

be a probability space where

Ω

is a sample space,

F

is a

σ

-field on

Ω

and

P

is a probability measure. Let

(E, E)

is the measurable space on the admissible set

E

where each

f \in E

is defined on a bounded domain

D \subseteq R^{d}, d = 1, 2, 3

. A random field is defined as a measurable mapping

ω : (Ω, F, P) \to (E, E)

. To explicitly represent this abstract mapping, an auxiliary measurable space

(R^{n}, B^{n}), n \in N^{+}

is needed and the corresponding two mappings

Θ : (Ω, F, P) \to (R^{n}, B^{n})

and

H : (R^{n}, B^{n}) \to (E, E)

are to be found such that

ω = H \circ Θ

. The task is to capture the major part of information in

ω

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)


$x (ξ)$	$= n_{1} \sum i = 1 R_{i}^{p} (ξ) B_{i},$	$ξ \in [0, 1]$	$(c u r v e s)$	(1a)
$x (ξ, η)$	$= n_{1} \sum i = 1 n_{2} \sum j = 1 R_{i, j}^{p, q} (ξ, η) B_{i, j},$	$(ξ, η) \in [0, 1]^{2}$	$(s u r f a c e s)$	(1b)
$x (ξ, η, ζ)$	$= n_{1} \sum i = 1 n_{2} \sum j = 1 n_{3} \sum k = 1 R_{i, j, k}^{p, q, r} (ξ, η, ζ) B_{i, j, k},$	$(ξ, η, ζ) \in [0, 1]^{3}$	$(s o l i d s)$	(1c)

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),


$R_{i}^{p} (ξ)$	$= \frac{N_{i}^{p} (ξ) w_{i}}{n_{1} \sum i = 1 N_{i}^{p} (ξ) w_{i}}$	(2a)
$R_{i, j}^{p, q} (ξ, η)$	$= \frac{N_{i}^{p} (ξ) N_{j}^{q} (η) w_{i, j}}{n_{1} \sum i = 1 n_{2} \sum j = 1 N_{i}^{p} (ξ) N_{j}^{q} (η) w_{i, j}}$	(2b)
$R_{i, j, k}^{p, q, r} (ξ, η, ζ)$	$= \frac{N_{i}^{p} (ξ) N_{j}^{q} (η) N_{k}^{r} (ζ) w_{i, j, k}}{n_{1} \sum i = 1 n_{2} \sum j = 1 n_{3} \sum k = 1 N_{i}^{p} (ξ) N_{j}^{q} (η) N_{k}^{r} (ζ) w_{i, j, k}}$	(2c)

$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

Ξ = (ξ_{1}, \dots, ξ_{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)


$N_{i}^{0}$	$= {\begin{matrix} 1 ξ \in [ξ_{i}, ξ_{i + 1}) 0 o t h e r w i s e \end{matrix}$	(3a)
$N_{i}^{p} (ξ)$	$= \frac{ξ - ξ_{i}}{ξ_{i + p} - ξ_{i}} N_{i}^{p - 1} (ξ) + \frac{ξ_{i + p + 1} - ξ}{ξ_{i + p + 1} - ξ_{i + 1}} N_{i + 1}^{p - 1} (ξ)$	(3b)

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).

x (ξ) = \sum I R_{I}^{p} (ξ) B_{I}

(4)

where $x \in D \subseteq R^{d}, ξ \in [0, 1]^{m}$ and $m ⩽ d$ .

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

After representing the physical coordinates with parametric ones, the original random field $ω (x, θ)$ is transformed to the parametric space. Denote $α (ξ, θ) = ω (x (ξ), θ)$ and the corresponding covariance function $~ C (ξ, ξ^{'}) = C (x (ξ), x (ξ^{'}))$ . We seek to represent $α (ξ, θ)$ with a generalized K-L expansion by improving the hierarchical SVD method in Zentner2016 . Taking the case of $d$ =2 as an example, the first step is to split the $ξ$ -modes $f (ξ)$ = $(f_{i} (ξ))_{1 \times n_{1}}$ from $α (ξ, η; θ)$ by computing the eigenfunctions of the kernel in Eq.(5).

{~ C}_{f} (ξ, ξ^{'}) = \int_{0}^{1} ~ C (ξ, η; ξ^{'}, η) d η

(5)

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

α (ξ, η; θ) \approx f (ξ) α_{- ξ} (η, θ)

(6)

where

α_{- ξ} (η, θ) = \int_{0}^{1} f (ξ)^{T} α (ξ, η; θ) d ξ

(7)

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

{~ C}_{- ξ} (η, η^{'})_{n_{1} \times n_{1}} = \int_{0}^{1} \int_{0}^{1} f (ξ)^{T} ~ C (ξ, η; ξ^{'}, η^{'}) f (ξ^{'}) d ξ d ξ^{'}

(8)

Then, rather than decomposing each component of $α_{- ξ} (η, θ)$ separately as in Zentner2016 , $α_{- ξ} (η, θ)$ is treated as whole object by regarding the index as a new coordinate. Next, the $η$ -modes of $α (ξ, η; θ)$ (also the $η$ -modes of $α_{- ξ} (η, θ)$ ), denoting $g (η)_{n_{1} \times n_{2}}$ , are split by computing the eigenpairs of ${~ C}_{- ξ} (η, η^{'})$ , see Eq.(9).

\int_{0}^{1} {~ C}_{- ξ} (η, η^{'}) g_{:, j} (η^{'}) d η^{'} = μ_{j} g_{:, j} (η)

(9)

Finally, we get the tensor train decomposition of $α (ξ, η; θ)$ as Eq.(10)

α (ξ, η; θ) \approx f (ξ) g (η) γ (θ)

(10)

where

γ (θ) = \int_{0}^{1} \int_{0}^{1} g (η)^{T} f (ξ)^{T} α (ξ, η; θ) d ξ d η

(11)

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


$⟨ f_{i}, f_{j} ⟩$	$= \int_{0}^{1} f_{i} (ξ) f_{j} (ξ) d ξ = δ_{i j}$	(12a)
$⟨ g_{:, i}, g_{:, j} ⟩$	$= \int_{0}^{1} g_{:, i} (η)^{T} g_{:, j} (η) d η = δ_{i j}$	(12b)
$⟨ γ_{i}, γ_{j} ⟩$	$= E [h_{i} (θ) h_{j} (θ)] = μ_{i} δ_{i j}$	(12c)

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

α (ξ, η, ζ; θ) \approx f (ξ)_{1 \times n_{1}} g (η)_{n_{1} \times n_{2}} h (ζ)_{n_{2} \times n_{3}} γ (θ)_{n_{3} \times 1}

(13)

where each column of $f (ξ)$ is an eigenfunction of the covariance kernel ${~ C}_{f} (ξ, ξ^{'})$ in Eq.(14),

{~ C}_{f} (η, η^{'}) = \int_{0}^{1} \int_{0}^{1} ~ C (ξ, η, ζ; ξ^{'}, η, ζ) d η d ζ

(14)

each column of $g (η)$ is an eigenfunction of ${~ C}_{- ξ} (η, η^{'})$ in Eq.(15)

{~ C}_{- ξ} (η, η^{'})_{n_{1} \times n_{1}} = \int_{0}^{1} \int_{0}^{1} f (ξ)^{T} \int_{0}^{1} ~ C (ξ, η, ζ; ξ^{'}, η^{'}, ζ) d ζ f (ξ^{'}) d ξ d ξ^{'}

(15)

and each column of $h (ζ)$ is an eigenfunction of ${~ C}_{- ξ η} (ζ, ζ^{'})$ in Eq.(16).

\begin{matrix} {~ C}_{- ξ η} (ζ, ζ^{'})_{n_{2} \times n_{2}} = & \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} g (η)^{T} f (ξ)^{T} ~ C (ξ, η, ζ; ξ^{'}, η^{'}, ζ^{'}) f (ξ^{'}) g (η^{'})) d ξ d η d ξ^{'} d η^{'} \end{matrix}

(16)

(17)


$⟨ f_{i}, f_{j} ⟩$	$= \int_{0}^{1} f_{i} (ξ) f_{j} (ξ) d ξ = δ_{i j}$	(18a)
$⟨ g_{:, i}, g_{:, j} ⟩$	$= \int_{0}^{1} g_{:, i} (η)^{T} g_{:, j} (η) d η = δ_{i j}$	(18b)
$⟨ h_{:, i}, h_{:, j} ⟩$	$= \int_{0}^{1} h_{:, i} (η)^{T} h_{:, j} (η) d η = δ_{i j}$	(18c)
$⟨ γ_{i}, γ_{j} ⟩$	$= E [γ_{i} (θ) γ_{j} (θ)] = ν_{i} δ_{i j}$	(18d)

where $ν_{i}$ is the $i$ th eigenvalue of the covariance kernel ${~ C}_{- ξ η} (ζ, ζ^{'})$ 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)} (ξ, η)}_{k = 1}^{\infty}$ and $Λ$ = ${λ_{k}}_{k = 1}^{\infty}$ as the sets of 2D modes and eigenvalues in the trivial PCA and $E_{2}$ = ${F_{k}^{(2)} (ξ, η) = f (ξ) g_{:, k} (η)}_{k = 1}^{\infty}$ and $M$ = ${μ_{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 $Λ$ = $M$ .

First, we prove $E_{1} \subseteq E_{2}$ and $Λ \in M$ . For each $k = 1, 2, \dots$ , by projecting $F_{k} (ξ, η)$ on $f (ξ)$ ,

F_{k} (ξ, η) = f (ξ) {^g}_{:, k} (η)

we derive

\begin{matrix} \int_{0}^{1} \int_{0}^{1} ~ C (ξ, η; ξ^{'}, η^{'}) f (ξ^{'}) {^g}_{:, k} (η^{'}) d ξ^{'} d η^{'} = λ_{k} f (ξ) {^g}_{:, k} (η) \Rightarrow & \int_{0}^{1} f (ξ)^{T} \int_{0}^{1} \int_{0}^{1} ~ C (ξ, η; ξ^{'}, η^{'}) f (ξ^{'}) {^g}_{:, k} (η^{'}) d ξ^{'} d η^{'} d ξ = λ_{k} \int_{0}^{1} f (ξ)^{T} f (ξ) d ξ {^g}_{:, k} (η) \Rightarrow & \int_{0}^{1} {~ C}_{- ξ} (η, η^{'}) {^g}_{:, k} (η^{'}) d η^{'} = λ_{k} {^g}_{:, k} (η) \Rightarrow & {^g}_{:, k} (η) i s a c o l u m n o f g (η), λ_{k} \in M \Rightarrow & E_{1} \subseteq E_{2}, Λ \in M . \end{matrix}

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

\begin{matrix} ε_{k} (ξ, η) = \int_{0}^{1} \int_{0}^{1} ~ C (ξ, η; ξ^{'}, η^{'}) f (ξ^{'}) g_{:, k} (η^{'}) d ξ^{'} d η^{'} - μ_{k} f (ξ) g_{:, k} (η) \Rightarrow & \int_{0}^{1} f (ξ)^{T} ε_{k} (ξ, η) d ξ = \int_{0}^{1} {~ C}_{- ξ} (η, η^{'}) g_{:, k} (η^{'}) d η^{'} - μ_{k} g_{:, k} (η) = 0 \Rightarrow & ε_{k} (ξ, η) \equiv 0 \Rightarrow & \int_{0}^{1} \int_{0}^{1} ~ C (ξ, η; ξ^{'}, η^{'}) f (ξ^{'}) g_{:, k} (η^{'}) d ξ^{'} d η^{'} = μ_{k} f (ξ) g_{:, k} (η) \Rightarrow & E_{2} \subseteq E_{1}, M \in Λ . \end{matrix}

Finally, we get $E_{1}$ = $E_{2}$ and $Λ$ = $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 $Ξ = (ξ_{1}, \dots, ξ_{n})$ in parametric space and arbitrary $u = (u_{1}, \dots, u_{n}) \in R^{n}$

, according to the A-type Gram-Charlier series, the corresponding marginal probability density function (PDF) can be represented as Eq.(

19)

\begin{matrix} f_{Ξ} (u; κ) = & f_{G} (u) ⎡ ⎣ 1 + \frac{1}{3!} \sum i, j, k κ^{i, j, k} h_{i j k} (u) + \frac{1}{4!} κ^{i, j, k, l} h_{i j k l} (u) + \frac{1}{5!} κ^{i, j, k, l, m} h_{i j k l m} (u) + \frac{1}{6!} (κ^{i, j, k, l, m, n} + 10 κ^{i, j, k} κ^{l, m, n}) h_{i j k l m n} (u) + \dots] \end{matrix}

(19)

where $f_{G}$

is the PDF of the Gaussian distribution with mean zero and covariance function

~ C (ξ, ξ^{'})

h

with subscripts are Hermite polynomials and

κ

is the family of cumulant functions consists of elements

κ^{i, j, k}

{~ C}_{3} (ξ_{i}, ξ_{j}, ξ_{k})

(the third-order cumulant function) and so on. Meanwhile,

κ

is defined by the family of marginal PDFs by definition. Thus,

κ

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 (

⩾

3) cumulant functions into the representation of

α (ξ; θ)

, 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).

\begin{matrix} {˘ C}_{3} (i, j, k) = & \int_{[0, 1]^{6}} {~ C}_{3} (ξ, η; ξ^{'}, η^{'}; ξ^{''}, η^{''}) (g (η)_{:, i}^{T} f (ξ)^{T}) (g (η^{'})_{:, j}^{T} f (ξ^{'})^{T}) (g (η^{''})_{:, k}^{T} f (ξ^{''})^{T}) d ξ d η d ξ^{'} d η^{'} d ξ^{''} d η^{''} \end{matrix}

(20)

This relationship can be directly extended to the cases of other values of $d$ and cumulant orders. Unfortunately, ${˘ C}_{3} (i, j, k)$ is generally not orthogonal decomposable Robeva2016 , 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 $ω (x; θ)$ 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 Piegl2012 .

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 $α (ξ; θ)$ . 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 Savostyanov2014 , we propose a tensor train decomposition algorithm to obtain a low-rank separable approximation to a $K$ th ( $K ⩾ 2$ ) cumulant function ${~ C}_{K} (ξ_{1}, \dots, ξ_{K})$ .

For the sake of simplicity, ${~ C}_{K} (ξ_{1}, \dots, ξ_{K})$ is redefined as an auxiliary function $G : [0, 1]^{a} \to R (a = m K)$ where $u = (u_{1}, \dots, 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)

\begin{matrix} G (u_{1}, \dots, u_{a}) \approx & G (u_{1}, I^{> 1}) [G (I^{⩽ 1}, I^{> 1})]^{- 1} G (I^{⩽ 1}, u_{2}, I^{> 2}) [G (I^{⩽ 2}, I^{> 2})]^{- 1} \dots G (I^{⩽ a - 1}, u_{a}) \end{matrix}

(21)

where $I^{⩽ k}$ and $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).

{I^{⩽ k *}, I^{> k *}} = arg max \begin{matrix} I^{⩽ k}, I^{> k} \end{matrix} | det G (I^{⩽ k}, I^{> k}) |

(22)

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 ${I^{⩽ k}, I^{> k}}_{k = 1}^{a - 1}$ and derive an adaptive tensor train decomposition of a cumulant function, see Algorithm 1.

1:Auxiliary function

G : [0, 1]^{a} \to R (a = m K)

; maximum number of sweeps

m a x s w p

; stopping tolerance

t o l

2:Tensor train decomposition of

G

, denoting

G_{T T}

; interpolation sets

I

; a matrix

e r r d m

containing iteration errors of each dimension during each iteration

u^{(0)} \leftarrow arg {max}_{\begin{matrix} u \end{matrix}} | G (u) |

from

m_{0}

quasi-random samples;

▹

Find the initial pivot

4:for all

k = 1 : a - 1

{I^{⩽ k}, I^{> k}} \leftarrow {{u^{(0) ⩽ k}}, {u^{(0) > k}}}

;

▹

Initialize the interpolation sets

6:end for

I \leftarrow {I^{⩽ k}, I^{> k}}_{k = 1}^{a - 1}

;

8:for all

k = 1 : a - 1

e r r d (k) = 1

;

▹

Initialize the iteration error in each dimension

10:end for

11:

S \leftarrow

▹

Initialize the sweep number

12:while

S < m a x s w p

13:

S \leftarrow S + 1

;

▹

Begin left-to-right sweep

14: for

k = 1 : a - 1

15: if

e r r d (k) < t o l

then

16: continue

17: else

18:

[I, e r r m a x] \leftarrow

ISE(

k

I

G

a - 1

)

▹

use Algorithm 2 to expand interpolation sets

19:

e r r d (k) \leftarrow e r r m a x

;

▹

update iteration errors

20: end if

21: end for

22: if

S = 1

then

▹

Record iteration errors of all dimensions during each iteration

23:

e r r d m \leftarrow e r r d

24: else

25:

e r r d m \leftarrow [e r r d m, e r r d]

26: end if

27: if

max (e r r d) < t o l

then

28: break

29: end if

30:

S \leftarrow S + 1

;

▹

Begin right-to-left sweep

31: for

k = a : - 1 : 2

32: if

e r r d (k - 1) < t o l

then

33: continue

34: else

35:

[I, e r r m a x] \leftarrow

ISE(

k - 1

I

G

a - 1

)

▹

use Algorithm 2 to expand interpolation sets

36:

e r r d (k - 1) \leftarrow e r r m a x

;

▹

update iteration errors

37: end if

38: end for

39:

e r r d m \leftarrow [e r r d m, e r r d]

;

▹

Record iteration errors of all dimensions during each iteration

Algorithm 1 Adaptive tensor train decomposition of a

K

th cumulant function on a

m

-dimensional parametric space

40: if

max (e r r d) < t o l

then

41: break

42: end if

43:end while

44:for all

k = 1 : a

▹

Output tensor train decomposition by adaptive reconstruction of each fiber

45: if

k = 1

then

46:

^G (u_{1})_{1 \times r_{1}} \leftarrow

chebmatrix(

G (u_{1}, I^{> 1})

)

[G (I^{⩽ 1}, I^{> 1})]^{- 1}

;

47: else if

k < d

then

48:

^G (u_{k})_{r_{k - 1} \times r_{k}} \leftarrow

chebmatrix(

G (I^{⩽ k - 1}, u_{k}, I^{> k})

)

[G (I^{⩽ k}, I^{> k})]^{- 1}

;

49: else

50:

^G (u_{a})_{r_{a - 1} \times 1} \leftarrow

chebmatrix(

G (I^{⩽ a - 1}, u_{a})

);

51: end if

52:end for

53:

G_{T T} (u) = {^G}_{1} (u_{1}) {^G}_{2} (u_{2}) \dots {^G}_{a} (u_{a})

(23)

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, \dots, a - 1}$ , $I^{⩽ k + 1} \subseteq I^{⩽ k} \times [0, 1]$ , $I^{> k} \subseteq [0, 1] \times I^{> k + 1}$ . Then, each fiber in Eq.(21) is adaptively represented in the chebfun format Driscoll2014 . Next, fibers of each two neighboring dimensions form a chebmatrix as in lines 46, 48, 50.

1:current dimension

k

; current interpolation sets

I

; auxiliary function

G

; number of variables of

G

, denoting

a

;

m_{k}

for solving Eq.(24); iteration error tolerance

t o l

2:Expanded interpolation sets

I

; iteration error

e r r m a x

3:if

k

=1 then

4: Join the supports of the matrix-valued function

G (u_{k}, u_{k + 1}, I^{> k + 1})

in the second dimension and define a new function

G_{J s t} : [0, 1] \times [0, r_{k + 1}] \to R

;

5: Redefine the matrix-valued function

G (u_{k}, I^{> k})

as a new function

G_{J s} : [0, 1] \to R^{1 \times r_{k}}

;

6: Join the supports of each row of the matrix-valued function

G (I^{⩽ k}, u_{k + 1}, I^{> k + 1})

and define a new function

G_{J t} : [0, r_{k + 1}] \to R^{r_{k} \times 1}

;

7:else if

k < d - 1

then

8: Join the supports of the matrix-valued function

G (I^{⩽ k - 1}, u_{k}, u_{k + 1}, I^{> k + 1})

and define a new function

G_{J s t} : [0, r_{k - 1}] \times [0, r_{k + 1}] \to R

;

9: Join the supports of each column of the matrix-valued function

G (I^{⩽ k - 1}, u_{k}, I^{> k})

and define a new function

G_{J s} : [0, r_{k - 1}] \to R^{1 \times r_{k}}

;

10: Join the supports of each row of the matrix-valued function

G (I^{⩽ k}, u_{k + 1}, I^{> k + 1})

and define a new function

G_{J t} : [0, r_{k + 1}] \to R^{r_{k} \times 1}

;

11:else

12: Join the supports of the matrix-valued function

G (I^{⩽ k - 1}, u_{k}, u_{k + 1})

in the first dimension and define a new function

G_{J s t} : [0, r_{k - 1}] \times [0, 1] \to R

;

13: Join the supports of each column of the matrix-valued function

G (I^{⩽ k - 1}, u_{k}, I^{> k})

and define a new function

G_{J s} : [0, r_{k - 1}] \to R^{1 \times r_{k}}

;

14: Redefine the matrix-valued function

G (I^{⩽ k}, u_{k + 1})

as a new function

G_{J t} : [0, 1] \to R^{r_{k} \times 1}

;

15:end if

16:Compute

[G (I^{⩽ k}, I^{> k})]^{- 1}

;

17:Solve Eq.(24)

(s^{*}, t^{*}) = arg max \begin{matrix} s, t \end{matrix} | G_{J s t} (s, t) - G_{J s} (s) [G (I^{⩽ k}, I^{> k})]^{- 1} G_{J t} (t) |

(24)

with

m_{k}

quasi-random samples;

18:

e r r m a x \leftarrow | G_{J s t} (s^{*}, t^{*}) - G_{J s} (s^{*}) [G (I^{⩽ k}, I^{> k})]^{- 1} G_{J t} (t^{*}) |

;

▹

Iteration error

19:if

e r r m a x < t o l

then

20: break

21:end if

22:

i^{*} \leftarrow [s^{*}]

j^{*} \leftarrow [t^{*}]

▹

[x]

denotes the largest integer no larger than

x

23:

u^{(S)}

is obtained with

u^{(S) ⩽ k} \leftarrow [I^{⩽ k - 1} (i^{*}, :), s^{*} - i^{*}]

and

u^{(S) > k} \leftarrow [t^{*} - j^{*}, I^{> k + 1} (j^{*}, :)]

;

▹

The new pivot

24:

I^{⩽ k} \leftarrow I^{⩽ k} \cup {u^{(S) ⩽ k}}

I^{> k} \leftarrow I^{> k} \cup {u^{(S) > k}}

;

▹

The expanded interpolation set of dimension

k

Algorithm 2

I S E

: Interpolation set expansion algorithm

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

1:Auxiliary function

G : [0, 1]^{2} \to R

; maximum number of sweeps

m a x s w p

; stopping tolerance

t o l

2:Tensor train decomposition of

G

, denoting

G_{T T}

; interpolation sets

I

; a matrix

e r r d m

containing iteration errors during each iteration

u^{(0)} \leftarrow arg {max}_{\begin{matrix} u \end{matrix}} | G (u) |

from

m_{0}

quasi-random samples;

▹

Find the initial pivot

{I^{⩽ 1}, I^{> 1}} \leftarrow {{u^{(0) ⩽ 1}}, {u^{(0) > 1}}}

;

▹

Initialize the interpolation sets

I \leftarrow {I^{⩽ 1}, I^{> 1}}

;

e r r m a x = 1

;

▹

Initialize the iteration error

7:for

S = 1 : m a x s w p

8: Compute

[G (I^{⩽ 1}, I^{> 1})]^{- 1}

;

9: Solve Eq.(25)

(s^{*}, t^{*}) = arg max \begin{matrix} s, t \end{matrix} | G (s, t) - G (s, I^{> 1}) [G (I^{⩽ 1}, I^{> 1})]^{- 1} G (I^{⩽ k}, t) |

(25)

with

m_{1}

quasi-random samples;

10:

e r r m a x \leftarrow | G (s^{*}, t^{*}) - G (s^{*}, I^{> 1}) [G (I^{⩽ 1}, I^{> 1})]^{- 1} G (I^{⩽ k}, t^{*}) |

;

▹

Iteration error

11: if

S = 1

then

▹

Record iteration errors of all dimensions during each iteration

12:

e r r d m \leftarrow e r r m a x

13: else

14:

e r r d m \leftarrow [e r r d m, e r r m a x]

15: end if

16: if

e r r m a x < t o l

then

17: break

18: end if

19:

u^{(S)}

is obtained with

u^{(S) ⩽ 1} \leftarrow s^{*}

and

u^{(S) > k} \leftarrow t^{*}

;

▹

The new pivot

20:

I^{⩽ 1} \leftarrow I^{⩽ 1} \cup {u^{(S) ⩽ 1}}

I^{> 1} \leftarrow I^{> 1} \cup {u^{(S) > 1}}

;

▹

The expanded interpolation set

21:end for

22:

^G (u_{1})_{1 \times r_{1}} \leftarrow

chebmatrix(

G (u_{1}, I^{> 1})

)

[G (I^{⩽ 1}, I^{> 1})]^{- 1}

;

▹

Output tensor train decomposition by adaptive reconstruction of each fiber

23:

^G (u_{2})_{r_{1} \times 1} \leftarrow

chebmatrix(

G (I^{⩽ 1}, u_{2})

);

24:

Tensor train-Karhunen-Loève expansion for continuous-indexed random fields using higher-order cumulant functions

Authors

A Galerkin Isogeometric Method for Karhunen-Loeve Approximation of Random Fields

A matrix-free isogeometric Galerkin method for Karhunen-Loève approximation of random fields using tensor product splines, tensor contraction and interpolation based quadratur

Certified and fast computations with shallow covariance kernels

Forming A Random Field via Stochastic Cliques: From Random Graphs to Fully Connected Random Fields

On the Numerical Approximation of the Karhunen-Loève Expansion for Lognormal Random Fields

A Higher Order Unscented Transform

Modelling Distributed Shape Priors by Gibbs Random Fields of Second Order

1 Introduction

2 Theoretical framework of the proposed method

2.1 Space transformation

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

Theorem 1.

Proof.

2.3 Representing higher-order cumulants of latent factors

3 Computational procedure

3.1 Space transformation

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

3.2.1 Tensor train decomposition of cumulant functions

Tensor train-Karhunen-Loève expansion for continuous-indexed random fields using higher-order cumulant functions

Authors

A Galerkin Isogeometric Method for Karhunen-Loeve Approximation of Random Fields

A matrix-free isogeometric Galerkin method for Karhunen-Loève approximation of random fields using tensor product splines, tensor contraction and interpolation based quadratur

Certified and fast computations with shallow covariance kernels

Forming A Random Field via Stochastic Cliques: From Random Graphs to Fully Connected Random Fields

On the Numerical Approximation of the Karhunen-Loève Expansion for Lognormal Random Fields

A Higher Order Unscented Transform

Modelling Distributed Shape Priors by Gibbs Random Fields of Second Order

This week in AI

1 Introduction

2 Theoretical framework of the proposed method

2.1 Space transformation

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

Theorem 1.

Proof.

2.3 Representing higher-order cumulants of latent factors

3 Computational procedure

3.1 Space transformation

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

3.2.1 Tensor train decomposition of cumulant functions