Modeling longitudinal data using matrix completion - PMC

2.1. Direct approach

If the set of observed values for each individual is dense, elementary interpolation using a continuous basis can be sufficient for approximating the entire trajectory. Let $\left(b_i:i\in N\right)$ be a basis of $L_2\left(\left(t_{min},t_{max}\right)\right)$, i.e. the space of functions $f:\left(t_{min},t_{max}\right)\to R$ such that the integral of $f^2$ on $\left(t_{min},t_{max}\right)$ is finite, such as splines or polynomials. In practice, we truncate the basis to a finite set of the first $K\in N$ basis elements. Let $b\left(t\right)={\left(b_1\left(t\right),b_2\left(t\right),\dots ,b_K\left(t\right)\right)}^{\prime }$ be a vector of $K$ basis elements evaluated at a timepoint $t\in \left(t_{min},t_{max}\right)$. Throughout this article we use the word basis to refer to some truncated basis of $K$ elements. The choice of the basis is relevant and should be done carefully. We refer the reader to the vast literature on functional data analysis for recommendations on this topic (Ramsay, 2006). Throughout this manuscript, we assume the basis is fixed upfront.

To find an individual trajectory, for each subject $i\in \{1,\dots ,N\}$, we might use least squares and estimate a set of coefficients $w_i\in R^K$, minimizing squared Euclidean distance to observed points

This direct approach has two main drawbacks. First, it ignores correlations within the curves, which could potentially improve the fit. Second, if the number of observations $n_i$ for an individual $i$ is smaller or equal to the size of the basis $K$, we can fit a curve with no error leading to overfitting and an unreliable estimator of the variance.

To remedy the first issue, it is common to estimate the covariance operator, compute principal functional components across all individuals and represent curves in the space spanned by the first few of them. Such representation, is referred to as a Karhunen-Loève expansion (Watanabe, 1965; Kosambi, 2016), has became a foundation of many functional data analysis workflows (Ramsay and Dalzell, 1991; Yao et al., 2005; Cnaan et al., 1997; Laird, 1988; Horváth and Kokoszka, 2012; Besse et al., 1997).

A basic idea to remedy the second issue is to estimate both the basis coefficients and the covariance structure simultaneously. Linear mixed-effect models provide a convenient solution to this problem.

2.2. Linear mixed-effect models

A common approach to modeling longitudinal observation $\left(t_i,y_i\right)$ is to assume that data come from a linear mixed-effect model (LMM) (Verbeke, 1997; Zeger et al., 1988). We operate in a functionals pace with a basis $b\left(t\right)$ of $K$ elements and we assume there exists a fixed effect $\mu \left(t\right)=m^{\prime }b\left(t\right)$, where $m={\left(m_1,\dots ,m_K\right)}^{\prime }$ for $m_k\in R$ for $1\le k\le K$. We model the individual random effect as a vector of basis coefficients. In the simplest form, we assume

where $\Sigma$ is a $K\times K$ covariance matrix. We model individual observations as

where ${\mu }_i={\left(\mu \left(t_{i,1}\right),\mu \left(t_{i,2}\right),\dots ,\mu \left(t_{i,n_i}\right)\right)}^{\prime },\sigma$ is the standard deviation of observation error and $B_i={\left(b\left(t_{i,1}\right),\dots ,b\left(t_{i,n_i}\right)\right)}^{\prime }$ is the basis evaluated at timepoints defined in $t_i$. Estimation of the model parameters is typically accomplished by the expectation-maximization (EM) algorithm (Laird and Ware, 1982). For predicting coefficients $w_i$ one can use the best unbiased linear predictor (BLUP) (Henderson, 1950; Robinson, 1991).

Since the LMM estimates the covariance structure and the individual fit simultaneously, it reduces the problem of overfitting of $w_i$, present in the direct approach. However, this model is only applicable if we observe a relatively large number of observations per subject since we attempt to estimate $K$ coefficients for every subject.

To model trajectories from a small number of observations, practitioners further constrain the covariance structure. If we knew the functions which contribute the most to the random effect, we could fit an LMM in a smaller space spanned by these functions. We explore possibilities to learn the basics from the data using low-rank approximations of the covariance matrix.

2.3. Low-rank approximations

There are multiple ways to constrain the covariance structure. We can use cross-validation or information criteria to choose the best basis, the number of elements, or the positions of spline knots (Rice and Wu, 2001; Bigelow and Dunson, 2009). Alternatively, we can place a prior on the covariance matrix (MacLehose and Dunson, 2009).

Another solution is to restrict the latent space to $q<K$ dimensions and learn from the data a mapping $A\in R^{K\times q}$ between the latent space and the basis. In the simplest scenario with Gaussian noise, observations can be modeled as

following the notation from (3).

James et al. (2000) propose an EM algorithm for finding model parameters and predicting latent variables $w_i\in R^q$ in (4). In the expectation stage, they compute the conditional mean of $w_i$ given $y_i$ and the model parameters, while in the maximization stage, with $w_i$ assumed observed, they maximize the likelihood with respect to $\{\mu ,A,\sigma \}$. The likelihood, given $w_i$, takes the form

Another approach to estimating parameters of (4) is to optimize over $w_i$ and marginalize $A$ (Lawrence, 2004). This approach allows modifying the distance measure in the latent space, using the kernel trick (Schulam and Arora, 2016).

