LEARNING HIGH-DIMENSIONAL PROBABILITY
DISTRIBUTIONS USING TREE TENSOR
NETWORKS

Erwan Grelier, Anthony Nouy, Regis Lebrun

► To cite this version:

Erwan Grelier, Anthony Nouy, Regis Lebrun. LEARNING HIGH-DIMENSIONAL PROBABILITY DISTRIBUTIONS USING TREE TENSOR NETWORKS. International Journal for Uncertainty Quantification, 2022, 12 (5), pp.47-69. 10.1615/Int.J.UncertaintyQuantification.2022034340. hal-03919932

HAL Id: hal-03919932

https://hal.science/hal-03919932v1

Submitted on 7 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Learning high-dimensional probability distributions using tree tensor networks

Erwan Grelier*

Anthony Nouy†

Régis Lebrun‡

## Abstract

We consider the problem of the estimation of a high-dimensional probability distribution from i.i.d. samples of the distribution using model classes of functions in tree-based tensor formats, a particular case of tensor networks associated with a dimension partition tree. The distribution is assumed to admit a density with respect to a product measure, possibly discrete for handling the case of discrete random variables. After discussing the representation of classical model classes in tree-based tensor formats, we present learning algorithms based on empirical risk minimization using a $L^2$ contrast. These algorithms exploit the multilinear parametrization of the formats to recast the nonlinear minimization problem into a sequence of empirical risk minimization problems with linear models. A suitable parametrization of the tensor in tree-based tensor format allows to obtain a linear model with orthogonal bases, so that each problem admits an explicit expression of the solution and cross-validation risk estimates. These estimations of the risk enable the model selection, for instance when exploiting sparsity in the coefficients of the representation. A strategy for the adaptation of the tensor format (dimension tree and tree-based ranks) is provided, which allows to discover and exploit some specific structures of high-dimensional probability distributions such as independence or conditional independence. We illustrate the performances of the proposed algorithms for the approximation of classical probabilistic models (such as Gaussian distribution, graphical models, Markov chain).

**Keywords:** density estimation, high-dimensional approximation, tensor networks, hierarchical tensor format, tensor train format, neural networks.

## 1 Introduction

The approximation of high-dimensional functions is a typical task in statistics and machine learning. We here consider the problem of the approximation of the probability distribution

<footnote>*Centrale Nantes, Laboratoire de Mathématiques Jean Leray, CNRS UMR 6629, Joint Laboratory of Marine Technology between Naval Group and Centrale Nantes</footnote>

<footnote>†Centrale Nantes, Laboratoire de Mathématiques Jean Leray, CNRS UMR 6629</footnote>

<footnote>‡Airbus Central Research & Technology, Virtual Product Engineering/XRV</footnote>
of a high-dimensional random vector $X = (X_1, \dots, X_d)$, characterized by its density $f(x)$ with respect to a measure $\mu$ (here assumed to be a product measure, e.g. the Lebesgue measure or another probability measure) or the function $f(x) = P(X = x)$ for discrete random variables (considered as a density with respect to a discrete measure). We assume that we are given independent and identically distributed samples of $X$, which is a classical setting in learning or statistics.

The approximation of a multivariate function $f(x_1, \dots, x_d)$ is a challenging problem when the dimension $d$ is large or when the available information on the distribution (evaluations of the density function or samples from the distribution) is limited. This requires to introduce model classes (or hypothesis sets) that exploit low-dimensional structures of the function. Typical model classes for high dimensional density approximation include

* multiplicative models $g^1(x_1) \cdots g^d(x_d)$, which correspond to the hypothesis that the components of $X$ are independent,

* generalized multiplicative models $\prod_{\alpha \in T} g^\alpha(x_\alpha)$, where $T$ is a collection of subsets $\alpha \subset \{1, \dots, d\}$ and $x_\alpha$ denotes the corresponding group of variables. These include bayesian networks or more general graphical models,

* mixture models $\sum_{k=1}^K \gamma_k g_k(x)$ with $\sum_{k=1}^K \gamma_k = 1$ and the $g_k$ in suitable model classes (possibly of different types). For example, a mixture of multiplicative models takes the form

$$ \sum_{k=1}^{K} \gamma_k g_k^1(x_1) \cdots g_k^d(x_d). \qquad (1) $$

Here, we consider the model classes of rank-structured functions, widely used in data analysis, signal and image processing and numerical analysis. The mixture of multiplicative models (1) is a particular case of rank-structured functions associated with the canonical notion of rank. Other classes of rank-structured approximations that have better approximation power and that are more amenable for numerical computations can be defined by considering different notions of rank. For any subset $\alpha$ of $\{1, \dots, d\} := D$, a natural notion of $\alpha$-rank of a function $g(x)$, denoted by $\text{rank}_\alpha(g)$, can be defined as the minimal integer $r_\alpha$ such that

$$ g(x) = \sum_{k=1}^{r_\alpha} g_k^\alpha(x_\alpha) g_k^{\alpha^c}(x_{\alpha^c}) $$

for some functions of two complementary groups of variables $x_\alpha$ and $x_{\alpha^c}$, $\alpha^c := D \setminus \alpha$. By considering a collection $T$ of subsets of $D$ and a tuple $r = (r_\alpha)_{\alpha \in T}$, a model class $\mathcal{T}_r^T(\mathcal{H})$ of rank-structured functions can then be defined as

$$ \mathcal{T}_r^T(\mathcal{H}) = \{g \in \mathcal{H} : \text{rank}_\alpha(g) \le r_\alpha, \alpha \in T\}, $$

where $\mathcal{H}$ is some function space. The particular case where $T$ is a dimension partition tree over $D$ corresponds to the model class of functions in tree-based tensor format [1], which
includes the Tucker format for a trivial tree, the hierarchical tensor format [2] for a balanced
binary tree, and the tensor train format [3, 4] for a linear tree. These model classes are
also known as tree tensor networks. This particular choice for $T$ provides $\mathcal{T}_r^T(\mathcal{H})$ with nice
topological properties [1] and geometrical properties [5, 6, 7, 8]. Also, and more importantly
from a practical point of view, elements of $\mathcal{T}_r^T(\mathcal{H})$ admit explicit and numerically stable
representations. The complexity of these representations is linear in $d$ and polynomial in
the ranks, with a polynomial degree depending on the arity of the tree. Model classes of
tree tensor networks have a high approximation power [9, 10, 11]. For more details on tree
tensor networks and their applications, the reader is referred to the monograph [12] and
surveys [13, 14, 15, 16, 17, 18].

The outline of the article is as follows. In Section 2, we present the model class of
multivariate functions in tree-based tensor formats, a particular case of tensor networks
associated with a dimension partition tree. In Section 3, we discuss the representation of
probabilistic models in tensor formats, with a focus on classical models such as Markov
processes, graphical models and mixtures models. Section 4 describes learning algorithms
with tree-based tensor formats in a classical empirical risk minimization setting. In such a
setting, selecting optimal tree-based ranks and optimal dimension trees (in the sense that it
gives the best performance in terms of the complexity or number of samples) are challenging
combinatorial problems for which several heuristic approaches have been proposed. Here,
we present strategies for the adaptation of ranks and tree, and a possible exploitation of
the sparsity in the tensors of the tensor networks. Finally, Section 5 presents numerical
experiments that show the performances of the proposed algorithms.

Note that tree-based tensor formats have been recently considered for the approximation
of high-dimensional probability distributions in different settings where the probability
density function is known (possibly up to some normalizing constant in Bayesian inference)
[19, 20] or solution of a variational problem [21]. Here, we consider the problem of estimating
a distribution in a different setting where i.i.d. samples from the distribution are available.
When one is interested in moments of a distribution $\mu$ or the expectation $\int fd\mu$ of some
function $f$, there is in general no significant benefit of computing these quantities from
an estimation of the distribution compared to a direct Monte-Carlo estimation of these
quantities. However, obtaining an explicit estimation of the distribution that can be easily
manipulated and simulated (that is the case with tree-based tensor formats) can be useful
for many purposes such as the generation of new (synthetic) data from the distribution
(an alternative to generative adversarial networks), the computation of level-sets of the
density (e.g., for anomaly detection) or the computation of distances or divergences between
measures.
## 2 Tree-based tensor formats

We consider the space $L_{\mu}^{2}(\mathcal{X})$ of square integrable functions defined on a product set $\mathcal{X} = \mathcal{X}_{1} \times \cdots \times \mathcal{X}_{d}$ equipped with a product measure $\mu = \mu_{1} \otimes \cdots \otimes \mu_{d}$. For $1 \le \nu \le d$, we consider spaces $\mathcal{H}_{\nu}$ of functions in $L_{\mu_{\nu}}^{2}(\mathcal{X}_{\nu})$ and the algebraic tensor space $\mathcal{H} = \mathcal{H}_{1} \otimes \cdots \otimes \mathcal{H}_{d}$ composed by functions of the form

$$ \sum_{k=1}^{r} g_{k}^{1}(x_{1}) \cdots g_{k}^{d}(x_{d}) \qquad (2) $$

for some $r \in \mathbb{N}$ and functions $g_{k}^{\nu} \in \mathcal{H}_{\nu}$. A function $g^{1}(x_{1}) \cdots g^{d}(x_{d}) := (g^{1} \otimes \cdots \otimes g^{d})(x)$ is called an elementary tensor. When spaces $\mathcal{H}_{\nu}$ are infinite dimensional, a tensor Banach space is obtained by the completion of $\mathcal{H}$ with respect to a certain norm. If $\mathcal{H}_{\nu} = L_{\mu_{\nu}}^{2}(\mathcal{X}_{\nu})$, then we obtain the space $L_{\mu}^{2}(\mathcal{X})$ after completion with respect to its natural norm. Hereafter, we assume that spaces $\mathcal{H}_{\nu}$ are equipped with the canonical norm in $L_{\mu_{\nu}}^{2}(\mathcal{X})$ and that $\mathcal{H}$ is equipped with the canonical norm in $L_{\mu}^{2}(\mathcal{X})$.

The canonical rank of a function $g \in \mathcal{H}$ is the minimal integer $r$ such that $g$ can be written in the form (2). The set of functions in $\mathcal{H}$ with canonical rank bounded by $r$ is denoted by $\mathcal{R}_{r}(\mathcal{H})$. An approximation in $\mathcal{R}_{r}(\mathcal{H})$ is called an approximation in canonical tensor format. For an order-two tensor ($d = 2$), the canonical rank coincides with the classical and unique notion of rank. For higher-order tensors ($d \ge 3$), there exist different and natural notions of rank that will be introduced for defining tree-based tensor formats.

### 2.1 $\alpha$-ranks

For a non-empty subset $\alpha$ in $\{1, \dots, d\} := D$ and its complementary subset $\alpha^{c} = D \setminus \alpha$, a tensor $g \in \mathcal{H}$ can be identified with an element $\mathcal{M}_{\alpha}(g)$ of the space of order-two tensors $\mathcal{H}_{\alpha} \otimes \mathcal{H}_{\alpha^{c}}$, where $\mathcal{H}_{\alpha} = \bigotimes_{\nu \in \alpha} \mathcal{H}_{\nu}$. This is equivalent to identifying the function $g(x)$ with a bivariate function of the complementary groups of variables $x_{\alpha} = (x_{\nu})_{\nu \in \alpha}$ and $x_{\alpha^{c}} = (x_{\nu})_{\nu \in \alpha^{c}}$ in $x$. The operator $\mathcal{M}_{\alpha}$ is called the $\alpha$-matricization operator. If $\mathrm{rank}_{\alpha}(g) = r_{\alpha}$, then $g$ admits the representation

