Higher-order principal component analysis for the approximation of tensors in tree-based low-rank formats

TL;DR

This paper introduces a hierarchical algorithm for tensor approximation in tree-based formats using only point evaluations, combining higher-order SVD and empirical PCA to achieve near-optimal accuracy efficiently.

Contribution

It presents a novel algorithm that constructs tensor approximations from point evaluations, with theoretical guarantees for near-optimality and applicability to high-dimensional functions.

Findings

Algorithm achieves approximation with a number of evaluations comparable to storage complexity.

Provides quasi-optimal or prescribed error approximations under certain assumptions.

Numerical examples demonstrate effectiveness for high-dimensional tensor formats.

Abstract

This paper is concerned with the approximation of tensors using tree-based tensor formats, which are tensor networks whose graphs are dimension partition trees. We consider Hilbert tensor spaces of multivariate functions defined on a product set equipped with a probability measure. This includes the case of multidimensional arrays corresponding to finite product sets. We propose and analyse an algorithm for the construction of an approximation using only point evaluations of a multivariate function, or evaluations of some entries of a multidimensional array. The algorithm is a variant of higher-order singular value decomposition which constructs a hierarchy of subspaces associated with the different nodes of the tree and a corresponding hierarchy of interpolation operators. Optimal subspaces are estimated using empirical principal component analysis of interpolations of partial random evaluations of the function. The algorithm is able to provide an approximation in any tree-based format with either a prescribed rank or a prescribed relative error, with a number of evaluations of the order of the storage complexity of the approximation format. Under some assumptions on the estimation of principal components, we prove that the algorithm provides either a quasi-optimal approximation with a given rank, or an approximation satisfying the prescribed relative error, up to constants depending on the tree and the properties of interpolation operators. The analysis takes into account the discretization errors for the approximation of infinite-dimensional tensors. Several numerical examples illustrate the main results and the behavior of the algorithm for the approximation of high-dimensional functions using hierarchical Tucker or tensor train tensor formats, and the approximation of univariate functions using tensorization.

Equations

$$v(x_1,\dots ,x_d)=\sum _{k=1}^r{v}_k^1(x_1)\dots v_k^d(x_d)$$

$$v(x_1,\dots ,x_d)=\sum _{k=1}^{{\mathrm{rank}}_{\alpha }(v)}{v}_k^{\alpha }(x_{\alpha })v_k^{{\alpha }^c}(x_{{\alpha }^c})$$

$${\mathcal{T}}_r^A=\{v:{\mathrm{rank}}_{\alpha }(v)\le r_{\alpha },\alpha \in A\}=\bigcap _{\alpha \in A}{\mathcal{T}}_{r_{\alpha }}^{\{\alpha \}},$$

$$\Vert u-u^{\star }\Vert \le c\underset{v\in {\mathcal{T}}_r^A}{\min }\Vert u-v\Vert +\rho ,$$

$$\Vert u-u^{\star }\Vert \le \tilde{c}ϵ\Vert u\Vert$$

$$\Vert v-P_{U}v{\Vert }_V=\underset{u\in U}{\min }\Vert v-u{\Vert }_V,$$

$$\dim (W)=\dim (U),\text{ and}$$

$$\{u\in U:⟨w,u⟩=0\text{ for all }w\in W\}=\{0\},$$

$$⟨w,P_U^{W}v-v⟩=0,\forall w\in W.$$

$$\alpha =\underset{0\ne u\in U}{\min }\underset{0\ne w\in W}{\max }\frac{⟨w,u⟩}{\Vert u{\Vert }_V\Vert w{\Vert }_W}.$$

$$\Vert P_U^{W}v{\Vert }_V\le {\alpha }^{-1}\underset{0\ne w\in W}{\max }\frac{⟨w,P_U^{W}v⟩}{\Vert w{\Vert }_W}={\alpha }^{-1}\underset{0\ne w\in W}{\max }\frac{⟨w,v⟩}{\Vert w{\Vert }_W}\le {\alpha }^{-1}\Vert v{\Vert }_V,$$

$$\tilde{P}P=P.$$

$$\tilde{P}P=P\tilde{P}=P.$$

$$\Vert P_U^{W}v-P_{U}v{\Vert }_V\le \Vert P_U^W-P_U{\Vert }_{\tilde{U}\to V}\Vert v-P_{U}v{\Vert }_V,$$

$$\Vert P_U^W-P_U{\Vert }_{\tilde{U}\to V}=\Vert P_U^W{\Vert }_{(id-P_U)\tilde{U}\to V}\le \Vert P_U^W{\Vert }_{\tilde{U}\to V}.$$

$$\Vert v-P_U^{W}v{\Vert }_V^2\le (1+\Vert P_U^W-P_U{\Vert }_{\tilde{U}\to V}^2)\Vert v-P_{U}v{\Vert }_V^2.$$

$$⟨{\delta }_x,I_{U}v-v⟩=I_{U}v(x)-v(x)=0\forall x\in \Gamma .$$

$$I_U{I}_{\tilde{U}}=I_{\tilde{U}}{I}_U=I_U.$$

$$|{\phi }_{i_1}(x^1)|=\underset{x\in {\Gamma }_{\star }}{\max }\underset{1\le i\le n}{\max }|{\phi }_i(x)|.$$

$$|{\psi }_{i_{k+1}}^{(k)}(x^{k+1})|=\underset{x\in {\Gamma }_{\star }}{\max }\underset{1\le i\le n}{\max }|{\psi }_i^{(k)}(x)|.$$

$$(u,v-Q_{U}v{)}_{\Gamma }=\sum _{x\in \Gamma }u(x)⟨{\delta }_x,v-Q_{U}v⟩=0\forall u\in U,$$

$$Q_U{Q}_{\tilde{U}}=Q_{\tilde{U}}{Q}_U=Q_U.$$

$$(v^1\otimes \dots \otimes v^d,w^1\otimes \dots w^d)=(v^1,w^1{)}_{{\mathcal{H}}_1}\dots (v^d,w^d{)}_{{\mathcal{H}}_d},$$

$$\mathcal{H}={\mathcal{H}}_1\otimes \dots \otimes {\mathcal{H}}_d=\underset{\nu =1}{\overset{d}{⨂}}{\mathcal{H}}_{\nu }.$$

$$\Vert \underset{\nu \in \alpha }{⨂}{v}^{\nu }{\Vert }_{{\mathcal{H}}_{\alpha }}=\prod _{\nu \in \alpha }\Vert v^{\nu }{\Vert }_{{\mathcal{H}}_{\nu }}$$

$${\mathcal{M}}_{\alpha }(v)\in {\mathcal{H}}_{\alpha }\otimes {\mathcal{H}}_{{\alpha }^c},$$

$${\mathcal{M}}_{\alpha }(v)=\sum _{k=1}^{{\mathrm{rank}}_{\alpha }(v)}{v}_k^{\alpha }\otimes v_k^{{\alpha }^c},\text{or equivalently}v(x)=\sum _{k=1}^{{\mathrm{rank}}_{\alpha }(v)}{v}_k^{\alpha }(x_{\alpha })v_k^{{\alpha }^c}(x_{{\alpha }^c}),$$

$$v\in U_{\alpha }^{min}(v)\otimes {\mathcal{H}}_{{\alpha }^c},$$

$$U_{\alpha }^{min}(v)\subset \underset{\beta \in S(\alpha )}{⨂}{U}_{\beta }^{min}(v).$$

$$v\in \underset{\beta \in S(D)}{⨂}{U}_{\beta }^{min}(v).$$

$$(A_{\alpha }\otimes {id}_{{\alpha }^c})(v^{\alpha }\otimes v^{{\alpha }^c})=(A_{\alpha }{v}^{\alpha })\otimes v^{{\alpha }^c},$$

Full text

Keywords: high-dimensional approximation, tree-based tensor formats, deep tensor networks, higher-order singular value decomposition, higher-order principal component analysis, interpolation.

2010 AMS Subject Classifications: 15A69, 41A05, 41A63, 65D15, 65J99

1 Introduction

The approximation of high-dimensional functions is one of the most challenging tasks in computational science. Such high-dimensional problems arise in many domains of physics, chemistry, biology or finance, where the functions are the solutions of high-dimensional partial differential equations (PDEs). Such problems also typically arise in statistics or machine learning, for the estimation of high-dimensional probability density functions, or the approximation of the relation between a certain random variable and some predictive variables, the typical task of supervised learning. The approximation of high-dimensional functions is also required in optimization or uncertainty quantification problems, where the functions represent the response of a system (or model) in terms of some parameters. These problems require many evaluations of the functions and are usually intractable when one evaluation requires a specific experimental set-up or one run of a complex numerical code.

The approximation of high-dimensional functions from a limited number of information on the functions requires exploiting low-dimensional structures of functions. This usually call for nonlinear approximation tools [11, 42]. A prominent approach consists of exploiting the sparsity of functions relatively to a basis, a frame, or a more general dictionary of functions [43, 4, 7]. Another approach consists of exploiting low-rank structures of multivariate functions, interpreted as elements of tensor spaces, which is related to notions of sparsity in (uncountably infinite) dictionaries of separable functions. For a multivariate function $v(x_{1},\ldots,x_{d})$ defined on a product set $\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{d}$, which is here identified with a tensor of order $d$, a natural notion of rank is the canonical rank, which is the minimal integer $r$ such that

[TABLE]

for some univariate functions $v^{\nu}_{k}$ defined on $\mathcal{X}_{\nu}$. For $d=2$, this corresponds to the unique notion of rank, which coincides with the matrix rank when the variables take values in finite index sets and $v$ is identified with a matrix. A function with low canonical rank $r$ has a number of parameters which scales only linearly with $r$ and $d$. However, it turns out that this format has several drawbacks when $d>2$ (see, e.g., [10, 21]), which makes it unsuitable for approximation. Then, other notions of rank have been introduced. For a subset of dimensions $\alpha$ in $\{1,\ldots,d\}$, the $\alpha$-rank of a function $v$ is the minimal integer $\mathrm{rank}_{\alpha}(v)$ such that

[TABLE]

for some functions $v^{\alpha}_{k}$ and $v^{\alpha^{c}}_{k}$ of complementary groups of variables $x_{\alpha}=(x_{\nu})_{\nu\in\alpha}\in\mathcal{X}_{\alpha}$ and $x_{\alpha^{c}}=(x_{\nu})_{\nu\in\alpha^{c}}\in\mathcal{X}_{\alpha^{c}}$, with $\alpha^{c}$ the complementary subset of $\alpha$ in $\{1,\ldots,d\}.$ Approximation formats can then be defined by imposing $\alpha$-ranks for a collection of subsets $\alpha$. More precisely, if $A$ is a collection of subsets in $\{1,\ldots,d\}$, we define an approximation format

[TABLE]

where $r=(r_{\alpha})_{\alpha\in A}$ is a tuple of integers. When $A$ is a tree-structured collection of subsets (a subset of a dimension partition tree), $\mathcal{T}_{r}^{A}$ is a tree-based tensor format whose elements admit a hierarchical and data-sparse representation. Tree-based tensor formats are tree tensor networks, i.e. tensor networks with tree-structured graphs [35]. They include the hierarchical Tucker (HT) format [20] and the tensor-train (TT) format [37]. Tree-based formats have many favorable properties that make them favorable for numerical use. As an intersection of subsets of tensors with bounded $\alpha$-rank, $\alpha\in A$, these formats inherit most of the nice properties of the low-rank approximation format for order-two tensors. In particular, under suitable assumptions on tensor norms, best approximation problems in the set $\mathcal{T}_{r}^{A}$ are well-posed [14, 15]. Also, the $\alpha$-rank of a tensor can be computed through singular value decomposition, and the notion of singular value decomposition can be extended (in different ways) to these formats [9, 17, 36]. Another interesting property, which is not exploited in the present paper, is the fact that the set $\mathcal{T}_{r}^{A}$ is a differentiable manifold [22, 44, 15, 16], which has interesting consequences in optimization or model order reduction of dynamical systems in tensor spaces [29]. There are only a few results available on the approximation properties of tree-based formats [41]. However, it has been observed in many domains of applications that tree-based formats have a high approximation power (or expressive power). Hierarchical tensor formats have been recently identified with deep neural networks with a particular architecture [8].

The reader is referred to the monograph [19] and surveys [27, 25, 18, 34, 33, 1] for an introduction to tensor numerical methods and an overview of recent developments in the field.

This paper is concerned with the problem of computing an approximation of a function $u(x_{1},\ldots,x_{d})$ using point evaluations of this function, where evaluations can be selected adaptively. This includes problems where the function represents the output of a black-box numerical code, a system or a physical experiment for a given value of the input variables $(x_{1},\ldots,x_{d})$. This also includes the solution of high-dimensional PDEs with a probabilistic interpretation, where Monte-Carlo methods can be used to obtain point evaluations of their solutions. This excludes problems where evaluations of the functions come as an unstructured data set. A multivariate function $u(x_{1},\ldots,x_{d})$ is here considered as an element of a Hilbert tensor space $\mathcal{H}_{1}\otimes\ldots\otimes\mathcal{H}_{d}$ of real-valued functions defined on a product set $\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{d}$ equipped with a probability measure. This includes the case of multidimensional arrays when the variables $x_{\nu}$ take values in finite sets $\mathcal{X}_{\nu}$. In this case, a point evaluation corresponds to the evaluation of an entry of the tensor.

Several algorithms have been proposed for the construction of approximations in tree-based formats using point evaluations of functions or entries of tensors. Let us mention algorithms that use adaptive and structured evaluations of tensors [38, 2] and statistical learning approaches that use unstructured (random) evaluations of functions [13, 12, 6, 5]. Let us also mention the recent work [30] for the approximation in Tucker format, with an approach similar to the one proposed in the present paper.

In the present paper, we propose and analyse a new algorithm which is based on a particular extension of the singular value decomposition for the tree-based format $\mathcal{T}_{r}^{A}$ which allows us to construct an approximation using only evaluations of a function (or entries of a tensor). The proposed algorithm constructs a hierarchy of subspaces $U_{\alpha}$ of functions of groups of variables $x_{\alpha}$, for all $\alpha\in A$, and associated interpolation operators $I_{U_{\alpha}}$ which are oblique projections onto $U_{\alpha}$. For the construction of $U_{\alpha}$ for a particular node $\alpha\in A$, we interpret the function $u$ as a random variable $u(\cdot,x_{\alpha^{c}})$ depending on a set of random variables $x_{\alpha^{c}}$ with values in the space of functions of the variables $x_{\alpha}$. Then $U_{\alpha}$ is obtained by estimating the principal components of this function-valued random variable using random samples $u(\cdot,x^{k}_{\alpha^{c}})$. In practice, we estimate the principal components from interpolations $I_{{V_{\alpha}}}u(\cdot,x^{k}_{\alpha^{c}})$ of these samples on a subspace ${V_{\alpha}}$ which is a certain approximation space when $\alpha$ is a leaf of the tree, or the tensor product of subspaces $\{U_{\beta}\}_{\beta\in S(\alpha)}$ associated with the sons $S(\alpha)$ of the node $\alpha$ when $\alpha$ is not a leaf of the tree. This construction only requires evaluations of $u$ on a product set of points which is the product of an interpolation grid in $\mathcal{X}_{\alpha}$ (unisolvent for the space ${V_{\alpha}}$), and a random set of points in $\mathcal{X}_{\alpha^{c}}$. It is a sequential construction going from the leaves to the root of the tree.

