# Functional Generalized Canonical Correlation Analysis for studying multiple
longitudinal variables

Lucas Sort  Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des
Signaux et Systèmes, Gif-sur-Yvette, France Laurent Le Brusquet  Université
Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des Signaux et Systèmes, Gif-
sur-Yvette, France Arthur Tenenhaus  Université Paris-Saclay, CNRS,
CentraleSupélec, Laboratoire des Signaux et Systèmes, Gif-sur-Yvette, France

###### Abstract

In this paper, we introduce Functional Generalized Canonical Correlation
Analysis (FGCCA), a new framework for exploring associations between multiple
random processes observed jointly. The framework is based on the multiblock
Regularized Generalized Canonical Correlation Analysis (RGCCA) framework. It
is robust to sparsely and irregularly observed data, making it applicable in
many settings. We establish the monotonic property of the solving procedure
and introduce a Bayesian approach for estimating canonical components. We
propose an extension of the framework that allows the integration of a
univariate or multivariate response into the analysis, paving the way for
predictive applications. We evaluate the method’s efficiency in simulation
studies and present a use case on a longitudinal dataset.

Keywords: Longitudinal data, Functional Data, Generalized Canonical
Correlation Analysis

##  1 Introduction

Measuring multiple biomarkers jointly over time is common in observational
studies and clinical trials. As they characterize various biological processes
which are often interdependent, those biomarkers are usually correlated.
Hence, separately analyzing those longitudinal variables may hide parts of the
biological mechanisms at stake and give redundant information. Furthermore, as
subjects often miss one or more visits, the biomarkers may be observed
sparsely and irregularly. Therefore, along with the complex time-dependent
continuous structure of the data, the statistical analysis of multiple
longitudinal variables requires using specific methodologies for efficiently
harvesting information and integrating the interaction between the variables.

In the multivariate setting, the analysis of data coming from multiple
sources, usually represented by multiple sets of variables, is often referred
to as ”multiblock” or ”multi-set” analysis. Canonical Correlation Analysis
(CCA) (Hotelling (1936)) is one of the most notorious approaches, but it is
limited to exploring associations between only two sets of variables.
Therefore, various methods were proposed to generalize the CCA problem to more
than two sets of variables (Horst (1961); Carroll (1968); Kettenring (1971)).
More recently, Tenenhaus and Tenenhaus (2011) introduced Regularized
Generalized Canonical Correlation Analysis (RGCCA), a regularized and more
flexible framework for studying multiple sets of variables, giving birth to a
new generation of methods (Tenenhaus et al. (2014, 2015); Singh et al.
(2019)). Various extensions of RGCCA were proposed to handle emerging data
types, such as tensor data (Girka et al. (2023)). However, most are still
limited to finite-dimensional Euclidian spaces and are not designed to handle
large amounts of missing observations.

In the longitudinal literature, approaches based on linear mixed-effects
models have been widely employed over the past decades for studying
longitudinal biomarkers (Rizopoulos (2011); German et al. (2021)). However,
for adapting multivariate data analysis methods to the longitudinal setting,
functional approaches have been preferred as functional spaces are easy to
handle and allow to describe the underlying smooth structure of the time-
dependent variables. In this context, adaptations of Principal Component
Analysis (PCA) to the longitudinal setting have flourished (Rice and Silverman
(1991)). Most notably, Yao et al. (2005) proposed an adaptation using a
covariance-based procedure and a Bayesian approach to estimate the principal
components. The method is thus robust to sparse and irregular data, making it
applicable to numerous problems.

Multiple extensions of CCA were proposed to explore associations between two
longitudinal variables (Leurgans et al. (1993); He et al. (2003); Zhou et al.
(2008); Shin and Lee (2015)). Regularization is crucial in this infinite-
dimensional context as CCA requires inverting covariance matrices. Inspired by
Yao et al. (2005) approach, Yang et al. (2011) introduced the Functional
Singular Value Decomposition (FSVD), which moves the CCA criterion to a
covariance criterion and also uses a Bayesian approach to estimate the
canonical components. However, as in the multivariate setting, those
adaptations are limited to a pair of longitudinal variables. Few methods can
go beyond this limitation. To our knowledge, Hwang et al. (2011) proposed the
first approach to find associations between any number of longitudinal
variables using a homogenous components criterion. More recently, Górecki et
al. (2020) proposed to adapt Horst (1961) approach to functional spaces using
basis decomposition. Although it is not designed to explore association among
several longitudinal variables, it is worth mentioning Multivariate Functional
Principal Component Analysis (MFPCA) (Happ and Greven (2018)), designed to
retrieve the principal modes of variation on multivariate longitudinal data.
The method can also handle sparse and irregular data.

In this context, we propose Functional Generalized Canonical Correlation
Analysis (FGCCA), a framework based on RGCCA that allows exploring and
studying associations between several longitudinal markers in a flexible way.
The method proposed is robust to sparse and irregular longitudinal data.
Furthermore, it is designed so it can integrate a multivariate block in the
analysis to perform, for instance, supervised learning. As RGCCA, the
framework provided by FGCCA is so vast that it encompasses many existing
methods, notably Yao et al. (2005) FPCA, FSVD, or Functional Partial Least
Squares (FPLS) as presented by Preda et al. (2007).

The paper is organized as follows. First, in Section 2.1, we recall the
Regularized Generalized Canonical Correlation Analysis (RGCCA) framework.
Then, in Section 2.2, we introduce the Functional Generalized Canonical
Correlation Analysis (FGCCA) context and optimization problem. We present the
solving procedure of FGCCA in Section 2.3 and introduce a scheme for
retrieving higher-order functions and estimating components in Section 2.4 and
2.5 respectively. We validate our method on simulation studies in Section 3.1
and propose an application to a real dataset in Section 3.2. Finally, our
approach’s limitations and possible extensions are discussed in Section 4.

Proofs of propositions are given in Supplementary Materials. The code used to
run the experiments and the R implementation of FGCCA are freely available on
github: https://github.com/Sort-L/FGCCA-Code.

##  2 Method

###  2.1 Regularized Generalized Canonical Correlation Analysis (RGCCA)

Regularized Generalized Canonical Correlation Analysis (RGCCA) (Tenenhaus et
al. (2017)) is an optimization and statistical framework for studying
associations between multiple sets of random variables. Denoting
$\operatorname{x}_{1},\dots,\operatorname{x}_{J}$ the sets of respectively
$p_{1},\dots,p_{J}$ random variables, and $\boldsymbol{\Sigma}_{jj^{\prime}}$
the $p_{j}\times p_{j^{\prime}}$ matrix of (cross-)covariance between
$\operatorname{x}_{j}$ and $\operatorname{x}_{j^{\prime}}$, the RGCCA
optimization problem can be expressed as:

$\operatorname*{argmax}_{\operatorname{a}_{1},\dots,\operatorname{a}_{J}\in\Omega_{1}\times\dots\times\Omega_{J}}\sum_{j\neq
j^{\prime}}c_{j,j^{\prime}}g(\operatorname{a}_{j}^{\top}\boldsymbol{\Sigma}_{jj^{\prime}}\operatorname{a}_{j^{\prime}})\\\
$ (1)

where $\Omega_{j}$ is defined as
$\Omega_{j}=\left\\{\operatorname{a}_{j}\in\mathbb{R}^{p_{j}}\mid\operatorname{a}_{j}^{\top}\boldsymbol{M}_{j}\operatorname{a}_{j}=1\right\\}$,
with $\boldsymbol{M}_{j}$ being a symmetric-positive matrix, $g$ is a convex
differentiable function and the matrix $C=(c_{j,j^{\prime}})$ is a $J\times J$
symmetric matrix with positive elements specifying the desired connection
design to study associations between the blocks. Classically we set
$c_{j,j^{\prime}}=1$ if we want to consider the interaction between the blocks
$j$ and $j^{\prime}$ and $c_{j,j^{\prime}}=0$ otherwise. Additionally, we
often consider
$\boldsymbol{M}_{j}=\tau_{j}\boldsymbol{I}_{p_{j}}+(1-\tau_{j})\boldsymbol{\Sigma}_{jj}$
with $\tau_{j}\in[0,1]$ to interpolate smoothly the criterion between a
correlation criterion, when $\tau_{j}=0$, and a covariance criterion, when
$\tau_{j}=1$. In this context, the goal of RGCCA is to retrieve block weight
vectors $\operatorname{a}_{j}$ giving block components
$\operatorname{y}_{j}=\operatorname{x}_{j}^{\top}\operatorname{a}_{j}$, which
are a compromise of the information from each set of variables and the
information shared with the other sets of variables.

Finally, two strategies are often used to compute higher-level weight vectors
for each set of variables. The first strategy leads to new weight vectors
associated with components uncorrelated to the previous ones. The second
strategy yields new weight vectors orthogonal to the previous ones. Both
strategies require transforming the original sets of variables
$\operatorname{x}_{j}$ into new sets of variables
$\operatorname{x}_{j}^{\prime}$ called ”deflated” vectors. The first
transformation, more often used, consists in regressing out from each set
$\operatorname{x}_{j}$ its associated component
$\operatorname{y}_{j}=\operatorname{x}_{j}^{\top}a_{j}$: the transformation
can be written
$\operatorname{x}_{j}^{\prime}=\operatorname{x}_{j}-(a_{j}^{\top}\boldsymbol{\Sigma}_{jj}a_{j})^{-1}\boldsymbol{\Sigma}_{jj}a_{j}a_{j}^{\top}\operatorname{x}_{j}$.
The second transformation consists in projecting each set
$\operatorname{x}_{j}$ onto the orthogonal of the space spanned by the
previous weight vectors: it is defined as
$\operatorname{x}_{j}^{\prime}=\operatorname{x}_{j}-\operatorname{a}_{j}\operatorname{a}_{j}^{\top}\operatorname{x}_{j}$.
To retrieve new weight vectors and new components with the desired properties,
the solving procedure is rerun, replacing original vectors
$\operatorname{x}_{j}$ with deflated vectors $\operatorname{x}_{j}^{\prime}$.
The transformation equations, often called ”deflation” equations, can be used
repeatedly to retrieve multiple weight vectors and components.