Estimation of the covariance structure of processes is central to the estimation of individual trajectories. Descary and Panaretos (2016) propose a method where the estimate of the covariance matrix is obtained through matrix completion.

Methods based on low-rank approximations are widely adopted and applied in practice (Berkey and Kent, 1983; Yan et al., 2017; Hall et al., 2006; Besse and Ramsay, 1986; Yao and Lee, 2006; Greven et al., 2011). However, due to their probabilistic formulation and reliance on the distribution assumptions, these models usually need to be carefully fine-tuned for specific situations and existing implementations tend to be slow. This shortcoming motivates us to develop an elementary optimization framework, using existing, extensively studied, and well–optimized tools for matrix algebra.

3.1. Notation

For each individual $i\in \{1,2,\dots ,N\}$ we observe $n_i$ measurements $\left({\dot{y}}_{i,1},\dots ,{\dot{y}}_{i,n_i}\right)$ at time-points $\left(t_{i,1},t_{i,2},\dots ,t_{i,n_i}\right)$. Unlike in the prior work introduced in Section 2 operating on entire curves, here we discretize time to $T$ time-points $G=\left({\tau }_1,{\tau }_2,\dots ,{\tau }_T\right)$, where $t_{min}={\tau }_1,t_{max}={\tau }_T$ and $T$ is arbitrarily large. Each individual $i$ is expressed as a partially observed vector $r_i\in R^T$. For each time-point $t_{i,j}$ for $1\le j\le n_i$ we find a corresponding grid-point $g_i\left(j\right)=\underset{1\le k\le T}{argmin}\left({\tau }_k-t_{i,j}\right)$. We assign $y_{i,g_i\left(j\right)}={\dot{y}}_{i,j}$. Let $O_i=\left(g_i\left(j\right):1\le j\le n_i\right)$ be a set of grid indices corresponding to observed grid points for an individual $i$. All elements outside $O_i$, i.e. $\left(y_{i,j}:j\notin O_i\right)$ are considered missing.

For $T$ sufficiently large, our observations can be uniquely represented as a $N\times T$ matrix $Y$ with missing values. We denote the set of all observed elements by pairs of indices as $\Omega =\left((i,j):i\in \{1,2,\dots ,N\},j\in O_i\right)$. Let $P_{\Omega }\left(Y\right)$ be the projection onto observed indices, i.e. $P_{\Omega }\left(Y\right)=\tilde{Y}$, such that ${\tilde{Y}}_{i,j}=Y_{i,j}$ for $(i,j)\in \Omega$ and ${\tilde{Y}}_{i,j}=0$ otherwise. We define $P_{\Omega }^{\perp }\left(Y\right)$ as the projection on the complement of $\Omega$, i.e. $P_{\Omega }^{\perp }\left(Y\right)=Y-P_{\Omega }\left(Y\right)$. We use $\parallel \cdot {\parallel }_F$ to denote the Frobenius norm, i.e. the square root of the sum of squares of matrix elements, and $\parallel \cdot {\parallel }_{*}$ to denote the nuclear norm, i.e. the sum of singular values.

As in Section 2 we imposed smoothness by using a continuous basis $b\left(t\right)={\left(b_1\left(t\right),b_2\left(t\right),\dots ,b_K\left(t\right)\right)}^{\prime }$. When we evaluate the basis on a grid $G$ we get a $T\times K$ matrix $B={\left(b\left({\tau }_1\right),b\left({\tau }_2\right),\dots ,b\left({\tau }_T\right)\right)}^{\prime }$.

In our algorithms, we use diagonal thresholding operators defined as follows. Let $D=diag\left(d_1,\dots ,d_p\right)$ be a diagonal matrix. We define soft-thresholding as

where $(x{)}_{+}=max(x,0)$, and hard-thresholding as

where $q=\underset{k}{argmin}\left(d_k<\lambda \right)$.

3.2. Direct approach

The optimization problem (1) of the direct approach described in Section 2.1 can be approximated in the matrix completion setting. Note that the bias introduced by the grid is negligible if the grid is sufficiently large, i.e. if each observed timepoint is relatively close to some grid point. We defined the basis as a set of continuous functions, so the finer the grid the closer we are to the actual value of the basis at a given timepoint. Eventually, the difference becomes negligible. Similarly, for the observed points, the coordinate corresponding to time is more accurate when the grid is finer, but again, due to continuity assumption this difference becomes negligible. The analysis of the properties of the grid is beyond the scope of the paper. As a rule of thumb, we recommend choosing the grid based on a practical indication. For example, if we expect patients to visit a clinic 1–2 times per year, then a grid corresponding to 12 points per year might be sufficient.

We approximate each observation $y_{i,j}$ with a point on the selected grid $y_{i,g\left(j\right)}~y_{i,j}$. Next, we rewrite the optimization problem (1) as a matrix completion problem

where by optimization over $\left(w_i\right)$ we mean optimization over all $\left(w_i:i\in \{1,2,\dots ,N\}\right)$ and $W$ is an $N\times K$ matrix composed of vectors $\left(w_i^{\prime }\right)$ stacked vertically.