The proposed algorithm can be interpreted as an extension of principal component analysis for tree-based tensors which provides a statistical estimation of low-dimensional subspaces of functions of groups of variables for the representation of a multivariate function. It is able to provide an approximation $u^{\star}$ in any tree-based format $\mathcal{T}_{r}^{A}$ with either a prescribed rank $r$ or a prescribed relative error (by adapting the rank $r$). For a given $r$, it has the remarkable property that it is able to provide an approximation in $\mathcal{T}_{r}^{A}$ with a number of evaluations equal to the storage complexity of the resulting approximation. Under some assumptions on the estimation of principal components, we prove that the algorithm, up to some discretization error $\rho$, provides with high probability a quasi-optimal approximation with a prescribed rank, i.e.

[TABLE]

where the constant $c$ depends on the set $A$ and the properties of orthogonal projections and interpolation operators associated with principal subspaces. Also, under some assumptions on the estimation of principal components and discretization error, we prove that the algorithm with prescribed tolerance $\epsilon$ is able to provide an approximation $u^{\star}$ such that

[TABLE]

holds with high probability, where the constant $\tilde{c}$ depends on the set $A$ and the properties of projections and interpolation operators. Sharp inequalities are obtained by considering the properties of projection and interpolation operators when restricted to minimal subspaces of tensors. The analysis takes into account the discretization errors for the approximation of infinite-dimensional tensors. For a tensor with finite and known rank in a tree-based format, and when there is no discretization error, the algorithm is able to recover the tensor in a stable way using a number of evaluations equal to the storage complexity of the representation of the tensor in this format. This algorithm may have important applications in the manipulation of big data, by providing a way to reconstruct a multidimensional array from a limited number of entries (tensor completion).

The outline of the paper is as follows. In section 2, we introduce some definitions and properties of projections in Hilbert spaces, with a particular attention on Hilbert spaces of functions and projections based on point evaluations. In section 3, we recall basic definitions on tensors and Hilbert tensor spaces of functions defined on measured product sets. Then we introduce some definitions and properties of operators on tensor spaces, with partial point evaluation functionals as a particular case. Finally, we introduce definitions and properties of projections on tensor spaces, with a particular attention on orthogonal projection and interpolation. In section 4, we introduce tree-based low-rank formats in a general setting including classical HT and TT formats. In section 5, we first introduce the notion of principal component analysis for multivariate functions and then propose an extension of principal component analysis to tree-based tensor format. This is based on a new variant of higher-order singular value decomposition of tensors in tree-based format. In section 6, we present and analyse a modified version of the algorithm presented in section 5 which only requires point evaluations of functions, and which is based on empirical principal component analyses and interpolations. In section 7, the behavior of the proposed algorithm is illustrated and analysed in several numerical experiments.

2 Projections

For two vector spaces $V$ and $W$ equipped with norms $\|\cdot\|_{V}$ and $\|\cdot\|_{W}$ respectively, we denote by $L(V,W)$ the space of linear operators from $V$ to $W$. We denote by $\mathcal{L}(V,W)$ the space of linear and continuous operators from $V$ to $W$, with bounded operator norm $\|A\|_{V\to W}=\max_{\|v\|_{V}=1}\|Av\|_{W}$. We denote by $V^{*}=L(V,\mathbb{R})$ the algebraic dual of $V$ and by $V^{\prime}=\mathcal{L}(V,\mathbb{R})$ the topological dual of $V$, and we let $\|\cdot\|_{V\to\mathbb{R}}=\|\cdot\|_{V^{\prime}}$. We denote by $\langle\cdot,\cdot\rangle$ the duality pairing between a space and its dual. We let $L(V):=L(V,V)$ and $\mathcal{L}(V):=\mathcal{L}(V,V)$, and we replace the notation $\|\cdot\|_{V\to V}$ by $\|\cdot\|_{V}$, where the latter notation also stands for the norm on $V$.

2.1 Projections

Let $V$ be a Hilbert space and $U$ be a finite-dimensional subspace of $V$. An operator $P$ is a projection onto a subspace $U$ if $\mathrm{Im}(P)=U$ and $Pu=u$ for all $u\in U$.

The orthogonal projection $P_{U}$ onto $U$ is a linear and continuous operator which associates to $v\in V$ the unique solution $P_{U}v\in U$ of

[TABLE]

or equivalently $(u,P_{U}v-v)=0,\;\forall u\in U.$ The orthogonal projection $P_{U}$ has operator norm $\|P_{U}\|_{V}=1$.

Let $W$ be a finite-dimensional subspace of $V^{*}$ such that

[TABLE]

where the latter condition is equivalent to $U\cap{}^{\perp}W=\{0\}$, with ${}^{\perp}W$ the annihilator of $W$ in $V$ (see [32, Definition 1.10.4]). Under the above assumptions, we have that for any $v\in V$, there exists a unique $u\in U$ such that $\langle w,u-v\rangle=0$ for all $w\in W$.111Uniqueness comes from (1b) while existence comes from (1a) and (1b). This allows to define the projection $P_{U}^{W}$ onto $U$ along $W$ which is the linear operator on $V$ which associates to $v\in V$ the unique solution $P_{U}^{W}v\in U$ of

[TABLE]

For $W=R_{V}U$, where $R_{V}:V\to V^{\prime}$ is the Riesz map, the projection $P_{U}^{W}$ coincides with the orthogonal projection $P_{U}$. A non orthogonal projection is called an oblique projection.

If $W\subset V^{\prime},$ then $P_{U}^{W}$ is a projection from $V$ onto $U$ parallel to $Ker(P_{U}^{W})=Z^{\perp}$, where $Z=R_{V}^{-1}W$. If $W\subset\tilde{U}^{\prime}$, with $\tilde{U}$ a closed subspace of $V$, then $P_{U}^{W}|_{\tilde{U}}$ is a projection from $\tilde{U}$ onto $U$ parallel to $Ker(P_{U}^{W})\cap\tilde{U}=Z^{\perp}\cap\tilde{U}$, where $Z=R_{\tilde{U}}^{-1}W$, with $R_{\tilde{U}}$ the Riesz map from $\tilde{U}$ to $\tilde{U}^{\prime}$.

Proposition 2.1**.**

*Let $\tilde{U}$ be a closed subspace of $V$ and assume that $U\subset\tilde{U}$ and $W\subset\tilde{U}^{\prime}.$222Note that $V^{\prime}\subset\tilde{U}^{\prime}$ and we may have $W\not\subset V^{\prime}.$ Then $P_{U}^{W}$ is a continuous operator from $\tilde{U}$ to $V$. *

Proof.

Let us equip $W$ with the norm $\|w\|_{W}=\|w\|_{\tilde{U}^{\prime}}=\max_{v\in\tilde{U}}\langle w,v\rangle/\|v\|_{V},$ such that for all $v\in\tilde{U}$, $\langle w,v\rangle\leq\|w\|_{W}\|v\|_{V}.$ Let

[TABLE]

Assumption (1b) implies that $\alpha>0$. Then for all $v\in\tilde{U}$, we have

[TABLE]

which ends the proof. ∎

Proposition 2.2**.**

Let $P$ and $\tilde{P}$ be projections onto subspaces $U$ and $\tilde{U}$ respectively and assume $U\subset\tilde{U}$. Then

[TABLE]

Moreover, if $P$ and $\tilde{P}$ are projections along $W$ and $\tilde{W}$ respectively, with $W\subset\tilde{W}$, then

[TABLE]

Proof.

For all $v\in V$, $Pv\in U\subset\tilde{U}$, and therefore $\tilde{P}Pv=Pv$, which proves the first statement. For the second statement, by definition of the projection $P$, we have that $\langle\phi,P\tilde{P}v-\tilde{P}v\rangle=0$ for all $\phi\in W$. Since $W\subset\tilde{W}$ and by definition of $\tilde{P}$, this implies that $\langle\phi,P\tilde{P}v-v\rangle=0$ for all $\phi\in W$. By definition of $Pv$ and since $P\tilde{P}v\in U,$ this implies $P\tilde{P}v=Pv=\tilde{P}Pv$. ∎

Proposition 2.3**.**

Let $U$ and $\tilde{U}$ be two closed subspaces of $V$, with $U$ of finite dimension. Let $P_{U}$ be the orthogonal projection onto $U$ and let $P_{U}^{W}$ be the projection onto $U$ along $W\subset\tilde{U}^{\prime}$. For all $v\in\tilde{U}$,

[TABLE]

with

[TABLE]

Also, for all $v\in\tilde{U}$,

[TABLE]

Proof.

For $v\in\tilde{U}$, $\|P_{U}^{W}v-P_{U}v\|_{V}=\|P_{U}^{W}(v-P_{U}v)\|_{V}=\|(P_{U}^{W}-P_{U})(v-P_{U}v)\|_{V}\leq\|P_{U}^{W}-P_{U}\|_{(id-P_{U})\tilde{U}\to V}\|v-P_{U}v\|_{V}$, with $\|P_{U}^{W}-P_{U}\|_{(id-P_{U})\tilde{U}\to V}=\|P_{U}^{W}-P_{U}\|_{\tilde{U}\to V}=\|P_{U}^{W}\|_{(id-P_{U})\tilde{U}\to V}$. This proves the first statement. The second statement directly follows from $\|v-P_{U}^{W}v\|_{V}^{2}=\|v-P_{U}v\|_{V}^{2}+\|P_{U}v-P_{U}^{W}v\|_{V}^{2}$. ∎

2.2 Projection of functions using point evaluations

Let $V$ be a Hilbert space of functions defined on a set $X$. For $x\in X$, the point evaluation functional $\delta_{x}\in V^{*}$ is defined by $\langle\delta_{x},v\rangle=v(x)$.

2.2.1 Interpolation

Let $U$ be a $n$-dimensional subspace of $V$ and let $\Gamma=\{x^{k}\}_{k=1}^{n}$ be a set of $n$ interpolation points in $X$. The set of interpolation points $\Gamma$ is assumed to be unisolvent for $U$, i.e. for any $(a_{k})_{k=1}^{n}\in\mathbb{R}^{n}$, there exists a unique $u\in U$ such that $u(x^{k})=a_{k}$ for all $1\leq k\leq n$. The interpolation operator $I_{U}$ associated with $\Gamma$ is a linear operator from $V$ to $U$ such that for $v\in V$, $I_{U}v$ is the unique element of $U$ such that

[TABLE]

The interpolation operator $I_{U}$ is an oblique projection $P_{U}^{W}$ onto $U$ along $W=\mathrm{span}\{\delta_{x}:x\in\Gamma\}$. Note that the condition that $\Gamma$ is unisolvent for $U$ is equivalent to the condition $\eqref{condition_UWb}$ on $U$ and $W$, which ensures that $I_{U}$ is well defined. From Proposition 2.2, we deduce the following property.

Proposition 2.4**.**

Let $U$ and $\tilde{U}$ be two subspaces associated with sets of interpolation points $\Gamma$ and $\tilde{\Gamma}$ respectively. If $U\subset\tilde{U}$ and $\Gamma\subset\tilde{\Gamma}$, then

[TABLE]

Magic points.

For a given basis $\{\varphi_{i}\}_{i=1}^{n}$ of $U$, a set of interpolation points $\Gamma=\{x^{k}\}_{k=1}^{n}$, called magic points, can be determined with a greedy algorithm proposed in [31, Remark 2]. The procedure for selecting the set $\Gamma$ in a subset $\Gamma_{\star}$ in $X$ is as follows. We first determine a point $x^{1}\in\Gamma_{\star}$ and an index $i_{1}$ such that

[TABLE]

Then for $k\geq 1$, we define $\psi_{i}^{(k)}(x)=\varphi_{i}(x)-\sum_{m=1}^{k}\sum_{p=1}^{k}\varphi_{i_{m}}(x)a^{(k)}_{m,p}\varphi_{i}(x^{p})$, with the matrix $(a^{(k)}_{m,p})_{1\leq m,p\leq k}$ being the inverse of the matrix

$(\varphi_{i_{m}}(x^{p}))_{1\leq p\leq k,1\leq m\leq k}$, such that $\psi_{i_{m}}^{(k)}(x)=0$ for all $1\leq m\leq k$ and $x\in X$, and $\psi_{i}^{(k)}(x^{p})=0$ for all $1\leq p\leq k$ and $1\leq i\leq n$. Then, we determine the point $x^{k+1}\in\Gamma_{\star}$ and an index $i_{k+1}$ such that

[TABLE]

2.2.2 Discrete least-squares projection

Let $U$ be a $n$-dimensional subspace of $V$ and let $\Gamma=\{x^{k}\}_{k=1}^{m}$ be a set of $m$ points in $X$, $m\geq n$, such that $\|v\|_{\Gamma}=(\sum_{x\in\Gamma}v(x)^{2})^{1/2}$ defines a norm on $U$. The discrete least-squares projection $Q_{U}$ is the linear operator from $V$ to $U$ such that for $v\in V$, $Q_{U}v$ is the unique element in $U$ which minimizes $\|v-u\|_{\Gamma}^{2}$ over all $u\in U$, or equivalently

[TABLE]

where $(\cdot,\cdot)_{\Gamma}$ is the inner product associated with the norm $\|\cdot\|_{\Gamma}$ on $U$. The discrete least-squares projection $Q_{U}$ is an oblique projection onto $U$ along $W=\{\sum_{x\in\Gamma}u(x)\delta_{x}:u\in U\}$. If $\#\Gamma=\dim(U)$ and $\Gamma$ is unisolvent for $U$, then $Q_{U}$ coincides with the interpolation operator $I_{U}$.

Proposition 2.5**.**

Let $U$ and $\tilde{U}$ be two finite-dimensional subspaces such that $U\subset\tilde{U}$. Let $Q_{U}$ be the discrete least-squares projection onto $U$ associated with a set of points $\Gamma$ in $X$, and let $Q_{\tilde{U}}$ be the discrete least-squares projection onto $\tilde{U}$ associated with a set of points $\tilde{\Gamma}$ in $X$. If either $\Gamma=\tilde{\Gamma}$ or $\Gamma\subset\tilde{\Gamma}$ and $\tilde{\Gamma}$ is unisolvent for $\tilde{U}$, then

[TABLE]

Proof.

