Clustering multivariate functional data in group-specific functional subspaces

Abstract

With the emergence of numerical sensors in many aspects of everyday life, there is an increasing need in analyzing multivariate functional data. This work focuses on the clustering of such functional data, in order to ease their modeling and understanding. To this end, a novel clustering technique for multivariate functional data is presented. This method is based on a functional latent mixture model which fits the data into group-specific functional subspaces through a multivariate functional principal component analysis. A family of parsimonious models is obtained by constraining model parameters within and between groups. An Expectation Maximization algorithm is proposed for model inference and the choice of hyper-parameters is addressed through model selection. Numerical experiments on simulated datasets highlight the good performance of the proposed methodology compared to existing works. This algorithm is then applied to the analysis of the pollution in French cities for 1 year.

Appendix: Proofs

1.1 Proof of Proposition 1

$$\begin{aligned} l(\theta )=\sum _{i=1}^n\sum _{k=1}^Kz_{ki}log(\pi _k f(x_i,\theta _k)), \end{aligned}$$

where $z_{ki}{=}1$ if $x_i$ belongs to the cluster k and $z_{ki}=0$ otherwise. $f(x_i,\theta _k)$ is a Gaussian density, with parameters $\theta _k=\{\mu _k,\Sigma _k\}$. So the complete log-likelihood is written:

$$\begin{aligned} l(\theta )= & {} \sum _{i=1}^n\sum _{k=1}^Kz_{ki}log[\pi _k\frac{1}{(2\pi )^{R/2}|\Sigma _k|^{1/2}} exp(\frac{-1}{2}(x_i-\mu _k)^t \Sigma _k^{-1}(x_i-\mu _k))] \\= & {} \sum _{i=1}^n\sum _{k=1}^Kz_{ki}[log(\pi _k)- \frac{1}{2}log|\Sigma _k|-\frac{1}{2}(x_i-\mu _k)^t\Sigma _k^{-1}(x_i-\mu _k)-\frac{R}{2}log(2\pi )]. \end{aligned}$$

For the $[a_{kj}b_kQ_kd_k]$ model, we have:

$$\begin{aligned} l(\theta )= & {} \frac{1}{2}\sum _{i=1}^{n}\sum _{k=1}^Kz_{ki}[-2log(\pi _k)+log(\prod _{j=1}^{d_k}a_{kj}\prod _{j=d_k+1}^Rb_k)+(x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)]\\&\quad -\,\frac{nR}{2}log(2\pi ) \\= & {} \frac{1}{2}\sum _{i=1}^{n}\sum _{k=1}^Kz_{ki}[-2log(\pi _k)+\sum _{j=1}^{d_k} log(a_{kj})+\sum _{j=d_k+1}^Rlog(b_k)]\\&\quad +\,(x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)]-\frac{nR}{2}log(2\pi ). \end{aligned}$$

Let $n_k=\sum _{i=1}^nz_{ki}$ be the number of curves within cluster k, the complete log-likelihood is then written:

$$\begin{aligned} l(\theta )= & {} -\frac{1}{2}\sum _{k=1}^Kn_k[-2log(\pi _k)+\sum _{j=1}^{d_k}log(a_{kj})+\sum _{j={d_k+1}}^Rlog(b_k) \\&\quad +\,\frac{1}{n_k}\sum _{i=1}^nz_{kj}(x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)]-\frac{nR}{2}log(2\pi ). \end{aligned}$$

The quantity $(x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)$ is a scalar, so it is equal to it trace:

$$\begin{aligned} \frac{1}{n_k}\sum _{i=1}^nz_{ki}(x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)=\frac{1}{n_k}\sum _{i=1}^nz_{ki}tr((x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)). \end{aligned}$$

Well $tr([(x_i-\mu _k)^tQ_k]\times [\Delta _k^{-1}Q_k^t(x_i-\mu _k)])=tr([\Delta _k^{-1}Q_k^t(x_i-\mu _k)]\times [(x_i-\mu _k)^tQ_k])$, consequently:

$$\begin{aligned} \frac{1}{n_k}\sum _{i=1}^nz_{ki}(x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)= & {} \frac{1}{n_k}\sum _{i=1}^nz_{ki}tr(\Delta _k^{-1}Q_k^t(x_i-\mu _k)(x_i-\mu _k)^tQ_k) \\= & {} tr(\Delta _k^{-1}Q_k^t[\frac{1}{n_k}\sum _{i=1}^nz_{ki}(x_i-\mu _k)^t(x_i-\mu _k)]Q_k)\\= & {} tr(\Delta _k^{-1}Q_k^tC_kQ_k), \end{aligned}$$

where $C_k=\frac{1}{n_k}\sum _{i=1}^{n}z_{ki}(x_i-\mu _k)^t(x_i-\mu _k)$ is the empirical covariance matrix of the kth element of the mixture model. The $\Delta _k$ matrix is diagonal, so we can write:

$$\begin{aligned} \frac{1}{n_k}\sum _{i=1}^nz_{ki}(x_i-\mu _k)^tQ_k\Delta _k^{-1}Q_k^t(x_i-\mu _k)= & {} \sum _{j=1}^{d_k}\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{a_{kj}} \\&\quad +\,\sum _{j={d_k+1}}^R\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{b_k}, \end{aligned}$$

where $q_{kj}$ is jth column of $Q_k$. Finally,

$$\begin{aligned} l(\theta )= & {} -\frac{1}{2}\sum _{k=1}^Kn_k[-2log(\pi _k)+\sum _{j=1}^{d_k}log(a_{kj})+\sum _{j={d_k+1}}^Rlog(b_k)+\sum _{j=1}^{d_k}\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{a_{kj}} \\&\quad +\,\sum _{j={d_k+1}}^R\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{b_k}]+\frac{nR}{2}log(2\pi ). \end{aligned}$$

1.2 Proof of Proposition 2

$$\begin{aligned} H_k(x)= & {} -2log(\pi _k f(x,\theta _k)) \\= & {} -2log(\pi _k) -2log(f(x,\theta _k)) \\= & {} -2log(\pi _k) - 2log(\frac{1}{(2\pi )^{R/2}|\Sigma _k|^{1/2}} exp(\frac{-1}{2}(x-\mu _k)^t \Sigma _k^{-1}(x-\mu _k)) \\= & {} -2log(\pi _k) - 2log(\frac{1}{(2\pi )^{R/2}|\Sigma _k|^{1/2})}) + (x-\mu _k)^t \Sigma _k^{-1}(x-\mu _k) \\= & {} -2log(\pi _k) + Rlog(2\pi ) +log|\Sigma _k|+(x-\mu _k)^t\Sigma _k^{-1}(x-\mu _k). \end{aligned}$$

But, $\Sigma _k=Q_k\Delta _k Q_k^t$ and $Q_k^tQ_k=I_R$, hence:

$$\begin{aligned} H_k(x) = -2log(\pi _k)+Rlog(2\pi )+log|\Sigma _k|+(x_k-\mu _k)^t(Q_k\Delta _k Q_k^t)^{-1}(x_k-\mu _k). \end{aligned}$$

Let $Q_k=\tilde{Q_k}+\bar{Q_k}$ where $\tilde{Q_k}$ is the $R\times R$ matrix containing the $d_k$ first columns of $Q_k$ completed by zeros and where $\bar{Q_k}=Q_k-\tilde{Q_k}$. Notice that $\tilde{Q_k}\Delta _k^{-1}\bar{Q_k^t}=\bar{Q_k}\Delta _k^{-1}\tilde{Q_k^t}=O_p$ where $O_p$ is the null matrix. So,

$$\begin{aligned} Q_k\Delta _k^{-1}Q_k^t= (\tilde{Q_k}+\bar{Q_k})\Delta _k^{-1}(\tilde{Q_k}+\bar{Q_k})= \tilde{Q_k}\Delta _k^{-1}\tilde{Q_k}+\bar{Q_k}\Delta _k^{-1}\bar{Q_k}. \end{aligned}$$

Hence,

$$\begin{aligned} H_k(x)= & {} -2log(\pi _k)+Rlog(2\pi )+log|\Sigma _k|+(x-\mu _k)^t\tilde{Q_k}\Delta _k^{-1}\tilde{Q_k^t}(x-\mu _k) \\&\quad +\,(x-\mu _k)^t\bar{Q_k}\Delta _k^{-1}\bar{Q_k^t}(x-\mu _k). \end{aligned}$$