$$ g(x) = \sum_{i=1}^{r_{\alpha}} g_{i}^{\alpha}(x_{\alpha}) g_{i}^{\alpha^{c}}(x_{\alpha^{c}}), \qquad (3) $$

for some functions $g_{i}^{\alpha} \in \mathcal{H}_{\alpha}$ and $g_{i}^{\alpha^{c}} \in \mathcal{H}_{\alpha^{c}}$. A non-zero function $g$ is such that $\mathrm{rank}_{D}(g) = 1$. The above definition of $\alpha$-ranks yields the following properties.

**Proposition 2.1.** If $\alpha = \bigcup_{i \in I} \beta_{i}$ with $\{\beta_{i} : i \in I\}$ a collection of disjoint subsets of $D$, then for any function $g$,

$$ \mathrm{rank}_{\alpha}(g) \le \prod_{i \in I} \mathrm{rank}_{\beta_{i}}(g). $$
**Proposition 2.2.** For two functions g and h, and for any $\alpha \subset D$,

* $\mathrm{rank}_\alpha(g+h) \le \mathrm{rank}_\alpha(g) + \mathrm{rank}_\alpha(h)$,

* $\mathrm{rank}_\alpha(gh) \le \mathrm{rank}_\alpha(g)\mathrm{rank}_\alpha(h)$.

**Example 2.3.**

* $g(x) = g^1(x_1) \cdots g^d(x_d)$ can be written $g(x) = g^\alpha(x_\alpha)g^{\alpha^c}(x_{\alpha^c})$, where $g^\alpha(x_\alpha) = \prod_{\nu \in \alpha} g^\nu(x_\nu)$. Therefore $\mathrm{rank}_\alpha(g) \le 1$ for all $\alpha \subset D$.

* $g(x) = \sum_{k=1}^r g_k^1(x_1) \cdots g_k^d(x_d)$ can be written $\sum_{k=1}^r g_k^\alpha(x_\alpha)g_k^{\alpha^c}(x_{\alpha^c})$ with $g_k^\alpha(x_\alpha) = \prod_{\nu \in \alpha} g_k^\nu(x_\nu)$. Therefore, $\mathrm{rank}_\alpha(g) \le r$ for all $\alpha \subset D$.

* $g(x) = g^1(x_1) + \cdots + g^d(x_d)$ can be written $g(x) = g^\alpha(x_\alpha) + g^{\alpha^c}(x_{\alpha^c})$, with $g^\alpha(x_\alpha) = \sum_{\nu \in \alpha} g^\nu(x_\nu)$. Therefore, $\mathrm{rank}_\alpha(g) \le 2$ for all $\alpha \subset D$.

* $g(x) = \prod_{\alpha \in A} g^\alpha(x_\alpha)$ with A a collection of disjoint subsets is such that $\mathrm{rank}_\alpha(g) = 1$ for all $\alpha \in A$, and $\mathrm{rank}_\gamma(g) \le \prod_{\alpha \in A, \alpha \cap \gamma \ne \emptyset} \mathrm{rank}_\alpha(g^\alpha)$ for all $\gamma$.

## 2.2 Tree-based tensor formats

For a collection $T$ of non-empty subsets of $D$, we define the $T$-rank of $g$ as the tuple $\mathrm{rank}_T(g) = \{\mathrm{rank}_\alpha(g) : \alpha \in T\}$. Then, we define the set of tensors $\mathcal{T}_r^T(\mathcal{H})$ with $T$-rank bounded by $r = (r_\alpha)_{\alpha \in T}$ by

$$ \mathcal{T}_r^T(\mathcal{H}) = \{g \in \mathcal{H} : \mathrm{rank}_T(g) \le r\}. $$

A dimension partition tree $T$ is a tree such that (i) all nodes $\alpha \in T$ are non-empty subsets of $D$, (ii) $D$ is the root of $T$, (iii) every node $\alpha \in T$ with $\#\alpha \ge 2$ has at least two children and the set of children of $\alpha$, denoted by $S(\alpha)$, is a non-trivial partition of $\alpha$, and (iv) every node $\alpha$ with $\#\alpha = 1$ has no child and is called a leaf (see the example in Figure 1). The level of a node $\alpha$ is denoted by level($\alpha$). The levels are defined such that level($D$) = 0 and level($\beta$) = level($\alpha$) + 1 for $\beta \in S(\alpha)$. We let depth($T$) = $\max_{\alpha \in T}$ level($\alpha$) be the depth of $T$, and $\mathcal{L}(T)$ be the set of leaves of $T$, which are such that $S(\alpha) = \emptyset$ for all $\alpha \in \mathcal{L}(T)$.

When $T$ is a dimension partition tree, $\mathcal{T}_r^T(\mathcal{H})$ is the set of tensors with tree-based rank bounded by $r$, and an approximation in $\mathcal{T}_r^T(\mathcal{H})$ is called an approximation in tree-based (or hierarchical) tensor format [2, 1]. The trivial tree of depth 1 corresponds to the Tucker format. A balanced tree of depth $O(\log_2(d))$ corresponds to the classical hierarchical tensor format. A linear tree with depth $d-1$ corresponds to the tensor train format.

A tree-based rank $r$ is said admissible if $\mathcal{T}_r^T(\mathcal{H}) \ne \emptyset$. Necessary conditions of admissibility can be found in [12, Section 11.2.3]. In particular, $r_D$ has to be less than or equal to 1 for $\mathcal{T}_r^T(\mathcal{H})$ to be non empty, and $\mathcal{T}_r^T(\mathcal{H})$ is reduced to $\{0\}$ if $r_D = 0$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Example of a dimension partition tree over $D = \{1, \dots, 10\}$.</figcaption></figure>

## 2.3 Representation of tensors in tree-based format

We consider that the $\mathcal{H}_\nu$ are finite dimensional subspace of $L_{\mu\nu}^2(\mathcal{X}_\nu)$, $1 \le \nu \le d$ (e.g. a space of polynomials with bounded degree, a space of wavelets with finite resolution). Then, $\mathcal{H}$ is a finite-dimensional subspace of $L_{\mu}^2(\mathcal{X})$. Let $\{\phi_i^\nu : i \in I^\nu\}$ be a basis of $\mathcal{H}_\nu$, and $N_\nu = \dim(\mathcal{H}_\nu) = \#I^\nu$. For a multi-index $i = (i_1, \dots, i_d) \in I^1 \times \dots \times I^d := I$, we let $\phi_i(x) = \phi_{i_1}^1(x_1) \cdots \phi_{i_d}^d(x_d)$. Then the set of functions $\{\phi_i : i \in I\}$ is a basis of $\mathcal{H}$ and any function $g \in \mathcal{H}$ can be written

$$g(x) = \sum_{i_1 \in I^1} \cdots \sum_{i_d \in I^d} g_{i_1, \dots, i_d} \phi_{i_1}^1(x_1) \cdots \phi_{i_d}^d(x_d) = \sum_{i \in I} g_i \phi_i(x)$$

and identified with the set of coefficients $g \in \mathbb{R}^{I^1 \times \dots \times I^d} = \mathbb{R}^I$. This defines a linear bijection $G : \mathbb{R}^I \to \mathcal{H}$ such that $g = G(g)$. Also, we use the notation

$$g(x) = \langle \Phi(x), g \rangle,$$

where $\langle \cdot, \cdot \rangle$ is the canonical inner product on $\ell^2(\mathbb{R}^I)$ and $\Phi(x) \in \mathbb{R}^I$ is defined by

$$\Phi(x) = \Phi^1(x_1) \otimes \cdots \otimes \Phi^d(x_d),$$

with $\Phi^\nu(x_\nu) = (\phi_{i_\nu}^\nu(x_\nu))_{i_\nu \in I^\nu} \in \mathbb{R}^{I^\nu}$. If $\{\phi_i^\nu : i \in I^\nu\}$ is an orthonormal basis of $\mathcal{H}_\nu$, then $\{\phi_i : i \in I\}$ is an orthonormal basis of $\mathcal{H}$, and the $L_\mu^2$-norm of $g$ coincides with the canonical (Frobenius) norm $\|g\| = \langle g, g \rangle^{1/2}$ of the tensor $g$, i.e. $\|g\|_{L_\mu^2}^2 = \|g\|^2$. In this case, the map $G$ is a linear isometry.

A function $g \in \mathcal{T}_r^T(\mathcal{H})$ is identified with a tensor $g = G^{-1}(g)$ with entries

$$g_{i_1, \dots, i_d} = \sum_{1 \le k_\alpha \le r_\alpha} \prod_{\alpha \in T \setminus \mathcal{L}(T)} C_{(k_\beta)\beta \in S(\alpha), k_\alpha}^{\alpha} \prod_{\alpha \in \mathcal{L}(T)} C_{i_\alpha, k_\alpha}^{\alpha}$$
where

$$C^{\alpha} \in \mathbb{R}^{K^{\alpha}}, \quad K^{\alpha} := I^{\alpha} \times \{1, \dots, r_{\alpha}\},$$

with $I^{\alpha} = \times_{\beta \in S(\alpha)} \{1, \dots, r_{\beta}\}$ for $\alpha \notin \mathcal{L}(T)$.

The storage complexity (number of parameters) of a function $g$ in $\mathcal{T}_r^T(\mathcal{H})$ is

$$C(T, r) = \sum_{\alpha \in T \setminus \mathcal{L}(T)} r_{\alpha} \prod_{\beta \in S(\alpha)} r_{\beta} + \sum_{\alpha \in \mathcal{L}(T)} \dim(\mathcal{H}_{\alpha}) r_{\alpha}.$$

If $r_{\alpha} = O(R)$ and $\dim(\mathcal{H}_{\alpha}) = O(N)$, then since $\#T = O(d)$, $C(T, r) = O(dNR + (\#T - d - 1)R^{s+1} + R^s)$, where $s = \max_{\alpha \in T \setminus \mathcal{L}(T)} \#S(\alpha)$ is the arity of the tree. For a binary tree, $s = 2$ and $\#T = 2d - 1$, so that $C(T, r) = O(dNR + (d - 2)R^3 + R^2)$.

**Remark 2.4** (Discrete set $\mathcal{X}$). If $\mathcal{X} = \mathcal{X}_1 \times \dots \times \mathcal{X}_d$ is a finite or countable set and $\mu = \sum_{x \in \mathcal{X}} \delta_x = \mu_1 \otimes \dots \otimes \mu_d$ with $\mu_\nu = \sum_{x_\nu \in \mathcal{X}_\nu} \delta_{x_\nu}$, then the space $\mathcal{H}_\nu = L_{\mu_\nu}^2(\mathcal{X}_\nu)$ is identified with $\ell^2(\mathcal{X}_\nu)$. If $\mathcal{X}_\nu = \{x_\nu^{i_\nu} : i_\nu \in I^\nu\}$, the canonical basis $\phi_{i_\nu}^\nu(x_\nu^{j_\nu}) = \mathbf{1}_{i_\nu=j_\nu}$ (which is equal to 1 if $i_\nu = j_\nu$ and 0 otherwise) is orthonormal in $\mathcal{H}_\nu$. A function $g(x) = \sum_{i \in I} g_i \phi_i(x)$ is then isometrically identified with the set of coefficients $g_{i_1, \dots, i_d} = g(x^{i_1}, \dots, x^{i_d})$.