As demonstrated in Tenenhaus et al. (2017), the framework of RGCCA is very
general and subsumes many notorious data analysis methods such as Principal
Component Analysis (PCA), Canonical Correlation Analysis (CCA) (Hotelling
(1936)), Partial Least Squares (PLS) regression (Wold et al. (2001)), and
Generalized Canonical Correlation Analysis (GCCA) (Carroll (1968), Horst
(1961), Kettenring (1971)), to name a few. Many extensions and adaptations
have been proposed to tackle a wide variety of problems, but, to our
knowledge, none exists to integrate the time continuous structure of
longitudinal data and, especially, to handle highly sparse and irregular
observations.

###  2.2 Functional Generalized Canonical Correlation Analysis (FGCCA)

We now consider multiple time-dependent variables, such as longitudinal
biomarkers. We propose to adapt the previous framework to functional spaces,
and more precisely square-integrable random processes, since random processes
can represent the time continuous structure of the data. Our approach, named
Functional Generalized Canonical Correlation Analysis (FGCCA) is now
introduced.

####  2.2.1 Notations and definitions

From now on, we consider $\operatorname{X}_{1},\dots,\operatorname{X}_{J}$,
$J$ square-integrable random processes defined on compact intervals of
$\mathbb{R}$, $I_{1},\dots,I_{J}$ respectively. Note that the random objects
are thus part of infinite-dimensional Hilbert spaces
$L^{2}(I_{1}),\dots,L^{2}(I_{J})$. In this context, we define for the process
$j$ the mean function $\mu_{j}$ as
$\mu_{j}(t)=\mathbb{E}(\operatorname{X}_{j}(t))$ for $t\in I_{j}$.
Additionally, the (cross-)covariance function (or ”surface”)
$\Sigma_{jj^{\prime}}$ between the processes $j$ and $j^{\prime}$ is defined
as
$\Sigma_{jj^{\prime}}(s,t)=\mathbb{E}((\operatorname{X}_{j}(s)-\mu_{j}(s))(\operatorname{X}_{j^{\prime}}(t)-\mu_{j^{\prime}}(t)))$
with $s\in I_{j}$ and $t\in I_{j^{\prime}}$. Finally, the (cross-)covariance
operator between the processes $j$ and $j^{\prime}$,
$\boldsymbol{\Sigma}_{jj^{\prime}}$, is defined as :

$\displaystyle\boldsymbol{\Sigma}_{jj^{\prime}}$
$\displaystyle:L^{2}(I_{j^{\prime}})\rightarrow L^{2}(I_{j}),\;f\mapsto
g,\;g(s)=\int_{I_{j^{\prime}}}\Sigma_{jj^{\prime}}(s,t)f(t)\text{dt}$

####  2.2.2 Model

Following the optimization problem (1), defining RGCCA, we define Functional
Generalized Canonical Correlation Analysis (FGCCA) optimization problem by
moving from the multivariate setting of $\mathbb{R}^{p_{j}}$ spaces to the
functional setting of $L^{2}(I_{j})$ spaces, replacing the sets of variables
$\operatorname{x}_{j}$ by the random processes $\operatorname{X}_{j}$ and,
therefore, the euclidean dot product $a^{\top}b$ by the functional scalar
product defined by $\langle f,g\rangle_{L^{2}}=\int fg$. The FGCCA
optimization problem can therefore be written:

$\operatorname*{argmax}_{f_{1},\dots,f_{J}\in\Omega_{1}\times\dots\times\Omega_{J}}\sum_{j\neq
j^{\prime}}c_{j,j^{\prime}}g(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})$ (2)

where $\Omega_{j}$ is defined as $\Omega_{j}=\left\\{f_{j}\in
L^{2}(I_{j})\mid\langle
f_{j},\boldsymbol{M}_{j}f_{j}\rangle_{L^{2}}=1\right\\}$ with
$\boldsymbol{M}_{j}$ being a symmetric positive-definite operator, and where
$g$ and $\mathbb{C}=(c_{j,j^{\prime}})$ are defined similarly as before.
Similarly to RGCCA, we suggest setting
$\boldsymbol{M}_{j}=\tau_{j}\boldsymbol{I}_{I_{j}}+(1-\tau_{j})\boldsymbol{\Sigma}_{jj}$.
However, to ensure the positive definiteness of the operator
$\boldsymbol{M}_{j}$, regularization parameters $\tau_{j}$ must be strictly
superior to $0$ (and thus, lie in $]0,1]$) as covariance operators
$\boldsymbol{\Sigma}_{jj}$ are not necessarily definite in the infinite-
dimensional setting. Functions $f_{j}$ and components
$\operatorname{y}_{j}=\langle\operatorname{X}_{j},f_{j},\rangle_{L^{2}}$,
allow capturing information for each process which, depending on the model
parameters, is a summary of both the information from each process and the
information shared with the others.

####  2.2.3 (Cross-)Covariance estimation

In the multivariate setting, the most straightforward estimation for
(cross-)covariance matrices $\boldsymbol{\Sigma}_{jj^{\prime}}$ is the sample
covariance matrix, which is fast and easy to compute. However, in the
functional setting, it is usually preferable to use alternative strategies
that integrate the data’s time-continuous structure.

Various methods have been proposed over the past decades to integrate this
structure. In the functional data analysis literature, (cross-)covariance
operators $\boldsymbol{\Sigma}_{jj^{\prime}}$ are often discretized and
estimated on dense & regular grids using kernel smoothing methods (Yao et al.
(2005), Yang et al. (2011)) or Generative Additive Models (GAM), which are
easy to implement. In this context, estimation methods often have several
hyperparameters that must be set. Usually, those hyperparameters are manually
specified prior to the analysis. However, cross-validation procedures, such as
leave-one-out cross-validation or criterion-based procedures, can be used for
selecting them (Leurgans et al. (1993)). Additionally, due to approximation
errors, estimated operators are rarely positive in practice, making the choice
of the regularization parameters $\tau_{j}$ in FGCCA crucial, as the interval
they are defined on may not be clearly identified. Consequently, we advise
setting regularization parameters to $1$ as it will always prevent the
optimization problem from being ill-posed.

Finally, we recommend normalizing the different processes before estimating
(cross-)covariance operators. Like in the multivariate setting, considerable
differences in the variance of the processes may lead to biased results. For
this purpose, we suggest using the normalization presented in a similar
context by Happ and Greven (2018), enforcing the integrated variance to be the
same for each process. To achieve this, each process $j$ is multiplied by the
following normalization quantity:

$w_{j}=\left(\int_{{I_{j}}}\operatorname*{var}(\operatorname{X}_{j}(t))\text{dt}\right)^{-1/2}$

###  2.3 Resolution

We now introduce a procedure to retrieve solutions to the FGCCA optimization
problem. Convergence properties are given, ensuring the stability of the
solving procedure.

####  2.3.1 Procedure

Let $\Psi$ be the objective function:

$\Psi(f_{1},\dots,f_{J})=\Psi(\boldsymbol{f})=\sum_{j\neq
j^{\prime}}c_{j,j^{\prime}}g(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})$

As $\Psi$ is a function of multiple arguments, we suggest using a block
coordinate ascent (BCA) strategy (de Leeuw (1994)) for finding solutions to
the maximization problem. This strategy consists in maximizing $\Psi$ argument
by argument until convergence is reached. The properties of $g$ implies that
the objective function is differentiable and multi-convex, meaning that it is
convex with respect to each argument $f_{j}$ when all the others are fixed.
Note also that we consider here the ”functional differentiability” since
$\Psi$ has functional arguments. A proof of the definition of the gradient is
given in Supplementary Materials. From these properties we can derive the
following inequality for $\tilde{f}_{j}\in\Omega_{j}$:

$\displaystyle\Psi(f_{1},\dots,\tilde{f}_{j},$
$\displaystyle\dots,f_{J})\geq\Psi(\boldsymbol{f})+\langle\nabla_{j}\Psi(\boldsymbol{f}),\tilde{f}_{j}-f_{j}\rangle=m_{j}(\boldsymbol{f},\tilde{f}_{j})$
(3)

where $\nabla_{j}\Psi(\boldsymbol{f})$ is the functional partial derivative of
$\Psi$ with respect to the $j$th function. With this expression, we notice
that maximizing $\Psi$ for the $j$th argument can be achieved by maximizing
the minorizing function $m_{j}$. In this minorizing function, only the term
$\langle\nabla_{j}\Psi(\boldsymbol{f}),\tilde{f}_{j}\rangle$ is relevant since
all the others are fixed. Therefore, the maximum of $m_{j}$ under the
constraint that $\tilde{f}_{j}\in\Omega_{j}$ is reached for:

$\hat{f}_{j}=\operatorname*{argmax}_{\tilde{f}_{j}\in\Omega_{j}}\langle\nabla_{j}\Psi(\boldsymbol{f}),\tilde{f}_{j}\rangle=\frac{\boldsymbol{M}_{j}^{-1}\nabla_{j}\Psi(\boldsymbol{f})}{||\boldsymbol{M}_{j}^{-1/2}\nabla_{j}\Psi(\boldsymbol{f})||}:=r_{j}(\boldsymbol{f})$
(4)

where the partial derivative can be expressed as:

$\nabla_{j}\Psi(\boldsymbol{f})=2\sum_{\begin{subarray}{c}j^{\prime}=1\\\
j^{\prime}\neq j\end{subarray}}^{J}c_{j,j^{\prime}}g^{\prime}(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle)\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}$
(5)

From this, we propose the Algorithm (1) for retrieving solutions to
optimization problem (13).

Algorithm 1 FGCCA algorithm

Data: $(\boldsymbol{\Sigma}_{jj^{\prime}})_{1\leq j,j^{\prime}\leq
J},\,g,C,\epsilon,\boldsymbol{f}^{0}$

Result: $f_{1}^{s+1},\dots,f_{J}^{s+1}$

repeat

for  _$j=1$ to $J$_ do

$\tilde{f}_{j}^{s+1}=\frac{\boldsymbol{M}_{j}^{-1}\nabla_{l}\Psi(f_{1}^{s+1},\dots,f_{j-1}^{s+1},f_{j}^{s},f_{j+1}^{s},\dots,f_{J}^{s})}{||\boldsymbol{M}_{j}^{-1/2}\nabla_{l}\Psi(f_{1}^{s+1},\dots,f_{j-1}^{s+1},f_{j}^{s},f_{j+1}^{s},\dots,f_{J}^{s})||}$