With definitions $\tilde{Q_k}[\tilde{Q_k^t}\tilde{Q_k}]=\tilde{Q_k}$ and $\bar{Q_k}[\bar{Q_k^t}\bar{Q_k}]=\bar{Q_k}$, we can rephrase $H_k(x)$ as:

$$\begin{aligned} H_k(x)= & {} -2log(\pi _k)+Rlog(2\pi )+log|\Sigma _k|+(x-\mu _k)^t\tilde{Q_k}\tilde{Q_k^t}\tilde{Q_k}\Delta _k^{-1}\tilde{Q_k^t}\tilde{Q_k}\tilde{Q_k^t}(x-\mu _k) \\&\quad +\,(x-\mu _k)^t\bar{Q_k}\bar{Q_k^t}\bar{Q_k}\Delta _k^{-1}\bar{Q_k^t}\bar{Q_k}\bar{Q_k^t}(x-\mu _k) \\= & {} -2log(\pi _k)+Rlog(2\pi )+log|\Sigma _k|+[\tilde{Q_k}\tilde{Q_k^t}(x-\mu _k)]^t\tilde{Q_k}\Delta _k^{-1}\tilde{Q_k}^t[\tilde{Q_k}\tilde{Q_k^t}(x-\mu _k)] \\&\quad + \,[\bar{Q_k}\bar{Q_k^t}(x-\mu _k)]^t\bar{Q_k}\Delta _k^{-1}\bar{Q_k}^t[\bar{Q_k}\bar{Q_k^t}(x-\mu _k)]. \end{aligned}$$

We define $\mathcal {D}_k=\tilde{Q_k}\Delta _k^{-1}\tilde{Q_k^t}$ and the norm $||.||_{\mathcal {D}_k}$ on $\mathbb {E}_k$ such as $||x||_{\mathcal {D}_k}=x^t\mathcal {D}_kx$. So, on one hand:

$$\begin{aligned}{}[\tilde{Q_k}\tilde{Q_k^t}(x-\mu _k)]^t\tilde{Q_k}\Delta _k^{-1}\tilde{Q_k^t}[\tilde{Q_k}\tilde{Q_k^t}(x-\mu _k)]= ||\tilde{Q_k}\tilde{Q_k^t}(x-\mu _k)||_{\mathcal {D}_k}^2. \end{aligned}$$

On the other hand:

$$\begin{aligned}{}[\bar{Q_k}\bar{Q_k^t}(x-\mu _k)]^t\bar{Q_k}\Delta _k^{-1}\bar{Q_k^t}[\bar{Q_k}\bar{Q_k^t}(x-\mu _k)]=\frac{1}{b_k}||\bar{Q_k}\bar{Q_k^t}(x-\mu _k)||^2. \end{aligned}$$

Consequently,

$$\begin{aligned} H_k(x)=-2log(\pi _k)+Rlog(2\pi )+log|\Sigma _k|+||\tilde{Q_k}\tilde{Q_k^t}(x-\mu _k)||_{\mathcal {D}_k}^2+\frac{1}{b_k}||\bar{Q_k}\bar{Q_k^t}(x-\mu _k)||^2. \end{aligned}$$

Knowing $P_k$, $P_k^{\bot }$ and $||\mu _k-P_k^{\bot }||^2=||x-P_k(x)||^2$, we have:

$$\begin{aligned} H_k(x)= ||\mu _k-P_k(x)||^2_{\mathcal {D}_k}+\frac{1}{b_k}||x-P_k(x)||^2+log|\Sigma _k|-2log(\pi _k)+Rlog(2\pi ). \end{aligned}$$

Moreover, $log|\Sigma _k|=\sum _{j=1}^{d_k}log(a_{kj})+(R-d_k)log(b_k)$. Finally,

$$\begin{aligned} H_k(x)= & {} ||\mu _k-P_k(x)||^2_{\mathcal {D}_k}+\frac{1}{b_k}||x-P_k(x)||^2+\sum _{j=1}^{d_k}log(a_{kj})+(R-d_k)log(b_k)-2log(\pi _k) \\&\quad \,+Rlog(2\pi ). \end{aligned}$$

1.3 Proof of Proposition 3