# 3 Representation of probabilistic models in tensor formats

In this section, we first discuss different types of representation of a probability distribution using tree-based tensor formats, and provide some results on the relations between the ranks of these representations. Then we provide several examples of standard probabilistic models, and discuss their representation in tree-based tensor format.

## 3.1 Representation of a probability distribution

The probability distribution of the random variable $X = (X_1, \dots, X_d)$ is characterized by its cumulative distribution function $F(x) = \mathbb{P}(X \le x)$. In the following we assume that the distribution admits a density $f(x)$ with respect to a product measure $\mu = \mu_1 \otimes \dots \otimes \mu_d$ on $\mathbb{R}^d$ (e.g. the Lebesgue measure), such that

$$F(x) = \int_{\{t \le x\}} f(t) d\mu(t).$$

This includes the case of a discrete random variable taking values in a finite or countable set $\mathcal{X} = \mathcal{X}_1 \times \dots \times \mathcal{X}_d$, with measure $\rho = \sum_{x \in \mathcal{X}} \mathbb{P}(X = x) \delta_x$, by letting $f(x) := \mathbb{P}(X = x)$ and $\mu := \sum_{x \in \mathcal{X}} \delta_x$. In this case, $f$ is identified with an element of $\mathbb{R}^{\mathcal{X}} = \mathbb{R}^{\mathcal{X}_1 \times \dots \times \mathcal{X}_d}$.

**Proposition 3.1.** Assume that the distribution $F$ admits a density $f$ with respect to a product measure $\mu$. Then for any $\alpha \subset D$,

$$\mathrm{rank}_{\alpha}(F) \le \mathrm{rank}_{\alpha}(f).$$
Moreover, if $\mu$ is the Lebesgue measure,

$$ \mathrm{rank}_{\alpha}(F) = \mathrm{rank}_{\alpha}(f). $$