The matrix formulation in equation (5) and the classic approach in Section 2.1 share multiple characteristics. In both cases, if data is dense, we may find an accurate representation of the underlying process simply by fitting least-squares estimates of $W$ or $\left(w_i\right)$. Conversely, if the data is too sparse, the problem becomes ill-posed, and the least-squares estimates can overfit.

However, the represenational difference between functional representation (1) and matrix representation (5) leads to different computational approaches. The matrix representation enables us to use the matrix completion framework and, in particular, in Section 3.4 we introduce convex optimization algorithms for solving (5).

Some low-rank constraints on the random effect from the mixed-effect model introduced in Section 2.2 can be expressed in terms of constraints on $W$. In particular, in Section 3.3 we show that the linear mixed-effect model can be expressed by constraining the $rank\left(W\right)$.

3.3. Low-rank matrix approximation

In the low-rank approach described in Section 2.3 we assume that individual trajectories can be represented in a low-dimensional space, by constraining the rank of $W$.

We use a similar approach in solving (5). One difficulty comes from the fact that optimization with a rank constraint on $W$ turns the original least-squares problem into a nonconvex problem. Motivated by the matrix completion literature, we relax the rank constraint in (5) to a nuclear norm constraint and we attempt to solve

for some parameter $\lambda >0$.

From the practical standpoint, low-rank representation corresponds to describing curves on a basis comprised of only a few functions which are the best functions for obtaining low reconstruction errors. In Figure 3 we illustrate two functions for describing processes introduced in Section 1 and Figure 1. Such low-dimensional representation of processes can also be used for clustering based on the type of their progression as we discuss in the data study (Section 6).

Cai et al. (2010) propose a first-order singular value thresholding algorithm (SVT), for solving a general class of problems involving a nuclear norm penalty and a linear form of $Y$, which includes (6). Their algorithm can handle large datasets and has strong guarantees of convergence, but it requires tuning the step size parameter, which can decrease computational efficiency as different steps have to be tried in practice. This limitation was addressed by Ma et al. (2011); Mazumder et al. (2010); Hastie et al. (2015a) who introduced iterative procedures which do not depend on such tuning.

In the last decade, machine learning and statistical learning communities have introduced multiple algorithms for matrix completion. Many of them are suitable for solving (6). However, in this article, we focus on analyzing the benefits of framing the trajectory prediction problem (1) in the matrix completion framework, rather than on benchmarking possible solutions.

3.4. Low-rank approximation with singular value thresholding

The low-rank probabilistic approach, introduced in Section 2.3, is based on the observation that the underlying processes for each subject can be approximated in a low-dimensional space. Here, we exploit the same characteristic using matrix-factorization techniques for solving the optimization problem (6).

For the sake of clarity and simplicity, we choose to solve the problem (6) with an extended version of the SOFT-IMPUTE algorithm designed by Hastie et al. (2015a); Mazumder et al. (2010). We find it easy to implement because it only uses matrix multiplication and Singular Value Decomposition (SVD). Using standard matrix operations also reduces the computational cost of our method, since implementations of these operations leverage well-optimized algorithms. However, as discussed in Section 3.3, many other convex optimization algorithms can be applied.

The key component to the solution is the following property linking problem (6) with the SVD. Consider the fully observed case of (6). Using Theorem 2.1 in Cai et al. (2010), one can show that the optimization problem

with $B^{\prime }B=I$, has a unique solution $W={𝒮}_{\lambda }\left(YB\right)$, where ${𝒮}_{\lambda }\left(X\right)=U{D}_{\lambda }{V}^{\prime }$ and $X=UD{V}^{\prime }$ is the SVD of $X$, i.e., $U$ is an $n\times n$ unitary matrix, $D$ is an $n\times K$ rectangular diagonal matrix, with non-negative real numbers on the diagonal in decreasing order, while $V$ is a $K\times K$ matrix. We refer to ${𝒮}_{\lambda }\left(X\right)$ as the singular value thresholding (SVT) of $X$. Note that, the decomposition of $YB$ is identifiable since the solution of SVD with constrained $D$ is identifiable.