end for

$s=s+1$

until  _$\Psi(f_{1}^{s+1},\dots,f_{J}^{s+1})-\Psi(f_{1}^{s},\dots,f_{J}^{s})
<\epsilon$_;

####  2.3.2 Monotone convergence

Denoting $\Omega=\Omega_{1}\times\dots\times\Omega_{J}$ we define
$c_{j}:\Omega\rightarrow\Omega$ as the operator
$c_{j}(\boldsymbol{f})=(f_{1},\dots;f_{j-1},r_{j}(\boldsymbol{f}),f_{j+1},\dots,f_{J})$
with $r_{j}(\boldsymbol{f})$ being the update function for the $j$th function
of the solving procedure in Section 2.3.1. We also define
$c:\Omega\rightarrow\Omega$ as the operator $c=c_{J}\circ\dots\circ c_{1}$

We consider the sequence
$\\{\boldsymbol{f}^{s}=(f_{1}^{s},\dots,f_{J}^{s})\\}$ generated by
$\boldsymbol{f}^{s+1}=c(\boldsymbol{f}^{s})$. The following proposition states
the monotone convergence of the generated sequence
$\\{\boldsymbol{f}^{s}\\}^{\infty}_{s=0}$ and holds as long as the update
$r_{j}(\boldsymbol{f})$ exists, is unique and $\Omega$ is bounded:

######  Proposition 1.

Considering any sequence $\\{\boldsymbol{f}^{s}\\}_{s=0}^{\infty}$ generated
recursively by the relation $\boldsymbol{f}^{s+1}=c(\boldsymbol{f}^{s})$ with
$\boldsymbol{f}^{0}\in\Omega$. The sequence $\\{\Psi(\boldsymbol{f}^{s})\\}$
is monotonically increasing and therefore convergent as $\Psi$ is bounded on
$\Omega$, implying the convergence of the FGCCA algorithm.

Using this proposition, since $\Omega$ is bounded and $r_{j}(\boldsymbol{f})$
properly and uniquely defined for all $j$ with the functional gradient, we can
conclude that the Algorithm (1) is monotone and convergent.

###  2.4 Retrieving higher-order orthogonal functions

Solving the optimization problem as described previously only yields one
function per block. However, it is often preferable to retrieve multiple
functions leading to multiple components. For this purpose, we suggest, as in
the RGCCA framework, using a deflation strategy.

First, we propose considering a deflation strategy for retrieving orthogonal
vectors. As detailed in Section 2.1, the deflation equation associated with
this strategy is
$\operatorname{x}_{j}^{\prime}=\operatorname{x}_{j}-\operatorname{a}_{j}\operatorname{a}_{j}^{\top}\operatorname{x}_{j}$.
This expression may be unusable in the functional setting, with sparse and
irregular observations, as processes $\operatorname{X}_{j}$ are not fully
observed. Moreover, only the (cross-)covariance operators are involved in the
solving procedure, making the deflation of the blocks appear as an unnecessary
step in the procedure. Therefore, we establish the following proposition,
allowing us to deflate the (cross-)covariance operators directly without
involving the possibly ill-defined blocks:

######  Proposition 2.

Denoting $\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}$ the deflated
(cross-)covariance operator of $\boldsymbol{\Sigma}_{jj^{\prime}}$, obtained
after projecting the processes onto the orthogonal of the space spanned by
their associated vectors. The following equality holds:

$\displaystyle\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}$
$\displaystyle=(\boldsymbol{I}_{I_{j}}-\boldsymbol{\Phi}_{j})\boldsymbol{\Sigma}_{jj^{\prime}}(\boldsymbol{I}_{I_{j^{\prime}}}-\boldsymbol{\Phi}_{j^{\prime}})$
(6)

where $\boldsymbol{\Phi}_{j}:L^{2}(I_{j})\rightarrow L^{2}(I_{j})$ is the
operator defined by :

$(\boldsymbol{\Phi}_{j})(f)=(f_{j}\otimes f_{j})(f)=\langle f_{j},f\rangle
f_{j}$

To retrieve new orthogonal functions, deflated operators are plugged into the
optimization problem and Algorithm 1 is run again. The deflation procedure can
be repeated multiple times, allowing to retrieve a set of canonical functions
$\\{f_{j}^{m}\\}_{1\leq m\leq M}$ for each random process. The number of
canonical functions to retrieve is often manually set but it can be chosen
using cross-validation or criterion-based approaches. Finally, in this
context, the set of components $\\{\operatorname{y}_{j}^{m}\\}_{1\leq m\leq
M}$ obtained for each process can be directly estimated from the original
processes $\operatorname{X}_{j}$ without using the deflated processes
$\operatorname{X}_{j}^{m}$, a desirable property in our setting as it may be
difficult to compute and manipulate the deflated processes. Indeed, we can
demonstrate easily that
$\langle\operatorname{X}_{j},f_{j}^{m}\rangle_{L^{2}}=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle_{L^{2}}=\operatorname{y}_{j}^{m}$,
where $\operatorname{X}_{j}^{m}$ stands for the $m$th deflation of process
$j$.

###  2.5 Estimating components

Computing the components
$\operatorname{y}_{j}^{m}=\langle\operatorname{X}_{j},f_{j}^{m}\rangle_{L^{2}}$
may be difficult in the sparse and irregular setting, as the numerical
estimation of the $L^{2}$ scalar product can be unstable and untracktable,
particularly if the number of observations is small. In this context, inspired
from Yao et al. (2005) and Yang et al. (2011), we propose to estimate the
components using a Bayesian approach.

####  2.5.1 Notations

In the following, subscripts $i$, $j$, $k$ denote respectively the subject
number, the process number, and the observation number. We denote $n_{ij}$ the
number of observations,
$\mathbb{U}_{ij}=(U_{ij1},\dots,U_{ijn_{ij}})^{\top}\in\mathbb{R}^{n_{ij}\times
1}$ the observations, and $t_{ij}=(t_{ij1},\dots,t_{ijn_{ij}})$ the
observation time points. Finally the observations are modeled as:

$U_{ijk}=X_{ij}(t_{ijk})+\varepsilon_{ijk}$ (7)

where $X_{ij}$ is the realization for the subject $i$ of the random process
$\operatorname{X}_{j}$ and $\varepsilon_{ijk}$ is a measurement error. The
measurement errors are supposed i.i.d and following a normal distribution
$\mathcal{N}(0,\sigma_{j}^{2})$.

####  2.5.2 Process modeling

As previously stated, considering the set of orthonormal canonical functions
$\\{f_{j}^{m}\\}_{1\leq m\leq M}$, each process $j$ from any subject $i$ can
be decomposed as:

$X_{ij}(t)=\mu_{j}(t)+\sum_{m=1}^{M}\xi_{ij}^{m}f^{m}_{j}(t)$ (8)

where the coefficients $\xi_{ij}^{m}$ are the basis coefficients associated
with the basis $\\{f_{j}^{m}\\}_{1\leq m\leq M}$. Therefore, at the sample
level we have:

$U_{ijk}=\mu_{j}(t_{ijk})+\sum_{m=1}^{M}\xi_{ij}^{m}f^{m}_{j}(t_{ijk})+\varepsilon_{ijk}$
(9)

This formulation allows to see the basis decomposition as a linear mixed-
effects model with the fixed-effects part being the mean term and the random-
effects part being the decomposition term. Moreover, as $\xi_{ij}^{m}=\langle
X_{ij},f_{j}^{m}\rangle_{L^{2}}$ and since the deflation strategy introduced
in Section 2.4 leads as previously stated to
$\langle\operatorname{X}_{j},f_{j}^{m}\rangle_{L^{2}}=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle_{L^{2}}=\operatorname{y}_{j}^{m}$
we have that $\xi_{ij}^{m}=\operatorname{y}_{ij}^{m}$. Therefore, estimating
basis coefficients is equivalent to estimating components.

####  2.5.3 Estimation

For simplifying expressions, denoting $N_{i}=\sum_{j}n_{ij}$, we write the
vector of observations
$\mathbb{U}_{i}=(\mathbb{U}_{i1}^{\top},\dots,\mathbb{U}_{iJ}^{\top})^{\top}\in\mathbb{R}^{N_{i}\times
1}$, and the mean function vector (at the observation time points)
$\boldsymbol{\mu}_{i}=(\boldsymbol{\mu}_{i,1}^{\top},\dots,\boldsymbol{\mu}_{i,J}^{\top})^{\top}\in\mathbb{R}^{N_{i}\times
1}$ with
$\boldsymbol{\mu}_{i,j}=(\mu_{i}(t_{ij1}),\dots,\mu_{i}(t_{ijn_{ij}}))^{\top}\in\mathbb{R}^{n_{ij}\times
1}$. We also write
$\mathbb{F}_{ij}^{m}=(f_{j}^{m}(t_{ij1}),\dots,f_{j}^{m}(t_{ijn_{ij}}))^{\top}\in\mathbb{R}^{n_{ij}\times
1}$ and
$\mathbb{F}_{ij}=(\mathbb{F}_{ij}^{1},\dots,\mathbb{F}_{ij}^{M})^{\top}\in\mathbb{R}^{n_{ij}\times
M}$, the matrix of the $M$ canonical functions at the observation time points
for subject $i$, process $j$, and finally
$\boldsymbol{\xi}_{j}=(\xi_{j}^{1},\dots,\xi_{j}^{M})^{\top}\in\mathbb{R}^{M\times
1}$,
$\boldsymbol{\xi}=(\boldsymbol{\xi}_{1}^{\top},\dots,\boldsymbol{\xi}_{J}^{\top})^{\top}\in\mathbb{R}^{MJ\times
1}$, the vector of basis coefficients. Considering that the basis coefficients
and the measurement errors are centered and jointly Gaussian, we establish the
following proposition, allowing estimating coefficients:

######  Proposition 3.

Denoting
$\mathbb{F}_{i}=\operatorname*{diag}(\mathbb{F}_{i1},\dots,\mathbb{F}_{iJ})$,
$\mathbb{\Sigma}=\mathbb{E}[\boldsymbol{\xi}\boldsymbol{\xi}^{\top}]$ and
$\boldsymbol{\sigma}_{i}=\operatorname*{diag}(\sigma_{1}^{2}\mathbb{I}_{n_{i1}},\dots,\sigma_{J}^{2}\mathbb{I}_{n_{iJ}})$,
the best linear unbiaised predictor (BLUP) for $\boldsymbol{\xi}_{i}$ is given
by