**Proof.** If $f(x) = \sum_{k=1}^{r} f_k^{\alpha}(x_{\alpha}) f_k^{\alpha^c}(x_{\alpha^c})$, then $F(x) = \sum_{k=1}^{r} F_k^{\alpha}(x_{\alpha}) F_k^{\alpha^c}(x_{\alpha^c})$ with $F_k^{\beta}(x_{\beta}) = \int_{\{t_{\beta} \le x_{\beta}\}} f_k^{\beta}(t_{\beta}) d\mu_{\beta}(t_{\beta})$ for $\beta = \alpha$ and $\alpha^c$. This implies $\mathrm{rank}_{\alpha}(F) \le \mathrm{rank}_{\alpha}(f)$. If $\mu$ is the Lebesgue measure and $F(x) = \sum_{k=1}^{r} F_k^{\alpha}(x_{\alpha}) F_k^{\alpha^c}(x_{\alpha^c})$, then almost everywhere, $f(x) = \sum_{k=1}^{r} f_k^{\alpha}(x_{\alpha}) f_k^{\alpha^c}(x_{\alpha^c})$ with $f_k^{\beta}(x_{\beta}) = \partial_{x_{\nu_1}} \cdots \partial_{x_{\nu_{\#\beta}}} F_k^{\beta}(x_{\beta})$, $\beta = \{\nu_1, \dots, \nu_{\#\beta}\}$. This implies $\mathrm{rank}_{\alpha}(F) \ge \mathrm{rank}_{\alpha}(f)$. $\square$

**Remark 3.2.** Note that the above framework and results can be extended to the case where a random variable $X_{\nu}$ is either continuous or discrete, by letting $\mu_{\nu}$ be either the Lebesgue measure or a discrete measure.

For $1 \le \nu \le d$, let us denote by $F_{\nu} : \mathcal{X}_{\nu} \to [0, 1]$ the marginal cumulative distribution function of $X_{\nu}$. By Sklar's theorem, there exists a copula $C : [0, 1]^d \to [0, 1]$ such that

$$ F(x) = C(F_1(x_1), \dots, F_d(x_d)). $$

**Proposition 3.3.** For all $\alpha \subset D$, if $C$ is a copula of $X$,

$$ \mathrm{rank}_{\alpha}(F) \le \mathrm{rank}_{\alpha}(C). $$

If $F$ admits a density $f$ with respect to the Lebesgue measure, then $X$ admits a unique copula $C$ with density $c$ and

$$ \mathrm{rank}_{\alpha}(F) = \mathrm{rank}_{\alpha}(C) = \mathrm{rank}_{\alpha}(f) = \mathrm{rank}_{\alpha}(c). $$

**Proof.** If $C(u) = \sum_{k=1}^{r} C_k^{\alpha}(u_{\alpha}) C_k^{\alpha^c}(u_{\alpha^c})$, then $F(x) = \sum_{k=1}^{r} C_k^{\alpha}(u_{\alpha}) C_k^{\alpha^c}(u_{\alpha^c})$ with $u_{\nu} = F_{\nu}(x_{\nu})$. This implies $\mathrm{rank}_{\alpha}(F) \le \mathrm{rank}_{\alpha}(C)$. If $F$ admits a density with respect to the Lebesgue measure, then $C(u) = F(F_1^{-1}(u_1), \dots, F_d^{-1}(u_d))$, and if $F(x) = \sum_{k=1}^{r} F_k^{\alpha}(x_{\alpha}) F_k^{\alpha^c}(x_{\alpha^c})$, then $C(u) = \sum_{k=1}^{r} F_k^{\alpha}(x_{\alpha}) F_k^{\alpha^c}(x_{\alpha^c})$ with $x_{\nu} = F_{\nu}^{-1}(u_{\nu})$. This implies $\mathrm{rank}_{\alpha}(C) \le \mathrm{rank}_{\alpha}(F)$, and therefore $\mathrm{rank}_{\alpha}(F) = \mathrm{rank}_{\alpha}(C)$. The other equalities are deduced from proposition 3.1. $\square$

## 3.2 Encoding particular distributions in tensor formats

### 3.2.1 Mixtures

Consider a random variable $X = (X_1, \dots, X_d)$ which is a mixture of $m$ random variables $Z^i = (Z_1^i, \dots, Z_d^i)$ with weights $\gamma_i$, $1 \le i \le m$, such that $\sum_{i=1}^m \gamma_i = 1$. Let $f$ and $f^i$ denote the densities with respect to a product measure $\mu$ of the probability distributions of $X$ and $Z^i$ respectively. We have

$$ f(x) = \sum_{i=1}^{m} \gamma_i f^i(x). $$
From Proposition 2.2, we know that for any $\alpha \subset D$, $\text{rank}_\alpha(f) \le \sum_{i=1}^m \text{rank}_\alpha(f^i)$, therefore, for any tree $T$,

$$ \text{rank}_T(f) \le \sum_{i=1}^m \text{rank}_T(f^i). $$

Assuming that $Z^i$ has independent components $Z_k^i$ with densities $f_k^i$, we have $f^i(x) = f_1^i(x_1) \cdots f_d^i(x_d)$ with $\text{rank}_\alpha(f^i) = 1$ for any $\alpha$, and therefore, for any tree $T$, $\text{rank}_T(f) \le m$. Assume now that the function $f^i$ is represented in a tree based format with tree $T^i$. For any $\alpha \subset D$, there exists a subset $T_\alpha^i$ of $T^i$ which forms a partition of $\alpha$, and from Proposition 2.1, we have $\text{rank}_\alpha(f^i) \le \prod_{\beta \in T_\alpha^i} \text{rank}_\beta(f^i)$, and therefore

$$ \text{rank}_\alpha(f) \le \sum_{i=1}^m \prod_{\beta \in T_\alpha^i} \text{rank}_\beta(f^i) := R_\alpha. $$

Then a dimension tree for the representation of $f$ could be chosen to minimize the complexity $C(T, R)$ using the above upper bound $R = (R_\alpha)_{\alpha \in T}$ of the $T$-rank of $f$.

### 3.2.2 Markov processes

Consider a discrete time Markov process $X = (X_1, \dots, X_d)$ whose density is given by

$$ f(x) = f_{d|d-1}(x_d|x_{d-1}) \cdots f_{2|1}(x_2|x_1)f_1(x_1), $$

where $f_1$ is the density of $X_1$ and $f_{i|i-1}(\cdot|x_{i-1})$ is the density of $X_i$ knowing $X_{i-1} = x_{i-1}$. Let $m_i$ be the rank of the bivariate function $(t, s) \mapsto f_{i|i-1}(t|s)$, $i = 2, \dots, d$.

Let

$$ T = \{\{1, \dots, d\}, \{1\}, \dots, \{d\}, \{1, 2\}, \dots, \{1, \dots, d-1\}\} $$

be the linear tree of Figure 2a. We note that $\text{rank}_{\{1\}}(f) = \text{rank}(f_{2|1}) = m_2$, $\text{rank}_{\{d\}}(f) = \text{rank}(f_{d|d-1}) = m_d$ and for $2 \le \nu \le d-1$, $\text{rank}_{\{\nu\}}(f) \le \text{rank}(f_{\nu|\nu-1})\text{rank}(f_{\nu+1|\nu}) = m_\nu m_{\nu+1}$. Also, for $1 < \nu < d$, we have that $\text{rank}_{\{1, \dots, \nu\}}(f) = \text{rank}(f_{\nu+1|\nu})$. Letting $m = \max_i m_i$, we deduce that $f$ has a representation in tree-based format with complexity in $O(m^4)$. Note that the choice of tree is here crucial. Indeed, a different ordering of variables may lead to ranks growing exponentially with the dimension $d$. For instance, consider the tree $\tilde{T}$ represented in Figure 2b, with $\tilde{T} = \{\sigma(\alpha) : \alpha \in T\}$ with the permutation

$$ \sigma = (1, 3, \dots, 2 \lfloor \frac{d+1}{2} \rfloor - 1, 2, 4, \dots, 2 \lfloor \frac{d}{2} \rfloor). $$

For $\alpha = \{1, 3, \dots, 2k+1\}$, with $k \le \lfloor \frac{d+1}{2} \rfloor - 1$, we have $\text{rank}_\alpha(f) \le m_2 m_3 \cdots m_{2k+2} \le m^{2k+1}$ if $2k+1 < d$, and $\text{rank}_\alpha(f) \le m_2 m_3 \cdots m_{2k+1} \le m^{2k}$ if $2k+1 = d$. Therefore, the representation in the corresponding tree-based format has a complexity in $O(m^{2d-2})$. Example 3.4 presents a Markov process for which the tree-based rank exhibits such a behavior.
Therefore, when the structure of the Markov process is not known, a procedure for finding a suitable tree should be used (see Section 4.4 for the description of a tree adaptation algorithm).

<figure><img src="image_3.png" /><figcaption>Figure 2: Examples of linear dimension trees.</figcaption></figure>

**Example 3.4 (Discrete state space Markov process).** We consider the discrete time discrete state space Markov process $X = (X_1, \dots, X_8)$, where each random variable $X_\nu$ takes values in $\mathcal{X}_\nu = \{1, \dots, 5\}$. The distribution of $X$ writes

$$f(i_1, \dots, i_8) := \mathbb{P}(X_1 = i_1, \dots, X_8 = i_8) = f_{8|7}(i_8|i_7) \cdots f_{2|1}(i_2|i_1)f_1(i_1)$$

with $f_1(i_1) = 1/5$ for all $i_1 \in \mathcal{X}_1$, and for $\nu = 1, \dots, d-1$, $f_{\nu+1|\nu}(i_{\nu+1}|i_\nu) = P_{i_\nu, i_{\nu+1}}^\nu$ the $(i_\nu, i_{\nu+1})$ component of a randomly chosen rank-2 transition matrix $P^\nu$. We then have $\text{rank}(f_{\nu+1|\nu}) = m = 2$ for $\nu = 1, \dots, d-1$.

We first compute a representation of $f$ in tree-based format with the tree $T$ depicted in Figure 2a (using Algorithm 3 in [22] at precision $10^{-13}$), the obtained $\alpha$-ranks are shown in Figure 3a. We then compute a representation of $f$ in tree-based format with the tree $\tilde{T}$ depicted in Figure 2b, to obtain the $\alpha$-ranks shown in Figure 3b. We see that $\max_{\alpha \in T} \text{rank}_\alpha(f) = 4 = m^2$ whereas $\max_{\alpha \in \tilde{T}} \text{rank}_\alpha(f) = 128 = 2^7 = m^{d-1}$. As a consequence, the storage complexity of the representation is equal to 240 with $T$, and to 35088 with $\tilde{T}$, more than 146 times larger.

### 3.2.3 Graphical models

Let us consider a graphical model with a density of the form

$$f(x) = \prod_{\beta \in C} g_{\beta}(x_{\beta})$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Obtained α-ranks when representing the Markov process of Example 3.4 in tree-based format with two different linear dimension trees.</figcaption></figure>

where $\mathcal{C} \subset 2^D$ represents the cliques of a graph $G$ with nodes $\{1\}, \dots, \{d\}$.

Consider $\alpha \subset D$. First note that if $\alpha \in \mathcal{C}$, then $\mathrm{rank}_\alpha(g_\alpha) = 1$. Also, for a clique $\beta$ such that either $\beta \subset \alpha$ or $\beta \subset \alpha^c$, $\mathrm{rank}_\alpha(g_\beta) = 1$. Then let $\mathcal{C}_\alpha$ be the set of cliques that intersect both $\alpha$ and $\alpha^c$,

$$ \mathcal{C}_{\alpha} = \{\beta \in \mathcal{C} : \beta \cap \alpha \neq \emptyset, \beta \cap \alpha^c \neq \emptyset\}. $$

Since $\mathcal{C} \setminus \mathcal{C}_\alpha = \{\beta \in \mathcal{C} : \beta \subset \alpha^c \text{ or } \beta \subset \alpha\}$, and from Proposition 2.2, we have

$$ \mathrm{rank}_{\alpha}(f) = \mathrm{rank}_{\alpha}(\prod_{\beta \in \mathcal{C}_{\alpha}} g_{\beta}) \leq \prod_{\beta \in \mathcal{C}_{\alpha}} \mathrm{rank}_{\alpha}(g_{\beta}). $$

Assuming that the α-ranks of all functions $g_\beta$ are bounded by $m$, we have

$$ \mathrm{rank}_{\alpha}(f) \leq m^{\#\mathcal{C}_{\alpha}} := R_{\alpha}. $$

For the representation of $f$ in tree-based tensor format, a tree $T$ could be chosen such that it minimizes the complexity $C(T, R)$, with $R = (R_\alpha)_{\alpha \in T}$ the above upper bound of the $T$-rank of $f$.

**Example 3.5.** An example of graphical model is provided in Figure 4. The dimension is $d = 10$. Here we consider discrete random variables $X_\nu$ taking $N = 5$ possible instances, so that $f(x) = \mathbb{P}(X = x)$ is identified with a tensor of size $N^d = 5^{10} = 9765625$. It has the form

$$ f(x) = g_{1,2,3,7}(x_1, x_2, x_3, x_7)g_{3,4,5,6}(x_3, x_4, x_5, x_6)g_{4,8}(x_4, x_8)g_{8,9,10}(x_8, x_9, x_{10}). $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Example of graphical model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Representation in tree-based format with an initial random tree. The storage complexity is 10595875.</figcaption></figure>

We first consider the random binary tree of Figure 5 and compute a representation of the graphical model (using Algorithm 3 in [22] at precision $10^{-13}$) in the corresponding tree-based format. We observe a storage complexity of 10595875, higher than the storage complexity of the full tensor. After the tree optimization (with Algorithm 8 provided in [22]), we obtain the tree in Figure 6 with a storage complexity of 3275.

# 4 Learning with tree-based tensor format

In this section, we describe a learning algorithm for estimating a probability distribution of a random vector $X$ from independent samples of this distribution. We assume that the distribution of $X$ has a density $f$ with respect to a measure $\mu$ over $\mathcal{X}$ (possibly discrete).

We introduce a contrast (or loss) function $\gamma: L_{\mu}^{0}(\mathcal{X}) \times \mathcal{X} \to \mathbb{R}$ and the associated risk functional $\mathcal{R}: L_{\mu}^{0}(\mathcal{X}) \to \mathbb{R}$ defined by

$$ \mathcal{R}(g) = \mathbb{E}(\gamma(g, X)), \qquad (4) $$

such that the minimizer of $\mathcal{R}$ over the set of $\mu$-measurable functions is the density $f$.
<figure><img src="image_7.png" /><figcaption>Figure 6: Representation in tree-based format after tree optimization. The storage complexity is 3275.</figcaption></figure>

Given independent samples $\{x_i\}_{i=1}^n := S$ of $X$ (training sample), an approximation $g_n^M$ of the density is then obtained by minimizing the empirical risk

$$ \mathcal{R}_n(g) = \frac{1}{n} \sum_{i=1}^{n} \gamma(g, x_i) $$

over a certain model class $M$, here a class of functions in tree-based tensor format. For a theoretical analysis of the performance of these model classes in this setting, the reader is referred to [23].

Choosing the contrast function as

$$ \gamma(g, x) = \|g\|_{L_{\mu}^{2}}^{2} - 2g(x), \qquad (5) $$

with $\| \cdot \|_{L_{\mu}^2}$ the natural norm in $L_{\mu}^2$, leads to

$$ \mathcal{R}(g) = \mathcal{R}(f) + \|f - g\|_{L_{\mu}^{2}}^{2}, $$

so that the minimization of $\mathcal{R}(g)$ is equivalent to the minimization of the distance (in $L_{\mu}^2$ norm) between $g$ and the density $f$. For the minimization of the empirical risk over the class of functions in tree-based formats, with an adaptation of the ranks and possibly of the dimension tree, we will adapt the algorithms proposed in [22].

**Remark 4.1 (Maximum likelihood estimation).** Choosing the contrast function equal to $\gamma(g, x) = -\log(g(x))$ leads to

$$ \mathcal{R}(g) = - \int_{\mathcal{X}} \log(g(x)) f(x) d\mu(x) = \mathcal{R}(f) + D_{KL}(f \| g) $$
with $D_{KL}(f||g)$ the Kullback-Leibler divergence between $f$ and $g$ (with reference measure $\mu$). The empirical risk minimization then corresponds to the maximum likelihood estimation.

**Remark 4.2.** We recall that the present framework applies to the case of a discrete random variable $X$ taking values in a set $\mathcal{X} = \mathcal{X}_1 \times \cdots \times \mathcal{X}_d$, with measure $\rho = \sum_{x \in \mathcal{X}} \mathbb{P}(X = x)\delta_x$, by letting $f(x) := \mathbb{P}(X = x)$ and $\mu := \sum_{x \in \mathcal{X}} \delta_x$. A function $g \in L_{\mu}^2(\mathcal{X})$ is identified with an algebraic tensor $g \in \mathbb{R}^{\mathcal{X}}$ with a norm

$$ \|g\|_{L_{\mu}^2}^2 = \sum_{x \in \mathcal{X}} g(x)^2 = \|g\|_{\ell^2(\mathcal{X})}^2, $$

where $\| \cdot \|_{\ell^2(\mathcal{X})}$ is the canonical (Frobenius) norm on $\mathbb{R}^{\mathcal{X}}$.

## 4.1 Empirical risk minimization with a fixed tree-based format

We here consider the problem of density estimation in a tensor format $\mathcal{T}_r^T(\mathcal{H})$ with fixed dimension tree $T$ and tree-based rank $r$.

As seen in Section 2.3, a function $g \in \mathcal{T}_r^T(\mathcal{H})$ admits a parametrized representation in terms of tensors $C^\alpha \in \mathbb{R}^{K^\alpha}$, $K^\alpha = I^\alpha \times \{1, \dots, r_\alpha\}$, $\alpha \in T$. Therefore, for a given $x \in \mathcal{X}$, $g(x)$ admits the following parametrized representation

$$ g(x) = \Psi(x)((C^{\alpha})_{\alpha \in T}), $$

where $\Psi(x) : \times_{\alpha \in T} \mathbb{R}^{K^\alpha} \to \mathbb{R}$ is a multilinear map which associates to $(C^\alpha)_{\alpha \in T}$ the value $g(x) = \langle \Phi(x), g \rangle$, with a tensor $g$ given by (2.3). For a given $\alpha \in T$, $g(x)$ can be written $g(x) = \Psi^\alpha(x)(C^\alpha)$ with $\Psi^\alpha(x) : \mathbb{R}^{K^\alpha} \to \mathbb{R}$ a linear map depending on $(C^\beta)_{\beta \in T \setminus \{\alpha\}}$, which associates to a tensor $C^\alpha$ the value $\Psi(x)((C^\beta)_{\beta \in T})$. Therefore $g(x)$ can be written

$$ g(x) = \Psi^{\alpha}(x)(C^{\alpha}) = \langle \Psi^{\alpha}(x), C^{\alpha} \rangle_{\ell^2} = \sum_{i_{\alpha} \in K^{\alpha}} \Psi_{i_{\alpha}}^{\alpha}(x) C_{i_{\alpha}}^{\alpha}, \quad (6) $$

where $\langle \cdot, \cdot \rangle_{\ell^2}$ is the canonical inner product on $\mathbb{R}^{K^\alpha}$ and $\Psi^\alpha(x) = (\Psi_{i_\alpha}^\alpha(x))_{i_\alpha \in K^\alpha}$ is a tensor of the same size as $C^\alpha$. The explicit expression of $\Psi^\alpha(x)$ is given in [22, section 2.3]. The parametrization of a tensor is not unique and it is possible to obtain a representation (6) where the set of functions $(\Psi_{i_\alpha}^\alpha)_{i_\alpha \in K^\alpha}$ forms an orthonormal system in $L_{\mu}^2$, so that

$$ \|g\|_{L_{\mu}^2}^2 = \|C_{\alpha}\|_{\ell^2}^2 := \sum_{i_{\alpha} \in K^{\alpha}} (C_{i_{\alpha}}^{\alpha})^2. \quad (7) $$

This orthonormality property improves the numerical stability of learning algorithms. The algorithm for obtaining such a representation is given in [22, section 2.4].

The empirical risk minimization problem writes

$$ \min_{(C^\alpha)_{\alpha \in T}} \frac{1}{n} \sum_{i=1}^{n} \gamma(\Psi(\cdot)((C^\alpha)_{\alpha \in T}), x_i) \quad (8) $$
and can be solved with an alternating minimization algorithm which consists in minimization alternatively on each parameter $C^{\alpha}$, $\alpha \in T$. For a given $\alpha \in T$, for fixed parameters $C^{\beta}$, $\beta \in T \setminus \{\alpha\}$, we optimize over $C^{\alpha}$ by solving

$$ \min_{C^{\alpha} \in \mathbb{R}^{K^{\alpha}}} \frac{1}{n} \sum_{i=1}^{n} \gamma(\Psi^{\alpha}(\cdot)(C^{\alpha}), x_i). \quad (9) $$

Using the contrast function (5) and (7), this problem becomes

$$ \min_{C^{\alpha} \in \mathbb{R}^{K^{\alpha}}} \|C^{\alpha}\|_{\ell^2}^2 - \frac{2}{n} \sum_{i=1}^{n} \langle \Psi^{\alpha}(x_i), C^{\alpha} \rangle_{\ell^2} $$

and admits the solution

$$ C^{\alpha} = \frac{1}{n} \sum_{i=1}^{n} \Psi^{\alpha}(x_i), $$

where $\Psi^{\alpha}$ depends on the parameters $C^{\beta}$, $\beta \in T \setminus \{\alpha\}$.

## 4.2 Sparsity exploitation and model selection

Exploiting sparsity in the tensors $(C^{\alpha})_{\alpha \in T}$ of the representation of a function $g \in \mathcal{T}_r^T(\mathcal{H})$ may be relevant.

For a node $\alpha \in T$, we let $\mathcal{J}$ be a given subset of $K^{\alpha}$ (a sparsity pattern). Using the contrast function (5), the solution of the problem

$$ \min_{\substack{C^{\alpha} \in \mathbb{R}^{K^{\alpha}} \\ C_{i_{\alpha}}^{\alpha} = 0, \forall i_{\alpha} \notin \mathcal{J}}} \|C_{\alpha}\|_{\ell^2}^2 - \frac{2}{n} \sum_{i=1}^{n} \langle \Psi^{\alpha}(x_i), C^{\alpha} \rangle_{\ell^2} \quad (10) $$

is

$$ C_{i_{\alpha}}^{\alpha, \mathcal{J}} = \frac{1}{n} \sum_{k=1}^{n} \Psi_{i_{\alpha}}^{\alpha}(x_k) \mathbf{1}_{i_{\alpha} \in \mathcal{J}}, $$

with $\mathbf{1}_{i_{\alpha} \in \mathcal{J}}$ equal to 1 when $i_{\alpha} \in \mathcal{J}$, and 0 otherwise. This means that the solution of (10) is simply the solution of the unconstrained problem (9) where the components with indices $i_{\alpha} \notin \mathcal{J}$ are set to zero.

We introduce sparsity in the tensor $C^{\alpha}$ using two strategies (working set or thresholding), both yielding $L$ solutions $C^{\alpha,1}, \dots, C^{\alpha,L}$ associated with different sparsity patterns $\mathcal{J}_1, \dots, \mathcal{J}_L$.

* (Working-set strategy) We assume that there exists a natural increasing sequence $I_1^{\alpha} \subset \dots \subset I_L^{\alpha} = I^{\alpha}$, and we take $\mathcal{J}_l = \{(i_{\alpha}, k) : i_{\alpha} \in I_l^{\alpha}, 1 \le k \le r_{\alpha}\}$ for $1 \le l \le L$.
• (Thresholding strategy) Denoting by $a_1 > \dots > a_L$ the ordered sequence of the values taken by $|C_{i_\alpha}^\alpha|$, $i_\alpha \in K^\alpha$, the sparsity patterns are defined by $\mathcal{J}_l = \{i_\alpha \in K^\alpha : |C_{i_\alpha}^\alpha| \ge a_l\}$, for $1 \le l \le L$.

Different sparsity patterns $\mathcal{J}_l$ correspond to different models and associated solutions $C^{\alpha,l}$. We wish to select, among the different solutions $C^{\alpha,l}$ associated with the sparsity patterns $\mathcal{J}_l$, the one, denoted by $C^{\alpha,\star}$, that minimizes the risk $\mathcal{R}(\Psi^\alpha(\cdot)(C^{\alpha,l}))$. The risk can be estimated, using a sample $V$ independent of the training sample $S$, by

$$ \mathcal{R}_V(\Psi^\alpha(\cdot)(C^{\alpha,l})) = \frac{1}{\#V} \sum_{x \in V} \gamma(\Psi^\alpha(\cdot)(C^{\alpha,l}), x). $$

However, one usually does not wish to use a sample $V$ for the sole purpose of model selection, that could have been otherwise used for the learning process. Cross-validation offers a way of estimating the risk without the need of an independent sample. A leave-one-out estimator of the risk writes

$$ \mathcal{R}_n^{\text{loo}}(\Psi^\alpha(\cdot)(C^{\alpha,l})) = \frac{1}{n} \sum_{i=1}^{n} \mathcal{R}_{\{x_i\}}(\Psi^\alpha(\cdot)(C^{\alpha,l,-i})), $$

where $C^{\alpha,l,-i}$ is the solution of the empirical risk minimization problem using the training set $S \setminus \{x_i\}$, that is

$$ \min_{\substack{C^\alpha \in \mathbb{R}^{K^\alpha} \\ C_{i_\alpha}^\alpha = 0, \forall i_\alpha \notin \mathcal{J}}} \|C_\alpha\|_{\ell^2}^2 - \frac{2}{n-1} \sum_{\substack{j=1 \\ j \neq i}}^{n} \langle \Psi^\alpha(x_j), C^\alpha \rangle_{\ell^2}, $$

It can be written

$$ C_{i_\alpha}^{\alpha,l,-i} = \frac{1}{n-1} \sum_{\substack{k=1 \\ k \neq i}}^{n} \Psi_{i_\alpha}^{\alpha}(x_k) \mathbf{1}_{i_\alpha \in \mathcal{J}_l}. $$

The following proposition provides an explicit expression of the leave-one-out risk estimate as a function of $C^{\alpha,l}$. It is a special case for $p=1$ of the result of [26, Prop. 2.1] for the leave-$p$-out estimator.

**Proposition 4.3.** *The leave-one-out estimator of the risk can be expressed as*

$$ \begin{aligned} \mathcal{R}_n^{\text{loo}}(\Psi^\alpha(\cdot)(C^{\alpha,l})) &= \frac{-n^2}{(n-1)^2} \|C^{\alpha,l}\|_{\ell^2}^2 + \frac{2n-1}{n(n-1)^2} \sum_{i=1}^{n} \sum_{i_\alpha \in \mathcal{J}_l} \Psi_{i_\alpha}^\alpha(x_i)^2 \\ &= \frac{1}{n(n-1)} \sum_{i_\alpha \in \mathcal{J}_l} \left[ \sum_{i=1}^{n} \Psi_{i_\alpha}^\alpha(x_i)^2 - \frac{n}{n-1} \sum_{i \neq j} \Psi_{i_\alpha}^\alpha(x_i) \Psi_{i_\alpha}^\alpha(x_j) \right]. \end{aligned} $$
## 4.3 Tree-based rank adaptation

In Section 4.1, the computed approximation $g$ is in $\mathcal{T}_r^T(\mathcal{H})$, with both the tree-based rank $r$ and dimension tree $T$ given. However, when the dimension $d$ is large or the size $n$ of the training sample is small, the choice of the tree-based rank can be crucial. We use the rank adaptation strategy of [22, Algorithm 5] that incrementally increases ranks associated with a subset of nodes of the dimension tree, selected with a criterion based on truncation errors.

At step $m$, given an approximation $g^m$ in $\mathcal{T}_{r^m}^T(\mathcal{H})$ with rank $r^m$, the algorithm selects the ranks to increase by estimating the truncation errors

$$ \eta_{\alpha}(f, r_{\alpha}^{m}) = \min_{\text{rank}_{\alpha}(g) \le r_{\alpha}^{m}} \mathcal{R}(g) - \mathcal{R}(f), \quad (11) $$

and by increasing the ranks associated with the nodes $\alpha$ with high truncation errors. When using a $L^2$ risk,

$$ \eta_{\alpha}(f, r_{\alpha}^{m}) = \min_{\text{rank}_{\alpha}(g) \le r_{\alpha}^{m}} \|f - g\|_{L_{\mu}^{2}}^{2} = \sum_{k > r_{\alpha}^{m}} (\sigma_{k}^{\alpha}(f))^{2}, $$

where $\sigma_k^\alpha(f)$ denotes the $k$-th $\alpha$-singular value of $f$ (see [12, Section 8.3] for more details about the higher-order singular value decomposition of a tensor). The ranks associated with the nodes $\alpha$ yielding the highest truncation errors are then increased by one. In practice, we estimate the truncation errors $\eta_\alpha(f, r_\alpha^m)$ by computing the $\alpha$-singular values of an approximation $\tilde{g}^m$ of $f$ with $\alpha$-ranks $\tilde{r}_\alpha^m \ge r_\alpha^m$, which are such that

$$ \eta(\tilde{g}^m, r_{\alpha}^m)^2 = \sum_{k=r_{\alpha}^m+1}^{\tilde{r}_{\alpha}^m} (\sigma_k^\alpha(\tilde{g}^m))^2. $$

Then, choosing a parameter $\theta \in [0, 1]$, we increase by one the $\alpha$-ranks associated with the nodes $\alpha$ in

$$ T_{\theta} = \{\alpha \in \hat{T} : \eta_{\alpha}(\tilde{g}^m, r_{\alpha}^m) \ge \theta \max_{\beta \in T} \eta_{\beta}(\tilde{g}^m, r_{\beta}^m)\} $$

where $\hat{T} \subset T$ is the subset of nodes for which the $\alpha$-rank can be increased by one without violating the admissibility conditions. The approximation $\tilde{g}^m$ is computed using as initial values the parameters of the tensor $g^m + c$, with $c = \Psi(x)((\tilde{C}^\alpha)_{\alpha \in T}) \in \mathcal{T}_{\hat{T}}^T(\mathcal{H})$ a correction of $g^m$ obtained by solving the problem

$$ \min_{c \in \mathcal{T}_{\hat{T}}^T} \mathcal{R}_n(g^m + c) $$

which is equivalent to

$$ \min_{c \in \mathcal{T}_{\hat{T}}^T} \|c\|_{L_{\mu}^2}^2 - \frac{2}{n} \sum_{i=1}^{n} c(x_i) + 2 \int_{\mathcal{X}} c(x) g^m(x) d\mu(x). $$
As in the previous section, this problem is solved using an alternating minimization algorithm which consists in minimizing alternatively on each parameter $\tilde{C}^\alpha$, $\alpha \in T$. For a given $\alpha$, using a representation of $c(x) = \langle \Psi^\alpha(x), \tilde{C}^\alpha \rangle_{\ell^2}$ with orthonormal functions $\{\Psi_i^\alpha(x)\}_{i \in K^\alpha}$ depending on the fixed parameters $\tilde{C}^\beta$, $\beta \in T \setminus \{\alpha\}$, the minimization problem is

$$ \min_{\tilde{C}^\alpha \in \mathbb{R}^{K^\alpha}} \|\tilde{C}^\alpha\|_{\ell^2}^2 - \frac{2}{n} \sum_{i=1}^n \langle \Psi^\alpha(x_i), \tilde{C}^\alpha \rangle_{\ell^2} + 2 \int_X \langle \Psi^\alpha(x), \tilde{C}^\alpha \rangle_{\ell^2} g^m(x) d\mu(x). \quad (12) $$

The solution of (12) is

$$ \tilde{C}^\alpha = \frac{1}{n} \sum_{i=1}^{n} \Psi^\alpha(x_i) - S^\alpha, \quad (13) $$

with $S^\alpha$ such that

$$ \langle S^\alpha, \tilde{C}^\alpha \rangle_{\ell^2} = \int_X \langle \Psi^\alpha(x), \tilde{C}^\alpha \rangle_{\ell^2} g^m(x) d\mu(x), $$

that is

$$ S^\alpha = \sum_{\substack{i_\beta \in I^\beta \\ \beta \in \mathcal{L}(T) \\ \beta \notin S(\alpha) \cup \{\alpha\}}} \sum_{1 \le k_\beta \le \tilde{r}_\beta^m} \sum_{1 \le l_\beta \le r_\beta^m} \prod_{\beta \notin \mathcal{L}(T) \cup \{\alpha\}} \tilde{C}_{(k_\nu)_{\nu \in S(\beta)}, k_\beta}^{\beta} \prod_{\beta \notin \mathcal{L}(T)} C_{(l_\nu)_{\nu \in S(\beta)}, l_\beta}^{\beta} \prod_{\beta \in \mathcal{L}(T)} \tilde{C}_{i_\beta, k_\beta}^{\beta} C_{i_\beta, l_\beta}^{\beta} $$

if $\alpha \in T \setminus \mathcal{L}(T)$, and

$$ S^\alpha = \sum_{\substack{i_\beta \in I^\beta \\ \beta \in \mathcal{L}(T) \setminus \{\alpha\}}} \sum_{1 \le k_\beta \le \tilde{r}_\beta^m} \sum_{\beta \in T \setminus \{\alpha\}} \prod_{\beta \in \mathcal{L}(T)} \tilde{C}_{(k_\nu)_{\nu \in S(\beta)}, k_\beta}^{\beta} C_{(l_\nu)_{\nu \in S(\beta)}, l_\beta}^{\beta} \prod_{\beta \in \mathcal{L}(T) \setminus \{\alpha\}} \tilde{C}_{i_\beta, k_\beta}^{\beta} C_{i_\beta, l_\beta}^{\beta} C_{i_\alpha, l_\alpha}^{\alpha} $$

if $\alpha \in \mathcal{L}(T)$. In the above expressions, a summation over $\beta \notin J$ means over $\beta \in T \setminus J$.

In the numerical experiments, we will only consider a rank-one correction $c$, which means $\tilde{r}_\alpha = 1$ for all $\alpha$ and $\tilde{C}^\beta = 1$ for all $\beta \notin \mathcal{L}(T)$. Then the truncation errors are estimated with $\eta(\tilde{g}^m, r_\alpha^m)^2 = (\sigma_{r_\alpha^m+1}^\alpha(\tilde{g}^m))^2$.

## 4.4 Dimension tree adaptation

As seen in Section 3, the choice of dimension tree $T$ can have a strong influence on the $T$-rank $r$ of the representation of some probabilistic models, hence on their complexity $C(T, r)$. Furthermore, the dimension tree yielding the smallest complexity can carry information about the dependence structure of the represented probabilistic model.

Several approaches have been proposed for adapting dimension trees (or ordering of variables for a given tree topology) using different types of information on the function (see, e.g. [27, 28, 29].

Here, in order to try to recover the optimal tree in terms of storage complexity, we rely on the stochastic algorithm 8 in [22], which generates a Markov process in the set of
all possible trees. Starting from a tree $T$ and a representation of a function $f \in \mathcal{T}_r^T(\mathcal{H})$, it explores the set of trees of given arity (the maximal number of children among all the nodes of the tree) and returns a function $g \in \mathcal{T}_{r^*}^{T^*}(\mathcal{H})$ with $C(T^*, r^*) \le C(T, r)$, such that $\|f - g\|_{L_\mu^2} \le \varepsilon \|f\|_{L_\mu^2}$ for a given tolerance $\varepsilon$. This algorithm performs a sequence of permutations of pairs of nodes using efficient linear algebra (contractions of tensors and truncated singular value decompositions of matricizations with an accuracy depending on the sought tolerance $\varepsilon$).

In order to reduce the storage complexity of the representation of a function, one ought to perform permutations that remove the nodes associated with high ranks, and adding nodes associated with small ranks. This, as well as the complexity of the computations needed to permute two given nodes, are the criteria used by the stochastic algorithm to choose which nodes to permute.

More precisely, at a given iteration of the algorithm, given a tree $T$ and a current approximation $f \in \mathcal{T}_r^T(\mathcal{H})$, we draw a new tree $T^*$ from a suitable distribution over the set of all possible trees, compute an approximation $g$ of $f$ in $\mathcal{T}_{r^*}^{T^*}(\mathcal{H})$ such that $\|f - g\|_{L_\mu^2} \le \varepsilon \|f\|_{L_\mu^2}$ (using truncated higher order singular value decomposition), and accept this tree if $C(T^*, r^*) \le C(T, r)$. If accepted, we set $T \leftarrow T^*$ and $f \leftarrow g$. Given $T_0 := T$, a new proposal $T^*$ is defined as $T^* = \sigma_m \circ \dots \circ \sigma_1(T)$, where $\sigma_i$ is a function that permutes two nodes of the tree $T_{i-1} := \sigma_{i-1} \circ \dots \circ \sigma_1(T)$, and $m$ is a random number of permutations drawn according the distribution $\mathbb{P}(m=k) \propto k^{-\gamma_1}, k \in \mathbb{N}^*$, with $\gamma_1 > 0$ that gives a higher probability to low numbers of permutations. The function $\sigma_i$ is defined by drawing two nodes $\alpha_i$ and $\beta_i$ in $T_{i-1}$. The first node $\alpha_i$ in $T_{i-1}$ is drawn according the distribution $\mathbb{P}(\alpha_i = \alpha) \propto \text{rank}(P(\alpha))^{\gamma_2}$, where $P(\alpha)$ is the parent node of $\alpha$ in $T_i$ and $\gamma_2 > 0$. This gives a higher probability to select a node $\alpha$ with high parent's rank and obtain a next tree $T_i$ without this parent node. Given $\alpha_i$, the second node $\beta_i$ is drawn in $T_i$ according a distribution over $T_i$ such that $\mathbb{P}(\beta_i = \beta) = 0$ if $\beta$ is either an ascendant or a descendant of $\alpha_i$, and $\mathbb{P}(\beta_i = \beta) \propto d(\alpha_i, \beta)^{-\gamma_3}$ otherwise, with $\gamma_3 > 0$. $d(\alpha_i, \beta)$ measures the computational complexity for performing the permutation of $\alpha_i$ and $\beta$ and compute an approximation of $f$ in the format associated with the next tree $T_i$ (see [22, section 4.2.1] for a precise definition of $d$). This choice gives a higher probability to modifications of the tree inducing a low computational complexity.

**Remark 4.4 (Permutation of nodes).** Let consider two nodes $\alpha$ and $\beta$ of a tree $T$ such that $\alpha \cap \beta = \emptyset$ (i.e., one node is not an ascendant of the other). A tree $\sigma(T)$ obtained by permuting the nodes $\alpha$ and $\beta$ of $T$ is such that the brothers $\alpha$ (resp. $\beta$) in $\sigma(T)$ are the brothers of $\beta$ (resp. $\alpha$) in $T$, and the descendants of $\alpha$ (resp. $\beta$) are the same in $T$ and $\sigma(T)$.

**Remark 4.5 (Choice of parameters).** The parameters $\gamma_1, \gamma_2$ and $\gamma_3$ of the distributions involved in the definition of the proposal have an influence on the performance of the stochastic algorithm. In our experiments, we choose $\gamma_1 = \gamma_2 = \gamma_3 = 2$. An analysis of the impact of this choice would merit further studies.
**Example 4.6.** The tree adaptation algorithm allows the transition from the linear dimension tree $\mathcal{I}a$ to the balanced tree $\mathcal{I}e$ (of same arity) by performing a sequence of permutations of pairs of nodes displayed in Figures 7b to 7e. The transition between the trees $\mathcal{I}a$ and $\mathcal{I}b$ removed the node $\{1,2,3\}$ and added the node $\{3,4\}$.

<figure><img src="image_8.png" /><figcaption>Figure 7: Example of a sequence of nodes permutations for a transition from a linear tree to a balanced tree in dimension 8.</figcaption></figure>
# 5 Numerical examples

In this section, we illustrate the performance of the proposed algorithm for learning probability distributions. For a validation purpose, we here consider synthetic data generated from given distributions for which standard sampling strategies are available. The results were obtained with the MATLAB toolbox Approximation Toolbox [30] and can be reproduced by running the scripts provided with it. The algorithms are also available in the python package tensap [31].

**Contrast function and reference measure.** In all examples, we consider the $L^2$ contrast function $\gamma(g, x) = \|g\|_{L^2_\mu}^2 - 2g(x)$. The reference measure $\mu$ is always a product measure and $\mathcal{X} = \mathcal{X}_1 \times \dots \times \mathcal{X}_d$. For examples involving discrete random variables, we consider $\mu = \sum_{x \in \mathcal{X}} \delta_x$, so that $L^2_\mu(\mathcal{X})$ is identified with $\ell^2(\mathcal{X})$. For continuous random variables, $\mu$ is taken as the Lebesgue measure or a uniform probability measure on $\mathcal{X}$ when $\mathcal{X}$ is a compact set.

**Approximation spaces.** In the case of discrete random variables with a finite set $\mathcal{X}$, we let $\mathcal{H} = L^2_\mu(\mathcal{X})$ so that there is no discretization error, and we use a canonical basis (see Remark 2.4). In the case of continuous random variables, for each dimension $\nu = 1, \dots, d$, we introduce a finite dimensional space $\mathcal{H}_\nu$ in $L^2_{\mu_\nu}(\mathcal{X}_\nu)$ and use orthonormal bases of $\mathcal{H}_\nu$ (e.g. polynomials, wavelets). We exploit sparsity in the leaf tensors $(C^\alpha)_{\alpha \in \mathcal{L}(T)}$ by using the working-set strategy presented in Section 4.1. For polynomial bases, we use the natural sequence of candidate patterns associated with spaces of polynomials with increasing degree. For wavelets bases, we use a sequence of candidate patterns associated with wavelets spaces with increasing resolution.

**Tensor formats.** We only consider tensor formats associated with binary dimension partition trees. When using the algorithm with tree-based rank and tree adaptation, the starting dimension tree is always a linear tree (such as in Figure 2a), where the dimensions $\nu = 1, \dots, d$ are randomly assigned to the leaf nodes.

**Error measures.** The quality of the obtained approximation $g$ is assessed by estimating the risk by

$$ \mathcal{R}_{S_{\text{test}}}(g) = \|g\|_{L^2_\mu}^2 - \frac{2}{\#S_{\text{test}}} \sum_{x \in S_{\text{test}}} g(x), $$

with $S_{\text{test}}$ a sample of $X$, independent of $S$, as well as, when $f$ can be evaluated, by computing the relative error

$$ \varepsilon(g) = \left( \frac{\sum_{x \in S_{\varepsilon}} (f(x) - g(x))^2}{\sum_{x \in S_{\varepsilon}} f(x)^2} \right)^{1/2}, $$
with $S_ε$ a sample from $\mu$ if $\mu$ is a probability measure, or from $\frac{1}{\mu(\mathcal{X})}\mu$ when $\mu$ is a finite measure (e.g. when $\mu$ is the Lebesgue measure over a compact set $\mathcal{X}$). In the case of discrete random variables, a function $f$ in $\mathbb{R}^\mathcal{X}$ is identified with a multi-dimensional array, and $S_\epsilon$ corresponds to a sample of the entries of the array.

## 5.1 Truncated multivariate normal distribution

We first consider the estimation of the density of a random vector $X = (X_1, \dots, X_6)$ following a truncated normal distribution with zero mean and covariance matrix $\Sigma$. Its support is $\mathcal{X} = \times_{\nu=1}^6 [-5\sigma_\nu, 5\sigma_\nu]$, with $\sigma_\nu^2 = \Sigma_{\nu\nu}$. The reference measure $\mu$ is the Lebesgue measure on $\mathcal{X}$ and the density to approximate is such that

$$f(x) \propto \exp \left( -\frac{1}{2} x^T \Sigma^{-1} x \right) \mathbf{1}_{x \in \mathcal{X}}. \qquad (14)$$

In each dimension $\nu$, we use polynomials of maximal degree 50, orthonormal in $L^2(-5\sigma_\nu, 5\sigma_\nu)$.

### 5.1.1 Groups of independent random variables

We consider the following covariance matrix

$$\Sigma = \begin{pmatrix} 2 & 0 & 0.5 & 1 & 0 & 0.5 \\ 0 & 1 & 0 & 0 & 0.5 & 0 \\ 0.5 & 0 & 2 & 0 & 0 & 1 \\ 1 & 0 & 0 & 3 & 0 & 0 \\ 0 & 0.5 & 0 & 0 & 1 & 0 \\ 0.5 & 0 & 1 & 0 & 0 & 2 \end{pmatrix}. \qquad (15)$$

Up to a permutation $(3, 6, 1, 4, 2, 5)$ of its rows and columns, it can be written

$$\begin{pmatrix} 2 & 1 & 0.5 & 0 & 0 & 0 \\ 1 & 2 & 0.5 & 0 & 0 & 0 \\ 0.5 & 0.5 & 2 & 1 & 0 & 0 \\ 0 & 0 & 1 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0.5 \\ 0 & 0 & 0 & 0 & 0.5 & 1 \end{pmatrix}$$

so that one can see that the random variables $(X_1, X_3, X_4, X_6)$ and $(X_2, X_5)$ are independent, as well as $X_4$ and $(X_3, X_6)$. Therefore, the density has the form

$$f(x) = f_{1,3,4,6}(x_1, x_3, x_4, x_6) f_{2,5}(x_2, x_5) = f_{4|1}(x_4|x_1) f_{1,3,6}(x_1, x_3, x_6) f_{2,5}(x_2, x_5).$$

One can then expect that, when approximating the density of $X$ in tree-based format, a suitable dimension tree $T$ would contain the nodes $\{2, 5\}$ and $\{1, 3, 4, 6\}$, since $\text{rank}_{\{2,5\}}(f) =$
<figure><img src="image_9.png" /><figcaption>Figure 8: Example of expected tree T for the approximation of density (14) with covariance matrix (15) in tree-based tensor format with tree adaptation.</figcaption></figure>

$\text{rank}_{\{1,3,4,6\}}(f) = 1$. If we further assume that $f_{1,3,6}$ has low ranks, it would contain the nodes $\{3,6\}$ and $\{1,4\}$, since $\text{rank}_{\{3,6\}}(f) = \text{rank}_{\{3,6\}}(f_{1,3,6})$ and $\text{rank}_{\{1,4\}}(f) = \text{rank}_{\{1,4\}}(f_{1,3,6})$ (see a possible tree in Figure 8).

Table 1 shows the results obtained with the learning algorithm with different sizes of training set. We notice that, as expected, $\mathcal{R}_{\text{Stest}}(g)$ and $\varepsilon(g)$ decrease with $n$. Figure 10 shows that, on average, the convergence of $\varepsilon(g)$ is approximately in $n^{-1/2}$. This convergence is better than the one expected (and observed in practice) when using a second-order kernel estimator, in $n^{-1/5}$. On Figure 9, we show the obtained tree (associated with the smallest error over 10 trials) for different training sample sizes $n$. For $n \ge 10^4$, the algorithm yields a tree with the expected nodes (see discussion above).

<table><caption>Table 1: Obtained results for the learning of density (14) with covariance matrix (15), with different training sample sizes $n$. Expectation and ranges (in brackets) over 10 trials of $\mathcal{R}_{\text{Stest}}(g)$, $\varepsilon(g)$ and $C(T, r)$, and tree T associated with the smallest error over 10 trials.</caption><thead><tr><th>n</th><th>R<sub>Stest</sub>(g) × 10<sup>-2</sup></th><th>ε(g)</th><th>T</th><th>C(T, r)</th></tr></thead><tbody><tr><td>10<sup>2</sup></td><td>22.76 [-5.50, 119]</td><td>1.71 [0.53, 4.06]</td><td>Fig. 9a</td><td>311.0 [311, 311]</td></tr><tr><td>10<sup>3</sup></td><td>-6.35 [-7.29, -5.93]</td><td>0.41 [0.22, 0.47]</td><td>Fig. 9b</td><td>443.9 [311, 637]</td></tr><tr><td>10<sup>4</sup></td><td>-7.23 [-7.60, -6.85]</td><td>0.23 [0.11, 0.33]</td><td>Fig. 9c</td><td>622.9 [521, 911]</td></tr><tr><td>10<sup>5</sup></td><td>-7.67 [-7.68, -7.66]</td><td>0.05 [0.04, 0.07]</td><td>Fig. 9c</td><td>1081.3 [911, 1213]</td></tr><tr><td>10<sup>6</sup></td><td>-7.69 [-7.70, -7.69]</td><td>0.01 [0.008, 0.012]</td><td>Fig. 9c</td><td>1357.7 [1283, 1546]</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Figure 9: Dimension trees $T$ obtained after computing an approximation in tree-based tensor format of density (14) with covariance matrix (15), using different training sample sizes $n$. The displayed trees are associated with the smallest error over 10 trials.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Convergence (in expectation) with the training sample size $n$ for the approximation of density (14) with covariance matrix (15). Expectations are computed over 10 trials for each $n$.</figcaption></figure>
### 5.1.2 Band-diagonal covariance matrix

We now consider the following covariance matrix

$$ \Sigma = \begin{pmatrix} 2 & 1/5 & 0 & 0 & 1/4 & 0 \\ 1/5 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 1/3 & 1/2 \\ 0 & 0 & 0 & 2 & 0 & 1 \\ 1/4 & 0 & 1/3 & 0 & 2 & 0 \\ 0 & 0 & 1/2 & 1 & 0 & 2 \end{pmatrix} \quad (16) $$

which is, after applying the permutation $\sigma = (4, 6, 3, 5, 1, 2)$, a band diagonal matrix

$$ \begin{pmatrix} 2 & 1 & 0 & 0 & 0 & 0 \\ 1 & 2 & 1/2 & 0 & 0 & 0 \\ 0 & 1/2 & 2 & 1/3 & 0 & 0 \\ 0 & 0 & 1/3 & 2 & 1/4 & 0 \\ 0 & 0 & 0 & 1/4 & 2 & 1/5 \\ 0 & 0 & 0 & 0 & 1/5 & 2 \end{pmatrix}. $$

The vector $(X_{\sigma(1)}, \dots, X_{\sigma(6)})$ therefore represents a Markov process and the density $f$ has the following form

$$ f(x) = f_{2|1}(x_2|x_1)f_{1|5}(x_1|x_5)f_{5|3}(x_5|x_3)f_{3|6}(x_3|x_6)f_{6|4}(x_6|x_4). $$

Given this structure, one might expect the density of $X$ to be efficiently and accurately represented in tree-based tensor format with one of the linear trees of Figure 11 or any tree containing the same internal nodes.

Table 2 shows the results obtained when using the learning algorithm to approximate the density of $X$. We can see that the risk and error decrease when the size of the training sample increases, with a convergence rate of about $n^{-2/5}$. Figure 12 shows the obtained trees (associated with the smallest error over 10 trials) for different sizes of training sets. We observe that except for $n = 10^2$, the algorithm yields trees that contain most of the expected internal nodes.
<figure><img src="image_12.png" /><figcaption>Figure 11: Example of expected trees *T* for the approximation of density (14) with covariance matrix (16) in tree-based tensor format.</figcaption></figure>

<table><caption>Table 2: Obtained results for the learning of density (14) with covariance matrix (16), with different training sample sizes *n*. Expectation and ranges (in brackets) over 10 trials of $\mathcal{R}_{\text{Stest}}(g)$, $\varepsilon(g)$ and $C(T, r)$, and tree *T* associated with the smallest error over 10 trials.</caption><thead><tr><th>n</th><th>R<sub>Stest</sub>(g) × 10<sup>-2</sup></th><th>ε(g)</th><th>T</th><th>C(T, r)</th></tr></thead><tbody><tr><td>10<sup>2</sup></td><td>-2.04 [-4.21, 1.91]</td><td>0.80 [0.58, 1.11]</td><td>Fig. 12a</td><td>311 [311, 311]</td></tr><tr><td>10<sup>3</sup></td><td>-5.64 [-5.88, -5.13]</td><td>0.33 [0.26, 0.45]</td><td>Fig. 12b</td><td>454.4 [311, 579]</td></tr><tr><td>10<sup>4</sup></td><td>-6.18 [-6.27, -5.58]</td><td>0.13 [0.08, 0.34]</td><td>Fig. 12c</td><td>670.0 [311, 830]</td></tr><tr><td>10<sup>5</sup></td><td>-6.30 [-6.30, -6.27]</td><td>0.04 [0.03, 0.07]</td><td>Fig. 12c</td><td>964.7 [724, 1012]</td></tr><tr><td>10<sup>6</sup></td><td>-6.31 [-6.31, -6.31]</td><td>0.014 [0.01, 0.02]</td><td>Fig. 12c</td><td>1140.2 [961, 1359]</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 12: Dimension trees $T$ obtained after computing an approximation in tree-based tensor format of density (14) with covariance matrix (16) using different training sample sizes $n$. The displayed trees are associated with the smallest error over 10 trials.</figcaption></figure>
## 5.2 Markov chain process

In this section, we study the discrete time discrete state space Markov process of Example 3.4. We recall that $X = (X_1, \dots, X_8)$, where each random variable $X_\nu$ takes values in $\mathcal{X}_\nu = \{1, \dots, 5\}$. The distribution of $X$ writes

$$f(i_1, \dots, i_8) = \mathbb{P}(X_1 = i_1, \dots, X_8 = i_8) = f_{d|d-1}(i_8|i_7) \cdots f_{2|1}(i_2|i_1)f_1(i_1) \quad (17)$$

with $f_1(i_1) = 1/5$ for all $i_1 \in \mathcal{X}_1$, and for $\nu = 1, \dots, d-1$, $f_{\nu+1|\nu}(i_{\nu+1}|i_\nu) = P_{i_\nu, i_{\nu+1}}$ the $(i_\nu, i_{\nu+1})$ component of a randomly chosen rank-2 transition matrix $P$, independent of the dimension $\nu$.

As shown in Example 3.4, the choice of the tree has a great impact on the storage complexity of the representation of the Markov process. We then expect the adaptive learning algorithm to compute an approximation of $f$ with a tree containing the same internal nodes as in Figure 2a. Table 3 shows the ranges over 10 trials of the obtained results. One can notice that, even though the algorithm did not recover an optimal tree for the representation of the Markov chain, it is able, with a sample size high enough, to represent it with dimension trees including most of the internal nodes yielding the smallest $\alpha$-ranks, limiting the complexity of the representation.

<table><caption>Table 3: Obtained results for the learning of (17), with different training sample sizes $n$. Expectation and ranges (in brackets) over 10 trials of $\mathcal{R}_{S_{\text{test}}}(g)$, $\varepsilon(g)$ and $C(T, r)$, and tree $T$ associated with the smallest error over 10 trials.</caption><thead><tr><th>n</th><th>R<sub>S<sub>test</sub></sub>(g)</th><th>ε(g)</th><th>T</th><th>C(T, r)</th></tr></thead><tbody><tr><td>10<sup>3</sup></td><td>-1.86 [-2.22, -1.25]</td><td>0.59 [0.49, 0.75]</td><td>Fig. 13a</td><td>84 [47, 109]</td></tr><tr><td>10<sup>4</sup></td><td>-2.58 [-2.85, -2.04]</td><td>0.31 [0.16, 0.55]</td><td>Fig. 13b</td><td>183 [72, 298]</td></tr><tr><td>10<sup>5</sup></td><td>-2.92 [-2.93, -2.91]</td><td>0.06 [0.04, 0.08]</td><td>Fig. 13c</td><td>422 [294, 519]</td></tr><tr><td>10<sup>6</sup></td><td>-2.93 [-2.93, -2.93]</td><td>0.016 [0.01, 0.02]</td><td>Fig. 13d</td><td>591.6 [384, 1010]</td></tr></tbody></table>

## 5.3 Graphical model with discrete random variables

We consider the graphical model already studied in Example 3.5, in dimension $d = 10$, represented in Figure 4. The random variable $X_\nu$ takes values in $\mathcal{X}_\nu = \{1, 2, 3, 4, 5\}$, $1 \le \nu \le d$, so that $f(i) = \mathbb{P}(X = i)$ (the probability table) is defined by

$$f(i_1, \dots, i_{10}) = f_{1,2,3,7}(i_1, i_2, i_3, i_7)f_{3,4,5,6}(i_3, i_4, i_5, i_6)f_{4,8}(i_4, i_8)f_{8,9,10}(i_8, i_9, i_{10}). \quad (18)$$

The tensors $f_\alpha$ are randomly selected under the constraint that any of their matricization has a rank equal to 3. In this example, the probability table can be stored in memory, so that a standard pseudo-random number generator can be used to generate samples from $X$.
<figure><img src="image_14.png" /><figcaption>Figure 13: Dimension trees $T$ obtained after computing an approximation in tree-based tensor format of (17), using different training sample sizes $n$. The displayed trees are associated with the smallest error over 10 trials.</figcaption></figure>
Representing this function in tree-based tensor format with the binary tree in Figure 5 yields a storage complexity of 117027, whereas using the tree in Figure 6, which exhibits the dependence structure of the graphical model, leads to a representation with a storage complexity of 675. We then expect our algorithms to be able to learn *f* with a tree representing its dependence structure.

Table 4 shows the ranges over 10 trials of the obtained results. We observe that, even though the obtained errors are high, the algorithm is able to provide approximations of *f* with a tree that exhibits the dependence structure of the graphical model (for instance by containing the nodes {1, 2, 3, 7}, {4, 8} or {4, 8, 9, 10}, which are cliques of the graph of *f*).

<table><caption>Table 4: Obtained results for the learning of (18), with different training sample sizes *n*. Expectation and ranges (in brackets) over 10 trials of $\mathcal{R}_{S_{\text{test}}}(g)$, $\varepsilon(g)$ and $C(T, r)$, and tree *T* associated with the smallest error over 10 trials.</caption><thead><tr><th>n</th><th>R<sub>S<sub>test</sub></sub>(g)</th><th>ε(g)</th><th>T</th><th>C(T, r)</th></tr></thead><tbody><tr><td>10<sup>3</sup></td><td>-1.24 [-1.38, -1.13]</td><td>0.56 [0.49, 0.61]</td><td>Fig. 14a</td><td>64.8 [59, 78]</td></tr><tr><td>10<sup>4</sup></td><td>-1.40 [-1.59, -1.23]</td><td>0.46 [0.35, 0.56]</td><td>Fig. 14b</td><td>109.3 [59, 208]</td></tr><tr><td>10<sup>5</sup></td><td>-1.71 [-1.77, -1.52]</td><td>0.21 [0.13, 0.40]</td><td>Fig. 14c</td><td>518.2 [131, 703]</td></tr><tr><td>10<sup>6</sup></td><td>-1.75 [-1.80, -1.60]</td><td>0.16 [0.06, 0.33]</td><td>Fig. 14d</td><td>784.3 [363, 1148]</td></tr></tbody></table>
<figure><img src="image_15.png" /></figure>

(a) $n = 10^3$.

<figure><img src="image_16.png" /></figure>

(c) $n = 10^5$.

<figure><img src="image_17.png" /></figure>

(b) $n = 10^4$.

<figure><img src="image_18.png" /></figure>

(d) $n = 10^6$.

<caption>Figure 14: Dimension trees $T$ obtained after computing an approximation in tree-based tensor format of (18), using different training sample sizes $n$. The displayed trees are associated with the smallest error over 10 trials.</caption>
# 6 Conclusion

We have proposed algorithms for learning high-dimensional probability distributions in tree-based tensor formats. The tree-based rank and dimension tree adaptation strategies enable the computation of approximations with trees that exhibit the dependence structure of the distribution.

Different contrast (or loss) functions could be used for learning a probabilistic distribution. However, this calls for a modification of the proposed rank-adaptation strategy, since the rank selection procedure as well as the computation of a correction depend on the choice of the contrast function.

Also, the proposed algorithms should be modified to be able to work efficiently with large data sets, for instance by using subsampling methods as in stochastic gradient methods.

# Acknowledgements

The authors acknowledge AIRBUS Group for the financial support with the project AtRandom, and the Joint Laboratory of Marine Technology between Naval Group and Centrale Nantes for the financial support with the project Eval PI.

# References

[1] Falcó, A., Hackbusch, W., and Nouy, A., Tree-based tensor formats, *SeMA Journal*, 2018.

[2] Hackbusch, W. and Kuhn, S., A New Scheme for the Tensor Representation, *Journal of Fourier analysis and applications*, 15(5):706-722, 2009.

[3] Oseledets, I. and Tyrtyshnikov, E., Breaking the curse of dimensionality, or how to use svd in many dimensions, *SIAM Journal on Scientific Computing*, 31(5):3744-3759, 2009.

[4] Oseledets, I., Tensor-train decomposition, *SIAM J. Sci. Comput.*, 33(5):2295-2317, 2011.

[5] Uschmajew, A. and Vandereycken, B., The geometry of algorithms using hierarchical tensors, *Linear Algebra and its Applications*, 439(1):133-166, 2013.

[6] Falcó, A., Nouy, A., and Hackbusch, W., Geometry of tree-based tensor formats in tensor banach spaces, *arXiv preprint arXiv:2011.08466*, 2020.

[7] Falcó, A., Hackbusch, W., and Nouy, A., On the Dirac-Frenkel variational principle on tensor Banach spaces, *Foundations of Computational Mathematics*, 19(1):159-204, Feb 2019.
[8] Holtz, S., Rohwedder, T., and Schneider, R., On manifolds of tensors of fixed tt-rank, *Numerische Mathematik*, 120(4):701–731, 2012.

[9] Ali, M. and Nouy, A., Approximation with tensor networks. part i: Approximation spaces, *ArXiv*, abs/2007.00118, 2020.

[10] Ali, M. and Nouy, A., Approximation with tensor networks. part ii: Approximation rates for smoothness classes, *ArXiv*, abs/2007.00128, 2020.

[11] Ali, M. and Nouy, A., Approximation with tensor networks. part iii: Multivariate approximation, *arXiv preprint arXiv:2101.11932*, 2021.

[12] Hackbusch, W., *Tensor spaces and numerical tensor calculus*, Vol. 42 of Springer series in computational mathematics, Springer, Heidelberg, 2012.

[13] Kolda, T.G. and Bader, B.W., Tensor decompositions and applications, *SIAM Review*, 51(3):455–500, 2009.

[14] Khoromskij, B., Tensors-structured numerical methods in scientific computing: Survey on recent advances, *Chemometrics and Intelligent Laboratory Systems*, 110(1):1–19, 2012.

[15] Grasedyck, L., Kressner, D., and Tobler, C., A literature survey of low-rank tensor approximation techniques, *GAMM-Mitteilungen*, 36(1):53–78, 2013.

[16] Nouy, A., *Low-Rank Tensor Methods for Model Order Reduction*, pp. 857–882. Springer International Publishing, Cham, 2017.

[17] Nouy, A., Low-rank methods for high-dimensional approximation and model order reduction, *In Model Reduction and Approximation: Theory and Algorithms*, SIAM, volume 15:chapter 4, pages 171–226, 2017.

[18] Bachmayr, M., Schneider, R., and Uschmajew, A., Tensor networks and hierarchical tensors for the solution of high-dimensional partial differential equations, *Foundations of Computational Mathematics*, pp. 1–50, 2016.

[19] Eigel, M., Marschall, M., and Schneider, R., Sampling-free bayesian inversion with adaptive hierarchical tensor representations, *Inverse Problems*, 34(3):035010, feb 2018.

[20] Dolgov, S., Anaya-Izquierdo, K., Fox, C., and Scheichl, R., Approximation and sampling of multivariate probability distributions in the tensor train decomposition, *Statistics and Computing*, 30(3):603–625, 2020.

[21] Eigel, M., Gruhlke, R., and Marschall, M., Low-rank tensor reconstruction of concentrated densities with application to bayesian inversion, 2020.
[22] Grelier, E., Nouy, A., and Chevreuil, M., Learning with tree-based tensor formats, *arXiv e-prints*, p. arXiv:1811.04455, Nov 2018.

[23] Michel, B. and Nouy, A., Learning with tree tensor networks: complexity estimates and model selection, *arXiv e-prints*, p. arXiv:2007.01165, July 2020.

[24] Dick, J., Kuo, F.Y., and Sloan, I.H., High-dimensional integration: the quasi-monte carlo way, *Acta Numerica*, 22:133–288, 2013.

[25] Haberstich, C., Nouy, A., and Perrin, G., Active learning of tree tensor networks using optimal least-squares, *arXiv preprint arXiv:2104.13436*, 2021.

[26] Celisse, A., Optimal cross-validation in density estimation with the $L^2$-loss, *The Annals of Statistics*, 42(5):1879–1910, 2014.

[27] Ballani, J. and Grasedyck, L., Tree adaptive approximation in the hierarchical tensor format, *SIAM journal on scientific computing*, 36(4):A1415–A1431, 2014.

[28] Bebendorf, M. and Kuske, C., Separation of variables for function generated high-order tensors, *Journal of Scientific Computing*, 61(1):145–165, 2014.

[29] Barcza, G., Legeza, O., Marti, K.H., and Reiher, M., Quantum-information analysis of electronic states of different molecular structures, *Phys. Rev. A*, 83:012508, Jan 2011.

[30] Nouy, A., Grelier, E., and Giraldi, L., Approximation Toolbox, February 2020. https://anthony-nouy.github.io/ApproximationToolbox/.

[31] Nouy, A. and Grelier, E., anthony-nouy/tensap v1.1, June 2020. https://anthony-nouy.github.io/tensap/.

[32] Falcó, A. and Hackbusch, W., On minimal subspaces in tensor representations, *Foundations of Computational Mathematics*, 12:765–803, 2012.