$Q_{U}$ is the projection onto $U$ along $W=\{\sum_{x\in\Gamma}u(x)\delta_{x}:u\in U\}$, and $Q_{\tilde{U}}$ is the projection onto $\tilde{U}$ along $\tilde{W}=\{\sum_{x\in\tilde{\Gamma}}\tilde{u}(x)\delta_{x}:\tilde{u}\in\tilde{U}\}$. If we prove that $W\subset\tilde{W}$, then the result follows from Proposition 2.2. Let $w=\sum_{x\in\Gamma}u(x)\delta_{x}\in W$, with $u\in U$. If $\Gamma=\tilde{\Gamma}$, then since $u\in\tilde{U}$, we clearly have $w\in\tilde{W}$. If $\Gamma\subset\tilde{\Gamma}$ and $\tilde{\Gamma}$ is unisolvent for $\tilde{U}$, there exists a function $\tilde{u}\in\tilde{U}$ such that $\tilde{u}(x)=u(x)$ for all $x\in\Gamma$ and $\tilde{u}(x)=0$ for all $x\in\tilde{\Gamma}\setminus\Gamma$. Therefore, $w=\sum_{x\in\tilde{\Gamma}}\tilde{u}(x)\delta_{x}$ is an element of $\tilde{W}$, which ends the proof. ∎

3 Tensors

Let $\mathcal{H}_{\nu}$ be Hilbert spaces of real-valued functions defined on sets $\mathcal{X}_{\nu}$ equipped with probability measures $\mu_{\nu}$, $1\leq\nu\leq d$. We denote by $\|\cdot\|_{\mathcal{H}_{\nu}}$ the norm on $\mathcal{H}_{\nu}$ and by $(\cdot,\cdot)_{\mathcal{H}_{\nu}}$ the associated inner product. Let $\mathcal{X}=\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{d}$ and $\mu=\mu_{1}\otimes\ldots\otimes\mu_{d}$. The tensor product of $d$ functions $v^{\nu}\in\mathcal{H}_{\nu}$, $1\leq\nu\leq d$, denoted $v^{1}\otimes\ldots\otimes v^{d}$, is a multivariate function defined on $\mathcal{X}$ such that $(v^{1}\otimes\ldots\otimes v^{d})(x)=v^{1}(x_{1})\ldots v^{d}(x_{d})$ for $x=(x_{1},\ldots,x_{d})\in\mathcal{X}$. Such a function is called an elementary tensor. The algebraic tensor space $\mathcal{H}_{1}\otimes_{a}\ldots\otimes_{a}\mathcal{H}_{d}$ is defined as the linear span of all elementary tensors, which is a pre-Hilbert space when equipped with the canonical inner product $(\cdot,\cdot)$ defined for elementary tensors by

[TABLE]

and then extended by linearity to the whole algebraic tensor space. We denote by $\|\cdot\|$ the norm associated with inner product $(\cdot,\cdot)$. A Hilbert tensor space $\mathcal{H}=\overline{\mathcal{H}_{1}\otimes_{a}\ldots\otimes_{a}\mathcal{H}_{d}}^{\|\cdot\|}$ is then obtained by the completion of the algebraic tensor space, which we simply denote

[TABLE]

Example 3.1**.**

Consider finite sets $\mathcal{X}_{\nu}$ and $\mathcal{H}_{\nu}=\mathbb{R}^{\mathcal{X}_{\nu}}$ equipped with the norm $\|v\|_{\mathcal{H}_{\nu}}^{2}=\sum_{x_{\nu}\in\mathcal{X}_{\nu}}\mu_{\nu}(\{x_{\nu}\})|v(x_{\nu})|^{2}$. Then, $\mathcal{H}$ is the space of multidimensional arrays $\mathbb{R}^{\mathcal{X}_{1}}\otimes\ldots\otimes\mathbb{R}^{\mathcal{X}_{d}}$ and $\|v\|^{2}=\sum_{x\in\mathcal{X}}\mu(\{x\})|v(x)|^{2}$, where $\mu(\{x_{1},\ldots,x_{d}\})=\prod_{\nu=1}^{d}\mu_{\nu}(\{x_{\nu}\})$.

Example 3.2**.**

Consider $\mathcal{X}_{\nu}=\mathbb{R}$, $\mu_{\nu}$ a finite measure on $\mathbb{R}$, and $\mathcal{H}_{\nu}=L^{2}_{\mu_{\nu}}(\mathcal{X}_{\nu})$ equipped with the natural norm $\|v\|_{\mathcal{H}_{\nu}}^{2}=\int|v(x_{\nu})|^{2}\mu_{\nu}(dx_{\nu})$. Then $\mathcal{H}$ is identified with $L^{2}_{\mu}(\mathcal{X})$, where $\mu=\mu_{1}\otimes\ldots\otimes\mu_{d}$, and $\|v\|^{2}=\int|v(x)|^{2}\mu(dx).$

Example 3.3**.**

Consider for $\mathcal{H}_{\nu}$ a reproducing kernel Hilbert space (RKHS) with reproducing kernel $k_{\nu}:\mathcal{X}_{\nu}\times\mathcal{X}_{\nu}\to\mathbb{R}$. Then $\mathcal{H}$ is a RKHS with reproducing kernel $k(x,x^{\prime})=k_{1}(x_{1},x_{1}^{\prime})\ldots k_{d}(x_{d},x_{d}^{\prime})$.

For a non-empty subset $\alpha$ in $\{1,\ldots,d\}:=D$, we let $\mathcal{X}_{\alpha}$ be the set $\bigtimes_{\nu\in\alpha}\mathcal{X}_{\nu}$ equipped with the product measure $\mu_{\alpha}=\bigotimes_{\nu\in\alpha}\mu_{\nu}$. We denote by $\mathcal{H}_{\alpha}=\bigotimes_{\nu\in\alpha}\mathcal{H}_{\nu}$ the Hilbert tensor space of functions defined on $\mathcal{X}_{\alpha}$, equipped with the canonical norm $\|\cdot\|_{\mathcal{H}_{\alpha}}$ such that

[TABLE]

for $v^{\nu}\in\mathcal{H}_{\nu}$, $1\leq\nu\leq d$. We have $\mathcal{H}_{D}=\mathcal{H}$ and we use the convention $\mathcal{H}_{\emptyset}=\mathbb{R}$.

Matricisations and $\alpha$-ranks.

Let $\alpha\subset D$, with $\alpha\notin\{\emptyset,D\}$, and let $\alpha^{c}=D\setminus\alpha$ be its complement in $D$. For $x\in\mathcal{X}$, we denote by $x_{\alpha}$ the subset of variables $(x_{\nu})_{\nu\in\alpha}$. A tensor $v\in\mathcal{H}$ can be identified with an order-two tensor

[TABLE]

where $\mathcal{M}_{\alpha}$ is the matricisation operator associated with $\alpha$, which defines a linear isometry between $\mathcal{H}$ and $\mathcal{H}_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$. We use the conventions $\mathcal{M}_{\emptyset}(v)=\mathcal{M}_{D}(v)=v$ and $\mathcal{H}_{\emptyset}\otimes\mathcal{H}_{D}=\mathcal{H}_{D}\otimes\mathcal{H}_{\emptyset}=\mathcal{H}$. The $\alpha$-rank of a tensor $v\in\mathcal{H}$, denoted $\mathrm{rank}_{\alpha}(v)$, is defined as the rank of the order-two tensor $\mathcal{M}_{\alpha}(v)$, which is uniquely defined as the minimal integer such that

[TABLE]

for some functions $v^{\alpha}_{k}\in\mathcal{H}_{\alpha}$ and $v^{\alpha^{c}}_{k}\in\mathcal{H}_{\alpha^{c}}$ of complementary subsets of variables $x_{\alpha}$ and $x_{\alpha^{c}}$ respectively. By convention, we have $\mathrm{rank}_{\emptyset}(v)=\mathrm{rank}_{D}(v)=1.$ From now on, when there is no ambiguity, $\mathcal{M}_{\alpha}(v)$ and $\mathcal{H}_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ will be identified with $v$ and $\mathcal{H}$ respectively.

Minimal subspaces.

The minimal subspace $U^{min}_{\alpha}(v)$ of $v$ is defined as the smallest closed subspace in $\mathcal{H}_{\alpha}$ such that

[TABLE]

and we have $\mathrm{rank}_{\alpha}(v)=\dim(U^{min}_{\alpha}(v))$ (see [14]). If $v$ admits the representation (2), then $U^{min}_{\alpha}(v)$ is the closure of ${\mathrm{span}\{v^{\alpha}_{k}\}_{k=1}^{\mathrm{rank}_{\alpha}(v)}}.$ For any partition $S(\alpha)$ of $\alpha$, we have

[TABLE]

We have $U^{min}_{D}(v)=\mathbb{R}v$ and for any partition $S(D)$ of $D$,

[TABLE]

3.1 Operators on tensor spaces

Let consider the Hilbert tensor space $\mathcal{H}=\bigotimes_{\nu=1}^{d}\mathcal{H}_{\nu}$ equipped with the canonical norm $\|\cdot\|$. For linear operators from $\mathcal{H}$ to $\mathcal{H}$, we also denote by $\|\cdot\|$ the operator norm $\|\cdot\|_{\mathcal{H}\to\mathcal{H}}=\|\cdot\|_{\mathcal{H}}$.

We denote by ${id}$ the identity operator on $\mathcal{H}$. For a non-empty subset $\alpha\subset D$, we denote by ${id}_{\alpha}$ the identity operator on $\mathcal{H}_{\alpha}$. For $A_{\alpha}$ in $L(\mathcal{H}_{\alpha})$, we define the linear operator $A_{\alpha}\otimes{id}_{\alpha^{c}}$ such that for $v^{\alpha}\in\mathcal{H}_{\alpha}$ and $v^{\alpha^{c}}\in\mathcal{H}_{\alpha^{c}}$,

[TABLE]

and we extend this definition by linearity to the whole algebraic tensor space $\mathcal{H}_{\alpha}\otimes_{a}\mathcal{H}_{\alpha^{c}}$. For a finite dimensional tensor space $\mathcal{H}$, this completely characterizes a linear operator on $\mathcal{H}$. For an infinite dimensional tensor space $\mathcal{H}$, if $A_{\alpha}\in\mathcal{L}(U_{\alpha},\mathcal{H}_{\alpha})$, with $U_{\alpha}\subset\mathcal{H}_{\alpha}$, then $A_{\alpha}\otimes{id}_{\alpha^{c}}$ can be extended by continuity to $U_{\alpha}\otimes\mathcal{H}$.

We denote by $\mathcal{A}_{\alpha}$, using calligraphic font style, the linear operator in $L(\mathcal{H})$ associated with an operator $A_{\alpha}$ in $L(\mathcal{H}_{\alpha})$, defined by $\mathcal{A}_{\alpha}=\mathcal{M}_{\alpha}^{-1}(A_{\alpha}\otimes{id}_{\alpha^{c}})\mathcal{M}_{\alpha},$ and simply denoted

[TABLE]

when there is no ambiguity. If $A_{\alpha}\in\mathcal{L}(\mathcal{H}_{\alpha})$, then $\mathcal{A}_{\alpha}\in\mathcal{L}(\mathcal{H})$ and the two operators have the same operator norm $\|\mathcal{A}_{\alpha}\|=\|A_{\alpha}\|_{\mathcal{H}_{\alpha}}.$ Also, we have the following more general result.

Proposition 3.4**.**

If $A_{\alpha}\in\mathcal{L}(U_{\alpha},\mathcal{H}_{\alpha})$, with $U_{\alpha}\subset\mathcal{H}_{\alpha}$, then $\mathcal{A}_{\alpha}\in\mathcal{L}(U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}},\mathcal{H})$ and the two operators have the same operator norm

[TABLE]

Corollary 3.5**.**

For a tensor $v\in\mathcal{H}$ and an operator $A_{\alpha}\in\mathcal{L}(U^{min}_{\alpha}(v),\mathcal{H}_{\alpha})$,

[TABLE]

Let $S=\{\alpha_{1},\ldots,\alpha_{K}\}$ be a collection of disjoint subsets of $D$ and let $A_{\alpha}\in L(\mathcal{H}_{\alpha})$ be linear operators, $\alpha\in S$. Then we can define a linear operator $A_{\alpha_{1}}\otimes\ldots\otimes A_{\alpha_{K}}:=\bigotimes_{\alpha\in S}A_{\alpha}$ on $\mathcal{H}_{\alpha_{1}}\otimes_{a}\ldots\otimes_{a}\mathcal{H}_{\alpha_{K}}$ such that

[TABLE]

for $v^{\alpha}\in\mathcal{H}_{\alpha}$, $\alpha\in S$. The operator $\bigotimes_{\alpha\in S}A_{\alpha}$ can be identified with an operator

[TABLE]

defined on the algebraic tensor space $\mathcal{H}_{1}\otimes_{a}\ldots\otimes_{a}\mathcal{H}_{d}$. The definition of $\mathcal{A}$ is independent of the ordering of the elements of $S$. If the operators $A_{\alpha}$ are continuous, then $\mathcal{A}$ defines a continuous operator from $\mathcal{H}$ to $\mathcal{H}$ and since $\|\cdot\|$ is a uniform crossnorm (see [19, Proposition 4.127]), the operator $\mathcal{A}$ has for operator norm

[TABLE]

Also, we have the following more general result.

Proposition 3.6**.**

Let $S$ be a collection of disjoint subsets of $D$ and let $\beta\subset D$ such that $\beta\cup(\cup_{\alpha\in S}\alpha)=D$. Let $U_{\alpha}$ be a subspace of $\mathcal{H}_{\alpha}$ and $A_{\alpha}\in\mathcal{L}(U_{\alpha},\mathcal{H}_{\alpha})$, for $\alpha\in S$. Then $\mathcal{A}=\prod_{\alpha\in S}\mathcal{A}_{\alpha}$ is a continuous operator from $\mathcal{U}:=(\bigotimes_{\alpha\in S}U_{\alpha})\otimes\mathcal{H}_{\beta}$ to $\mathcal{H}$ such that

[TABLE]

Corollary 3.7**.**

Let $S$ be a collection of disjoint subsets of $D$. For a tensor $v\in\mathcal{H}$ and operators $A_{\alpha}$, $\alpha\in S$, such that $A_{\alpha}\in\mathcal{L}(U^{min}_{\alpha}(v),\mathcal{H}_{\alpha})$, the operator $\mathcal{A}=\prod_{\alpha\in S}\mathcal{A}_{\alpha}$ is such that

[TABLE]

3.2 Partial evaluations of tensors