$\mathbb{E}(\boldsymbol{\xi}_{i}|\boldsymbol{U}_{i})=\mathbb{\Sigma}\mathbb{F}_{i}^{\top}(\mathbb{F}_{i}\mathbb{\Sigma}\mathbb{F}_{i}^{\top}+\boldsymbol{\sigma}_{i})^{-1}(\boldsymbol{U}_{i}-\boldsymbol{\mu}_{i})$
(10)

In this expression, all the terms can be estimated. Notably, the canonical
functions can be interpolated at the observation time points of each subject
allowing estimating matrices $\mathbb{F}_{i}$. The noise standard deviation
$\sigma_{j}$ can be approximated for each process using the estimated
covariance surface of each process (a procedure is presented in Yao et al.
(2005)). Finally, the mean functions $\boldsymbol{\mu}_{i}$ are usually
estimated when estimating (cross-)covariance surfaces, often using smoothing
techniques or GAMs.

###  2.6 Retrieving higher-order uncorrelated components

The deflation strategy introduced in Section 2.4 allows recovering multiple
orthogonal functions for each process. However, as introduced in section 2.1,
retrieving uncorrelated components is often preferable. In the multivariate
setting, this deflation strategy is carried out using the following deflation
equation
$\operatorname{x}_{j}^{\prime}=\operatorname{x}_{j}-(a_{j}^{\top}\boldsymbol{\Sigma}_{jj}a_{j})^{-1}\boldsymbol{\Sigma}_{jj}a_{j}a_{j}^{\top}\operatorname{x}_{j}$.
As previously, we propose to adapt this strategy in the functional setting to
deflate the (cross-)covariance operators directly. For this purpose, we
establish the following proposition:

######  Proposition 4.

Denoting $\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}$ the deflated
(cross-)covariance operator of $\boldsymbol{\Sigma}_{jj^{\prime}}$, obtained
from regressing out the components from their associated block. The following
equality holds :

$\displaystyle\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}$
$\displaystyle=(\boldsymbol{I}_{I_{j}}-d_{j}\boldsymbol{\Sigma}_{jj}\boldsymbol{\Phi}_{j})\boldsymbol{\Sigma}_{jj^{\prime}}(\boldsymbol{I}_{I_{j^{\prime}}}-d_{j^{\prime}}\boldsymbol{\Phi}_{j^{\prime}}\boldsymbol{\Sigma}_{j^{\prime}j^{\prime}})$
(11)

where $d_{j}=(\operatorname{y}_{j}^{\top}\operatorname{y}_{j})^{-1}$ and, as
previously, $\boldsymbol{\Phi}_{j}:L^{2}(I_{j})\rightarrow L^{2}(I_{i})$ is
the operator defined by :

$(\boldsymbol{\Phi}_{j})(f)=(f_{j}\otimes f_{j})(f)=\langle f_{j},f\rangle
f_{j}$

As previously, new functions associated with uncorrelated components can be
obtained by replacing (cross-)covariance operators in the optimization problem
by deflated ones and running the solving procedure. However, to obtain
uncorrelated estimates of $\operatorname{y}_{j}^{m}$, additional steps are
required. Indeed, the equality
$\langle\operatorname{X}_{j},f_{j}^{m}\rangle_{L^{2}}=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle_{L^{2}}=\operatorname{y}_{j}^{m}$,
which was previously used to estimate the components in the mixed-effects
model, no longer holds in this setting. Nevertheless, the orthogonal property
of the retrieved functions is still holding by construction. Therefore,
canonical functions $f_{j}^{m}$ can still be used as a decomposition basis in
the mixed-effects framework presented previously. Furthermore, using deflation
equations, basis coefficients and components can be linked with the following
recursive equation:

$\operatorname{y}_{j}^{m+1}=\xi_{j}^{m+1}-\sum_{k=1}^{m}P_{j}^{k}\xi_{j}^{m+1}$
(12)

With
$P_{j}^{k}=({\operatorname{y}_{j}^{k}}^{\top}\operatorname{y}_{j}^{k})^{-1}\operatorname{y}_{j}^{k}{\operatorname{y}_{j}^{k}}^{\top}$
being the projection matrix from the regression of $\xi_{j}^{k}$ on
$\operatorname{y}_{j}^{k}$, starting with
$\xi_{j}^{1}=\operatorname{y}_{j}^{1}$. This equation can be seen as a
decorrelation procedure: for each process, each new basis coefficient estimate
is decorrelated from the previous components.

Finally, the choice of the deflation type depends on the desired usage of the
components and the canonical functions. For reconstructing trajectories,
orthogonal functions deflation may be more adapted as the orthonormal
functions retrieved are preferable for decomposition purposes. For doing
clustering or dimension reduction for further analysis, such as regression,
using uncorrelated components seem more suitable.

###  2.7 Integrating a multivariate response

Inspired from the PLS-framework, we propose to modify slightly the FGCCA
optimization problem to include a multivariate response
$\operatorname{Y}\in\mathbb{R}^{p}$ with $p\in\mathbb{N}^{*}$ for expanding
the possibilities given by the framework:

$\operatorname*{argmax}_{\begin{subarray}{c}f_{1},\dots,f_{J},\in\Omega_{1}\times\dots\times\Omega_{J}\\\
||a||_{2}=1\end{subarray}}\sum_{j\neq j^{\prime}}c_{j,j^{\prime}}g(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})+2\sum_{j}g(\langle
f_{j},\boldsymbol{\Sigma}_{j\operatorname{Y}}a\rangle_{L^{2}})$ (13)

where $\boldsymbol{\Sigma}_{j\operatorname{Y}}$ is the cross-covariance
operator between the process $j$ and the response $\operatorname{Y}$, defined
as:

$\displaystyle\boldsymbol{\Sigma}_{j\operatorname{Y}}$
$\displaystyle:\mathbb{R}^{p}\rightarrow L^{2}(I_{j}),\;a\mapsto
g,\;g(s)=\sum_{i=1}^{p}\mathbb{E}[\operatorname{X}_{j}(t)\operatorname{Y}_{i}]a_{i}$

As previously, this operator can be estimated with kernel smoothing methods or
GAMs. This additional interaction allows recovering canonical functions for
each process that explains best the process, interaction with the others
(depending on the design matrix $C$) and interaction with the response. This
design is particularly relevant in a predictive framework as it could be
interesting to use the components retrieved to predict the response
$\operatorname{Y}$.

Finally, note that the various procedures presented previously are not
significantly affected by this change and can be easily rewritten to integrate
the multivariate vector using the multivariate-functional cross-covariance
operator presented above.

##  3 Results

###  3.1 Simulations studies

####  3.1.1 Simulation 1: validating the Bayesian approach.

For $J=3$ processes, we propose to compare the components retrieved using the
Bayesian approach previously described in Section (2.5) to the components
obtained by computing the scalar product, as usually done. For this purpose,
we generate data according to Equation (8) so that true component values are
known. We choose the first $M=6$ Fourier basis functions in the $[0,1]$
interval as our orthonormal basis. The components $\xi_{ij}^{m}$ for each
subject are generated jointly with a centered Gaussian distribution of
covariance $\Sigma$ having a decreasing variance structure. For each subject
and each process, the time points are generated by sparsifying a grid of size
$50$ in the $[0,1]$ interval. Various sparsity levels are compared: Dense (100
% of observations retained), Low Sparsity (100 % to 80 % of observations
retained), Medium Sparsity (80 % to 40 % of observations retained), and High
Sparsity (40 % to 10 % of observations retained). We compare the two
approaches by computing the mean squared error on the canonical components for
100 simulations. Results are reported in Figure 1.

We observe a clear advantage of the Bayesian approach in the Medium and in the
High sparsity settings. The gap between the two approaches increases for
higher-order components. Furthermore, we notice a slight advantage of the
Bayesian approach in the Dense and in the Low sparsity cases, again,
especially for the higher-order components. Therefore, we can conclude that
the estimation error when using the Bayesian approach is smaller than the
estimation error due to the numerical approximation of the integral when using
the standard scalar product approach. Thus, we advise using the Bayesian
approach in all cases as it is also computationally inexpensive.

Figure 1: (top) Mean squared error (MSE) boxplots of the components
$\operatorname{y}_{j}^{m}$ estimated with the Bayesian approach or the scalar
product (integration) for $m=1,2,3,4,5,6$ obtained from $100$ simulations with
$N=100$ and averaged over the $J=3$ processes, $\sigma^{2}=1$ and various
sparsity settings.

####  3.1.2 Simulation 2: comparing results

To validate FGCCA further, we suggest to compare it to some of its subsumed
methods. We propose to consider two other approaches for $J=2$ processes:
PACE-based Functional Principal Component Analysis (FPCA) (Yao et al. (2005))
and Functional Singular Value Decomposition (FSVD) (Yang et al. (2011)). Both
approaches can handle sparse and irregular data. The first gives similar
results to FGCCA for specific component’s covariance designs. The latter is
based on an optimization problem equivalent to FGCCA in the 2 processes
setting when regularization parameters $\tau_{j}$ are set to 1. The simulation
setting previously described is used again to generate data. The covariance
matrix $\Sigma$ is designed so FPCA, FSVD and FGCCA recover components and
functions in the same order. Results are reported on Figure 2.

As FGCCA and FSVD are supposed to retrieve similar canonical functions, we
observe similar distributions over the mean squared errors for functions.
Components estimation rely on a slightly different formula for FSVD. Indeed,
the singular values retrieved from the analysis are used, which should, in
theory, increase estimation accuracy. However, it seems that in our case, the
estimations obtained from a FSVD are significantly worse than with FGCCA,
especially for high and medium sparsity settings and for high-order functions.
On another hand, FPCA gives slightly better function estimations both for
functions and components, especially in the sparse setting and for the first
functions and components. For high order components, FGCCA sometimes
outperforms FPCA. Additional results along with simulation details are
presented in Supplementary Materials.

