Functional regression with multivariate responses

We consider the functional regression model with multivariate response and functional predictors. Compared to fitting each individual response variable separately, taking advantage of the correlation between the response variables can improve the estimation and prediction accuracy. Using information in both functional predictors and multivariate response, we identify the optimal decomposition of the coefficient functions for prediction in population level. Then we propose methods to estimate this decomposition and fit the regression model for the situations of a small and a large number $p$ of functional predictors separately. For a large $p$, we propose a simultaneous smooth-sparse penalty which can both make curve selection and improve estimation and prediction accuracy. We provide the asymptotic results when both the sample size and the number of functional predictors go to infinity. Our method can be applied to models with thousands of functional predictors and has been implemented in the R package FRegSigCom.

💡 Research Summary

The paper addresses linear regression where the predictors are functional curves and the response is multivariate. While most existing work treats each response component separately or applies scalar‑on‑function methods (FPCR, FPLS) without exploiting inter‑response correlation, this study proposes a unified framework that leverages that correlation for both estimation and prediction.

The model is written as
(Y = \mu + \sum_{j=1}^{p}\int_{0}^{1} X_{j}(t),b_{j}(t),dt + \varepsilon),
with (Y\in\mathbb{R}^{m}) and (p) functional predictors (X_{j}(t)). The coefficient matrix (B(t)=(b_{1}(t),\dots,b_{p}(t))^{\top}) is a (p\times m) function.

A key contribution is the identification of an optimal decomposition of the coefficient functions of the form
(B(t)=\sum_{k=1}^{K}\alpha_{k}(t)W_{k}),
where each (\alpha_{k}(t)) is a scalar function common to all predictors and each (W_{k}) is a (p\times m) matrix. The optimal ({\alpha_{k},W_{k}}) are obtained by solving a generalized eigenvalue problem that maximizes the covariance between the projected predictor (\int X(t)\alpha(t)dt) and the multivariate response, subject to a unit‑variance constraint. The eigenvalues (\sigma_{k}^{2}) quantify the contribution of each component; the sum of the discarded eigenvalues gives the minimal mean‑squared error for a truncation at (K). Theorem 2.1 proves that this decomposition yields the smallest possible approximation error for the regression function and that at most (m) components are needed for a perfect representation.

For the case of a single predictor ((p=1)), the authors derive the eigenfunctions explicitly and compare them with those generated by FPCR (which uses the predictor covariance) and FPLS (which maximizes squared covariance with the response). Numerical illustrations show that the optimal components differ substantially and lead to markedly lower approximation error.

When (p>1), the same eigenfunctions are shared across predictors, providing a natural multivariate dimension reduction that respects both predictor and response structures.

Estimation strategies are split into two regimes.

• Small‑(p) regime: standard spline basis expansion and penalized least squares are employed.
• Large‑(p) regime: a simultaneous smooth‑sparse penalty is introduced:
  (\displaystyle \sum_{j=1}^{p}\big{\lambda_{1}\int (b_{j}’’(t))^{2}dt + \lambda_{2}|b_{j}|_{1}\big}).
  The first term enforces smoothness via a second‑derivative penalty, while the second term induces sparsity at the curve level (group‑Lasso style). An alternating coordinate‑descent algorithm (ADMM‑like) solves the penalized problem efficiently, scaling to thousands of functional predictors.

The authors develop asymptotic theory for the joint limit (n\to\infty) and (p\to\infty). They establish consistency of the estimated eigenfunctions and coefficient matrices, selection consistency (the probability of correctly identifying non‑zero predictor curves tends to one), and rates for the prediction risk that depend on the decay of the eigenvalues (\sigma_{k}^{2}). Stronger correlation among response components accelerates eigenvalue decay, allowing fewer components to achieve a given accuracy.

Extensive simulation studies explore various correlation levels ((\rho=0,0.4,0.8)), noise variances, and numbers of predictors (up to 1,000). The proposed method consistently outperforms FPCR and FPLS in terms of mean‑squared prediction error and variable‑selection metrics (precision, recall, F1).

A real‑data application uses near‑infrared (NIR) spectra of corn samples (thousands of wavelengths) to predict multiple nutritional outcomes (protein, fat, fiber, etc.). The multivariate functional regression with the smooth‑sparse penalty selects a parsimonious set of informative wavelengths and achieves a 15 % reduction in root‑mean‑square error compared with separate scalar‑on‑function models.

All methodology, including eigenfunction computation, penalized estimation, and cross‑validation for tuning parameters, is implemented in the R package FRegSigCom, which the authors make publicly available. The package provides user‑friendly functions for both small‑ and large‑(p) settings, making the approach readily applicable to high‑dimensional functional data.

In summary, the paper delivers a comprehensive solution for multivariate functional linear regression: it defines an optimal population‑level decomposition, supplies practical estimation procedures for both low‑ and high‑dimensional predictor spaces, proves rigorous asymptotic properties, and validates the approach through simulations and a substantive real‑world example. This work significantly advances the ability to exploit response correlation in functional data analysis and offers a scalable tool for modern high‑throughput applications.