Let $\alpha$ be a non-empty subset of $D$. For a linear form $\psi_{\alpha}\in\mathcal{H}_{\alpha}^{*}$, $\psi_{\alpha}\otimes id_{\alpha^{c}}$ is a linear operator from $\mathcal{H}_{\alpha}\otimes_{a}\mathcal{H}_{\alpha^{c}}$ to $\mathcal{H}_{\alpha^{c}}$ such that $(\psi_{\alpha}\otimes id_{\alpha^{c}})(v^{\alpha}\otimes v^{\alpha^{c}})=\psi_{\alpha}(v^{\alpha})v^{\alpha^{c}}.$ If $\psi_{\alpha}\in\mathcal{H}_{\alpha}^{\prime}$, the definition of $\psi_{\alpha}\otimes id_{\alpha^{c}}$ can be extended by continuity to $\mathcal{H}$. Then $\psi_{\alpha}\otimes id_{\alpha^{c}}$ is a continuous operator from $\mathcal{H}$ to $\mathcal{H}_{\alpha^{c}}$ with operator norm $\|\psi_{\alpha}\otimes id_{\alpha^{c}}\|_{\mathcal{H}\to\mathcal{H}_{\alpha^{c}}}=\|\psi_{\alpha}\|_{\mathcal{H}_{\alpha}^{\prime}}.$ Also, we have the following result.

Proposition 3.8**.**

If $\psi_{\alpha}\in U_{\alpha}^{\prime}$, with $U_{\alpha}$ a subspace of $\mathcal{H}_{\alpha}$, then $\psi_{\alpha}\otimes id_{\alpha^{c}}\in\mathcal{L}(U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}},\mathcal{H}_{\alpha^{c}})$ and

[TABLE]

Corollary 3.9**.**

For a tensor $v\in\mathcal{H}$ and $\psi_{\alpha}\in U^{min}_{\alpha}(v)^{\prime}$, we have

[TABLE]

For a point $x_{\alpha}\in\mathcal{X}_{\alpha}$, we denote by $\delta_{x_{\alpha}}\in\mathcal{H}_{\alpha}^{*}$ the point evaluation functional at $x_{\alpha}$, defined by $\langle\delta_{x_{\alpha}},v^{\alpha}\rangle=v^{\alpha}(x_{\alpha})$ for $v^{\alpha}\in\mathcal{H}_{\alpha}$. Then $\delta_{x_{\alpha}}\otimes id_{\alpha^{c}}$ defines a partial evaluation functional, which is a linear operator from $\mathcal{H}$ to $\mathcal{H}_{\alpha^{c}}$ such that

[TABLE]

From Corollary 3.9, we deduce that for a given tensor $v\in\mathcal{H}$, if $\delta_{x_{\alpha}}\in U^{min}_{\alpha}(v)^{\prime}$, then the definition of $\delta_{x_{\alpha}}\otimes id_{\alpha^{c}}$ can be extended by continuity to $U^{min}_{\alpha}(v)\otimes\mathcal{H}_{\alpha^{c}}$ and the partial evaluation

[TABLE]

is an element of $\mathcal{H}_{\alpha^{c}}$ such that

[TABLE]

3.3 Projection of tensors

Let $\alpha$ be a non-empty and strict subset of $D$ and let $U_{\alpha}$ be a finite-dimensional subspace of $\mathcal{H}_{\alpha}$. If $P_{\alpha}$ is a projection from $\mathcal{H}_{\alpha}$ onto $U_{\alpha}$, then $P_{\alpha}\otimes id_{\alpha^{c}}$ is a projection from $\mathcal{H}_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ onto $U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$.

Proposition 3.10**.**

Let $v\in\mathcal{H}$ and $\alpha,\beta\subset D$. Let $P_{\beta}$ be a projection from $\mathcal{H}_{\beta}$ to a subspace $U_{\beta}$ and let $\mathcal{P}_{\beta}$ be the corresponding projection onto $U_{\beta}\otimes\mathcal{H}_{\beta^{c}}$. If $\beta\subset\alpha$ or $\beta\subset D\setminus\alpha$, we have

[TABLE]

Proof.

A tensor $v$ admits a representation $v=\sum_{k=1}^{\mathrm{rank}_{\alpha}(v)}v^{\alpha}_{k}\otimes w^{\alpha^{c}}_{k}.$ If $\beta\subset\alpha$, then $\mathcal{P}_{\beta}=(P_{\beta}\otimes{id}_{\alpha\setminus\beta})\otimes{id}_{D\setminus\alpha}$ and $\mathcal{P}_{\beta}v=\sum_{k=1}^{\mathrm{rank}_{\alpha}(v)}((P_{\beta}\otimes{id}_{\alpha\setminus\beta})v^{\alpha}_{k})\otimes w^{\alpha^{c}}_{k}.$ If $\beta\subset D\setminus\alpha$, then $\mathcal{P}_{\beta}={id}_{\alpha}\otimes(P_{\beta}\otimes{id}_{D\setminus\{\alpha\cup\beta\}})$ and $\mathcal{P}_{\beta}v=\sum_{k=1}^{\mathrm{rank}_{\alpha}(v)}v^{\alpha}_{k}\otimes((P_{\beta}\otimes{id}_{D\setminus\{\alpha\cup\beta\}})w^{\alpha^{c}}_{k}).$ The result follows from the definition of the $\alpha$-rank. ∎

If $P_{U_{\alpha}}$ is the orthogonal projection from $\mathcal{H}_{\alpha}$ onto $U_{\alpha}$, then $P_{U_{\alpha}}\otimes id_{\alpha^{c}}$ coincides with the orthogonal projection $P_{U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}}$ from $\mathcal{H}_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ onto $U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$, and is identified with the orthogonal projection $\mathcal{P}_{U_{\alpha}}=P_{U_{\alpha}}\otimes id_{\alpha^{c}}$ in $\mathcal{L}(\mathcal{H})$. If $P_{U_{\alpha}}^{W_{\alpha}}$ is the oblique projection onto $U_{\alpha}$ along $W_{\alpha}\subset\mathcal{H}_{\alpha}^{*}$, then $\mathcal{P}_{U_{\alpha}}^{W_{\alpha}}:=P_{U_{\alpha}}^{W_{\alpha}}\otimes id_{\alpha^{c}}$ is the oblique projection from $\mathcal{H}_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ onto $U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ along $W_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}^{\prime}$. If $W_{\alpha}\subset\mathcal{H}_{\alpha}^{\prime}$, then $P_{U_{\alpha}}^{W_{\alpha}}$ and $\mathcal{P}_{U_{\alpha}}^{W_{\alpha}}$ are continuous operators with equal norms $\|\mathcal{P}_{U_{\alpha}}^{W_{\alpha}}\|=\|P_{U_{\alpha}}^{W_{\alpha}}\|_{\mathcal{H}_{\alpha}}.$

Proposition 3.11**.**

Let $U_{\alpha}$ be a finite-dimensional subspace of $\mathcal{H}_{\alpha}$ and let $P_{U_{\alpha}}^{W_{\alpha}}$ be the projection onto $U_{\alpha}$ along $W_{\alpha}$. For a tensor $v\in\mathcal{H}$ such that $W_{\alpha}\subset U^{min}_{\alpha}(v)^{\prime}$, $\mathcal{P}^{W_{\alpha}}_{U_{\alpha}}v$ is an element of $U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ such that

[TABLE]

and

[TABLE]

with

[TABLE]

Also,

[TABLE]

Proof.

We have $v\in U^{min}_{\alpha}(v)\otimes\mathcal{H}_{\alpha^{c}}$. Noting that $\|\mathcal{P}^{W_{\alpha}}_{U_{\alpha}}\|_{U^{min}_{\alpha}(v)\otimes\mathcal{H}_{\alpha^{c}}\to\mathcal{H}}=\|P^{W_{\alpha}}_{U_{\alpha}}\|_{U^{min}_{\alpha}(v)\to\mathcal{H}_{\alpha}}$ and $\|\mathcal{P}^{W_{\alpha}}_{U_{\alpha}}-\mathcal{P}_{U_{\alpha}}\|_{U^{min}_{\alpha}(v)\otimes\mathcal{H}_{\alpha^{c}}\to\mathcal{H}}=\|P^{W_{\alpha}}_{U_{\alpha}}-P_{U_{\alpha}}\|_{U^{min}_{\alpha}(v)\to\mathcal{H}_{\alpha}}$, the results directly follow from Proposition 2.3. ∎

Now, let $\alpha$ be a non-empty subset of $D$ and let $S(\alpha)$ be a partition of $\alpha$. Let $P_{U_{\beta}}^{W_{\beta}}$ be oblique projections onto subspaces $U_{\beta}$ of $\mathcal{H}_{\beta}$ along $W_{\beta}\subset\mathcal{H}_{\beta}^{*}$, $\beta\in S(\alpha)$. Then $\bigotimes_{\beta\in S(\alpha)}P_{U_{\beta}}^{W_{\beta}}:=P_{U_{S(\alpha)}}^{W_{S(\alpha)}}$ is the oblique projection from $\mathcal{H}_{S(\alpha)}=\bigotimes_{\beta\in S(\alpha)}\mathcal{H}_{\beta}$ onto $\bigotimes_{\beta\in S(\alpha)}U_{\beta}:=U_{S(\alpha)}$ along $\bigotimes_{\beta\in S(\alpha)}W_{\beta}:=W_{S(\alpha)}$, and $\mathcal{P}^{W_{S(\alpha)}}_{U_{S(\alpha)}}=P_{U_{S(\alpha)}}^{W_{S(\alpha)}}\otimes id_{\alpha^{c}}$ is the oblique projection from $\mathcal{H}_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ to $U_{S(\alpha)}\otimes\mathcal{H}_{\alpha^{c}}$ along $W_{S(\alpha)}\otimes\mathcal{H}_{\alpha^{c}}^{\prime}$. From Proposition 2.2, we directly obtain the following result.

Proposition 3.12**.**

If $U_{\alpha}\subset\bigotimes_{\beta\in S(\alpha)}U_{\beta}$ and $W_{\alpha}\subset\bigotimes_{\beta\in S(\alpha)}W_{\beta}$, then

[TABLE]

4 Tree-based tensor formats

Let $T\subset 2^{D}\setminus\emptyset$ be a dimension partition tree over $D$, with root $D$. The elements of $T$ are called the nodes of the tree. Every node $\alpha\in T$ with $\#\alpha\geq 2$ has a set of sons $S(\alpha)$ which form a partition of $\alpha$, i.e. $\bigcup_{\beta\in S(\alpha)}\beta=\alpha$. A node $\alpha\in T$ with $\#\alpha=1$ is such that $S(\alpha)=\emptyset$ and is called a leaf of the tree. The set of leaves of $T$ is denoted $\mathcal{L}(T)$ (see an example on Figure 1).

For $\alpha\in T$, we denote by $\mathrm{level}(\alpha)$ the level of $\alpha$ in $T$, such that $\mathrm{level}(D)=0$ and $\mathrm{level}(\beta)=\mathrm{level}(\alpha)+1$ if $\beta\in S(\alpha)$. We let $L=0pt(T)=\max_{\alpha\in T}\mathrm{level}(\alpha)$ be the depth of $T$, which is the maximum level of the nodes in $T$, and $T_{\ell}=\{\alpha\in T:\mathrm{level}(\alpha)=\ell\}$ be the subset of nodes with level $\ell$, $0\leq\ell\leq L$. We let $t_{\ell}=\bigcup_{\alpha\in T_{\ell}}\alpha.$ We have $t_{\ell+1}\subset t_{\ell}$ and $t_{\ell}\setminus t_{\ell+1}\subset\mathcal{L}(T)$ (see example on Figure 2).

We introduce a subset of active nodes $A\subset T\setminus\{D\}$ such that $T\setminus A\subset\{D\}\cup\mathcal{L}(T)$, which means that the set of non active nodes in $T\setminus\{D\}$ is a subset of the leaves (see Figure 3). A set $A$ is admissible if for any $\alpha\in A$, the parent node of $\alpha$ is in $A\cup\{D\}$. We let $\mathcal{L}(A)=A\cap\mathcal{L}(T)$, $A_{\ell}=A\cap T_{\ell}$ for $1\leq\ell\leq L$, and $a_{\ell}=\cup_{\alpha\in A_{\ell}}\alpha$. We define the $A$-rank of a tensor $v\in\mathcal{H}$ as the tuple $\mathrm{rank}_{A}(v)=\{\mathrm{rank}_{\alpha}(v)\}_{\alpha\in A}$.

Now we consider a tensor $v\in\mathcal{H}$ with $\mathrm{rank}_{A}(v)=(r_{\alpha})_{\alpha\in A}$. We let $r_{D}=\mathrm{rank}_{D}(v)=1$. For all $\alpha\in A\cup\{D\}$, we denote by $\{v^{\alpha}_{k_{\alpha}}\}_{k_{\alpha}=1}^{r_{\alpha}}$ a basis of the minimal subspace $U^{min}_{\alpha}(v)\subset\mathcal{H}_{\alpha}$, and we let $v^{D}_{1}=v$. For $\alpha\in A\cup\{D\}$ such that $\emptyset\neq S(\alpha)\subset A$, since $U^{min}_{\alpha}(v)\subset\bigotimes_{\beta\in S(\alpha)}U^{min}_{\beta}(v)$, the tensor $v^{\alpha}_{k_{\alpha}}$ admits a representation

[TABLE]

with a tensor of coefficients $C^{\alpha}\in\mathbb{R}^{r_{\alpha}\times\bigtimes_{\beta\in S(\alpha)}r_{\beta}}$. For $\alpha\in A\cup\{D\}$ such that $\emptyset\neq S(\alpha)\not\subset A$, we have $U^{min}_{\alpha}(v)\subset(\bigotimes_{\beta\in S(\alpha)\cap A}U^{min}_{\beta}(v))\otimes(\bigotimes_{\beta\in S(\alpha)\setminus A}\mathcal{H}_{\beta})$, and therefore the tensor $v^{\alpha}_{k_{\alpha}}$ admits a representation

[TABLE]

with $C^{\alpha}\in\mathbb{R}^{r_{\alpha}\times\bigtimes_{\beta\in S(\alpha)\cap A}r_{\beta}}\otimes(\bigotimes_{\beta\in S(\alpha)\setminus A}\mathcal{H}_{\beta}).$ Finally, a tensor $v$ such that $\mathrm{rank}_{A}(v)=(r_{\alpha})_{\alpha\in A}$ admits a representation

[TABLE]

For a tuple $r=(r_{\alpha})_{\alpha\in A}$, we define the subset $\mathcal{T}^{A}_{r}(\mathcal{H})$ of tensors in $\mathcal{H}$ with $A$-rank bounded by $r$,

[TABLE]

Remark 4.1**.**

A tensor $v\in\mathcal{T}^{A}_{r}(\mathcal{H})$ admits a representation as a composition of functions. For $\alpha\in A$, let $v^{\alpha}(x_{\alpha})=(v_{1}^{\alpha},\ldots,v_{r_{\alpha}}^{\alpha})\in\mathbb{R}^{r_{\alpha}}$. If $\emptyset\neq S(\alpha)\subset A$, the tensor $C^{\alpha}$ can be identified with a multilinear function $f^{\alpha}:\bigtimes_{\beta\in S(\alpha)}\mathbb{R}^{r_{\beta}}\to\mathbb{R}^{r_{\alpha}}$, and $v^{\alpha}(x_{\alpha})$ admits the representation

[TABLE]