Figure 2: Mean squared errors (MSE) of functions $f^{m}$ (top) and components
$\xi^{m}$ (bottom) for $m=1,2,3,4,5,6$ obtained from $100$ simulations with
$N=100$, $\sigma^{2}=1$ and various sparsity settings. Comparison between
FPCA, FSVD and FGCCA.

####  3.1.3 Simulation 3 : comparing reconstruction

Alternatively, we propose to compare the estimation quality of the
reconstructed trajectories between FGCCA with a fully-connected design and an
orthogonal deflation, and Multivariate Functional Principal Component Analysis
(MFPCA) (Happ and Greven (2018)). The reconstructed trajectories are obtained
for FGCCA using the decomposition equation (8) with the estimated canonical
functions and components, and for MFPCA using the multi-dimensional Karhunen-
Loeve decomposition, presented in the aforementioned paper. This time, using
the package funData (Happ-Kurz (2020)), 3 processes are generated based on the
first $M=6$ Fourier basis functions and with a linear decreasing variance over
the components.

The mean squared relative error is compared over 100 simulations in various
configurations. Results are presented on Figure 3. We can see that our
approach improves slightly the reconstruction of the processes, especially
when the number of subject $N$ is not too small. The gap between the two
methods appears to be stable among all sparsity settings. Further results,
available in the Supplementary Materials, show that the difference is bigger
as the noise $\sigma^{2}$ is smaller.

Figure 3: Mean relative squared errors (MRSE) of reconstructed trajectories,
using estimated canonical functions and components. Comparing FGCCA and MFPCA
with $M=6$, $\sigma^{2}=1$, for various number of subjects $N$ and various
sparsity settings. Statistical significance displayed : (***) $p<0.001$ (****)
$p<0.0001$

###  3.2 Application to Primary Biliary Cirrhosis dataset

The Primary Biliary Cirrhosis dataset (Murtaugh et al. (1994)) is a dataset
from Mayo Clinic containing the follow-up of various biomarkers extracted from
blood analyses of 312 patients who have been diagnosed with primary biliary
cirrhosis of the liver, a rare autoimmune disease. We propose to use this
multi-biomarker dataset to show various usages of FGCCA. For this purpose,
three biomarkers were considered: albumin, bilirubin, which is log
transformed, and prothrombin time, observed up to 10 years after the first
visit. Those biomarkers were chosen as they have been proven to be good
predictors of patient outcomes. Figure 4 represents an aggregated view of the
data.

Figure 4:  Longitudinal trajectories for albumin, bilirubin and prothrombin
time in the pbc2 dataset for all individuals. As usually done, the bilirubin
marker is log transformed.

####  3.2.1 Exploratory analysis

We first propose visualizing the canonical functions and components obtained
with FGCCA when using a fully connected design and deflation leading to
uncorrelated components. We compare the results to the principal functions and
components obtained using PACE-based FPCA (Yao et al. (2005)). For both
methods, the bandwidths are manually set to 1 for interpretability. The
results are displayed in Figure 5.

The first principal and canonical functions have a similar flat shape,
implying that the difference in the trajectories between subjects for all
biomarkers comes primarily from an overall shift around the mean. On the other
hand, the second canonical and principal functions are either decreasing or
increasing during the 10 years interval, indicating that the next source of
variation between subjects thus comes from the monotonicity of the
trajectories: for the different biomarkers, subjects have either an increasing
or decreasing trend. Additionally, we notice a more significant difference
between principal and canonical functions for prothrombin, suggesting that
this biomarker is particularly correlated to the others. The differences are,
however, difficult to interpret. Additionally, note that the functions
retrieved with FGCCA are slightly smoother, notably at the end of the
interval. We can explain this by the border effects when estimating the
(cross-)covariance operators. Indeed, for FGCCA, the functions are estimated
using information from multiple processes and not just one (as done in FPCA),
leading to more stable and reliable results.

Component plots allow us to see the differences between the two approaches
more clearly. First, we notice that the components are spread more evenly for
FGCCA, especially for prothrombin, thanks to the normalization. For the three
biomarkers, the components given by FGCCA seem to separate better the two
outcomes. This property will be confirmed in the predictive analysis.

Figure 5: (top) First 3 functional modes retrieved by FGCCA (canonical
functions) and FPCA (principal functions), for each biomarker. Functions were
flipped to minimize the differences between the two methods. (bottom) Biplots,
for each biomarker, of the first 2 components obtained with FGCCA and FPCA
coloured by final status. Ellipses represent the estimated Gaussian
distributions of the components for each outcome.

####  3.2.2 Prediction

Inspired by Singh et al. (2019), we propose using a multiblock functional PLS
framework to predict each patient’s outcome and demonstrate the ability of
FGCCA to integrate a multivariate or univariate response. The multiblock
functional PLS design is defined as a FGCCA design where only the associations
with the response are considered in the problem. It allows to recover
biomarker information correlated to the response, which is particularly useful
for predictive purposes. A simple logistic regression model is fitted per
biomarker to predict the response using the first component retrieved with
FGCCA. For a new subject, the components are predicted using the observed
trajectories. The final prediction is obtained by computing a weighted average
of the predicted outcomes, where each biomarker prediction is weighted by the
correlation of its component with the response. The predictive performances
obtained are compared to a similar model where the principal components from
FPCA are used instead.

The results, summarized in Figure 6, show that the FGCCA-based components
provide a significantly better outcome estimation. The apparent difference in
the canonical and principal functions implies that FGCCA has indeed retrieved
components highlighting the association with the outcome for each biomarker.
As the canonical functions have a monotonic trend, we can argue that survival
is mainly associated with a decreasing or increasing behavior of the various
biomarkers. More precisely, mortality seems to be associated with an increase
in bilirubin, prothrombin and a decrease in albumin. It is consistent with the
fact that a decreasing albumin or increasing bilirubin and prothrombin are
usually associated with a bad prognosis of the liver.

Figure 6: (left) First principal/canonical function retrieved with FPCA/FGCCA.
For FGCCA a multiblock-FPLS design is used, integrating only the interaction
between each biomarker and the response. Functions are flipped to ensure that
the first component is positively correlated to the outcome, to improve
interpretability. (right) Boxplot of the balanced accuracy computed on the
test set for 100 monte-carlo runs. p-value and significance level of the
difference between the two distributions (t-test) are given.

####  3.2.3 Reconstruction

Finally, we propose to evaluate the ability of FGCCA to reconstruct
trajectories and predict biomarker values at unobserved times. To this end, we
propose dividing the data into training and test datasets. The training
dataset is used to estimate (cross-)covariance operators and to run FGCCA. The
test dataset contains trajectories on which we seek to predict the last
observation, which has been removed. Each subject’s prediction is made from
reconstructed trajectories computed from previously obtained canonical
functions (from the training dataset) and components estimated using data
before the last observation. We propose to compare the results to an FPCA-
based approach, where the components and functions are estimated from an FPCA.
Furthermore, to evaluate the robustness of the two approaches, we sparsify the
test trajectories at various levels. The results are reported in Figure 7.

We observe a significant advantage of the FGCCA-based approach over the FPCA-
based approach except for bilirubin in the (1M) scenario, asserting that FGCCA
has integrated additional or more stable knowledge. We note that the gap
widens as the sparsity increases. These results pave the way to joint
modeling, as the components could be used both to predict the future
trajectories of biomarkers, as it is done here, and the survival of subjects.
This promising application is currently investigated.

Figure 7: (left) Reconstruction obtained with FGCCA and FPCA in 3 scenarios :
(1M) last observation removed, (2M) two last observations removed (3M) three
last observations removed. Crosses correspond to observations used to estimate
the components, circles correspond to observations we aim to predict. (right)
Boxplots of last observation prediction mean squared error obtained on 100
runs. Statistical significance displayed : (*) $p<0.05$ (****) $p<0.0001$

##  4 Discussion

We introduced Functional Generalized Canonical Correlation Analysis (FGCCA), a
flexible framework for exploring associations among multiple longitudinal
variables, by finding the main joint modes of variation. The method relies on
a monotone and globally convergent algorithm, which only requires
(cross-)covariance operators. We proposed a Bayesian approach for estimating
the components. Consequently, the method is robust to irregular and sparse
data, making it applicable to numerous settings. In addition, we allow
integrating a multivariate response in the analysis by slightly modifying the
optimization problem, paving the way to mixed data uses. Simulation studies
assess the validity of our approach and its underlying design. A wide variety
of usages are presented in the application.

As previously mentioned, the method relies significantly on the estimations of
(cross-)covariance operators. Therefore, studying more in-depth those
estimation procedures could be interesting as new methods have been proposed
recently (Xiao et al. (2018)). Computing confidence bands for the estimated
scores as it is done in sparse and irregular FPCA (Yao et al. (2005)) is
investigated. However, difficulties arise since the FGCCA algorithm does not
have a closed-form solution. Finally, an implementation allowing the user to
change the regularization parameter was developed. In this context, analyzing
the impact of the parameter on the algorithm and the results obtained could be
further investigated.

In numerous studies, longitudinal variables can be grouped in blocks
representing different modalities. For instance, in imaging genetics, multiple
longitudinal variables representing the evolution of several neuroimaging
features can be observed along multiple genetic features. In this context,
considering a block for each variable, as done with FGCCA, may be inefficient
as it would require a complex design and intensive computational resources.
Another approach, which is currently being investigated, would be to integrate
the multiple longitudinal variables in blocks. This approach was used by Happ
and Greven (2018) and can be referred to as multivariate functional data
modeling. Inspired by the multi-way/tensor literature, a reduced rank model
could also significantly help reduce the problem’s complexity. In this
context, numerous works have been proposed, notably for tensor regression
(Zhou et al. (2008)) and, as evoked in the introduction, for RGCCA (Girka et
al. (2023)).