In order to solve (7) with a sparsely observed $Y$, we modify the SOFT-IMPUTE algorithm to account for the basis $B$. Our Algorithm 1 iteratively constructs better approximations of the global solution for each $\lambda$ in some predefined set $\left({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k\right)$ for $k\in N$. For a given approximation of the solution $W^{old}$, we use $W^{old}{B}^{\prime }$ to impute unknown elements of $Y$ obtaining $\tilde{Y}$. Then, we construct the next approximation $W^{new}$ by computing SVT of $\tilde{Y}$.

As a stopping criterion, we compute the change between subsequent solutions, relative to the magnitude of the solution, following the methodology in Cai et al. (2010). We set a fixed threshold of $\epsilon >0$ for this criterion, depending on the desired accuracy.

A common bottleneck of the algorithms introduced by Cai et al. (2010); Mazumder et al. (2010); Ma et al. (2011) as well as other SVT-based approaches is the computation of the SVD of large matrices. This problem is particularly severe in standard matrix completion settings, such as the Netflix problem, where the matrix size is over 400, 000 × 20, 000. However, in our problem,

and $K~10$ in our motivating data example. While algorithms for large matrices are applicable here, the property (8) makes the computation of SVD feasible in practice with generic packages such as PROPACK (Larsen, 2004).

For estimating curves on new data, we need to find $w_i=V{a}_i$ minimizing the prediction error.

To this end, we regress observed values $Y_i$ on $B_{i}V$ with a penalty $\frac{\lambda }{2}{\left(a_i\right)}^2$ on the coeffcients, i.e. we optimize

where $B_i$ are basis functions evaluated at the corresponding timepoints.

While our algorithms for finding $W$ do not depend on the density of the grid, it is important to highlight that with the growing $T$ the norm of the solution $W_T$ grows to infinity and $W_T/\sqrt{T}\to W_0$, for some fixed $W_0$ (See Appendix A). In order to obtain comparable results for different $T$, parameters $\lambda$ need to be adjusted to this scale (for example, by modifying the penalty to $\lambda \sqrt{T}\parallel W{\parallel }_{*}$).

The algorithm converges at a rate $\frac{1}{t}$, where $t$ is the number of iterations (this can be established by an elementary extension of the proofs in Mazumder et al. (2010)) as we demonstrate in Appendix A. Each step of our algorithm requires matrix multiplications and an SVD. Thus the complexity of each step is $O\left(min\left(K{n}^2,n{K}^2\right)\right)$. For EM-based algorithms, the minimization step requires multiple matrix multiplications and therefore we expect a larger constant in each step. For practical use, we report computing times on synthetic data in Section 5.

3.5. $l_0$ regularization

While the nuclear norm relaxation (6) is motivated by making the problem convex, Mazumder et al. (2010) argue that in many cases it can also give better results than the rank constraint. They draw an analogy to the relation between best-subset regression $\left({\ell }_0\right)$ regularization) and LASSO (${\ell }_1$, regularization as in Tibshirani (1996); Friedman et al. (2001)). In LASSO, by shrinking model parameters, we allow more parameters to be included, which can potentially improve predictions if the true subset is larger than the one derived through ${\ell }_0$ regularization. In the case of (6), shrinking the nuclear norm allows us to include more dimensions of $W$ again potentially improving the prediction if the true dimension is high.

In theory, estimating too many components can lead to problems if the true dimension is small. In such cases, shrinking may cause the inclusion of unnecessary dimensions emerging from the noise. To remedy this issue, following the analogy with penalized linear regression, we may consider another class of penalties. In particular, we may consider coming back to the rank constraint by modifying the nuclear norm penalty (6) to ${\ell }_0$. We define $\parallel W{\parallel }_0=rank\left(W\right)$. The problem

has a solution $W={𝒮}_{\lambda }^H\left(YB\right)$, where ${𝒮}_{\lambda }^H\left(X\right)=U{D}_{\lambda }^H{V}^{\prime }$ and $X=UD{V}^{\prime }$ is the SVD of $X$. We refer to ${𝒮}_{\lambda }^H\left(X\right)$ as the hard singular value thresholding (hard SVT) of $X$. We use Algorithm 2 to find the hard SVT of $YB$.

The problem (9) is non-convex, however, by starting from the solutions of Algorithm 1 we explore the space around the global minimum of the ${\ell }_1$, version of the problem. Ge et al. (2016) show that this strategy is successful empirically.

In practice, the continuous parameterization via $\lambda$ leads to similar results as rank-truncation methods for different ranks. There are two advantages of the former. First, it allows for a continuum of warm starts and, as such, it is a natural post-processor for Soft-Impute solutions. Second, it can be a basis for further generalizations with penalties between ${\ell }_0$ and ${\ell }_1$ (Mazumder et al., 2010). However, if the true rank can be estimated with high confidence and the signal is strong, then a rank-truncated method can lead to better predictive performance.

3.6. Dense grid and proximal gradient descent

Growing grid size parameter $T$ improves the accuracy of approximation of the data on the grid, but benefits in approximation are negligible when $T$ becomes sufficiently large. However, we show that $T\to \infty$ also leads to a continuous version of the algorithm, equivalent to the proximal gradient descent method applied to the initial formulation of the problem (1) with nuclear norm penalty.

We start by observing that in Algorithm 1 the Step 2(a)(i) depends on the grid and requires multiplication with all evaluations of the basis on that grid. We can rewrite the expression in $S_{{\lambda }_i}$ such that only a sparse number of evaluations of the basis is required. Indeed, we write

and to compute $W^{old}{B}^{\prime }$ on $\Omega$ in (10) we only need $W^{old}$ and an evaluation of the basis on gridpoints corresponding to $\Omega$. Therefore, for computing $\left(P_{\Omega }\left(Y-W^{old}{B}^{\prime }\right)\right)B$ we also only need $B$ evaluated on gridpoints with a nonzero number of observations. Now, if we increase the number of grid points to infinity we derive a solution where the basis needs to be evaluated only at observed timepoints.

Since $B$ is orthogonal, when $T\to \infty$ the product $\left(P_{\Omega }\left(Y-W^{old}{B}^{\prime }\right)\right)B$ converges to 0 and the update in Step 2(a)(i) in Algorithm 1 also converges to 0 at the rate $1/T$. Therefore, in order to obtain comparable results for different $T$ we need to scale the update direction $\left(P_{\Omega }\left(Y-W^{old}{B}^{\prime }\right)\right)B$ with $\alpha /T$, for some $\alpha >0$. We refer to $\alpha$ as the step size.

We show that at the limit $T\to \infty$ our approach is equivalent to proximal gradient descent. Let us define the loss function as in (1), i.e.