For $\alpha\in A\cup\{D\}$ such that $\emptyset\neq S(\alpha)\not\subset A$, the tensor $C^{\alpha}((x_{\beta})_{\beta\in S(\alpha)\setminus A})$ can be identified with a multilinear function $f^{\alpha}(\cdot,(x_{\beta})_{\beta\in S(\alpha)\setminus A}):\bigtimes_{\beta\in S(\alpha)\cap A}\mathbb{R}^{r_{\beta}}\to\mathbb{R}^{r_{\alpha}}$, and $v^{\alpha}(x_{\alpha})$ admits the representation

[TABLE]

where the $f^{\alpha}$ is linear in the arguments associated with active nodes $\beta\in S(\alpha)\cap A$. As an example, for the case of Figure 3, the tensor $v$ admits the representation

[TABLE]

Proposition 4.2**.**

Let $V=V_{1}\otimes\ldots\otimes V_{d}\subset\mathcal{H}$, with $V_{\nu}$ a subspace of $\mathcal{H}_{\nu}$ with dimension $\dim(V_{\nu})=n_{\nu}$, $1\leq\nu\leq d$. The storage complexity of a tensor in $\mathcal{T}_{r}^{A}(\mathcal{H})\cap V=\mathcal{T}_{r}^{A}(V)$ is

[TABLE]

Example 4.3** (Tucker format).**

The Tucker format corresponds to a trivial tree $T=\{\{1,\ldots,d\},\{1\},\ldots,\{d\}\}$ with depth $L=1$, and $A=T\setminus\{D\}$ (see Figure 4). A tensor $v$ with $A$-rank bounded by $(r_{1},\ldots,r_{d})$ admits a representation of the form

[TABLE]

where $C\in\mathbb{R}^{r_{1}\times\ldots\times r_{d}}$, and $v^{\nu}_{k_{\nu}}\in\mathcal{H}_{\nu}$, $1\leq\nu\leq d$, or equivalently

[TABLE]

Example 4.4** (Degenerate Tucker format).**

A degenerate Tucker format corresponds to a trivial tree $T=\{\{1,\ldots,d\},\{1\},\ldots,\{d\}\}$ with depth $L=1$, and an active set of nodes $A$ strictly included in $T\setminus\{D\}$. Up to a permutation of dimensions, this corresponds to $A=\{\{1\},\ldots,\{p\}\}$, with $p<d$. A tensor $v$ with $A$-rank bounded by $(r_{1},\ldots,r_{p})$ admits a representation of the form

[TABLE]

where $C\in\mathbb{R}^{r_{1}\times\ldots\times r_{p}}\otimes\mathcal{H}_{\{p+1,\ldots,d\}}$, and $v^{\nu}_{k_{\nu}}\in\mathcal{H}_{\nu}$, $1\leq\nu\leq p$, or equivalently

[TABLE]

Example 4.5** (Tensor train format).**

*The tensor train (TT) format corresponds to a linear tree $T=\{\{1\},\{2\},\ldots,\{d\},\{1,2\},\ldots,\{1,\ldots,d\}\}$ and

$A=\{\{1\},\{1,2\},\ldots,\{1,\ldots,d-1\}\}$ (see Figure 5). Here, $A$ is a strict subset of $T\setminus\{D\}$. The nodes $\{2\},\ldots,\{d\}$ in $T$ are not active333Note that since $\mathrm{rank}_{\{d\}}(v)=\mathrm{rank}_{\{1,\ldots,d-1\}}(v)$, adding the node $\{d\}$ in the set of active nodes $A$ would yield an equivalent tensor format.. A tensor $v$ with $A$-rank bounded by $(r_{1},\ldots,r_{d-1})$ admits a representation of the form*

[TABLE]

where $v^{1}\in\mathbb{R}^{r_{1}}\otimes\mathcal{H}_{1}$, $C^{\nu}\in\mathbb{R}^{r_{\nu-1}\times r_{\nu}}\otimes\mathcal{H}_{\nu}$ for $2\leq\nu\leq d$, with the convention $r_{d}=1$. Here $L=d-1$, and for $1\leq\ell\leq L$, $T_{\ell}=\{\{1,\ldots,d-\ell\},\{d-\ell+1\}\}$, $t_{\ell}=\{1,\ldots,d-\ell+1\}$, $A_{\ell}=\{\{1,\ldots,d-\ell\}\}$ and $a_{\ell}=\{1,\ldots,d-\ell\}$. The tensor $v$ admits the equivalent representation

[TABLE]

Example 4.6** (Tensor train Tucker format).**

The tensor train Tucker (TTT) format corresponds to a linear tree $T=\{\{1\},\ldots,\{d\},\{1,2\},\ldots,\{1,\ldots,d\}\}$ and $A=T\setminus\{D\}$ (see Figure 6). A tensor $v$ having a $A$-rank bounded by $(r_{1},\ldots,r_{d},s_{2},\ldots,s_{d-1})$ admits a representation of the form (4) with a tensor $C\in\mathbb{R}^{r_{1}\times\ldots\times r_{d}}$ such that

[TABLE]

where $C^{2}\in\mathbb{R}^{r_{1}\times r_{2}\times s_{2}}$ and $C^{k}\in\mathbb{R}^{s_{k-1}\times r_{k}\times s_{k}}$ for $3\leq k\leq d$, with the convention $s_{d}=1$. Here $L=d-1$, $T_{\ell}=A_{\ell}=\{\{1,\ldots,d-\ell\},\{d-\ell+1\}\}$ and $t_{\ell}=a_{\ell}=\{1,\ldots,d-\ell+1\}$ for $1\leq\ell\leq L$. The tensor $v$ admits the equivalent representation

[TABLE]

5 Principal component analysis for tree-based tensor format

5.1 Principal component analysis of multivariate functions

Here we introduce the notion of principal component analysis for multivariate functions. We consider a given non-empty and strict subset $\alpha$ of $D$. Any tensor in $\mathcal{H}$ is identified (through the linear isometry $\mathcal{M}_{\alpha}$) with its $\alpha$-matricisation in $\mathcal{H}_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$. A tensor $u$ with $\alpha$-rank $\mathrm{rank}_{\alpha}(u)\in\mathbb{N}\cup\{+\infty\}$ admits a singular value decomposition (see [19, Section 4.4.3])

[TABLE]

where $\{u^{\alpha}_{k}\}_{k=1}^{\mathrm{rank}_{\alpha}(u)}$ and $\{u^{\alpha^{c}}_{k}\}_{k=1}^{\mathrm{rank}_{\alpha}(u)}$ are orthonormal vectors in $\mathcal{H}_{\alpha}$ and $\mathcal{H}_{\alpha^{c}}$ respectively, and where the $\sigma_{k}^{\alpha}$ are the $\alpha$-singular values of $u$ which are supposed to be arranged by decreasing values. The minimal subspace $U^{min}_{\alpha}(u)$ of $u$ is given by

[TABLE]

For $r_{\alpha}<\mathrm{rank}_{\alpha}(u)$, the truncated singular value decomposition

[TABLE]

is such that

[TABLE]

The functions $\{u^{\alpha}_{k}\}_{k=1}^{r_{\alpha}}$ are the $r_{\alpha}$ principal components of $u$ associated with dimensions $\alpha$, hereafter called $\alpha$-principal components. The corresponding subspace $U_{\alpha}^{\star}=\mathrm{span}\{u^{\alpha}_{k}\}_{k=1}^{r_{\alpha}}$, which is a subspace of $U^{min}_{\alpha}(u)$, is hereafter called a $\alpha$-principal subspace of dimension $r_{\alpha}$. Denoting $\mathcal{P}_{U_{\alpha}^{\star}}=P_{U_{\alpha}^{\star}}\otimes id_{\alpha^{c}}$ the orthogonal projection from $\mathcal{H}$ to $U_{\alpha}^{\star}\otimes\mathcal{H}_{\alpha^{c}}$, we have $u_{r_{\alpha}}=\mathcal{P}_{U_{\alpha}^{\star}}u,$444For all $m\geq r_{\alpha}$, we have $\mathcal{P}_{U_{\alpha}^{\star}}u_{m}=\sum_{k=1}^{m}\sigma_{k}^{\alpha}(P_{U_{\alpha}^{\star}}u_{k}^{\alpha})\otimes u_{k}^{\alpha^{c}}=\sum_{k=1}^{r_{\alpha}}\sigma_{k}^{\alpha}u_{k}^{\alpha}\otimes u_{k}^{\alpha^{c}}=u_{r_{\alpha}}$. Then using the continuity of $\mathcal{P}_{U_{\alpha}^{\star}}$ and taking the limit with $m$, we obtain $\mathcal{P}_{U_{\alpha}^{\star}}u=u_{r_{\alpha}}$. and

[TABLE]

Remark 5.1**.**

The optimization problem (7) over subspaces of dimension $r_{\alpha}$ in $\mathcal{H}_{\alpha}$ admits a unique solution $U_{\alpha}^{\star}$ if and only if $\sigma_{r_{\alpha}+1}^{\alpha}>\sigma_{r_{\alpha}}^{\alpha}$.

5.2 Principal component analysis for tree-based tensor format

Here, we propose and analyse an algorithm for the construction of an approximation $u^{\star}$ of a function $u$ in tree-based format $\mathcal{T}^{A}_{r}(\mathcal{H})$. It is based on the construction of a hierarchy of subspaces $U_{\alpha}$, $\alpha\in A$, from principal component analyses of approximations of $u$ in low-dimensional spaces in $\mathcal{H}_{\alpha}$. This is a variant of the leaves-to-root higher-order singular value decomposition method proposed in [17] (see also [19, Section 11.4.2.3]).

For each leaf node $\alpha\in\mathcal{L}(T)$, we introduce a finite-dimensional approximation space $V_{\alpha}\subset\mathcal{H}_{\alpha}$ with dimension $\dim(V_{\alpha})=n_{\alpha}$, and we let $V=\bigotimes_{\alpha\in\mathcal{L}(T)}V_{\alpha}\subset\mathcal{H}$. For each non active node $\alpha\in\mathcal{L}(T)\setminus A$, we let $U_{\alpha}=V_{\alpha}$. The algorithm then goes through all active nodes of the tree, going from the leaves to the root. For each $\alpha\in A$, we let

[TABLE]

where for $\alpha\notin\mathcal{L}(T)$, $V_{\alpha}$ is defined by

[TABLE]

where the $U_{\beta}$, $\beta\in S(\alpha)$, have been determined at a previous step.

Then we determine the $r_{\alpha}$-dimensional $\alpha$-principal subspace $U_{\alpha}$ of $u_{\alpha}$, which is solution of

[TABLE]

Finally, we define

[TABLE]

where $\mathcal{P}_{V_{D}}$ is the orthogonal projection from $\mathcal{H}$ onto $V_{D}=\bigotimes_{\beta\in S(D)}U_{\beta}.$

5.3 Analysis of the algorithm

Lemma 5.2**.**

For $\alpha\in\mathcal{L}(A)$, $U_{\alpha}\subset V_{\alpha}$. For $\alpha\in A\setminus\mathcal{L}(A)$,

[TABLE]

Proof.

For $\alpha\in A$, we have $U_{\alpha}\subset U^{min}_{\alpha}(u_{\alpha})$. If $\alpha\in\mathcal{L}(A)$, we have $U^{min}_{\alpha}(u_{\alpha})\subset V_{\alpha}$ since $u_{\alpha}=\mathcal{P}_{V_{\alpha}}u$. If $\alpha\in A\setminus\mathcal{L}(A)$, we have $U^{min}_{\alpha}(u_{\alpha})\subset\bigotimes_{\beta\in S(\alpha)}U^{min}_{\beta}(u_{\alpha})$, and $U^{min}_{\beta}(u_{\alpha})\subset U_{\beta}$ since $u_{\alpha}=\prod_{\beta\in S(\alpha)}\mathcal{P}_{U_{\beta}}u.$ ∎

Proposition 5.3**.**

The approximation $u^{\star}$ is an element of $\mathcal{T}_{r}^{A}(\mathcal{H})\cap V=\mathcal{T}^{A}_{r}(V)$.

Proof.

Since $u^{\star}={\mathcal{P}_{V_{D}}u}$, we have $u^{\star}\in\bigotimes_{\alpha\in S(D)}U_{\alpha}$. Then using Lemma 5.2, we prove by recursion that $u^{\star}\in\bigotimes_{\alpha\in\mathcal{L}(T)}V_{\alpha}=V$. Also, for any $\beta\in A$, Lemma 5.2 implies that $u^{\star}\in U_{\beta}\otimes\mathcal{H}_{\beta^{c}}$. Therefore, $U^{min}_{\beta}(u^{\star})\subset U_{\beta}$, and $\mathrm{rank}_{\beta}(u^{\star})\leq\dim(U_{\beta})=r_{\beta}$. This proves that $u^{\star}\in\mathcal{T}_{r}^{A}(\mathcal{H})$. ∎

For any level $\ell$, $1\leq\ell\leq L$, let $\mathcal{P}_{T_{\ell}}=\prod_{\alpha\in T_{\ell}}\mathcal{P}_{U_{\alpha}}$ be the orthogonal projection from $\mathcal{H}$ onto $U_{T_{\ell}}\otimes\mathcal{H}_{t_{\ell}^{c}}$, with $U_{T_{\ell}}=\bigotimes_{\alpha\in T_{\ell}}U_{\alpha},$ and let

[TABLE]

with the convention $u^{L+1}=u.$

Lemma 5.4**.**

For all $1\leq\ell<\ell^{\prime}\leq L$, we have

[TABLE]

Proof.

For $1\leq\ell<L$, we deduce from Lemma 5.2 that

[TABLE]

and then $U_{T_{\ell}}\otimes\mathcal{H}_{t_{\ell}^{c}}\subset U_{T_{\ell+1}}\otimes\mathcal{H}_{t_{\ell+1}^{c}}.$ Therefore, for $1\leq\ell<\ell^{\prime}\leq L$, we have $U_{T_{\ell}}\otimes\mathcal{H}_{t_{\ell}^{c}}\subset U_{T_{\ell^{\prime}}}\otimes\mathcal{H}_{t_{\ell^{\prime}}^{c}},$ and the result follows from Proposition 2.2. ∎

From Lemma 5.4, we have that

[TABLE]

for $1\leq\ell\leq L$, and

[TABLE]

We now state the two main results about the proposed algorithm.

Theorem 5.5**.**

For a given $r$, the approximation $u^{\star}\in\mathcal{T}_{r}^{A}(\mathcal{H})\cap V$ satisfies

[TABLE]

Proof.

We first note that for all $1\leq\ell<\ell^{\prime}\leq L$, $u^{\ell}-u^{\ell+1}$ is orthogonal to $u^{\ell^{\prime}}-u^{\ell^{\prime}+1}$. Indeed, using Lemma 5.4, we obtain that

[TABLE]

Then, we have

[TABLE]