Parameter$Q_{k}$ We have to maximise the log-likelihood under the constraint $q_{kj}^tq_{kj}=1$, which is equivalent to looking for a saddle point of the Lagrange function:

$$\begin{aligned} \mathcal {L}=-2l(\theta )-\sum _{j=1}^R \omega _{kj}(q_{kj}^tq_{kj}-1), \end{aligned}$$

where $\omega _{kj}$ are Lagrange multiplier. So we can write:

$$\begin{aligned} \mathcal {L}&=\sum _{k=1}^K \eta _k\left[ \sum _{j=1}^{d_k}\left( log(a_{kj})+\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{a_{kj}}\right) \right. \\&\left. \quad +\,\sum _{j=d_k+1}^R\left( log(b_k)+\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{b_k}\right) -2log(\pi _k)\right] +\frac{nR}{2}log(2\pi ) \\&\quad -\,\sum _{j=1}^R\omega _{kj}(q_{kj}^tq_{kj}-1). \end{aligned}$$

Therefore, the gradient of $\mathcal {L}$ in relation to $q_{kj}$ is:

$$\begin{aligned} \nabla _{q_{kj}}\mathcal {L}&=\nabla _{q_{kj}}\left( \sum _{k=1}^K\eta _k\left[ \sum _{j=1}^{d_k}\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{a_{kj}}+\sum _{j=d_k+1}^{R}\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{b_k}\right] \right. \\&\left. \quad - \,\sum _{j=1}^R\omega _{kj}(q_{kj}^tq_{kj}-1)\right) . \end{aligned}$$

As a reminder, when W is symmetric, then $\frac{\partial }{\partial x}(x-s)^TW(x-s)=2W(x-s)$ and $\frac{\partial }{\partial x}(x^Tx)=2x$ (cf. Petersen and Pedersen (2012)), so:

$$\begin{aligned} \nabla _{q_{kj}}\mathcal {L}= & {} \eta _k\left[ 2\frac{W^{1/2}C_kW^{1/2}}{\sigma _{kj}}q_{kj}\right] -2\omega _{kj}q_{kj} \end{aligned}$$

where $\sigma _{kj}$ is the jth diagonal term of matrix $\Delta _k$.

So,

$$\begin{aligned} q_{kj}^t\nabla _{q_{kj}}\mathcal {L}=0\Leftrightarrow & {} \omega _{kj}q_{kj}=\frac{\eta _k}{\sigma _{kj}}q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj} \\\Leftrightarrow & {} W^{1/2}C_kW^{1/2}q_{kj}=\frac{\omega _{kj}\sigma _{kj}}{\eta _k}q_{kj}. \end{aligned}$$

$q_{kj}$ is the eigenfunction of $W^{1/2}C_kW^{1/2}$ which match the eigenvalue $\lambda _{kj}=\frac{\omega _{kj}\sigma _{kj}}{\eta _k}=W^{1/2}C_kW^{1/2}$. We can write $q_{kj}^tq_{kl}=0$ if $j\ne l$. So the log-likelihood can be written:

$$\begin{aligned} -2l(\theta )=\sum _{k=1}^K\eta _k\left[ \sum _{j=1}^{d_k}\left( log(a_{kj})+\frac{\lambda _{kj}}{a_{kj}}\right) +\sum _{j=d_k+1}^R\left( log(b_k)+\frac{\lambda _{kj}}{b_k}\right) \right] +C^{te}, \end{aligned}$$

we substitute the equation $\sum _{j=d_k+1}^R\lambda _{kj}=tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}\lambda _{kj}$:

$$\begin{aligned} -2l(\theta )= & {} \sum _{k=1}^K\eta _k\left[ \sum _{j=1}^{d_k}log(a_{kj})+\sum _{j=1}^{d_k}\frac{\lambda _{kj}}{a_{kj}}+\sum _{j=d_i+1}^Rlog(b_k)+\frac{1}{b_k}\left( tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}\lambda _{kj}\right) \right] +C^{te} \\= & {} \sum _{k=1}^K\eta _k\left[ \sum _{j=1}^{d_k}log(a_{kj})+\sum _{j=1}^{d_k}\lambda _{kj}\left( \frac{1}{a_{kj}}-\frac{1}{b_k}\right) +\sum _{j=d_i+1}^Rlog(b_k)+\frac{1}{b_k}tr(W^{1/2}C_kW^{1/2})\right] +C^{te} \\= & {} \sum _{k=1}^K\eta _k\left[ \sum _{j=1}^{d_k}log(a_{kj})+\sum _{j=1}^{d_k}\lambda _{kj}\left( \frac{1}{a_{kj}}-\frac{1}{b_k}\right) +(p-d_k)log(b_k)+\frac{tr(W^{1/2}C_kW^{1/2})}{b_k}\right] +C^{te}. \end{aligned}$$