where ${\tau }_{i,j}$ are timepoints, $y_{i,j}$ are observations. We optimize the loss function $f$ with a nuclear norm penalty for some $\lambda >0$

where $W={\left(w_1^{\prime },w_2^{\prime },\dots ,w_N^{\prime }\right)}^{\prime }$. To solve problem (11) we define

and the gradient of $f$ takes form

Note that ${\nabla }_{W}f\left(W\right)$ is corresponding to $\left(P_{\Omega }\left(Y-W^{old}{B}^{\prime }\right)\right)B$ in (10) with infinitely dense grid. For the update step in the proximal gradient descent method, we use step size $\alpha >0$ and singular value thresholding with a threshold $\alpha \lambda$ as recommended in Hastie et al. (2015b). We describe the full procedure in Algorithm 3.

Nesterov (2013) showed sufficient conditions for convergence of Algorithm 3. If $\alpha \in (\mathrm{0,1}/L]$ where $L$ is the Lipschitz constant of ${\nabla }_{w}f$, then the algorithm converges at the rate $1/t$, where $t$ is the number of steps.

3.7. A link between the reduced-rank model and Soft-Longitudinal-Impute

Intuitively, we might expect similarity between the principal directions derived using the probabilistic approach (3) and their counterparts derived from the SVT-based approach. We investigate this relation by analyzing the behavior of SVT for matrices sampled from the probabilistic model given by (3).

For simplicity, let us assume that $\mu =0$ and the data is fully observed on a grid $G$ of $T$ timepoints. Assume that observations $i\in \{1,2,\dots ,N\}$ come from the mixed-effect model

where $\Sigma$ is an unknown covariance matrix of rank $q<K$ and variables $\left(w_i,y_i\right)$ are independent. By the spectral decomposition theorem, we decompose $\Sigma =V\Lambda V^{\prime }$, where $V$ is a $K\times K$ orthogonal and $\Lambda$ is a diagonal $K\times K$ matrix with $q$ positive coefficients in decreasing order. Since $y_i$ and $w_i$ are independent, the distribution (12) can be rewritten as

The model (13) is a factor model with uncorrelated errors and $q$ factors — the first $q$ columns of $BV$.

The following theorem constitutes a link between the mixed-effect model and SVT (Theorem 9.4.1 in Mardia et al. (1980)),

4.1. Illustrative example: Univariate regression

Suppose that the true processes are in a low-dimensional space of some continuous functions (e.g., splines) and that we observe them with noise. More precisely, let

for $1\le i\le N$, where $x_i,y_i,e_{x,i},e_{y,i}$ are $T\times 1$ vectors, $w_i,u_i$ are $K\times 1$ vectors and $B$ is a $T\times K$ matrix of $K$ splines evaluated on a grid of $T$ points. We assume zero-mean independent errors $e_{x,i},e_{y,i}$ with fixed covariance matrices ${\Sigma }_X,{\Sigma }_Y$ respectively, and that the true processes underlying the observed $x_i$ and $y_i$ follow a linear relation in the spline space, i.e.

where $A$ is an unknown $K\times K$ matrix.

Let $X,Y,U,W$ be matrices formed from $x_i^{\prime },y_i^{\prime },w_i^{\prime },u_i^{\prime }$ stacked vertically. From (15) we have $U=WA$, while (14) implies

where $E_X,E_Y$ are matrices of observation noise. Without loss of generality, we assume that $B$ is orthonormal. We have full freedom to choose the basis $B$, since any basis can be orthogonalized using, for example, the singular value decomposition.

Due to the orthogonality of $B$ and after multiplying both expressions in (16) by $B$ we can map the problem to the classical multivariate errors-in-variables models. Let

where ${\tilde{E}}_X=E_{X}B$ and ${\tilde{E}}_Y=E_{Y}B$. In errors-in-variables models it is assumed that the parameters $W$ and $A$ are unknown, and are to be estimated. Both regressors and responses are measured with noise (here ${\tilde{E}}_X$ and ${\tilde{E}}_Y$). The parameter $W$ can be interpreted as a latent representation of both $\tilde{X}$ and $\tilde{Y}$.

The problem of estimating parameters in (17) has been extensively studied in literature dating back to Adcock (1878). Two main methods for estimating parameters in (17) are maximum likelihood estimates (MLE) and generalized least squares estimators (GLSE). The estimators in MLE are derived under the assumption of certain distributions of errors. In GLSE, the only assumption about errors is that they are independent, zero-mean, and they have a common covariance matrix. Then, $\tilde{X}-W$ and $\tilde{Y}-WA$ are zero-mean and estimates for $W$ and $B$ can be found by optimizing some norm of these expressions. Gleser and Watson (1973) show that in the no-intercept model for $\tilde{X}$ and $\tilde{Y}$ of the same size (as in our case) and under the assumption of Gaussian errors, MLE and GLSE give the same estimates of $A$, if GLSE is derived for the Frobenius norm.

In this work, we focus on the GLSE approach as it can be directly solved in our matrix factorization framework and we find it easier to deploy and extend in practice. To account for different magnitudes of the noise in $X$ and $Y$, we consider the optimization problem with weights

where ${\sigma }_X,{\sigma }_Y>0$. Let $\gamma ={\sigma }_X^2/{\sigma }_Y^2$. Then (18) can be transformed to