From Lemma 5.4, we know that $u^{\ell+1}=\mathcal{P}_{T_{\ell+1}}u$, where we use the convention $\mathcal{P}_{T_{L+1}}=id$. For $\alpha\in\mathcal{L}(T_{\ell})$, since $\mathcal{P}_{U_{\alpha}}$ and $\mathcal{P}_{T_{\ell+1}}$ commute, $\|u^{\ell+1}-\mathcal{P}_{U_{\alpha}}u^{\ell+1}\|=\|\mathcal{P}_{T_{\ell+1}}u-\mathcal{P}_{U_{\alpha}}\mathcal{P}_{T_{\ell+1}}u\|=\|\mathcal{P}_{T_{\ell+1}}(u-\mathcal{P}_{U_{\alpha}}u)\|\leq\|u-\mathcal{P}_{U_{\alpha}}u\|.$ Therefore, for $\alpha\in\mathcal{L}(T_{\ell})\setminus\mathcal{L}(A)$, we have $\|u^{\ell+1}-\mathcal{P}_{U_{\alpha}}u^{\ell+1}\|\leq\|u-\mathcal{P}_{V_{\alpha}}u\|$, and for $\alpha\in\mathcal{L}(A_{\ell})$, we have $\|u^{\ell+1}-\mathcal{P}_{U_{\alpha}}u^{\ell+1}\|^{2}\leq\|u-\mathcal{P}_{V_{\alpha}}u\|^{2}+\|\mathcal{P}_{V_{\alpha}}u-\mathcal{P}_{U_{\alpha}}u\|^{2}=\|u-\mathcal{P}_{V_{\alpha}}u\|^{2}+\|u_{\alpha}-\mathcal{P}_{U_{\alpha}}u_{\alpha}\|^{2}$. For $\alpha\in A_{\ell}\setminus\mathcal{L}(A)$, we have

[TABLE]

so that

[TABLE]

Gathering the above results, we obtain

[TABLE]

For $\alpha\in A$, we let $U_{\alpha}^{\star}$ be the subspace in $\mathcal{H}_{\alpha}$ such that

[TABLE]

For $\alpha\in\mathcal{L}(A)$, we have $u_{\alpha}=\mathcal{P}_{V_{\alpha}}u$. From Proposition 3.10, we know that $\mathrm{rank}_{\alpha}(\mathcal{P}_{V_{\alpha}}\mathcal{P}_{U_{\alpha}^{\star}}u)\leq\mathrm{rank}_{\alpha}(\mathcal{P}_{U_{\alpha}^{\star}}u)\leq r_{\alpha}$. The optimality of $U_{\alpha}$ then implies that

[TABLE]

Now consider $\alpha\notin A\setminus\mathcal{L}(A)$. We know that $\mathrm{rank}_{\alpha}(\prod_{\beta\in S(\alpha)}\mathcal{P}_{U_{\beta}}\mathcal{P}_{U_{\alpha}^{\star}}u)\leq\mathrm{rank}_{\alpha}(\mathcal{P}_{U_{\alpha}^{\star}}u)\leq r_{\alpha}$ from Proposition 3.10. The optimality of $U_{\alpha}$ then implies that

[TABLE]

Finally, we obtain

[TABLE]

which ends the proof. ∎

Theorem 5.6**.**

For any $\epsilon\geq 0$, if for all $\alpha\in A$, the rank $r_{\alpha}$ is chosen such that

[TABLE]

the approximation $u^{\star}$ satisfies

[TABLE]

Proof.

Starting from (10), we obtain

[TABLE]

and the result follows from $\|u_{\alpha}\|=\|\prod_{\beta\in S(\alpha)}\mathcal{P}_{U_{\alpha}}u\|\leq\|u\|$ if $\alpha\notin A\setminus\mathcal{L}(A)$, and $\|u_{\alpha}\|=\|\mathcal{P}_{V_{\alpha}}u\|\leq\|u\|$ if $\alpha\in\mathcal{L}(A)$. ∎

6 Empirical principal component analysis for tree-based tensor format

6.1 Empirical principal component analysis of multivariate functions

Here we present the empirical principal component analysis for the statistical estimation of $\alpha$-principal subspaces of a multivariate function (see Section 5.1). We consider that $\mathcal{H}=L^{2}_{\mu}(\mathcal{X})$ or that $\mathcal{H}$ is a separable reproducing kernel Hilbert space compactly embedded in $L^{2}_{\mu}(\mathcal{X})$, equipped with the natural norm in $L^{2}_{\mu}(\mathcal{X})$. Let $(X_{\alpha},X_{\alpha^{c}})$ be the random vector with values in $\mathcal{X}_{\alpha}\times\mathcal{X}_{\alpha^{c}}$ with probability law $\mu_{\alpha}\otimes\mu_{\alpha^{c}}$. The tensor $u$ can be identified with a random variable defined on $\mathcal{X}_{\alpha^{c}}$ with values in $\mathcal{H}_{\alpha}$ which associates to $x_{\alpha^{c}}\in\mathcal{X}_{\alpha^{c}}$ the function $u(\cdot,x_{\alpha^{c}})=(id_{\alpha}\otimes\delta_{x_{\alpha^{c}}})u,$ this random variable being an element of the Bochner space $L^{2}_{\mu_{\alpha^{c}}}(\mathcal{X}_{\alpha^{c}};\mathcal{H}_{\alpha})$. Then problem (7) is equivalent to find a $r_{\alpha}$-dimensional subspace in $\mathcal{H}_{\alpha}$ solution of

[TABLE]

Given a set $\{x_{\alpha^{c}}^{k}\}_{k=1}^{m_{\alpha}}$ of $m_{\alpha}$ samples of $X_{\alpha^{c}}$, the $\alpha$-principal subspace can be estimated by an empirical $\alpha$-principal subspace $\widehat{U}_{\alpha}$ solution of

[TABLE]

where

[TABLE]

The problem is equivalent to finding the $r_{\alpha}$ left principal components of $\{u(\cdot,x_{\alpha^{c}}^{k})\}_{k=1}^{m_{\alpha}}$, which is identified with an order-two tensor in $\mathcal{H}_{\alpha}\otimes\mathbb{R}^{m_{\alpha}}$. We note that the number of samples $m_{\alpha}$ must be such that $m_{\alpha}\geq r_{\alpha}$ in order to estimate $r_{\alpha}$ principal components. In the case of i.i.d. samples, the semi-norm $\|\cdot\|_{\alpha,m_{\alpha}}$ on $\mathcal{H}$ is the natural statistical estimation of the Bochner norm $\|\cdot\|_{\alpha}$ in $L^{2}_{\mu_{\alpha^{c}}}(\mathcal{X}_{\alpha^{c}};\mathcal{H}_{\alpha})$, defined by $\|v\|_{\alpha}^{2}=\mathbb{E}(\|v(X_{\alpha^{c}})\|_{\mathcal{H}_{\alpha}})$. This norm $\|\cdot\|_{\alpha}$ coincides with the norm $\|\cdot\|$ on $\mathcal{H}$ when $\mathcal{H}$ is equipped with the $L^{2}_{\mu}(\mathcal{X})$-norm.555Note that when $\mathcal{H}$ is equipped with a norm stronger than the norm in $L^{2}_{\mu}(\mathcal{X})$, then $\|\cdot\|_{\alpha}$ does not coincides with the norm $\|\cdot\|$ on $\mathcal{H}$, so that the subspaces solutions of (7) and (11) will be different in general.

For some results on the comparison between $\|u-\mathcal{P}_{\widehat{U}_{\alpha}}u\|$ and the best approximation error $\|u-\mathcal{P}_{U_{\alpha}^{\star}}u\|,$ see [3, 40, 23, 24]. Under suitable assumptions on $u$ (e.g., $u$ uniformly bounded), for any $\eta>0$ and $\epsilon>0$, there exists a $m_{\alpha}$ sufficiently large (depending on $\eta,\epsilon,r_{\alpha}$ and $u$) such that

[TABLE]

holds with probability higher than $1-\eta$. Then, for any $\tau>0$, there exists a $m_{\alpha}$ sufficiently large (depending on $\eta,\tau,r_{\alpha}$ and $u$) such that

[TABLE]

holds with probability higher than $1-\eta$.

6.2 Empirical principal component analysis for tree-based format

Now we propose a modification of the algorithm proposed in Section 5.2 using only evaluations of the function $u$ at some selected points in $\mathcal{X}$. It is based on the construction of a hierarchy of subspaces $\{U_{\alpha}\}_{\alpha\in A}$, from empirical principal component analysis, and a corresponding hierarchy of commuting interpolation operators associated with nested sets of points.

For each leaf node $\alpha\in\mathcal{L}(T)$, we introduce a finite-dimensional approximation space $V_{\alpha}\subset\mathcal{H}_{\alpha}$ with dimension $\dim(V_{\alpha})=n_{\alpha}$, we introduce a set $\Gamma_{V_{\alpha}}$ of points in $\mathcal{X}_{\alpha}$ which is unisolvent for $V_{\alpha}$, we denote by $I_{V_{\alpha}}$ the corresponding interpolation operator from $\mathcal{H}_{\alpha}$ to $V_{\alpha}$, and we let $\mathcal{I}_{V_{\alpha}}=I_{V_{\alpha}}\otimes id_{\alpha^{c}}$ be the corresponding oblique projection from $\mathcal{H}$ to $V_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$. We let $V=\bigotimes_{\alpha\in\mathcal{L}(T)}V_{\alpha}\subset\mathcal{H}$. For each non active $\alpha\in\mathcal{L}(T)\setminus A$, we let $U_{\alpha}=V_{\alpha}$ and $\Gamma_{U_{\alpha}}=\Gamma_{V_{\alpha}}$.

The algorithm then goes through all active nodes of the tree, going from the leaves to the root.

For each active node $\alpha\in A$, we let

[TABLE]

where for $\alpha\notin\mathcal{L}(A)$, the space $V_{\alpha}$ is defined by

[TABLE]

where the $U_{\beta}$, $\beta\in S(\alpha)$, have been determined at a previous step. For $\alpha\notin\mathcal{L}(A)$, $\mathcal{I}_{V_{\alpha}}=I_{V_{\alpha}}\otimes id_{\alpha^{c}}$, where $I_{V_{\alpha}}$ is the interpolation operator onto $V_{\alpha}=\bigotimes_{\beta\in S(\alpha)}U_{\beta}$ associated with the product grid $\Gamma_{V_{\alpha}}=\bigtimes_{\beta\in S(\alpha)}\Gamma_{U_{\beta}}$, where each $\Gamma_{V_{\beta}}$ have been determined at a previous step. Then we determine a $r_{\alpha}$-dimensional empirical $\alpha$-principal subspace $U_{\alpha}$ of $u_{\alpha}$, which is solution of

[TABLE]

where

[TABLE]

and where $\{x_{\alpha^{c}}^{k}\}_{k=1}^{m_{\alpha}}$ are $m_{\alpha}$ random samples of $X_{\alpha^{c}}$, with $m_{\alpha}\geq r_{\alpha}$. The problem is equivalent to finding the $r_{\alpha}$ left principal components of $\{u_{\alpha}(\cdot,x_{\alpha^{c}}^{k})\}_{k=1}^{m_{\alpha}}$, which is identified with an order two tensor in $V_{\alpha}\otimes\mathbb{R}^{m_{\alpha}}$. The number of evaluations of the function $u$ for computing $U_{\alpha}$ is $m_{\alpha}\times{\dim(V_{\alpha})}$. We let $\{\varphi^{\alpha}_{k}\}_{k=1}^{r_{\alpha}}$ be the set of principal components, such that $U_{\alpha}=\mathrm{span}\{\varphi^{\alpha}_{k}\}_{k=1}^{r_{\alpha}}$. We then construct a set of points $\Gamma_{U_{\alpha}}$ which is unisolvent for $U_{\alpha}$, and such that

[TABLE]

For the practical construction of the set $\Gamma_{U_{\alpha}}$, we use the procedure described in Section 2.2.1. We denote by $I_{U_{\alpha}}$ the interpolation operator from $\mathcal{H}_{\alpha}$ onto $U_{\alpha}$ associated with the grid $\Gamma_{U_{\alpha}}$, and we let $\mathcal{I}_{U_{\alpha}}=I_{U_{\alpha}}\otimes id_{\alpha^{c}}$ be the corresponding projection from $\mathcal{H}$ onto $U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$.

Finally, we compute

[TABLE]

where $\mathcal{I}_{V_{D}}=\bigotimes_{\beta\in S(D)}I_{U_{\beta}}$ is the interpolation operator from $\mathcal{H}$ onto $V_{D}=\bigotimes_{\beta\in S(D)}U_{\beta},$ associated with the product grid $\Gamma_{V_{D}}=\bigtimes_{\beta\in S(D)}\Gamma_{U_{\beta}}$.

6.3 Analysis of the algorithm

Let us first prove that the algorithm produces an approximation $u^{\star}$ in the desired tensor format.

Lemma 6.1**.**

For $\alpha\in\mathcal{L}(T)\setminus A$, $U_{\alpha}=V_{\alpha}$. For $\alpha\in\mathcal{L}(A)$, $U_{\alpha}\subset V_{\alpha}$. For $\alpha\in A\setminus\mathcal{L}(A)$,

[TABLE]

Proof.

For $\alpha\in A$, we have $U_{\alpha}\subset U^{min}_{\alpha}(u_{\alpha})$. If $\alpha\in\mathcal{L}(A)$, we have $U^{min}_{\alpha}(u_{\alpha})\subset V_{\alpha}$ since $u_{\alpha}=\mathcal{I}_{V_{\alpha}}u$. If $\alpha\in A\setminus\mathcal{L}(A)$, we have $U^{min}_{\alpha}(u_{\alpha})\subset\bigotimes_{\beta\in S(\alpha)}U^{min}_{\beta}(u_{\alpha})$, and $U^{min}_{\beta}(u_{\alpha})\subset U_{\beta}$ since $u_{\alpha}=\prod_{\beta\in S(\alpha)}\mathcal{I}_{U_{\beta}}u.$ ∎

Proposition 6.2**.**

The algorithm produces an approximation

[TABLE]

Proof.

Since $u^{\star}={\mathcal{I}_{V_{D}}u}$, we have $u^{\star}\in{V_{D}=}\bigotimes_{\alpha\in S(D)}U_{\alpha}$. Then using Lemma 6.1, we prove by recursion that $u^{\star}\in\bigotimes_{\alpha\in\mathcal{L}(T)}V_{\alpha}=V$. Also, for any $\alpha\in A$, Lemma 6.1 implies that $u^{\star}\in U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$. Therefore, $U^{min}_{\alpha}(u^{\star})\subset U_{\alpha}$, and $\mathrm{rank}_{\alpha}(u^{\star})\leq\dim(U_{\alpha})=r_{\alpha}$. This proves that $u^{\star}\in\mathcal{T}_{r}^{A}(\mathcal{H})$. ∎