## References

  * Carroll [1968] J. D. Carroll.  Generalization of canonical correlation analysis to three of more sets of variables.  1968\.
  * de Leeuw [1994] J. de Leeuw.  Block-relaxation algorithms in statistics.  In _Studies in Classification, Data Analysis, and Knowledge Organization_ , pages 308–324. Springer Berlin Heidelberg, 1994.
  * German et al. [2021] C. A. German, J. S. Sinsheimer, J. Zhou, and H. Zhou.  WiSER: Robust and scalable estimation and inference of within-subject variances from intensive longitudinal data.  _Biometrics_ , 78(4):1313–1327, Aug. 2021.
  * Girka et al. [2023] F. Girka, A. Gloaguen, L. L. Brusquet, V. Zujovic, and A. Tenenhaus.  Tensor generalized canonical correlation analysis, 2023.  URL https://arxiv.org/abs/2302.05277.
  * Górecki et al. [2020] T. Górecki, M. Krzyśko, and W. Wołyński.  Generalized canonical correlation analysis for functional data.  _Biometrical Letters_ , 57(1):1–12, 2020.
  * Happ and Greven [2018] C. Happ and S. Greven.  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  _Journal of the American Statistical Association_ , 113(522):649–659, Feb. 2018.
  * Happ-Kurz [2020] C. Happ-Kurz.  Object-oriented software for functional data.  _Journal of Statistical Software_ , 93(5), 2020.
  * He et al. [2003] G. He, H. G. Müller, and J. L. Wang.  Functional canonical analysis for square integrable stochastic processes.  _Journal of Multivariate Analysis_ , 85(1):54–77, apr 2003.  ISSN 0047-259X.
  * Horst [1961] P. Horst.  Relations amongm sets of measures.  _Psychometrika_ , 26(2):129–149, June 1961.
  * Hotelling [1936] H. Hotelling.  Relations Between Two Sets of Variates.  _Biometrika_ , 28(3/4):321, dec 1936.  ISSN 00063444.
  * Hwang et al. [2011] H. Hwang, K. Jung, Y. Takane, and T. S. Woodward.  Functional Multiple-Set Canonical Correlation Analysis.  _Psychometrika 2011 77:1_ , 77(1):48–64, oct 2011\.  ISSN 1860-0980.
  * Kettenring [1971] J. R. Kettenring.  Canonical analysis of several sets of variables.  _Biometrika_ , 58(3):433–451, 1971.
  * Leurgans et al. [1993] S. E. Leurgans, R. A. Moyeed, and B. W. Silverman.  Canonical Correlation Analysis When the Data are Curves.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 55(3):725–740, jul 1993.  ISSN 2517-6161.
  * Murtaugh et al. [1994] P. A. Murtaugh, E. R. Dickson, G. M. Van Dam, M. Malinchoc, P. M. Grambsch, A. L. Langworthy, and C. H. Gips.  Primary biliary cirrhosis: Prediction of short‐term survival based on repeated patient visits.  _Hepatology_ , 20(1):126–134, 1994.  ISSN 15273350.
  * Preda et al. [2007] C. Preda, G. Saporta, and C. Lévéder.  PLS classification of functional data.  _Computational Statistics_ , 22(2):223–235, Feb. 2007.
  * Rasmussen and Williams [2006] C. E. Rasmussen and C. K. I. Williams.  Gaussian processes for machine learning.  2006\.
  * Rice and Silverman [1991] J. A. Rice and B. W. Silverman.  Estimating the Mean and Covariance Structure Nonparametrically When the Data are Curves.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 53(1):233–243, sep 1991.
  * Rizopoulos [2011] D. Rizopoulos.  Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data.  _Biometrics_ , 67(3):819–829, Feb. 2011.
  * Shin and Lee [2015] H. Shin and S. Lee.  Canonical correlation analysis for irregularly and sparsely observed functional data.  _Journal of Multivariate Analysis_ , 134:1–18, Feb. 2015\.
  * Singh et al. [2019] A. Singh, C. P. Shannon, B. Gautier, F. Rohart, M. Vacher, S. J. Tebbutt, and K.-A. L. Cao.  DIABLO: an integrative approach for identifying key molecular drivers from multi-omics assays.  _Bioinformatics_ , 35(17):3055–3062, Jan. 2019\.
  * Tenenhaus and Tenenhaus [2011] A. Tenenhaus and M. Tenenhaus.  Regularized Generalized Canonical Correlation Analysis.  _Psychometrika_ , 76(2):257–284, apr 2011.  ISSN 00333123.
  * Tenenhaus et al. [2014] A. Tenenhaus, C. Philippe, V. Guillemot, K.-A. L. Cao, J. Grill, and V. Frouin.  Variable selection for generalized canonical correlation analysis.  _Biostatistics_ , 15(3):569–583, Feb. 2014.
  * Tenenhaus et al. [2015] A. Tenenhaus, C. Philippe, and V. Frouin.  Kernel generalized canonical correlation analysis.  _Computational Statistics & Data Analysis_, 90:114–131, Oct. 2015.
  * Tenenhaus et al. [2017] M. Tenenhaus, A. Tenenhaus, and P. J. Groenen.  Regularized Generalized Canonical Correlation Analysis: A Framework for Sequential Multiblock Component Methods.  _Psychometrika_ , 82(3):737–777, sep 2017.  ISSN 1860-0980.
  * Wold et al. [2001] S. Wold, M. Sjöström, and L. Eriksson.  PLS-regression: a basic tool of chemometrics.  _Chemometrics and Intelligent Laboratory Systems_ , 58(2):109–130, Oct. 2001.  doi: 10.1016/s0169-7439(01)00155-1.
  * Xiao et al. [2018] L. Xiao, C. Li, W. Checkley, and C. Crainiceanu.  Fast covariance estimation for sparse functional data.  _Statistics and Computing_ , 28(3):511–522, may 2018.  ISSN 15731375.
  * Yang et al. [2011] W. Yang, H. G. Müller, and U. Stadtmüller.  Functional singular component analysis.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 73(3):303–324, jun 2011.  ISSN 1467-9868.
  * Yao et al. [2005] F. Yao, H. G. Müller, and J. L. Wang.  Functional data analysis for sparse longitudinal data.  _Journal of the American Statistical Association_ , 100(470):577–590, jun 2005.  ISSN 01621459.
  * Zhou et al. [2008] L. Zhou, J. Z. Huang, and R. J. Carroll.  Joint modelling of paired sparse functional data using principal components.  _Biometrika_ , 95(3):601–619, sep 2008.  ISSN 0006-3444.

##  Appendix A Proofs

###  A.1 Definition of the functional gradient

Let $\Psi$ be the objective function:

$\Psi(f_{1},\dots,f_{J})=\Psi(\boldsymbol{f})=\sum_{j\neq
j^{\prime}}c_{j,j^{\prime}}g(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})$

Let $j\in\\{1,\ldots,J\\}$, $f_{j^{\prime}}\,\forall j^{\prime}\neq j$. Let

$\begin{array}[]{lcl}\phi_{j}:L^{2}(I_{j})&\to&\mathbb{R}\\\
f_{j}&\mapsto&\Psi(f_{1},\ldots,f_{j},\ldots,f_{J})\\\ \end{array}$

Proposition. $\phi_{j}$ is Gâteaux-differentiable and its gradient is:

$\nabla_{j}\phi_{j}(f_{j})=2\sum_{j^{\prime}\neq
j}c_{j,j^{\prime}}g^{\prime}(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle)\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}$

Proof. Let $f_{j}\in L^{2}(I_{j})$. Let’s show that it exists a continuous
(bounded) linear operator $\phi_{j}^{\prime}(f_{j})$ such that

$\forall h\in
L^{2}(I_{j}),\,\phi_{j}^{\prime}(f_{j})h=\lim_{\begin{subarray}{c}\alpha\to
0\\\ \alpha>0\end{subarray}}\dfrac{\phi_{j}(f_{j}+\alpha
h)-\phi_{j}(f_{j})}{\alpha}$ (14)

First notice that according to Fubini theorem, expression of $\phi_{j}$ may be
simplified:

$\begin{array}[]{lcl}\phi_{j}(f_{j})&=&\displaystyle\sum_{j^{\prime}\neq
j}c_{j,j^{\prime}}g(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})+\sum_{j^{\prime}\neq
j}c_{j^{\prime},j}g(\langle
f_{j}^{\prime},\boldsymbol{\Sigma}_{j^{\prime}j}f_{j}\rangle_{L^{2}})\\\
&=&\displaystyle 2\sum_{j^{\prime}\neq j}c_{j,j^{\prime}}g(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})\\\
\end{array}$

Then, let’s show that the Gâteaux-differential of $\phi_{j}$ at $f_{j}$ in the
direction $h$ is defined by studying the limit in (14):

$\begin{array}[]{lcl}\dfrac{\phi_{j}(f_{j}+\alpha
h)-\phi_{j}(f_{j})}{\alpha}&=&2\displaystyle\sum_{j^{\prime}\neq
j}c_{j,j^{\prime}}\dfrac{g(\langle f_{j}+\alpha
h,\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})-g(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})}{\alpha}\\\
&=&2\displaystyle\sum_{j^{\prime}\neq j}c_{j,j^{\prime}}g^{\prime}(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})\langle
h,\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}}+o(\alpha)\\\
\end{array}$

Thus, limit in (14) exists and leads to a linear operator:

$\phi_{j}^{\prime}(f_{j})h=\left\langle h,2\displaystyle\sum_{j^{\prime}\neq
j}c_{j,j^{\prime}}g^{\prime}(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}})\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\right\rangle_{L^{2}}$

Applying Cauchy-Schwarz and the monotonicity of $g^{\prime}$, it comes that
$\phi_{j}^{\prime}(f_{j})$ is bounded (and thus continuous) and so that
$\phi_{j}(f_{j})$ is Gâteaux-differentiable at $f_{j}$. By definition of the
gradient:

$\nabla_{j}\phi_{j}(f_{j})=2\sum_{j^{\prime}\neq
j}c_{j,j^{\prime}}g^{\prime}\left(\langle
f_{j},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}\rangle_{L^{2}}\right)\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}$

###  A.2 Monotone convergence of the solving procedure (Proposition 1)

The proof of Proposition 1 is based on the proof given in Tenenhaus et al.
[2017]. In the following, we use the notations of the paper. We consider the
following lemma:

######  Lemma 1.

Consider the set $\Omega=\Omega_{1}\times\dots\times\Omega_{J}$, the function
$\Psi:\Omega\rightarrow\mathbb{R}$ and the operator
$c:\Omega\rightarrow\Omega$ defined in the paper. The following properties
hold:

  1. (a)

$\Omega$ is a bounded set

  2. (b)

$c$ is a continuous operator

  3. (c)