where $(\cdot :\cdot )$ operator stacks horizontally matrices with the same number of rows. To solve (19), we show that the SVD of $(XB:\gamma YB)$ truncated to the first $K$ dimensions, can be decomposed to $W(I:\gamma A)$. Let $US{V}^{\prime }$ be the SVD of $(XB:\gamma YB)$, with

where each $S_{ij}$ and $V_{ij}$ is a $K\times K$ matrix for $1\le i,j\le 2$ and each $U_i$ is a $N\times K$ matrix for $1\le i\le 2$. By Lemma 2.1 and Lemma 2.2 in (Gleser, 1981) matrix $V_{11}$ is almost surely nonsingular. Therefore, $V_{11}^{-1}$ almost surely exists and we can transform the decomposition such that $(I:\gamma A)={\left(V_{11}^{\prime }\right)}^{-1}\left(\begin{array}{cc}{V}_{11}^{\prime }& V_{21}^{\prime }\end{array}\right)$ and $W=U_1{S}_{11}{V}_{11}^{\prime }$, solving (19).

For partially observed data, if they are very sparse, it might be essential to constrain the rank of the solution. The partially-observed and rank-constrained version of the problem (19) takes the form

$\text{subject}\text{to}\text{rank}\left(W\left(B^{\prime }:\gamma A{B}^{\prime }\right)\right)=k$, where $k$ is the desired rank of the solution and $P_{\tilde{\Omega }}$ is a projection on

As previously, for an unknown $k$ we can introduce a rank penalty using the nuclear norm

The algorithm in the general case of multiple processes is derived in Section 4.2.

Although we motivate the problem as a regression of $Y$ on $X,X$ and $Y$ are symmetric in (19). The low-rank matrix $W$ is, therefore, a joint low-rank representation of matrices $X$ and $Y$ and thus our method can be seen as a dimensionality reduction technique or as a latent space model. In Section 4.2 we extended this idea to a larger number of variables. In Section 4.3 we discuss how this approach can be used for regression.

The linear structure of (16) allows us to draw analogy not only to the errors-in-variables models but also to the vast literature on canonical correlation analysis (CCA), partial least squares (PLS), factor analysis (FA), and linear functional equation (LFE) models. Borga et al. (1997) show that solutions of CCA and PLS can also be derived from the SVD of stacked matrices, as we did with $(XB:\gamma YB)$ in (19). Gleser (1981) thoroughly discusses the relation between errors-in-variables, FA, and LFE.

Finally, the method of using SVD for stacked matrices has also been directly applied in the recommender systems context. Condli et al. (1999) showed that for improving the prediction of unknown entries in some partially observed matrix $Q$ one might consider a low-rank decomposition of ($Q:R$), where $R$ are additional covariates for each row, e.g., demographic features of individuals.

4.2. Dimensionality reduction

Suppose that for every individual we observe multiple variables varying in time (e.g. results of multiple medical tests at different times in a clinic) and we want to find a projection on $R^d$ maximizing the variance explained for some $d\in N$. This projection would correspond to characterizing patients by their progression trends of multiple metrics simultaneously.

We extend the equation (19) to account for a larger number of covariates, and we do not impose decomposition of the solution yet. We formulate the following optimization problem

where $X_i$ are some $N\times T$ matrices corresponding to the processes measured on a grid of $T$ points, $B$ is a basis evaluated on the same grid and orthogonalized (a $T\times K$ matrix), and $W$ is a $N\times pK$ matrix.

Let $B=I_p\otimes B$ be the Kronecker product of $p\times p$ identity matrix and $B$, i.e. a $pT\times pK$ matrix with $B$ stacked $p$ times on the diagonal, and let $X=\left(X_1:X_2:\dots :X_p\right)$. Matrix $B$ is orthogonal and therefore results developed in Section 2.3 apply here. In particular, singular value thresholding solves

and we can use Algorithm 1 for solving

where $P_{\Omega }$ is the projection on observed indices $\Omega$.

The optimization problem (21) can be further extended. For example, we can choose a different basis for each process or scale individual contributions of processes using scaling factors. By modifying Equation (22) correspondingly, we can extend our solution to these cases.

4.3. Regression

In practice, we are often interested in the relationship between the progression of an individual parameter (e.g., cancer growth) and some individual features constant over time (e.g., demographics) or progressions of other covariates (e.g., blood tests, vitals, biopsy results). We show how our framework can be used to exploit relations between curves in order to improve the estimation of the progression of the main parameter of interest.

We start with a regression problem with fixed covariates and sparsely observed response trajectories. Assume that for each subject $1\le i\le N$ we observe covariates $x_i\in R^d$ and samples from trajectories $y_i\in R^{n_i}$ following the distribution

where $B_i$ is the matrix of evaluations of the basis on the same evaluation points as $y_i$ (See Section 2.2). We map the problem to the matrix completion framework as in Section 3. Let $X$ be a $N\times d$ matrix of observed covariates, $Y$ be a sparsely observed $N\times T$ matrix of trajectories, and $B$ be a $T\times K$ matrix representing a basis of $K$ splines evaluated on a grid of $T$ points. We consider the optimization problem

where $A$ is a $d\times K$ matrix and $P_{\Omega }$ is a projection on the observed indices $\Omega$. To solve (23) we propose an iterative Algorithm 4.