For all $\alpha\in T$, the operator $\mathcal{I}_{V_{\alpha}}=I_{V_{\alpha}}\otimes id_{\alpha^{c}}$ is a projection from $\mathcal{H}$ onto $V_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ along $W_{\alpha}^{\star}\otimes\mathcal{H}^{*}_{\alpha^{c}},$ with $W_{\alpha}^{\star}=\mathrm{span}\{\delta_{x}:x\in\Gamma_{V_{\alpha}}\}$. For all $\alpha\in T\setminus\{D\}$, the operator $\mathcal{I}_{U_{\alpha}}=I_{U_{\alpha}}\otimes id_{\alpha^{c}}$ is an oblique projection from $\mathcal{H}$ onto $U_{\alpha}\otimes\mathcal{H}_{\alpha^{c}}$ along $W_{\alpha}\otimes\mathcal{H}^{*}_{\alpha^{c}},$ with $W_{\alpha}=\mathrm{span}\{\delta_{x}:x\in\Gamma_{U_{\alpha}}\}$. From the property (14) of the grids, we deduce the following result.

Lemma 6.3**.**

For $\alpha\in\mathcal{L}(T)\setminus A$, $W_{\alpha}=W_{\alpha}^{\star}$. For $\alpha\in\mathcal{L}(A)$, $W_{\alpha}\subset W_{\alpha}^{\star}$. For $\alpha\in A\setminus\mathcal{L}(A)$,

[TABLE]

Remark 6.4**.**

Note that interpolation operators $I_{U_{\alpha}}$, $\alpha\in A$, could be replaced by oblique projections $P^{W_{\alpha}}_{U_{\alpha}}$ onto $U_{\alpha}$ along subspaces $W_{\alpha}$ in $\mathcal{H}_{\alpha}^{*}$, with subspaces $W_{\alpha}$ satisfying for $\alpha\notin\mathcal{L}(T)$, $W_{\alpha}\subset\bigotimes_{\beta\in S(\alpha)}W_{\beta}$. Under this condition, all results of this section remain valid.

For any level $\ell$, $1\leq\ell\leq L$, let

[TABLE]

where $I_{U_{T_{\ell}}}=\bigotimes_{\alpha\in T_{\ell}}I_{U_{\alpha}}$ is the interpolation operator from $\mathcal{H}_{t_{\ell}}$ to $U_{T_{\ell}}=\bigotimes_{\alpha\in T_{\ell}}U_{\alpha}$ associated with the tensor product grid $\Gamma^{T_{\ell}}=\bigtimes_{\alpha\in T_{\ell}}\Gamma^{\alpha},$ and let

[TABLE]

with the convention $u^{L+1}=u.$ We then prove that operators $\mathcal{I}_{T_{\ell}}$, $1\leq\ell\leq L$, are commuting oblique projections.

Lemma 6.5**.**

For all $1\leq\ell\leq L$, the operator $\mathcal{I}_{T_{\ell}}$ is an oblique projection from $\mathcal{H}$ to $\mathcal{U}_{\ell}:=U_{T_{\ell}}\otimes\mathcal{H}_{t_{\ell}^{c}}$ along $\mathcal{W}_{\ell}:=W_{T_{\ell}}\otimes\mathcal{H}_{t_{\ell}^{c}}^{*}$. For all $1\leq\ell<\ell^{\prime}\leq L$, we have $\mathcal{U}_{\ell}\subset\mathcal{U}_{\ell^{\prime}}$ and $\mathcal{W}_{\ell}\subset\mathcal{W}_{\ell^{\prime}}$, and therefore

[TABLE]

Proof.

For $1\leq\ell<L$, we have

[TABLE]

From Lemma 6.1, we know that $\bigotimes_{\alpha\in T_{\ell}\setminus\mathcal{L}(T)}U_{\alpha}$ is a subspace of $\bigotimes_{\beta\in T_{\ell+1}}U_{\beta}\subset\mathcal{H}_{t_{\ell+1}}.$ Therefore, we obtain $\mathcal{U}_{\ell}\subset\mathcal{U}_{\ell+1}$. In the same way, using Lemma 6.3, we obtain that $\mathcal{W}_{\ell}\subset\mathcal{W}_{\ell+1}$. We then deduce $\mathcal{I}_{T_{\ell}}\mathcal{I}_{T_{\ell+1}}=\mathcal{I}_{T_{\ell+1}}\mathcal{I}_{T_{\ell}}=\mathcal{I}_{T_{\ell}}$ from Proposition 2.2, which ends the proof. ∎

Lemma 6.6**.**

The approximation $u^{\star}$ satisfies

[TABLE]

where $\delta=\max_{\ell}\delta_{T_{\ell}}$ and

[TABLE]

with $t_{\ell}=\cup_{\alpha\in T_{\ell}}\alpha$. If $u\in V$, then

[TABLE]

with $a_{\ell}=\cup_{\alpha\in A_{\ell}}\alpha.$

Proof.

Since $u-u^{\star}=\sum_{\ell=1}^{L}(u^{\ell+1}-u^{\ell})$, we have

[TABLE]

For $\ell^{\prime}<\ell$, since $\mathcal{P}_{T_{\ell}}(u^{\ell^{\prime}+1}-u^{\ell^{\prime}})=u^{\ell^{\prime}+1}-u^{\ell^{\prime}}$, we have

[TABLE]

where we have used the fact that $\mathcal{P}_{T_{\ell}}\mathcal{I}_{T_{\ell}}=\mathcal{I}_{T_{\ell}}$ and $(\mathcal{P}_{T_{\ell}}-\mathcal{I}_{T_{\ell}})u^{\ell}=0.$ Since $\mathcal{P}_{T_{\ell}}-\mathcal{I}_{T_{\ell}}=(P_{U_{T_{\ell}}}-I_{U_{T_{\ell}}})\otimes id_{t_{\ell}^{c}}$ and $u^{\ell+1}-u^{\ell}=u^{\ell+1}-\mathcal{I}_{T_{\ell}}u^{\ell+1}\subset U^{min}_{T_{\ell}}(u^{\ell+1})\otimes\mathcal{H}_{t_{\ell}^{c}}$, we obtain from Proposition 3.11 that

[TABLE]

for $\ell^{\prime}<\ell$. We deduce that

[TABLE]

where the matrix $B\in\mathbb{R}^{L\times L}$ is such that $B_{\ell,\ell}=1$ and $B_{\ell,\ell^{\prime}}=\delta_{T_{\max\{\ell,\ell^{\prime}\}}}$ if $\ell\neq\ell^{\prime}$. Using the theorem of Gerschgorin, we have that

[TABLE]

with $\delta=\max_{\ell}\delta_{T_{\ell}}$.

Finally, when $u\in V$, we have $U^{min}_{\alpha}(u^{\ell+1})\subset U^{min}_{\alpha}(u)\subset V_{\alpha}$ for all $\alpha\in\mathcal{L}(T).$ Therefore, $I_{V_{\alpha}}v=P_{V_{\alpha}}v$ for all $v\in U^{min}_{\alpha}(u^{\ell+1})$ and $\alpha\in\mathcal{L}(T)$, and

[TABLE]

∎

Lemma 6.7**.**

For $1\leq\ell\leq L$,

[TABLE]

where for $S\subset T$,

[TABLE]

[TABLE]

and

[TABLE]

for $\alpha\in\mathcal{L}(T)$. Moreover, if $u\in V$, then $\delta_{\alpha}=0$ for all $\alpha\in\mathcal{L}(T)$, and $a_{\alpha}=0$ for all $\alpha\in T$.

Proof.

For all $1\leq\ell\leq L$, we have

[TABLE]

For $\alpha\in T_{\ell}\setminus\mathcal{L}(T)=A_{\ell}\setminus\mathcal{L}(T)$,

[TABLE]

and since $\mathcal{P}_{U_{\alpha}}$ and $\prod_{\delta\in T_{\ell+1}\setminus S(\alpha)}\mathcal{I}_{U_{\delta}}$ commute, we have

[TABLE]

Now for $\alpha\in T_{\ell}\cap\mathcal{L}(T)$, we have that $\mathcal{P}_{U_{\alpha}}$ and $\mathcal{I}_{T_{\ell+1}}$ commute, and therefore

[TABLE]

If $\alpha\in T_{\ell}\setminus A_{\ell}$, we have $U_{\alpha}=V_{\alpha}$. If $\alpha\in A_{\ell}\cap\mathcal{L}(T)$, we have

[TABLE]

so that if $\delta_{\alpha}=\|I_{V_{\alpha}}-P_{V_{\alpha}}\|_{U^{min}_{\alpha}(u)\to\mathcal{H}_{\alpha}}=0$, we have $\mathcal{P}_{V_{\alpha}}u=\mathcal{I}_{V_{\alpha}}u=u_{\alpha}$ and

[TABLE]

and if $\delta_{\alpha}\neq 0$, we have

[TABLE]

where we have used Proposition 3.11. We conclude from (18) and (19) that if $\alpha\in A_{\ell}\cap\mathcal{L}(T)$,

[TABLE]

Gathering the above results, we obtain

[TABLE]

which ends the proof.

∎

We now state the two main results about the proposed algorithm.

Theorem 6.8**.**

Assume that for all $\alpha\in A$, the subspace $U_{\alpha}$ is such that

[TABLE]

holds with probability higher than $1-\eta$, for some $\tau\geq 1$. Then the approximation $u^{\star}\in\mathcal{T}_{r}^{A}(\mathcal{H})\cap V$ is such that

[TABLE]

holds with probability higher than $1-\#A\eta$, where $C$ is defined by

[TABLE]

with

[TABLE]

and $a_{\alpha}$ and $\delta_{\alpha}$ defined by (16) and (17) respectively, and where $D_{\alpha}$ is defined by

[TABLE]

for $\alpha\in\mathcal{L}(T)\cap T_{\ell}$.

Proof.

For $\alpha\in A$, let $\widehat{U}_{\alpha}$ be a subspace such that

[TABLE]

and let $U_{\alpha}^{\star}\subset U^{min}_{\alpha}(u)$ be a subspace such that

[TABLE]

For $\alpha\in\mathcal{L}(A)$, we have $u_{\alpha}=\mathcal{I}_{V_{\alpha}}u$. We know that $\mathrm{rank}_{\alpha}(\mathcal{I}_{V_{\alpha}}\mathcal{P}_{U_{\alpha}^{\star}}u)\leq r_{\alpha}$ from Proposition 3.10. By the optimality of $\widehat{U}_{\alpha}$, we obtain

[TABLE]

Now consider $\alpha\in A\setminus\mathcal{L}(A)$. We know that $\mathrm{rank}_{\alpha}(\prod_{\beta\in S(\alpha)}\mathcal{I}_{U_{\beta}}\mathcal{P}_{U_{\alpha}^{\star}}u)\leq\mathrm{rank}_{\alpha}(\mathcal{P}_{U_{\alpha}^{\star}}u)\leq r_{\alpha}$ from Proposition 3.10. By the optimality of $\widehat{U}_{\alpha}$, we obtain

[TABLE]

Then, using Lemma 6.7 and assumption (20), we obtain

[TABLE]

Then, using Lemma 6.6, we obtain (21). ∎

Remark 6.9**.**

Assume $u\in V$ (no discretization). Then $\delta_{\alpha}=0$ and $\|I_{V_{\alpha}}\|_{U^{min}_{\alpha}(u)\to\mathcal{H}_{\alpha}}=1$ for all $\alpha\in\mathcal{L}(T)$, $a_{\alpha}=0$ and $\lambda_{\alpha}=1$ for all $\alpha\in T$, $\Lambda_{T_{\ell}}=\Lambda_{A_{\ell}}$ and $\delta_{T_{\ell}}=\delta_{A_{\ell}}$ for all $\ell$. Also, the constant $C$ defined by (22) is such that

[TABLE]

Moreover, if $U_{\alpha}=U^{min}_{\alpha}(u)$ for all $\alpha$, then $\Lambda_{T_{\ell}}=\Lambda_{A_{\ell}}=1$ and $\delta_{T_{\ell}}=\delta_{A_{\ell}}=0$ for all $\ell$, which implies

[TABLE]

Theorem 6.10**.**

Let $\epsilon,{\tilde{\epsilon}}\geq 0$. Assume that for all $\alpha\in A$, the subspace $U_{\alpha}$ is such that

[TABLE]

holds with probability higher than $1-\eta$, and further assume that the subspaces $V_{\alpha}$, $\alpha\in\mathcal{L}(T)$, are such that

[TABLE]

Then the approximation $u^{\star}$ is such that

[TABLE]

holds with probability higher than $1-\#A\eta$, where $C$ is defined by (22) and where $D^{2}=\sum_{\alpha\in\mathcal{L}(T)}D_{\alpha}^{2}$, with $D_{\alpha}$ defined by (24), is such that

[TABLE]

Proof.

We first note that for $\alpha\in A\setminus\mathcal{L}(A),$ we have $\|u_{\alpha}\|\leq\Lambda_{S(\alpha)}\|u\|.$ Also, for $\alpha\in\mathcal{L}(T)$, we have $\|u_{\alpha}\|\leq\lambda_{\alpha}\|u\|,$ with $\lambda_{\alpha}$ defined in (23). Using Lemma 6.7 and assumptions (27) and (28), we then obtain

[TABLE]

Finally, we obtain the desired result by using Lemma 6.6. ∎

Example 6.11**.**

For the Tucker format described in Example 4.3, the constants $C$ and $D$ are given by

[TABLE]

[TABLE]

with

[TABLE]

If $u\in V$, then

[TABLE]

Example 6.12**.**

For the tensor train format described in Example 4.5, the constant $C$ and $D$ are given by

[TABLE]

[TABLE]

with

[TABLE]

If $u\in V$, then

[TABLE]

Example 6.13**.**

For the tensor train Tucker format described in Example 4.6, the constant $C$ and $D$ are given by

[TABLE]

[TABLE]

If $u\in V$, then

[TABLE]

6.4 Complexity

Here we analyse the complexity of the algorithm in terms of the number of evaluations of the function. Evaluations of the function $u$ are required (i) for the computation of the subspaces $\{U_{\alpha}\}_{\alpha\in A}$ through empirical principal component analysis of the ${V_{\alpha}}$-valued functions $u_{\alpha}(\cdot,X_{\alpha^{c}})$, with $V_{\alpha}$ a given approximation space if $\alpha\in\mathcal{L}(A)$ or ${V_{\alpha}}=\bigotimes_{\beta\in S(\alpha)}U_{\beta}$ if $\alpha\in A\setminus\mathcal{L}(A)$, and (ii) for the computation of the final interpolation ${\mathcal{I}_{V_{D}}}u$.

We then obtain the following result about the number of evaluations of the function required by the algorithm

Proposition 6.14**.**

The total number of evaluations of $u$ required by the algorithm for computing an approximation $u^{\star}$ in the tensor format $\mathcal{T}_{r}^{A}(V)$ is

[TABLE]

where $n=(n_{\alpha})_{\alpha\in\mathcal{L}(T)}$, with $n_{\alpha}=\dim(V_{\alpha}),$ and ${m}=(m_{\alpha})_{\alpha\in A}$, with $m_{\alpha}$ the number of samples of the $Z_{\alpha}$-valued random variable $u_{\alpha}(\cdot,X_{\alpha^{c}})$ used for computing $U_{\alpha}$.