In order to minimize $-2l(\theta )$ compared to $q_{kj}$, we minimize the quantity $\sum _{k=1}^K\eta _k\sum _{j=1}^{d_k}\lambda _{kj}(\frac{1}{a_{kj}}-\frac{1}{b_k})$ compared to $\lambda _{kj}$. Knowing that $(\frac{1}{a_{kj}}-\frac{1}{b_k})\le 0, \forall j=1,\ldots ,d_k,$$\lambda _{kj}$ has to be as high as feasible. So, the jth column $q_{kj}$ of matrix Q is estimated by the eigenfunction associated to the jth highest eigenvalue of $W^{1/2}C_kW^{1/2}$.

Parameter$a_{kj}$ As a reminder $(ln(x))'=\frac{x'}{x}$ and $(\frac{1}{x})'=-\frac{1}{x^2}$. The partial derivative of $l(\theta )$ in relation to $a_{kj}$ is:

$$\begin{aligned} -2\frac{\partial l(\theta )}{\partial a_{kj}}=\eta _k\left( \frac{1}{a_{kj}}-\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{a_{kj}^2}\right) \end{aligned}$$

The condition $\frac{\partial l(\theta )}{\partial a_{kj}}=0$ is equivalent to:

$$\begin{aligned}&\eta _k\left( \frac{1}{a_{kj}}-\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{a_{kj}^2}\right) =0 \\&\quad \Leftrightarrow \frac{1}{a_{kj}}=\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{a_{kj}^2} \\&\quad \Leftrightarrow a_{kj}=q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj} \\&\quad \Leftrightarrow a_{kj}= \lambda _{kj} \end{aligned}$$

Parameter$b_k$ The partial derivative of $l(\theta )$ in relation to $b_{k}$ is:

$$\begin{aligned} -2\frac{\partial l(\theta )}{\partial b_{k}}= & {} \eta _k\sum _{j=d_k+1}^{R}\left( \frac{1}{b_{k}}-\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{b_{k}^2}\right) \\= & {} \eta _k\left( \frac{R-d_k}{b_k}-\sum _{j={d_k+1}}^R\frac{q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}}{b_k^2}\right) \end{aligned}$$

But,

$$\begin{aligned} \sum _{j=d_k+1}^R q_{j}^tW^{1/2}C_kW^{1/2}q_j=tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}q_j^tW^{1/2}C_kW^{1/2}q_j, \end{aligned}$$

so:

$$\begin{aligned} -2\frac{\partial l(\theta )}{\partial b_{k}}= & {} \eta _k\frac{(R-d_k)}{b_k}-\frac{\eta _k}{b_k^2}\left( tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}q_{kj}^tW^{1/2}C_kW^{1/2}q_{kj}\right) \\= & {} \eta _k\frac{(R-d_k)}{b_k}-\frac{\eta _k}{b_k^2}\left( tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}\lambda _{kj}\right) \end{aligned}$$

The condition $\frac{\partial l(\theta )}{\partial b_{k}}=0$ is equivalent to:

$$\begin{aligned}&\eta _k\frac{(R-d_k)}{b_k}-\frac{\eta _k}{b_k^2}\left( tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}\lambda _{kj}\right) =0 \\&\quad \Leftrightarrow \eta _k\frac{(R-d_k)}{b_k}=\frac{\eta _k}{b_k^2}\left( tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}\lambda _{kj}\right) \\&\quad \Leftrightarrow b_k=\frac{\eta _k}{\eta _k(R-d_k)}\left( tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}\lambda _{kj}\right) \\&\quad \Leftrightarrow b_k=\frac{1}{(R-d_k)}\left( tr(W^{1/2}C_kW^{1/2})-\sum _{j=1}^{d_k}\lambda _{kj}\right) \end{aligned}$$