Suppose we want to incorporate other variables varying in time for the prediction of the response process. We can directly apply the method proposed in Section 4.2 and model the response and regressors together. However, it might be suboptimal for prediction, as it optimizes the least-squares distance in all variables rather than only the response. This difference is analogous to the difference between the regression line of some univariate $y$ on independent variables $x_1,\dots ,x_p$ and the first principal component analysis of $\left(y,x_1,\dots ,x_p\right)$, used for the prediction of $y$. While the first one minimizes the distance to $y$, the latter minimizes the distance to $\left(y,x_1,\dots ,x_p\right)$, which is usually less efficient for predicting $y$.

Alternatively, we can use methodology from Section 4.2 only for covariates. The solution of (22) can be decomposed into $W=US{V}^{\prime }$, and we regress $Y$ on $U$ as in (23). Let $U$ be an $N\times d_2$ orthogonal matrix derived from $\left(X_1,\dots ,X_p\right)$, where $d_2$ is the numbers of latent components. We search for a $d_2\times K$ matrix $A$ solving

where $P_{\Omega }$ is a projection on a set of observed coefficients. We propose a two-step iterative procedure (Algorithm 5).

5.1. Data generation

We seek to simulate data similar to observations in Section 6. We want to parametrize the processes to draw insights depending on sample size and relation between variables. Let $G$ be a grid of $T$ equidistributed points and let $B$ be a basis of $K$ spline functions evaluated on the grid $G$. We simulate three $N\times K$ matrices using the same procedure $𝒢(r,K,N)$, where $r,\in \underset{+}{R^K}$:

1. Define the procedure ${ℳ}_{\left(r\right)}$ generating symmetric matrices $K\times K$ for a given vector $r\in R_{+}^K:$

   1. Simulate $K\times K$ matrix $S$

   2. Use SVD to decompose $S$ to $UDV\text{'}$

   3. Return $Q=Vdiag\left[r\right]V^{\prime }$, where $diag\left[r\right]$ is a diagonal matrix with $r$ on the diagonal

2. Let $\Sigma ={ℳ}_{\left(r\right)}$ and $\mu ={𝒩}_{\left(0,I_K\right)}$.

3. Draw

   $N$

   vectors

   $v_i$

   from the distribution

4. Return a matrix with rows ${\left(v_i\right)}_{1\le i\le N}$.

We chose values of $r$ close to values on the diagonal of the spectral decomposition of covariance structure estimated in Section 6. We define

Let $X_1,X_2,Z$ be generated following the procedure $𝒢(r,K,N)$ and let $Y=Z+X_1+X_2$. We consider $X_1,X_2$, and $Y$ as coefficients in a spline space $B$. We derive corresponding functions by multiplying these matrices by $B^{\prime }$, i.e. $X_{f,1}=X_1{B}^{\prime },X_{f,2}=X_2{B}^{\prime }$ and $Y_f=Y{B}^{\prime }$. We set $N=100,K=7$ and $T=31$.

In order to reflect a realistic sampling of curves, we use timepoints at which subjects in the data study were observed (see Section 6). As we are particularly interested in predictive performance, we select subjects with at least 4 observed timepoints to enable leave-one-out estimation. For each simulated curve we randomly select a subject and use grid-aligned timepoints at which this subject visited the clinic. We acknowledge that specifics of sampling might affect relative performance of methods, however a thorough analysis of sampling is beyond the scope of our work. Here we only focus on the case directly derived from our data study.

Each observed element of each matrix $X_{f,1},X_{f,2}$ and $Y$ is drawn with Gaussian noise with mean 0 and standard deviation 0.25. Our goal is to recover $Y$ from such sparsely sampled noisy data.

5.2. Methods

We compare Soft-Longitudinal-Impute (SLI) defined in Algorithm 1, Proximal-Gradient (PG) method defined in Algorithm 3, the fPCA procedure (James et al., 2000), implemented by Peng and Paul (2009), fdapace (Zhou et al., 2022), and face (Xiao et al., 2018). The first three algorithms require a set of basis functions as one of the inputs. In all cases, we use a polynomial $B$. In SLI and PG we also need to specify the tuning parameter $\lambda$, while in fPCA we need to choose the rank $R$. We use cross-validation to choose $\lambda$ and $R$ by optimizing for the prediction error on held-out (validation) observations.

We divide the observed coefficients into training (81%), validation (9%), and test (10%) sets. We choose the best parameters of the three models on the validation set and then retrain on the entire training and validation sets combined. We compute the error by taking the mean squared Frobenius distance between $Y$ and estimated $\hat{Y}$, i.e.

on the test set $S$.

We train the algorithms with all combinations of parameters: regularization parameter for SLI and PG procedures $\lambda \in \{10,15,20,\dots ,50\}$ and the rank for fPCA procedure $d\in \{2,3,4\}$. We compare the five methods fPCA, fdapace, face, PG, and SLI, to the baseline null model which we define as the population means across all visits.

For testing our multivariate approach, we model $\left(Y,X_1,X_2\right)$ together and we compare our Sparse-Longitudinal-Regression (SLR) with MFACE package (Li et al., 2020). We compare compute time and predictive performance in terms of recovering of $Y$.

5.3. Results and discussion

The SLI and PG methods achieve superior performance at a good running time as presented in Figure 4. Fdapace (Zhou et al., 2022) outperforms all methods in terms of speed, but it performs worse in terms of the mean squared error.

Results indicate that both our methods, SLI and PG, should be considered on par with fdapace and face for practical applications since each method has its own advantages. A solution derived from fdapace might also be used as a warm start for other algorithms to speed them up. We also observed that the performance gap decreases with the growing number of curves.

One additional advantage of SLI and PG is that they return a smooth space of solutions, parametrized with $\lambda$. In Figure 5, we present the estimated rank and cross-validation error of one of the simulation runs. Compared to fixed-rank alternatives, our methods allow for tuning the predictive accuracy and the rank simultaneously in a smooth manner.

In the multivariate case, both mface and our SLR successfully recovered addional information from variables $X_1$ and $X_2$ as their NMSE in modeling $Y$ was better compared to univariate counterpart methods (SLI and face respectively). SLR outperformed mface in terms of NMSE, but required more time to converge. We conclude that SLR should also be considered a viable candidate method for multivariate sparsely sampled functions on function regression.

6.1. Materials and methods

We analyze a dataset of Gillette Children’s Hospital patients visiting the clinic between 1994 and 2014, ages ranging from 4 to 19 years, mostly diagnosed with Cerebral Palsy. The dataset contains 84 visits of 36 patients without gait disorders and 6066 visits of 2898 patients with gait pathologies.

Motion capture data was collected at 120Hz and joint angles in time were extracted. These joint angles were then normalized in time to the gait cycle. Points from these curves were then subsampled (51 equidistributed points) for the downstream analysis established in the clinic.

In the dataset which we received from the hospital, for each patient we know their birthday and disease subtype. From each visit, we observe the following variables: patient ID, time of the visit, GDI of the left leg, GDI of the right leg, walking speed, and O₂ expenditure. Other clinical variables that we received were not included in this study. Walking speed is related to information we lose during the normalization of the gait cycle in time. O₂ expenditure is a measure of a subject’s energy expenditure during walking. Pathological gait is often energy inefficient and reduction of O₂ expenditure is one of the objectives of treatments. Finally, GDI is computed for two legs while in many cases the neurological disorder affects only one limb. To simplify the analysis, we focus on the more impaired limb by analyzing the minimum of the left and the right GDI.

Our objective is to model individual progression curves. We test two methods: functional principal components (fPCA) and Soft-Longitudinal-Impute (SLI). We compare the results to the null model – the population mean across all visits (mean). In SLR, we approximate GDI using latent variables of sparsely observed covariates O₂ expenditure and walking speed, following the methodology from Section 4.3.

Let us denote the test set as $\Omega \subset \{1,2,\dots ,N\}\times \{1,2,\dots ,T\}$. We validate each model $M$ on heldout indices by computing the mean squared error as defined in (25). We select the parameters of each of the three methods using cross-validation, using the same validation set.

In our evaluation procedure, for the test set, we randomly select 5% of observations of patients who visited the clinic at least 4 times. We chose the cutoff at 4 in order to have at least 3 observations to fit the curve, and keep 1 observation for testing. We verified with the collaborators at the hospital that three visits are reasonable for the practical use of our methods. We split the remaining 95% of observations into training and validation sets in 90: 10 proportion. We also ran the procedures using a smaller test set and did not observe any relevant difference in results, likely because the training set is already relatively large for the task. We train the algorithms with the following combinations of parameters: the regularization parameter for SLI and SLR procedures $\lambda \in \{0,0.1,0.2,\dots ,2.0\}$ and the rank for fPCA procedure $d\in \{2,3,4,\dots ,K\}$. We define the grid of $T=51$ points. We repeat the entire evaluation procedure 20 times.

6.2. Results

Compared to the null model, fPCA and SLI explain around ~30% of the variance. We present detailed results in Table 1. SLR using additional predictors, O₂ expenditure and walking speed, yielded a mean MSE of 0.68 with a standard deviation of 0.8. We conclude that O₂ expenditure and walking speed provide additional information for the prediction of GDI progression.

Both fPCA and Sparse-Longitudinal-Impute provide latent representations of patients’ progression curves. We analyze the singular value vectors from our SVD solution which we refer to as principal components. In the left plot in Figure 3 we show the first two estimated principal components. We found that the first component estimates the change between GDI before and after age of 20. The second component model changes around the age of 10 and around the age of 18. In the right plot in Figure 3, by adding a principle component to the population mean curve, we illustrate how differences in the first component are reflected in the patient’s trajectory. By visual investigation of curves returned by our Sparse-Longitudinal-Impute and by fPCA we found similar trends in the first two components.

Since our SVD decomposition defines a low-rank representation of the progression trends, we can also use it to gain insights into progression in different groups of patients. In Cerebral Palsy, we divide paralysis into subtypes depending on which limbs are affected: monoplegia (one leg), diplegia (two legs), hemiplegia (one side of the body), triplegia (three limbs), quadriplegia (four limbs). Hemiplegia is the most prevalent in our population and it might be divided depending on severity, from type I (weak muscles, drop foot) to type IV (severe spasticity). We find differences between trends of progression for different subtypes of paralysis of patients $\left(F_{\mathrm{6,541}}=17.17,p<{10}^{-15}\right)$. We illustrate these distributions in Figure 7.