Proof.

For $\alpha\in A$, the function $u_{\alpha}$ is an interpolation of $u$ in $Z_{\alpha}=V_{\alpha}$ if $\alpha\in\mathcal{L}(A)$, or in $Z_{\alpha}=\bigotimes_{\beta\in S(\alpha)}U_{\beta}=\left(\bigotimes_{\beta\in S(\alpha)\cap A}U_{\beta}\right)\otimes\left(\bigotimes_{\beta\in S(\alpha)\setminus A}V_{\beta}\right)$ if $\alpha\notin\mathcal{L}(A).$ Therefore, computing $u_{\alpha}(\cdot,x_{\alpha^{c}}^{k})$ for one realization $x^{k}_{\alpha^{c}}$ of $X_{\alpha^{c}}$ requires $\dim(V_{\alpha})=n_{\alpha}$ evaluations of $u$ if $\alpha\in\mathcal{L}(A)$ or $\dim(\bigotimes_{\beta\in S(\alpha)}U_{\beta})=\prod_{\beta\in S(\alpha)\cap A}r_{\beta}\prod_{\beta\in S(\alpha)\setminus A}n_{\beta}$ if $\alpha\notin\mathcal{L}(A)$. Finally, the computation of the interpolation $\mathcal{I}_{T_{1}}u=\mathcal{I}_{S(D)}u$ requires $\dim(\bigotimes_{\alpha\in S(D)}U_{\alpha})=\prod_{\beta\in S(D)\cap A}r_{\beta}\prod_{\beta\in S(D)\setminus A}n_{\beta}$ evaluations of $u$. ∎

For computing a $r_{\alpha}$-dimensional subspace $U_{\alpha}$, the number of samples $m_{\alpha}$ of $u_{\alpha}(\cdot,X_{\alpha^{c}})$ has to be at least $r_{\alpha}$.

Corollary 6.15**.**

If the number of samples $m_{\alpha}=r_{\alpha}$ for all $\alpha\in A$, then the number of evaluations of the function required by the algorithm is

[TABLE]

The above result states that for a prescribed rank $r=(r_{\alpha})_{\alpha\in A}$, the algorithm is able to construct an approximation of $u$ using a number os samples equal to the storage complexity of the tensor format $\mathcal{T}_{r}^{A}(V)$.

When using the algorithm with a prescribed tolerance $\epsilon$, the rank $r_{\alpha}$ is not fixed a priori but defined as the minimal integer such that the condition (27) is satisfied. Since samples of $u_{\alpha}(\cdot,X_{\alpha^{c}})$ belongs to the subspace $U^{min}_{\alpha}(u_{\alpha})\subset Z_{\alpha}$ with dimension $\mathrm{rank}_{\alpha}(u_{\alpha})\leq\dim(Z_{\alpha})$, the selected rank $r_{\alpha}$ is at most $\dim(Z_{\alpha})$. Therefore, by taking $m_{\alpha}=\dim(Z_{\alpha})$ for all $\alpha\in A$, if we assume that the set of $m_{\alpha}$ samples of $u(\cdot,X_{\alpha^{c}})$ contains $\mathrm{rank}_{\alpha}(u_{\alpha})$ linearly independent functions in $Z_{\alpha}$, then the algorithm is able to produce an approximation with arbitrary small tolerance $\epsilon.$

Corollary 6.16**.**

If the number of samples $m_{\alpha}=\dim(Z_{\alpha})$ for all $\alpha\in A$, then

[TABLE]

Remark 6.17**.**

For numerical experiments, when working with prescribed tolerance, we will use $m_{\alpha}=\dim(Z_{\alpha})$ for al $\alpha\in A$.

7 Numerical examples

In all examples, we consider functions $u$ in the tensor space $L^{2}_{\mu}(\mathcal{X})$, with $\mathcal{X}\subset\mathbb{R}^{d}$, equipped with the natural norm $\|\cdot\|$ (see example 3.2)666For the last example, $\mathcal{X}$ is a finite product set equipped with the uniform measure and $L^{2}_{\mu}(\mathcal{X})$ then corresponds to the space of multidimensional arrays equipped with the canonical norm.. For an approximation $u^{\star}$ provided by the algorithm, we estimate the relative error $\varepsilon(u^{\star})=\|u-u^{\star}\|/\|u\|$ using Monte-Carlo integration. We denote by $M$ the total number of evaluations of the function $u$ required by the algorithm to provide an approximation $u^{\star}$, and by $S$ the storage complexity of the approximation $u^{\star}$. Since the algorithm uses random evaluations of the function $u$ (for the estimation of principal components), we run the algorithm several times and indicate confidence intervals of level $90\%$ for $\varepsilon(u^{\star})$, and also for $M$, $S$ and approximation ranks when these quantities are random.

For the approximation with a prescribed $A$-rank, we use $m_{\alpha}=\gamma r_{\alpha}$ samples for the estimation of principal subspaces $U_{\alpha}$, $\alpha\in A$. If $\gamma=1$, then $M=S$ (see corollary 6.15).

For the approximation with a prescribed tolerance $\epsilon$, we use $m_{\alpha}=\dim(Z_{\alpha})$ for all $\alpha\in A$ (see corollary 6.16 for the estimation of $M$).

In all examples except the last one, we use polynomial approximation spaces $V_{\nu}=\mathbb{P}_{p}(\mathcal{X}_{\nu})$ over $\mathcal{X}_{\nu}\subset\mathbb{R}$, $\nu\in D$, with the same polynomial degree $p$ in all dimensions. For each $\nu\in D$, we use an orthonormal polynomial basis of $V_{\nu}=\mathbb{P}_{p}(\mathcal{X}_{\nu})$ (Hermite polynomials for a Gaussian measure, Legendre polynomials for a uniform measure,…), and associated interpolation grids $\Gamma^{\star}_{\nu}$ selected in a set of 1000 random points (drawn from the measure $\mu_{\nu}$) by using the greedy algorithm described in Section 2.2.1.

7.1 Henon-Heiles potential

We consider $\mathcal{X}=\mathbb{R}^{d}$ equipped with the standard Gaussian measure $\mu$ and the modified Henon-Heiles potential [28]

[TABLE]

with $\sigma_{\star}=0.2$. We consider approximation in the tensor train format $\mathcal{T}_{r}^{A}(V)$ described in example 4.5. The function is such that $\mathrm{rank}_{\alpha}(u)=3$ for all $\alpha\in A$. We use a polynomial degree $p=4$, so that there is no discretization error, i.e. $u\in V$.

In Table 1, we observe that the algorithm with a prescribed rank $r=(3,\ldots,3)$ is able to recover the function at very high precision with high probability with a number of samples equal to the storage complexity of the approximation (when $\gamma=1$), with no deterioration when the dimension $d$ increases from $5$ to $100$. The accuracy is slightly improved when $\gamma=100$ but with a much higher number of evaluations of the function.

7.2 Sine of a sum

We consider $\mathcal{X}=[-1,1]^{d}$ equipped with the uniform measure and the function

[TABLE]

We consider approximation in the tensor train Tucker format $\mathcal{T}_{r}^{A}(V)$ described in example 4.6. The function is such that $\mathrm{rank}_{\alpha}(u)=2$ for all $\alpha\in A$. In Table 2, we observe the behavior of the algorithm with a prescribed rank $r=(2,\ldots,2)$ for different polynomial degrees $p$ and different values of $d$. We observe a linear dependence of the complexity with respect to $d$.

In Table 3, we observe the behavior of the algorithm with prescribed tolerance $\epsilon=10^{-12}$ and fixed polynomial degree $p=17$, for different values of $d$. For this value of $\epsilon$, the algorithm always provides an approximation with rank $(2,\ldots,2)$ with a fixed number of evaluations which is about ten times the storage complexity.

7.3 Sum of bivariate functions

We consider $\mathcal{X}=[-1,1]^{d}$ equipped with the uniform measure and the function

[TABLE]

where $g$ is a bivariate function, and $d=10$. We consider approximation in the tensor train Tucker format $\mathcal{T}_{r}^{A}(V)$ described in example 4.6. The function is such that $\mathrm{rank}_{\{\nu\}}(u)=\mathrm{rank}(g)+1$ for all $\nu\in D$, and $\mathrm{rank}_{\{1,\ldots,\nu\}}(u)=2$ if $\nu$ is even, or $\mathrm{rank}_{\{1,\ldots,\nu\}}(u)=\mathrm{rank}(g)+1$ if $\nu$ is odd. Here, we use the algorithm we a prescribed tolerance $\epsilon$.

We first consider the function $g(y,z)=\sum_{j=0}^{3}y^{j}z^{j}$ whose rank is $4$ and we use polynomial spaces of degree $p=5$, so that there is no discretization error. We observe in Table 4 the behavior of the algorithm for decreasing values of $\epsilon$. For $\epsilon=10^{-4}$, the algorithm always provides the solution at almost machine precision, with an exact recovery of the rank of the function $u$. We observe that increasing $\gamma$ (i.e. the number of evaluations for the estimation of principal components) allows us to obtain a more accurate approximation for a given prescribed tolerance but with a significant increase in the number of evaluations.

We now consider the function $g(y,z)=\exp^{-\frac{1}{8}(y-z)^{2}}$ with infinite rank. We observe in Tables 5 and 6 the behavior of the algorithm for decreasing values of $\epsilon$, and for a fixed polynomial degree $p=10$ in Table 5, and an adaptive polynomial degree $p(\epsilon)=\log_{10}(\epsilon^{-1})$ in Table 6. We observe that the relative error of the obtained approximation is below the prescribed tolerance with high probability. Also, we clearly see the interest of adapting the discretization to the desired precision, which yields a lower complexity for small or moderate $\epsilon$.

7.4 Borehole function

We here consider the function

[TABLE]

which models the water flow through a borehole as a function of $8$ independent random variables $Y_{1}\sim\mathcal{N}(0.1,0.0161812)$, $Y_{2}\sim\mathcal{N}(7.71,1.0056)$, $Y_{3}\sim\mathcal{U}(63070,115600)$, $Y_{4}\sim\mathcal{U}(990,1110)$, $Y_{5}\sim\mathcal{U}(63.1,116)$, $Y_{6}\sim\mathcal{U}(700,820)$, $Y_{7}\sim\mathcal{U}(1120,1680)$, $Y_{8}\sim\mathcal{U}(9855,12045)$. We then consider the function

[TABLE]

where $g_{\nu}$ are functions such that $Y_{\nu}=g_{\nu}(X_{\nu})$, with $X_{\nu}\sim\mathcal{N}(0,1)$ for $\nu\in\{1,2\}$, and $X_{\nu}\sim\mathcal{U}(-1,1)$ for $\nu\in\{3,\ldots,8\}.$ Function $u$ is then defined on $\mathcal{X}=\mathbb{R}^{2}\times[-1,1]^{6}.$ We use polynomial approximation spaces $V_{\nu}=\mathbb{P}_{p}(\mathcal{X}_{\nu})$, $\nu\in D$. We consider approximation in the tensor train Tucker format $\mathcal{T}_{r}^{A}(V)$ described in example 4.6.

In Table 7, we observe the behavior of the algorithm with prescribed ranks $(r,\ldots,r)$ and fixed degree $p=10$. We observe a very fast convergence of the approximation with the rank. Increasing $\gamma$ (i.e. the number of evaluations for the estimation of principal components) allows us to improve the accuracy for a given rank but it we look at the error as a function of the complexity $M$, $\gamma=1$ is much better than $\gamma=100$.

In Table 8, we observe the behavior of the algorithm for decreasing values of $\epsilon$, and for an adaptive polynomial degree $p(\epsilon)=\log_{10}(\epsilon^{-1})$. We observe that for all $\epsilon,$ the relative error of the obtained approximation is below $\epsilon$ with high probability. We note that the required number of evaluations $M$ is about 2 to 4 times the storage complexity.

7.5 Tensorization of a univariate function

We consider the approximation of the univariate function $f:[0,1]\to\mathbb{R}$ using tensorization of functions [26, 39]. We denote by $f_{N}$ the piecewise constant approximation of $f$ on a uniform partition $0=t_{0}\leq t_{1}\leq\ldots\leq t_{N}=1$ with $N=2^{d}$ elements, such that $f_{N}(ih)=f(ih)$ for $0\leq i\leq N$ and $h=N^{-1}=2^{-d}.$ We denote by $v\in\mathbb{R}^{N}$ the vector with components $v(i)=f(ih)$, $0\leq i\leq N-1$. The vector $v\in\mathbb{R}^{2^{d}}$ can be identified with an order-$d$ tensor $u\in\mathcal{H}=\mathbb{R}^{2}\otimes\ldots\otimes\mathbb{R}^{2}$ such that

[TABLE]

where $(i_{1},\ldots,i_{d})\in\{0,1\}^{d}=\mathcal{X}$ is the binary representation of the integer $i\in\{0,\ldots,2^{d}-1\}$. The set $\mathcal{X}$ is equipped with the uniform measure $\mu$. Then we consider approximation of the tensor $u$ in tensor train format. The algorithm evaluates the tensor $u$ at some selected entries $(i_{1},\ldots,i_{d})$, which corresponds to evaluating the function $f$ at some particular points $t_{i}$.

In this finite-dimensional setting, we consider $V=\mathcal{H}.$ In all examples, we consider $d=40,$ and $N=2^{d}\approx 10^{12}$. This corresponds to a storage complexity of one terabyte for the standard representation of $f_{N}$ as a vector $v$ of size $N$.

We observe in Tables 9 and 10 the behavior of the algorithm with prescribed tolerance $\epsilon$ applied to the functions $f(t)=t^{2}$ and $f(t)=t^{1/2}$ respectively. We indicate relative errors in $\ell^{2}$ and $\ell^{\infty}$ norms between the tensor $u$ and the approximation $u^{\star}$. Let us recall that for $f(t)=t^{\alpha}$, the approximation error $\|f-f_{N}\|_{L^{\infty}}=O(N^{-\beta})=O(2^{-d\beta})$ with $\beta=\min\{1,\alpha\}$, which is an exponential convergence with respect to $d$. For the function $f(t)=t^{2}$, we observe that the relative error in $\ell^{2}$ norm is below the prescribed tolerance with high probability. For the function $f(t)=t^{1/2}$, the probability of obtaining a relative error in $\ell^{2}$ norm below the prescribed tolerance decreases with $\epsilon$ but the ratio between the true relative error and the prescribed tolerance remains relatively small (below 100). We note that for $f(t)=t^{2}$, the approximation ranks are bounded by $3$, which is the effective rank of $f_{N}$. For $f(t)=t^{1/2}$, the approximation ranks slowly increase with $\epsilon^{-1}$.

In both cases, we observe a very good behavior of the algorithm, which requires a number of evaluations which scales as $\log(\epsilon^{-1})$.