2.1 Notations

2.1.1 Tensor notations

In the following, we denote a vector by a small bold letter $\mathbf {a}$ , a matrix by a capital bold letter $\mathbf {A}$ , and a tensor by a calligraphic letter $\mathcal {A}$ . Scalars are denoted so they do not fall in any of those categories. Additionally, we denote $[D] = \{1, \dots , D\}$ . Elements of a vector $\mathbf {a}$ are denoted by $a_j$ , of a matrix $\mathbf {A}$ by $a_{jj'}$ , and of a tensor $\mathcal {A}$ of order D by $a_{j_1 \dots j_D}$ , often abbreviated using the vector of indices $\mathbf {j} = (j_1, \dots , j_D)$ as $a_{\mathbf {j}}$ . The jth column of $\mathbf {A}$ is denoted by $\mathbf {a}_j$ . The Frobenius norm of a matrix is defined as $\|\mathbf {A}\|_F^2 = \operatorname {\mathrm {tr}}(\mathbf {A}^\top \mathbf {A}) = \sum _{j, j'} a_{jj'}^2$ . We extend this definition to tensors as $\|\mathcal {A}\|_F^2 = \sum _{\mathbf {j}} a_{\mathbf {j}}^2$ . We denote the mode-d matricization of a tensor $\mathcal {X}$ as the matrix $\mathbf {X}_{(d)} \in \mathbb {R}^{p_d \times p_{(-d)}}$ , where $p_{(-d)} = \prod _{d'\neq d} p_{d'}$ , which columns correspond to mode-d fibers of the tensor $\mathcal {X}$ . We denote the vectorization operator, stacking the columns of a matrix $\mathbf {A}$ into a column vector, as $\operatorname *{\mathrm {vec}}(\mathbf {A})$ . In this context, we define the mode-d vectorization as $\mathbf {x}_{(d)} = \operatorname *{\mathrm {vec}}(\mathbf {X}_{(d)}) \in \mathbb {R}^{P\times 1}$ , where $P = \prod _{d} p_{d}$ . The mode- $1$ vectorization of $\mathcal {X}$ is sometimes simply denoted $\mathbf {x} = \mathbf {x}_{(1)}$ . For two matrices $\mathbf {A} \in \mathbb {R}^{p\times d}$ and $\mathbf {B} \in \mathbb {R}^{p' \times d'}$ , the Kronecker product between $\mathbf {A}$ and $\mathbf {B}$ is defined as

$$ \begin{align*} \mathbf{A} \otimes \mathbf{B} = \begin{bmatrix} a_{11} \mathbf{B} & \dots & a_{1d} \mathbf{B} \\ \vdots & \ddots & \vdots \\ a_{p1} \mathbf{B} & \dots & a_{pd} \mathbf{B} \end{bmatrix} \in \mathbb{R}^{pp'\times dd'}. \end{align*} $$

The Khatri–Rao product between two matrices with a fixed number of columns $\mathbf {A} = [\mathbf {a}_1 \dots \mathbf {a}_d] \in \mathbb {R}^{p\times d}$ and $\mathbf {B} = [\mathbf {b}_1 \dots \mathbf {b}_d]\in \mathbb {R}^{p' \times d}$ is defined as a column-wise Kronecker product:

$$ \begin{align*} \mathbf{A} \odot \mathbf{B} = \begin{bmatrix} \mathbf{a}_{1} \otimes \mathbf{b}_{1} & \dots & \mathbf{a}_{d} \otimes \mathbf{b}_{d} \end{bmatrix} \in \mathbb{R}^{pp'\times d}. \end{align*} $$

Finally, the Hadamard product between two matrices with similar dimensions $\mathbf {A} \in \mathbb {R}^{p\times d}$ and $\mathbf {B} \in \mathbb {R}^{p\times d}$ is defined as the element-wise product of the two matrices:

$$ \begin{align*} \mathbf{A} * \mathbf{B} = \begin{bmatrix} a_{11}b_{11} &\dots &a_{1d}b_{1d}\\ \vdots &\ddots &\vdots\\ a_{p1}b_{p1} &\dots &a_{pd}b_{pd} \end{bmatrix} \in \mathbb{R}^{p\times d}. \end{align*} $$

2.1.2 Functional notations

Considering an interval $\mathcal {I} \subset \mathbb {R}$ , we consider now the space $L^2(\mathcal {I})$ of square integrable functions on $\mathcal {I}$ (i.e., $\int _{\mathcal {I}} f(s)^2 \mathrm {d}s < \infty $ ). The inner product in $L^2(\mathcal {I})$ is defined for any $f \in L^2(\mathcal {I})$ and $g \in L^2(\mathcal {I})$ by $\langle f, g \rangle _{L^2} = \int _{\mathcal {I}} f(s)g(s) \mathrm {d}s$ . In the multivariate functional data setting (see Happ & Greven, Reference Happ and Greven2018), we can define the inner product between two multivariate functions $\mathbf {f} \in L^2(\mathcal {I}){}^p$ and $\mathbf {g} \in L^2(\mathcal {I}){}^p$ as

$$ \begin{align*} \langle \mathbf{f}, \mathbf{g} \rangle = \sum_{j=1}^p \langle f_j, g_j \rangle_{L^2} = \sum_{j=1}^p \int_{\mathcal{I}} f_j(s) g_j(s)\mathrm{d} s, \end{align*} $$

which defines a Hilbert Space. Considering now the space $L^2(\mathcal {I}){}^{p_1 \times \dots \times p_D}$ of tensors of order D with entries in $L^2(\mathcal {I})$ . We define the inner product between two functional tensors $\mathcal {F} \in L^2(\mathcal {I}){}^{p_1 \times \dots \times p_D}$ and $\mathcal {G} \in L^2(\mathcal {I}){}^{p_1 \times \dots \times p_D}$ as

(1) $$ \begin{align} \langle \mathcal{F}, \mathcal{G} \rangle = \sum_{\mathbf{j}} \langle f_{\mathbf{j}}, g_{\mathbf{j}} \rangle_{L^2} = \sum_{\mathbf{j}} \int_{\mathcal{I}} {f_{\mathbf{j}}}(s) {g_{\mathbf{j}}}(s) \mathrm{d} s, \end{align} $$

where we denote $f_{\mathbf {j}}(s) = (f(s))_{\mathbf {j}}$ .

2.2 Existing methods

The PARAFAC decomposition Carroll & Chang, Reference Carroll and Chang1970; Harshman, Reference Harshman1970 is a powerful and well-known method in the tensor literature. Its purpose is to represent any tensor as a sum of rank- $1$ tensors. In $\mathbb {R}^{p_1 \times \dots \times p_D}$ , the rank-R PARAFAC decomposition is defined as

(2)

where $\mathbf {a}_{dr} \in \mathbb {R}^{p_d}$ and $\mathbf {A}_d = [\mathbf {a}_{d1}, \dots , \mathbf {a}_{dR}] \in \mathbb {R}^{p_d \times R}$ . Here, each vector $\mathbf {a}_{dr}$ captures latent information along the mode d of the tensor. Consequently, the PARAFAC decomposition can be effectively employed as a dimension reduction technique, particularly compelling in the tensor setting, as the number of variables grows exponentially with the order of the tensor.

To approximate any tensor $\mathcal {X} \in \mathbb {R}^{p_1 \times \dots \times p_D}$ using a PARAFAC decomposition, factor matrices $\mathbf {A}_d$ are estimated by minimizing the squared difference between the original tensor $\mathcal {X}$ and its decomposition approximation,

(3)

The original optimization process used an alternating least squares (ALS) procedure Carroll & Chang, Reference Carroll and Chang1970; Harshman, Reference Harshman1970, where, iteratively, each factor matrix is updated while keeping the others fixed until convergence is reached. The procedure belongs to the more general class of block relaxation algorithms, known to be locally and globally convergent under mild assumptions Lange, Reference Lange2013. Note that PARAFAC can be seen as a generalization of the factor analysis framework Bartholomew et al., Reference Bartholomew, Knott and Moustaki2011, which is commonly used in psychometrics. Indeed, the loading vector associated with the statistical/sample mode (which we will describe after) can be seen as factors, and all other parameters can be used to derive loading factors.

A good property of the PARAFAC decomposition is that aside from scaling and order indeterminacies, the decomposition is generally unique Sidiropoulos & Bro, Reference Sidiropoulos and Bro2000. The scaling indeterminacy stems from the fact that, for $r = 1, \dots , R$ , when considering any set of scalars $c_{1r}, \dots , c_{Dr}$ such that $\prod _d c_{dr} = 1$ , using $\tilde {\mathbf {a}}_{dr} = c_{dr} \mathbf {a}_{dr}$ leads to a similar tensor as when using $\mathbf {a}_{dr}$ . This indeterminacy is usually fixed by scaling vectors $\mathbf {a}_{dr}$ for $d = 1, \dots , D-1$ and $r = 1, \dots , R$ , such that they have unit-norm. On the other hand, the order indeterminacy is caused by the fact that permuting columns of matrices $\mathbf {A}_d$ jointly does not change the resulting tensor. This indeterminacy can be fixed by ordering matrix columns such that $\| \mathbf {a}_{D1} \| < \| \mathbf {a}_{D2} \| <\dots < \| \mathbf {a}_{DR}\| $

Applying the PARAFAC model straightforwardly to decompose a tensor with smooth functional entries may be inefficient as the method does not integrate the smooth property. The PARAFAC decomposition also assumes a consistent structure across tensor modes. Therefore, it cannot handle schemes where functional entries are sampled irregularly, e.g., at different time points. To answer this last problem, several alternatives to the PARAFAC model were introduced, especially the PARAFAC2 model Harshman, Reference Harshman1972; Kiers et al., Reference Kiers, ten Berge and Bro1999, which was widely employed but still does not integrate any smoothness structure into the reconstructed tensor.

Another limitation of the PARAFAC model in our setting is that it is not a statistical model. Consequently, it does not integrate the potential sampling structure of the data. In the tensor literature, several approaches are able to integrate this property (see Giordani et al. (Reference Giordani, Rocci and Bove2020) for references). The SupCP model of Lock and Li (Reference Lock and Li2018) or the SupSVD model Li et al. (Reference Li, Yang, Nobel and Shen2016) recently went even further in this direction by allowing the introduction of covariates in the sampling modeling.

In the functional setting, the Karhunen–Loève theorem is often used to retrieve a low-dimensional representation of a high (or infinite) dimensional random process. The theorem guarantees the existence of an orthonormal basis ${\Phi _k}$ and a sequence ${\xi _k}$ of independent random variables for any squared integrable random process X, such that

(4) $$ \begin{align} X(t) = \sum_k^\infty \xi_k \Phi_k(t). \end{align} $$

By truncating the above sum to a limited number of functions, we can recover a good enough representation of X. More recently, this decomposition was extended to multivariate functional data, allowing the representation of any multivariate random process as a sum of multivariate orthonormal functions and scores. By arranging the sequence of independent random variables ${\xi _k}$ in descending order of variance, this decomposition can be viewed as an extension of the well-known principal component analysis (PCA). Similarly to the PCA, estimating parameters in model (4) relies on the eigen decomposition of the covariance operator of the random process. As evoked previously, numerous approaches have been proposed to achieve this. The most notable of which, Yao et al. (Reference Yao, Müller and Wang2005), proposed a linear least squares smoothing procedure with asymptotic guarantees, allowing estimation of the covariance surface in sparse and irregular data sampling settings. Furthermore, using simple assumptions on the joint distribution of principal components and residual errors, they introduce an efficient approach for estimating principal components.

When dealing with multivariate or tensor-structured functional data, it appears natural to extend this decomposition. As evoked in Section 1, various approaches have been introduced to extend this decomposition to multivariate functional data. In the tensor-structured functional data setting, vectorizing the tensor appears as a straightforward way to proceed. However, similar to before, as the number of processes grows exponentially with the order of the tensor, this procedure would be computationally unfeasible for large tensors, as it would require estimating orthonormal functions for each functional feature. In addition, the vectorization of the tensor would lead to the loss of the tensor structure, and the resulting representation would still be high-dimensional, causing the approach to be inefficient.

2.3 Proposed model

We propose to extend the rank R PARAFAC decomposition introduced in (2) to functional tensors with elements in $L^2(\mathcal {I})$ by further multiplying each rank-1 sub-tensor $\mathbf {a}_{1r} \circ \dots \circ \mathbf {a}_{Dr}$ by a functional term $\phi _{r} \in \mathscr {H}$ , where $\mathscr {H}$ is a subset of $L^2(\mathcal {I})$ ,

(5)

denoting $\boldsymbol {\phi } = (\phi _{1}, \dots , \phi _{r}) \in \mathscr {H}^R$ . We propose to name this decomposition the functional PARAFAC (F-PARAFAC) decomposition.

As discussed earlier, it is common for a functional tensor to have a sample mode, i.e., a mode on which sub-tensors can be seen as samples drawn from an underlying tensor distribution. Similar to the structural PARAFAC model Giordani et al., Reference Giordani, Rocci and Bove2020, we propose to account for the inherent randomness of random functional tensors with elements in $L^2(\mathcal {I})$ by multiplying each rank- $1$ functional tensor $\phi _{r}(t) (\mathbf {a}_{1r} \circ \dots \circ \mathbf {a}_{Dr})$ by a random variable $u_r$ ,

(6)

denoting $\mathbf {u} = (u_1, \dots , u_R) \in \mathbb {R}^{R}$ . We propose naming this model the latent functional PARAFAC (LF-PARAFAC) decomposition. A visualization of the model is given in Figure 1. Assumptions on the distribution of $\mathbf {u}$ are further discussed in Section 3.

Figure 1

Latent functional PARAFAC decomposition for an order-3 functional tensor.

Figure 1 Long description

The diagram flows from left to right.

* On the far left is a three-dimensional cube representing a tensor. It is labeled with dimensions n for height, p sub c for depth, and l for width. The front face of the cube contains three horizontal wavy lines and an ellipsis in the center.

* To the right of the cube is an equals sign, followed by a series of additive components.

* The first component consists of a vertical rectangle labeled u sub 1, a diagonal rectangle labeled c sub 1, and a wavy line labeled Phi sub 1. An arrow points from the vertical rectangle u sub 1 to a bell curve above it.

* A plus sign and an ellipsis lead to the final component on the right, which mirrors the first. It features a vertical rectangle u sub R, a diagonal rectangle c sub R, and a wavy line Phi sub R. An arrow also points from u sub R to a second bell curve above it.

* A long horizontal double-headed arrow spans across the top of both bell curves, labeled with the normal distribution notation N(0, Lambda).

Given a random functional tensor $\mathcal {X}$ , which we assume to be centered throughout this section, we are interested in finding the best LF-PARAFAC decomposition approximating $\mathcal {X}$ :

(7)

In this optimization problem, $\boldsymbol {\psi }^{\ast}$ is the operator $\mathcal {X} \mapsto \boldsymbol {\psi }^{\ast}(\mathcal {X}) = \mathbf {u}^{\ast}$ outputting the optimal (over the expectation) random vector $\mathbf {u}^{\ast} \in \mathbb {R}^{1 \times R}$ for the random functional tensor distribution $\mathscr {D}_{\mathcal {X}}$ followed by $\mathcal {X}$ . It implies that for any sample $\mathcal {X}_n$ obtained from $\mathscr {D}_{\mathcal {X}}$ , the criterion is optimal, in average, using the sample mode vector $\mathbf {u}_n = \boldsymbol {\psi }^{\ast}(\mathcal {X}_n)$ in its decomposition along with the other optimal parameters $\boldsymbol {\phi }^{\ast}, \mathbf {A}_1^{\ast}, \dots , \mathbf {A}_D^{\ast}$ . The asymmetry between $\mathbf {u}$ , which is random, and $\boldsymbol {\phi }, \mathbf {A}_1, \dots , \mathbf {A}_D$ , which are fixed, allows us to rewrite the optimization problem as

(8)

Proposition 1 gives the expression of $\Psi ^{\ast}$ , corresponding to the optimal matrix $\mathbf {u}$ of the optimization sub-problem on the right-hand side of (8).

Proposition 1. Denoting $\mathbf {A}_{(D)} = (\mathbf {A}_D \odot \dots \odot \mathbf {A}_{1})$ , and fixing $\boldsymbol {\phi }, \mathbf {A}_1, \dots , \mathbf {A}_D$ , for any $\mathcal {X}$ , the optimal operator $\boldsymbol {\psi }^{\ast}$ is given by

(9) $$ \begin{align} \boldsymbol{\psi}^{\ast}(\mathcal{X}) = \left ( \int \mathbf{x}(t)^\top (\mathbf{A}_{(D)} \odot \boldsymbol{\phi}(t)) \mathrm{d} t \right ) \left ( \int (\mathbf{A}_{(D)} \odot \boldsymbol{\phi}(t))^\top (\mathbf{A}_{(D)} \odot \boldsymbol{\phi}(t)) \mathrm{d} t \right )^{-1}. \end{align} $$

Plugging the optimal operator $\Psi ^{\ast}$ in model (6) eliminates the dependency on $\mathbf {u}$ , and the problem can now be seen as depending only on $\boldsymbol {\phi }, \mathbf {A}_1, \dots , \mathbf {A}_D$ :

(10)

Denoting collapsed functional covariance matrices $\boldsymbol {\Sigma }_{[f]}(s, t) = \mathbb {E} [\mathbf {x}(s) (\mathbf {I}_{P} \otimes \mathbf {x}(t)^\top )]$ and $\boldsymbol {\Sigma }_{[d]}(s, t) = \mathbb {E} [\mathbf {X}_{(d)}(s) (\mathbf {I}_{p_{(-d)}} \otimes \mathbf {x}(t)^\top )],$ where $\mathbf {I}_{P}$ and $\mathbf {I}_{p_{(-d)}}$ denote the identity matrix of order P and $p_{(-d)}$ respectively, for which an interpretation of $\boldsymbol {\Sigma }_{[f]}(s, t)$ and $\boldsymbol {\Sigma }_{[d]}(s, t)$ is given in Section 3.2, we can derive Corollary 1 from Proposition 1.

Corollary 1. Denoting $\boldsymbol {\Lambda }^{\ast} = \mathbb {E}[\mathbf {u}^{\ast}{\mathbf {u}^{\ast}}^\top ]$ the covariance matrix of $\mathbf {u}^{\ast}$ and $\mathbf {K}(t) = (\mathbf {A}_{(D)} \odot \boldsymbol {\phi }(t)) \left ( \int (\mathbf {A}_{(D)} \odot \boldsymbol {\phi }(t))^\top (\mathbf {A}_{(D)} \odot \boldsymbol {\phi }(t)) \right )^{-1}$ , we have

(11) $$ \begin{align} \boldsymbol{\Lambda}^{\ast} = \int \int \mathbf{K}(s)^\top \boldsymbol{\Sigma}_{[f]}(s, t) \mathbf{K}(t) \mathrm{d} s \mathrm{d} t. \end{align} $$

Using Corollary 1, it is possible to show that the optimization problem only depends on the optimal covariance $\boldsymbol {\Lambda }^{\ast}$ and other operators derive via an operator C. Expressions of this operator C are derived in Lemma A1 in the Appendix.

For a given tensor, this operator C depends on the collapsed covariance $\boldsymbol {\Sigma }_{[f]}$ . This dependency can be replaced by a dependency on $\boldsymbol {\Sigma }_{[d]}$ for any $d \in [D]$ , which gives multiple expressions for C, introduced in Lemma A1. Lemma A2 in the Appendix allows us to derive simple stationary equations, which are introduced in Proposition 2.

Proposition 2. For the functional mode, for fixed $\boldsymbol {\phi }$ , we have

(12) $$ \begin{align} \boldsymbol{\phi}^{\ast}(s) = \left ( \int \boldsymbol{\Sigma}_{[f]}(s, t) (\mathbf{A}_{(D)} \odot \mathbf{K}(t)) \mathrm{d} t \right ) \left ( (\mathbf{A}_{(D)}^\top \mathbf{A}_{(D)}) \ast \boldsymbol{\Lambda}^{\ast} \right )^{-1}. \end{align} $$

Denoting $\mathbf {A}_{(-d)} = (\mathbf {A}_D \odot \dots \odot \mathbf {A}_{d+1} \odot \mathbf {A}_{d-1} \odot \dots \odot \mathbf {A}_{1})$ , and $\mathbf {A}_{(-d)}^{\boldsymbol {\phi }}(t) = \mathbf {A}_{(-d)} \odot \boldsymbol {\phi }(t),$ we have for any $d\in [D]$ tabular mode, for fixed $\mathbf {A}_d$ ,

(13) $$ \begin{align} \mathbf{A}_d^{\ast} = \left ( \int \int \boldsymbol{\Sigma}_{[d]}(s, t) (\mathbf{A}_{(-d)}^{\boldsymbol{\phi}} (s) \odot \mathbf{K}(t)) \mathrm{d} s \mathrm{d} t \right ) \left ( \int \mathbf{A}_{(-d)}^{\boldsymbol{\phi}}(t)^\top \mathbf{A}_{(-d)}^{\boldsymbol{\phi}} (t) \mathrm{d} t \ast \boldsymbol{\Lambda}^{\ast} \right )^{-1}. \end{align} $$

2.4 Identifiability

As said previously, a well-known advantage of the PARAFAC decomposition is the identifiability of model parameters under mild assumptions. In the probabilistic setting, we define identifiability as the fact that if for two sets of model parameters $\boldsymbol {\Theta } = (\boldsymbol {\Lambda }, \boldsymbol {\phi }, \mathbf {A}_1, \dots , \mathbf {A}_D)$ and $\tilde {\boldsymbol {\Theta }} = (\tilde {\boldsymbol {\Lambda }}, \tilde {\mathbf {A}}_f, \tilde {\mathbf {A}}_1, \dots , \tilde {\mathbf {A}}_D)$ we have, $L( \boldsymbol {\Theta }| \mathcal {X} ) = L( \tilde {\boldsymbol {\Theta }} | \mathcal {X} )$ , for any data $\mathcal {X}$ , where $L( \boldsymbol {\Theta }| \mathcal {X} )$ stands for the model likelihood, then $\boldsymbol {\Theta } = \tilde {\boldsymbol {\Theta }}$ Lock & Li, Reference Lock and Li2018. In this context, similar to the standard setting, two well-known indeterminacies preventing parameters from being identifiable can be noted: the scaling and order indeterminacies. To remove these, we propose to consider the following assumptions:

1. (C1) For any $d\in [D]$ and $r\in [R]$ , $\| \mathbf {a}_{dr} \|_2 = 1$ , and the first non-zero value of $\mathbf {a}_{dr}$ is positive. Similarly, $\| \phi _{r} \|_{\mathcal {H}} = 1$ and $\phi _{r}$ is positive at the lower bound of its support.

2. (C2) Function $\phi _{r}$ and vectors $\mathbf {a}_{dr}$ are ordered such that sample-mode random variables $u_r$ have decreasing variance: $\boldsymbol {\Lambda }_{rr} < \boldsymbol {\Lambda }_{r'r'}$ for any $r < r'$ (we suppose no equal variance).

Considering the scaling and order indeterminacies fixed, a simple condition on the k-rank of the feature matrices and R was derived Sidiropoulos & Bro, Reference Sidiropoulos and Bro2000 for tensors of any order to guarantee identifiability. In the probabilistic setting, it can be shown that the any-order structural PARAFAC decomposition model Giordani et al., Reference Giordani, Rocci and Bove2020, with just the probabilistic property considered, is identifiable (see Lemma A3 in the Appendix). We propose to extend this identifiability result in the LF-PARAFAC setting. Proposition 3 gives conditions for which model (6) is identifiable. We recall that we denote $\operatorname {\mathrm {kr}}(\mathbf {A}_d)$ the k-rank of the feature matrix $\mathbf {A}_d$ and define it as the maximum number k such that any k columns of $\mathbf {A}_d$ are linearly independent.

Proposition 3. Denoting $\mathscr {H}$ the space to which belong functions $\phi _r$ , assuming assumptions (C1) and (C2) verified, model (6) is identifiable if $\mathscr {H}$ is infinite-dimensional and if the following inequality is verified:

(14) $$ \begin{align} \sum_{d=1}^D \operatorname{\mathrm{kr}}(\mathbf{A}_d) \geq D - 1. \end{align} $$

Conversely, assuming $\mathscr {H}$ is finite-dimensional, and, considering $\{ \psi _k \}_{1 \leq k \leq p}$ any basis of $\mathscr {H}$ , we can denote $\mathbf {C}$ the $R\times p$ basis decomposition matrix for functions $\{\phi _r\}_{1 \leq r \leq R}$ (i.e., with entries corresponding to $\langle \psi _k, \phi _r\rangle _{\mathscr {H}}$ ). In this context, model (6) is identifiable if deleting any row of $\mathbf {C} \odot \mathbf {A}_1 \odot \dots \odot \mathbf {A}_D$ results in a matrix having two distinct full-rank submatrices and if we have:

(15) $$ \begin{align} \operatorname{\mathrm{kr}}(\mathbf{C}) + \sum_{d=1}^D \operatorname{\mathrm{kr}}(\mathbf{A}_d) \geq R + (D - 1). \end{align} $$

Proposition 3 ensures that for a small enough R, the LF-PARAFAC decomposition will be identifiable. This result also implies that using high values of R often leads to unidentifiable and, therefore, poorly conditioned decomposition. Consequently, we recommend in practice considering not too large rank values for the decomposition to ensure identifiability and therefore stability of the interpretation. Additionally, note that since we enforce unit-norm constraints (C2), the number of free parameters in the decomposition is therefore reduced by a factor $RD$ .

3.1 Observation model

We now consider N samples of a random functional tensor $\mathcal {X}$ of order D, denoted $\{\mathcal {X}_n\}_{1 \leq n \leq N}$ . For a given $\mathcal {X}_n$ , we assume that entries of $\mathcal {X}_n$ are observed at similar time points, e.g., all functional entries of $\mathcal {X}_n$ are observed at the same time points. Observation time points of $\mathcal {X}_n$ are denoted $t_{nk}$ with $ 1 \leq k \leq K_n$ . The number of observations $K_n$ is allowed to be different across samples. In general, observations are contaminated with noise. Denoting $\mathcal {E}$ the random functional tensor modeling this noise, and $\mathcal {E}_{nk}$ the measurement error associated with $\mathcal {X}_{n}(t_{nk})$ , for any $t \in \mathcal {I,}$ we consider that $\mathcal {E}(t)$ has i.i.d. entries drawn from a distribution $\mathcal {N}(0, \sigma ^2)$ . The kth observation of the sample i, denoted $\mathcal {Y}_{nk}$ , is thus modeled as

(16) $$ \begin{align} \mathcal{Y}_{nk} = \mathcal{X}_{n}(t_{nk}) + \mathcal{E}_{nk}. \end{align} $$

As we aim to decompose $\mathcal {X}$ , we model the random functional tensor using the mean functional tensor $\mathcal {M} = \mathbb {E}[\mathcal {X}]$ and the LF-PARAFAC decomposition model as introduced in (6),

(17)

In practice, as introduced in Section 2.3, we, therefore, seek parameters of the LF-PARAFAC decomposition model that best fit the data given by observations $\{\mathcal {Y}_{nk}\}_{i,k}$ . The noise term $\mathcal {E}$ can also be seen as the unexplained part of the LF-PARAFAC decomposition. Furthermore, under the normality assumption on the error term $\mathcal {E}_{nk}$ entries, the identifiability result of Section 2.4, derived in the maximum likelihood setting, is similar for estimates obtained in the least squares setting (as seen in Section 2.3).

3.2 Covariance terms

In (12) and (13), we must estimate terms $\boldsymbol {\Sigma }_{[f]}(s, t) = \mathbb {E} [\mathbf {x}(s) (\mathbf {I}_{P} \otimes \mathbf {x}(t)^\top )]$ and $\boldsymbol {\Sigma }_{[d]}(s, t) = \mathbb {E} [\mathbf {X}_{(d)}(s) (\mathbf {I}_{p_{(-d)}} \otimes \mathbf {x}(t)^\top )]$ to retrieve optimal functions $\boldsymbol {\phi }^{\ast}$ and factor matrices $\mathbf {A}_d^{\ast}$ . As evoke previously, for any $d \in [D]$ , the term $\boldsymbol {\Sigma }_{[d]}(s, t)$ can be seen as a transformation of the mode-d functional covariance matrix $\boldsymbol {\Sigma }_{(d)}(s, t) = \mathbb {E}[\mathbf {x}_{(d)}(s) \mathbf {x}_{(d)}(t)^\top ] \in \mathbb {R}^{P \times P}$ . Recalling that $P = \prod _{d'} p_{d'}$ and $p_{(-d)} = \prod _{d'\neq d} p_{d'}$ , the transformation consists of a row-space collapse of non-mode-d entries into the column space, resulting in a matrix of dimensions $p_d\times P \cdot p_{(-d)}$ .

To better grasp the idea, let’s consider an order- $2$ functional tensor of dimensions $p_1\times p_2$ . The mode-1 functional covariance matrix can be decomposed as

$$ \begin{align*} \boldsymbol{\Sigma}_{(1)}(s,t) = \begin{bmatrix} \boldsymbol{\Sigma}_{(1, 1)}(s,t) & \dots & \boldsymbol{\Sigma}_{(1, p_2)}(s,t)\\ \vdots & \ddots & \vdots \\ \boldsymbol{\Sigma}_{(p_2, 1)}(s,t) & \dots & \boldsymbol{\Sigma}_{(p_2, p_2)}(s,t)\\ \end{bmatrix} \in \mathbb{R}^{p_1p_2 \times p_1p_2}, \end{align*} $$

where each block functional matrix $\boldsymbol {\Sigma }_{(j, j')}(s, t) = \mathbb {E}[\mathcal {X}_{.j}(s) \mathcal {X}_{. j'}(t)^\top ] \in \mathbb {R}^{p_1\times p_1}$ , where $\mathcal {X}_{. j}$ is the jth column of $\mathcal {X}$ , corresponds to functional covariance of mode- $1$ variables with fixed mode- $2$ variables $j, j' \in [p_2]$ . The functional matrix $\boldsymbol {\Sigma }_{[1]}(s, t)$ can then be expressed as

$$ \begin{align*} \boldsymbol{\Sigma}_{[1]}(s,t) = \begin{bmatrix} [\boldsymbol{\Sigma}_{(1, 1)}(s,t) & \dots & \boldsymbol{\Sigma}_{(1, p_2)}(s,t)] \; \dots \; [\boldsymbol{\Sigma}_{(p_2, 1)}(s,t) & \dots & \boldsymbol{\Sigma}_{(p_2, p_2)}(s,t)] \end{bmatrix} \in \mathbb{R}^{p_1 \times p_1p_2p_2}. \end{align*} $$

For $\boldsymbol {\Sigma }_{[2]}(s, t) \in \mathbb {R}^{p_2\times p_1 p_2 p_1}$ , a similar reasoning can help see the structure of the matrix with respect to the functional block entries of $\boldsymbol {\Sigma }_{(2)}$ .

Since functional matrices $\boldsymbol {\Sigma }_{(d)}$ can be obtained by reordering $\boldsymbol {\Sigma }$ , only $\boldsymbol {\Sigma }$ must be estimated to obtain an estimate of $\boldsymbol {\Sigma }_{[d]}$ . For this purpose, we rely on a local linear smoothing procedure outlined in Fan and Gijbels (Reference Fan and Gijbels2018). First, assuming functional tensor entries are not centered ( $\mathcal {M} \neq 0$ ), a centering pre-processing step is carried out. The mean function $m_{\mathbf {j}} = \mathbb {E}[y_{\mathbf {j}}]$ , which is estimated by aggregating observations and using the smoothing procedure evoked earlier, is removed from each functional entry of $\mathcal {Y}$ . Next, for each pair of centered functional entries $\mathcal {Y}_{\mathbf {j}}$ and $\mathcal {Y}_{\mathbf {j}'}$ , we estimate the (cross-)covariance surface $\Sigma _{\mathbf {j}\mathbf {j}'}(s, t)$ using a similar approach. The raw-(cross-)covariance is computed for each pair of observations within each sample. Raw (cross-)covariances are then aggregated, and the values are smoothed on a regular two-dimensional grid using the same smoothing procedure as before.

Additionally, we can estimate $\sigma $ by averaging the residual variance of functional entry. This residual variance is obtained for each process by running a one-dimensional smoothing of the diagonal using the previously removed pairs of observations and then averaging the difference between this smoothed diagonal and the previously computed smoothed diagonal (without the noise). For more details on this procedure, we refer the reader to the seminal paper by Yao et al. (Reference Yao, Müller and Wang2005).

This traditional smoothing approach was chosen for its simplicity and theoretical guarantees Yao et al., Reference Yao, Müller and Wang2005. However, several other approaches exist to estimate (cross-)covariance surfaces. Most notably, different smoothing procedures can be used Eilers & Marx, Reference Eilers and Marx2003; Wood, Reference Wood2003. For a univariate process, the fast covariance estimation (FACE) proposed by Xiao et al. (Reference Xiao, Zipunnikov, Ruppert and Crainiceanu2014) allows estimating the covariance surface with linear complexity (with respect to the number of observations). An alternative FACE method, robust to sparsely observed functional data, was introduced in Xiao et al. (Reference Xiao, Li, Checkley and Crainiceanu2017). More recently, the method was adapted to multivariate functional data in Li et al. (Reference Li, Xiao and Luo2020).

3.3 Solving procedure

Solution equations introduced in Proposition 2 along the expression of sample mode optimal vector covariance matrix introduced in Corollary 1 emphasize the use of a block relaxation algorithm outlined in Algorithm 1 to estimate parameters in (6).

Long description

The algorithm begins with a Result line defining Theta hat as a set containing Lambda hat, phi hat, and A sub 1 hat through A sub D hat.

Initialization: Theta hat super 0 is set to the initial values of Lambda, phi, and A sub 1 through A sub D, with the counter i set to 0.

The repeat loop contains the following steps:

1. Lambda super i plus 1 is updated using a double integral of K of s transpose, Sigma hat of s comma t, and K of t with respect to d s d t.

2. phi super i plus 1 is updated as the argmin over phi of the cost function C, taking Lambda super i plus 1 and the previous values of A sub 1 through A sub D as arguments.

3. A nested for loop runs from d equals 1 to D. Inside, A sub d super i plus 1 is updated as the argmin over A sub d of the cost function C. The arguments include the updated Lambda and phi, updated A values from 1 to d minus 1, the current A sub d, and the previous A values from d plus 1 to D.

4. The counter i is incremented by 1.

The loop terminates when the difference between the cost function C of Theta hat super i plus 1 and C of Theta hat super i is less than epsilon.

This algorithm uses all non-mode-d collapsed functional covariance matrices $\boldsymbol {\Sigma }_{[d]}$ as input. It produces estimates $\hat {\boldsymbol {\Lambda }}, \hat {\boldsymbol {\phi }}, \hat {\mathbf {A}_1}, \dots , \hat {\mathbf {A}_D}$ of the parameters in model (6).

Initialization can significantly influence the algorithm’s output, particularly in high-dimensional settings. In practice, random initialization is often chosen for simplicity, but it comes with the cost of running the algorithm multiple times. In our setting, we suggest initializing feature functions and matrices with the ones obtained from a standard PARAFAC decomposition. Experience shows that this simple approach is stable and provides good results. Note that alternative approaches were proposed in Guan (Reference Guan2023) and Han et al. (Reference Han, Shi and Zhang2023), but were not studied here.

3.4 Sample-mode inference

In this section, we consider that $\mathbf {u}_n = (u_{n1}, \dots , u_{nR})^\top \in \mathbb {R}^R$ . For a given sample $\boldsymbol {\mathcal {Y}}_n$ , when dealing with sparse and irregular observations, the estimation of $\mathbf {u}_n$ given by (9) can be unstable (or untracktable), since it requires integrating $\mathbf {x}_n$ on possibly very few observation time points. Furthermore, this estimator can give inconsistent estimates as it does not incorporate the error term $\mathcal {E}$ . Therefore, similar to the FPCA PACE-based procedure of Yao et al. (Reference Yao, Müller and Wang2005), we propose using a conditional approach to predict the factors $\mathbf {u}_n$ . In this context, we assume a joint Gaussian distribution on vectors $\mathbf {u}_{i}$ and residual errors $\mathbf {e}_{i}$ . Then, considering $\mathbf {u}_n \sim \mathcal {N}(\mathbf {0}, \boldsymbol {\Lambda }),$ we propose using the conditional expectation as our predictor of factors $\mathbf {u}_n$ . Denoting $\mathbf {y}_i$ the vector of observations of sample i, and $\mathbf {m}_{i}$ the vectorized mean tensor of $\mathcal {Y}$ , at time points $\mathbf {t}_i$ , Proposition 4 gives an expression of the conditional expectation of $\mathbf {u}_n$ .

Proposition 4. Denoting $\mathbf {u}_n = (u_{n1}, \dots , u_{nR}) \in \mathbb {R}^{1 \times R} \mathbf {A}_{f,n} = [ \boldsymbol {\phi }(t_{n1})^\top , \dots , \boldsymbol {\phi }(t_{nK_n})^\top ]^\top \in \mathbb {R}^{K_n \times R}$ , $\mathbf {F}_n = \mathbf {A}_{(D)} \odot \mathbf {A}_{f,n} \in \mathbb {R}^{K_n P \times R}$ , assuming joint Gaussian distribution on $\mathbf {u}$ and residual errors, a conditional estimator of $\mathbf {u}_n$ is

(18) $$ \begin{align} \tilde{\mathbf{u}}_n = \mathbb{E}[\mathbf{u}_n | \boldsymbol{\mathcal{Y}}_n] = \boldsymbol{\Lambda} \mathbf{F}_n^\top (\mathbf{F}_n \boldsymbol{\Lambda} \mathbf{F}_n^\top + \sigma^2\mathbf{I})^{-1} (\boldsymbol{{y}}_n - \mathbf{m}_n). \end{align} $$

Coupled with estimations $\hat {\boldsymbol {\Lambda }}, \hat {\boldsymbol {\phi }}, \hat {\mathbf {A}_1}, \dots , \hat {\mathbf {A}_D}$ obtained from Algorithm 1, and estimations $\hat {\sigma }$ and $\hat {\mathbf {m}}_n$ of $\sigma $ and $\mathbf {m}_n$ , obtained as described in Section 3.2, Proposition 4 allows to derive the following empirical Bayes estimator for $\mathbf {u}_n$ :

(19) $$ \begin{align} \hat{\mathbf{u}}_n = \hat{\boldsymbol{\Lambda}} \hat{\mathbf{F}}_n^\top (\hat{\mathbf{F}}_n \hat{\boldsymbol{\Lambda}} \hat{\mathbf{F}}_n^\top + \hat{\sigma}^2 \mathbf{I})^{-1} (\boldsymbol{{y}}_n - \hat{\mathbf{m}}_n). \end{align} $$