$\Psi(\boldsymbol{f})\leq\Psi(c(\boldsymbol{f}))$ for any
$\boldsymbol{f}\in\Omega$

Proof of lemma 1 :

(a) $\Omega$ is bounded as it is the Cartesian product of $J$ bounded sets

(b) As $\Psi$ is a continuous differentiable operator, $r_{j}$ is continuous.
Thus, $c_{j}$ is continuous and $c=c_{1}\circ\dots c_{J}$ being the
composition of $J$ continuous operators is also continuous

(c) Let
$\operatorname*{\boldsymbol{f}}=({\operatorname*{\boldsymbol{f}}}_{1},\dots,{\operatorname*{\boldsymbol{f}}}_{J})\in\Omega$.
We want to find an update
$\hat{\operatorname*{\boldsymbol{f}}}_{j}\in\Omega_{j}$ of
${\operatorname*{\boldsymbol{f}}}_{j}$ such that
$\Psi(\operatorname*{\boldsymbol{f}})\leq\Psi({\operatorname*{\boldsymbol{f}}}_{1},\dots,{\operatorname*{\boldsymbol{f}}}_{j-1},\hat{\operatorname*{\boldsymbol{f}}}_{j},{\operatorname*{\boldsymbol{f}}}_{j+1},\dots,{\operatorname*{\boldsymbol{f}}}_{J})$.
For that purpose we use the fact that a convex function lies above its linear
approximation. Thus at ${\operatorname*{\boldsymbol{f}}}_{j}$ and for any
$\tilde{\operatorname*{\boldsymbol{f}}}_{j}\in\Omega_{j}$ we have

$\Psi({\operatorname*{\boldsymbol{f}}}_{1},\dots,{\operatorname*{\boldsymbol{f}}}_{j-1},\tilde{\operatorname*{\boldsymbol{f}}}_{j},{\operatorname*{\boldsymbol{f}}}_{j+1},\dots,{\operatorname*{\boldsymbol{f}}}_{J})\geq\Psi(\operatorname*{\boldsymbol{f}})+\langle\nabla_{j}\Psi(\operatorname*{\boldsymbol{f}}),\tilde{\operatorname*{\boldsymbol{f}}}_{j}-{\operatorname*{\boldsymbol{f}}}_{j}\rangle=l_{j}(\tilde{\operatorname*{\boldsymbol{f}}}_{j},\operatorname*{\boldsymbol{f}})$

Using the Cauchy-Schwartz inequality we obtain the unique maximizer
$\hat{\operatorname*{\boldsymbol{f}}}_{j}\in\Omega_{j}$ of
$l_{j}(\tilde{\operatorname*{\boldsymbol{f}}}_{j},\operatorname*{\boldsymbol{f}})$
w.r.t. $\tilde{\operatorname*{\boldsymbol{f}}}_{j}\in\Omega_{j}$

$\hat{\operatorname*{\boldsymbol{f}}}_{j}=\operatorname*{argmax}_{\tilde{\operatorname*{\boldsymbol{f}}}_{j}\in\Omega_{j}}l_{j}(\tilde{\operatorname*{\boldsymbol{f}}}_{j},\operatorname*{\boldsymbol{f}})=\frac{\nabla_{j}\Psi(\operatorname*{\boldsymbol{f}})}{||\nabla_{j}\Psi(\operatorname*{\boldsymbol{f}})||}=r_{j}(\operatorname*{\boldsymbol{f}})$

Thus we have,

$\Psi(\operatorname*{\boldsymbol{f}})=l_{j}({\operatorname*{\boldsymbol{f}}}_{j},\operatorname*{\boldsymbol{f}})\leq
l_{j}(r_{j}(\operatorname*{\boldsymbol{f}}),\operatorname*{\boldsymbol{f}})\leq\Psi({\operatorname*{\boldsymbol{f}}}_{1},\dots,{\operatorname*{\boldsymbol{f}}}_{j-1},r_{j}(\operatorname*{\boldsymbol{f}}),{\operatorname*{\boldsymbol{f}}}_{j+1},\dots,{\operatorname*{\boldsymbol{f}}}_{J})=\Psi(c_{j}(\operatorname*{\boldsymbol{f}}))$
(15)

And also,

$\Psi(c_{j-1}\circ\dots\circ
c_{1}(\operatorname*{\boldsymbol{f}}))\leq\Psi(c_{j}\circ\dots\circ
c_{1}(\operatorname*{\boldsymbol{f}}))$

Leading to the desired inequality for any
$\operatorname*{\boldsymbol{f}}\in\Omega$:

$\Psi(\operatorname*{\boldsymbol{f}})\leq\Psi(c_{1}(\operatorname*{\boldsymbol{f}}))\leq\Psi(c_{2}\circ
c_{1}(\operatorname*{\boldsymbol{f}}))\leq\dots\leq\Psi(c_{J}\circ\dots\circ
c_{1}(\operatorname*{\boldsymbol{f}}))=\Psi(c(\operatorname*{\boldsymbol{f}}))$

Now, using Lemma 1, we can prove Proposition 1. Indeed, Point (c) of Lemma 1
implies that the sequence $\\{\Psi(\operatorname*{\boldsymbol{f}}^{s})\\}$ is
monotonically increasing and, since $\Psi$ is bounded on $\Omega$ (Cauchy-
Schwartz), it is convergent.

###  A.3 Deflation equations (Proposition 2 and 4)

To retrieve orthogonal functions, we project each process
$\operatorname{X}_{j}$ onto the orthogonal of the space spanned by its
associated normalized canonical function $f_{j}$. The deflation equation
associated to this transformation, giving the transformed process
$\operatorname{X}_{j}^{\prime}$, is:

$\operatorname{X}_{j}^{\prime}=\operatorname{X}_{j}-\langle
f_{j},\operatorname{X}_{j}\rangle f_{j}$

Denoting $\boldsymbol{F}_{j}=f_{j}\otimes f_{j}$ where $\otimes$ denotes the
functional tensor product, defined as $(f\otimes g)(h)=\langle g,h\rangle f$,
we can rewrite the previous expression:

$\operatorname{X}_{j}^{\prime}=(\boldsymbol{I}_{I_{j}}-\boldsymbol{F}_{j})(\operatorname{X}_{j})$

Now, we consider the deflated cross-covariance between deflated processes
$\operatorname{X}_{j}^{\prime}$ and $\operatorname{X}_{j^{\prime}}^{\prime}$:

$\displaystyle(\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}f)(s)$
$\displaystyle=\int_{I_{j^{\prime}}}\Sigma_{jj^{\prime}}^{\prime}(s,t)f(t)$
$\displaystyle=\int_{I_{j^{\prime}}}\mathbb{E}[\operatorname{X}_{j}^{\prime}(s)\operatorname{X}_{j^{\prime}}^{\prime}(t)]f(t)$
$\displaystyle=\mathbb{E}[\int_{I_{j^{\prime}}}\operatorname{X}_{j}^{\prime}(s)\operatorname{X}_{j^{\prime}}^{\prime}(t)f(t)]$
$\displaystyle=\mathbb{E}[\langle\operatorname{X}_{j^{\prime}}^{\prime},f\rangle\operatorname{X}_{j}^{\prime}(s)]$

Allowing us to write :

$\displaystyle\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}f$
$\displaystyle=\mathbb{E}[\langle\operatorname{X}_{j^{\prime}}^{\prime},f\rangle\operatorname{X}_{j}^{\prime}]$
$\displaystyle=\mathbb{E}[\langle(\boldsymbol{I}_{I_{j^{\prime}}}-\boldsymbol{F}_{j^{\prime}})(\operatorname{X}_{j^{\prime}}),f\rangle(\boldsymbol{I}_{I_{j}}-\boldsymbol{F}_{j})(\operatorname{X}_{j})]$
$\displaystyle=(\boldsymbol{I}_{I_{j}}-\boldsymbol{F}_{j})\mathbb{E}[\langle\operatorname{X}_{j^{\prime}},(\boldsymbol{I}_{I_{j^{\prime}}}-\boldsymbol{F}_{j^{\prime}})(f)\rangle\operatorname{X}_{j}]$
$\displaystyle=(\boldsymbol{I}_{I_{j}}-\boldsymbol{F}_{j})\boldsymbol{\Sigma}_{jj^{\prime}}(\boldsymbol{I}_{I_{j^{\prime}}}-\boldsymbol{F}_{j^{\prime}})(f)$

Leading to desired expression of the cross covariance operators:

$\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}=(\boldsymbol{I}_{I_{j}}-\boldsymbol{F}_{j})\boldsymbol{\Sigma}_{jj^{\prime}}(\boldsymbol{I}_{I_{j^{\prime}}}-\boldsymbol{F}_{j^{\prime}})$

For retrieving uncorrelated components, $\operatorname{X}_{j}^{\prime}$ is
defined by regressing out from the process the previously retrieved components
$\operatorname{y}_{j}$. The deflation process can be written as:

$\operatorname{X}_{j}^{\prime}=(\boldsymbol{I}_{I_{j}}-d_{j}\boldsymbol{\Sigma}_{jj}\boldsymbol{F}_{j})(\operatorname{X}_{j})$
(16)

where $d_{j}=\langle f_{j},\boldsymbol{\Sigma}_{jj}f_{j}\rangle^{-1}$.
Following the same steps as before we obtain

$\boldsymbol{\Sigma}_{jj^{\prime}}^{\prime}=(\boldsymbol{I}_{I_{j}}-d_{j}\boldsymbol{\Sigma}_{jj}\boldsymbol{F}_{j})\boldsymbol{\Sigma}_{jj^{\prime}}(\boldsymbol{I}_{I_{j^{\prime}}}-d_{j^{\prime}}\boldsymbol{F}_{j^{\prime}}\boldsymbol{\Sigma}_{j^{\prime}j^{\prime}})$

###  A.4 Components conditional expectation (Proposition 3)

In the following we denote
$\boldsymbol{\xi}_{ij}=(\xi_{ij}^{1},\dots,\xi_{ij}^{M})^{\top}$,
$\boldsymbol{\epsilon}_{ij}=(\epsilon_{ij1},\dots,\epsilon_{ijn_{ij}})^{\top}$,
$\boldsymbol{\xi}_{i}=(\boldsymbol{\xi}_{i1}^{\top},\dots,\boldsymbol{\xi}_{iJ}^{\top})^{\top}$
and
$\boldsymbol{\epsilon}_{i}=(\boldsymbol{\epsilon}_{i1}^{\top},\dots,\boldsymbol{\epsilon}_{iJ}^{\top})^{\top}$.
We propose to write the vector of observations for the process $j$ and subject
$i$ with matrix notations:

$\boldsymbol{U}_{ij}=\boldsymbol{\mu}_{ij}+\mathbb{F}_{ij}\boldsymbol{\xi}_{ij}+\boldsymbol{\epsilon}_{ij}$

With $\mathbb{F}_{ij}$ corresponding to a matrix whose columns are the
canonical functions at the observed time points. Denoting
$\mathbb{F}_{i}=\operatorname*{diag}(\mathbb{F}_{i1},\dots,\mathbb{F}_{iJ})$,
we also suggest writing the vector of all the observations for the subject $i$
using matrix notations:

$\boldsymbol{U}_{i}=\boldsymbol{\mu}_{i}+\mathbb{F}_{ij}\boldsymbol{\xi}_{i}+\boldsymbol{\epsilon}_{i}$

From this, we propose to rewrite jointly the observations and the scores, as:

$\begin{bmatrix}\boldsymbol{U}_{i}\\\ \boldsymbol{\xi}_{i}\\\
\end{bmatrix}=\begin{bmatrix}\boldsymbol{\mu}_{i}\\\ 0\\\
\end{bmatrix}+\begin{bmatrix}\mathbb{F}_{i}&\mathbb{I}\\\
\mathbb{I}&0\end{bmatrix}\begin{bmatrix}\boldsymbol{\xi}_{i}\\\
\boldsymbol{\epsilon}_{i}\end{bmatrix}$

Using this expression, and the jointly Gaussian assumptions on $\xi_{i}$ and
$\epsilon_{i}$ we clearly see that $\xi_{i}$ and $\boldsymbol{U}_{i}$ are also
jointly Gaussian with joint law:

$\displaystyle\begin{bmatrix}\boldsymbol{U}_{i}\\\ \boldsymbol{\xi}_{i}\\\
\end{bmatrix}$ $\displaystyle\sim
N\left(\begin{bmatrix}\boldsymbol{\mu}_{i}\\\ 0\\\
\end{bmatrix},\begin{bmatrix}\mathbb{F}_{i}&\mathbb{I}\\\
\mathbb{I}&0\end{bmatrix}\begin{bmatrix}\mathbb{\Sigma}&0\\\
0&\sigma^{2}\mathbb{I}\end{bmatrix}\begin{bmatrix}\mathbb{F}_{i}^{\top}&\mathbb{I}\\\
\mathbb{I}&0\end{bmatrix}\right)$ $\displaystyle\sim
N\left(\begin{bmatrix}\boldsymbol{\mu}_{i}\\\ 0\\\
\end{bmatrix},\begin{bmatrix}\mathbb{F}_{i}\mathbb{\Sigma}\mathbb{F}_{i}^{\top}+\sigma^{2}\mathbb{I}&\mathbb{F}_{i}\mathbb{\Sigma}\\\
\mathbb{\Sigma}\mathbb{F}_{i}^{\top}&\mathbb{\Sigma}\end{bmatrix}\right)$

To obtain the best predictions for the scores, we now consider the conditional
distribution of $\xi_{i}$ and suggest using the conditional expectation as the
best predictor. Using the standard formulation of the Gaussian conditional
distribution (Rasmussen and Williams [2006]) we obtain :

$\mathbb{E}(\boldsymbol{\xi}_{i}|\boldsymbol{U}_{i})=\mathbb{\Sigma}\mathbb{F}_{i}^{\top}(\mathbb{F}_{i}\mathbb{\Sigma}\mathbb{F}_{i}^{\top}+\sigma^{2}\mathbb{I})^{-1}(\boldsymbol{U}_{i}-\boldsymbol{\mu}_{i})$

###  A.5 Properties of deflated processes and components

####  A.5.1 Orthogonal deflation: correspondence between $\xi_{j}^{m}$ and
$\operatorname{y}_{j}^{m}$

In the decomposition model (8), basis coefficients are defined as
$\langle\operatorname{X}_{j},f_{j}^{m}\rangle=\xi_{j}^{m}$, since $f_{j}^{m}$
are orthonormal functions. As in RGCCA, the components are defined in the
FGCCA framework as
$\operatorname{y}_{j}^{m}=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle$
where $\operatorname{X}_{j}^{m}$ stands for the $m$th deflated processes $j$,
with $\operatorname{X}_{j}^{0}=\operatorname{X}_{j}$, and $f_{j}^{0}$ and
$\operatorname{y}_{j}^{0}$ being respectively the first canonical function and
the first component retrieved. We can easily demonstrate recursively that we
have for $1\leq m\leq M$:

$\operatorname{X}_{j}^{m}=\operatorname{X}_{j}-\sum_{k=0}^{m-1}\langle\operatorname{X}_{j},f_{j}^{k}\rangle
f_{j}^{k}$

from this we can write

$\displaystyle\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle$
$\displaystyle=\langle\operatorname{X}_{j},f_{j}^{m}\rangle-\sum_{k=0}^{m-1}\langle\operatorname{X}_{j},f_{j}^{k}\rangle\langle
f_{j}^{k},f_{j}^{m}\rangle$
$\displaystyle=\langle\operatorname{X}_{j},f_{j}^{m}\rangle$

since, again, the set of canonical functions $\\{f_{j}^{m}\\}_{0\leq m\leq M}$
is a set of orthonormal functions. Therefore, we have proven that
$\operatorname{y}_{j}^{m}=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle=\langle\operatorname{X}_{j},f_{j}^{m}\rangle=\xi_{j}^{m}$

####  A.5.2 Uncorrelated deflation: orthogonal property of $f_{j}^{m}$

First, we demonstrate that the functions retrieved using an uncorrelated
components deflation strategy are orthogonal. The $(m+1)$th deflation equation
for the process $j$ is:

$\operatorname{X}_{j}^{m+1}=\operatorname{X}_{j}^{m}-d_{j}^{m}\boldsymbol{\Sigma}_{jj}^{m}\boldsymbol{F}_{j}^{m}\operatorname{X}_{j}^{m}$
(17)

where $d_{j}^{m}=\langle
f_{j}^{m},\boldsymbol{\Sigma}_{jj}^{m}f_{j}^{m}\rangle^{-1}$. We recall that
the gradient operator, giving at each iteration the updated block weight
function is:

$\langle
f_{j}^{m},\hat{f}_{j}^{m+1}\rangle=\nabla_{j}\Psi(\boldsymbol{f})=2\sum_{\begin{subarray}{c}j^{\prime}=1\\\
j^{\prime}\neq j\end{subarray}}^{J}c_{j,j^{\prime}}g^{\prime}(\langle
f_{j}^{m+1},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}^{m+1}\rangle)\langle
fm_{j}^{m},\boldsymbol{\Sigma}_{jj^{\prime}}^{m+1}f_{j^{\prime}}^{m+1}\rangle$
(18)

Focusing on the elements in the sum we have that

$\displaystyle\langle
f_{j}^{m},\boldsymbol{\Sigma}_{jj^{\prime}}^{m+1}f_{j^{\prime}}^{m+1}\rangle$
$\displaystyle=\langle\mathbb{E}[\langle\operatorname{X}_{j^{\prime}}^{m+1},f_{j^{\prime}}^{m+1}\rangle\operatorname{X}_{j}^{m+1}],f_{j}^{m}\rangle$
$\displaystyle=\mathbb{E}[\langle\operatorname{X}_{j^{\prime}}^{m+1},f_{j^{\prime}}^{m+1}\rangle\langle\operatorname{X}_{j}^{m+1},f_{j}^{m}\rangle]$

Additionally,

$\displaystyle\langle\operatorname{X}_{j}^{m+1},f_{j}^{m}\rangle$
$\displaystyle=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle-
d_{j}^{m}\langle(\boldsymbol{\Sigma}_{jj}^{m}\mathbb{F}_{j}^{m})(\operatorname{X}_{j}^{m}),f_{j}^{m}\rangle$
$\displaystyle=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle-
d_{j}^{m}\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle\langle\boldsymbol{\Sigma}_{jj}^{m}f_{j}^{m},f_{j}^{m}\rangle$
$\displaystyle=\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle-
d_{j}^{m}\langle\operatorname{X}_{j}^{m},f_{j}^{m}\rangle(d_{j}^{m})^{-1}$
$\displaystyle=0$

Thus we have $\langle
f_{j}^{m},\boldsymbol{\Sigma}_{jj^{\prime}}f_{j^{\prime}}^{m+1}\rangle=0$ and
finally $\langle f_{j}^{m},f_{j}^{m+1}\rangle=0$, since all the terms in the
sum of the update function are null.

##  Appendix B Additional experiments

###  B.1 Simulation 2: additional results

In this additional section, we investigate the role of $\sigma$ on the
estimation accuracy. Results are displayed in figure 8 and 9.

###  B.2 Simulation 3: additional results

As previously, we propose additional results obtained with different values of
noise $\sigma$. Results are displayed in figure 10.

Figure 8: Mean squared errors (MSE) of functions $f_{j}^{m}$ (top) and
components $\xi_{j}^{m}$ (bottom) for $m=1,2,3,4,5,6$ obtained from $100$
simulations with $N=100$, $\sigma^{2}=0$ and various sparsity settings.
Comparison between FPCA (Yao et al. [2005]), FSVD (Yang et al. [2011]) and
FGCCA.

Figure 9: Mean squared errors (MSE) of functions $f_{j}^{m}$ (top) and
components $\xi_{j}^{m}$ (bottom) for $m=1,2,3,4,5,6$ obtained from $100$
simulations with $N=100$, $\sigma^{2}=0.1$ and various sparsity settings.
Comparison between FPCA (Yao et al. [2005]), FSVD (Yang et al. [2011]) and
FGCCA.

Figure 10: Mean relative squared errors (MRSE) of reconstructed trajectories,
using estimated canonical functions and components. Comparing FGCCA and MFPCA
with $M=6$, for various number of subjects $N$, various sparsity settings and
various noise levels: (top) $\sigma_{2}=0$ and (bottom) $\sigma^{2}=0.1$.