Functional Partial Least-Squares: Adaptive Estimation and Inference

Functional Partial Least-Squares: Adaptive Estimation and Inference

\def\spacingset#1{ {#1}} \spacingset{1}

\if00 \fi

\if10 {

} \fi

{\it Keywords:} Functional Partial Least-Squares, Inference, Rate Optimal and Adaptive Estimation, Functional Linear Regression, Climate Science.

\spacingset{1.9}

Introduction

With the increasing availability of data, functional data analysis has become widely applied across fields such as chemometrics, climate science, and economics. In this paper, we study a linear functional regression model with a scalar response $Y$ and functional predictor $X$:

The primary objective is to estimate the functional slope $\beta $ to predict the response $Y$. Both $X$ and $\beta $ lie in an infinite-dimensional Hilbert space, and the high dimensionality of $\beta$ leads to an ill-posed inverse problem. Consistent estimation of the slope coefficient thus requires either a dimension-reduction technique or some form of regularization.

Two main approaches have been popular in the functional data analysis literature: penalization methods, e.g., cardot2003spline, li2007rates, crambes2009smoothing, yuan2010reproducing, cai2012minimax, florens2015instrumental, and functional principal component analysis (PCA), e.g., cardot1999functional, ramsay2002applied, cai2006prediction, hall2007methodology. The PCA-based approach approximates $\beta$ using a finite expansion over the leading principal components, which correspond to the eigenfunctions of the covariance operator of X. As noted by jolliffe1982note, this method performs well when the response Y is primarily correlated with the leading principal components. Moreover, accurate recovery of the principal components requires a sufficient separation between the eigenvalues of the covariance operator.

In this paper, we explore an alternative method, functional partial least squares (PLS). PLS is a widely used technique in the statistical learning, see friedman2009elements, frank1993; economics and finance, see carrasco2016sample and kelly2015three; chemometrics, see helland1988structure, wold1984collinearity. However, it has been somewhat less commonly employed in the empirical functional data analysis. The method constructs components as linear transformations of the predictors X designed to maximize their correlation with the response Y. As a result, fewer components are typically needed to achieve good predictive performance compared to PCA. It was first introduced in preda2005pls to solve the high-dimensional problems with multicolinearity associated with the scalar-on-function linear model. Several interesting papers compare PCA and various functional PLS estimators with data-driven choice of latent components in simulations; see reiss2007functional, Kraemer2008, baillo2009, aguilera2010basis, febrero2017overview, and saricam2022partial. This literature concludes that the prediction ability of functional PCA and PLS approaches is similar but functional PLS requires fewer components and provides a much more accurate estimation of the parameter function than PCA.

However, the existing literature lacks supporting theoretical results, which are challenging to obtain given that the PLS estimator depends non-linearly on the response variable and is computed iteratively. An important step toward the theoretical analysis of functional PLS was made by delaigle2012methodology, who introduced an alternative PLS (APLS). Nevertheless, it remains unknown whether: (1) functional PLS achieves (nearly) minimax-optimal convergence rates under weak conditions; (2) a rate-adaptive method to select the number of PLS components exists; and (3) inference for functional PLS can be conducted.\footnote{Deriving an approximation to the distribution of PLS is particularly challenging because PLS is a nonlinear estimator converging at a nonparametric rate.} This paper fills these gaps and provides a comprehensive, rigorous theoretical analysis of functional PLS.

The paper makes several original contributions. First, we derive the convergence rate of our estimator and the prediction error under a source condition. This condition measures the complexity of the problem through the so-called degree of ill-posedness and relates the slope coefficient $\beta$ to the spectral decomposition of the covariance operator. Our results do not require assuming a separation between adjacent eigenvalues, as is common for PCA (see, e.g., hall2007methodology), and hold even in the presence of repeated eigenvalues. Second, we establish a lower bound on the minimax convergence rate and show that our estimator is (nearly) minimax-optimal. Because the optimal number of PLS components depends on the unknown degree of ill-posedness, we propose an adaptive early stopping rule for selecting the number of components. We show that this single, early stoping rule yields a rate-optimal PLS estimator simultaneously for estimation and prediction errors with high probability. We also characterize how the selected number of components evolves with sample size under various scenarios.

Lastly, we develop new test and confidence sets for PLS. We show that the test can detect alternatives converging at a parametric rate. Interestingly, while early stopping is crucial for estimation and prediction to prevent overfitting, it should be avoided for inference. An efficient iterative algorithm to compute the estimator is provided. Our simulation results reveal that the estimator outperforms several alternative methods combined with cross-validation in terms of estimation error and remains competitive for prediction. We also find that our test has excellent power properties with samples as small as $n=100$ observations.

To establish our theoretical results, we rely on the inverse problem literature and exploit the close connection between PLS and the conjugate gradient method, as presented in hanke1995conjugate, engl1996regularization, and blanchard2016convergence. To the best of our knowledge, this connection has not previously been established in the functional data analysis literature. Recently and independently of our work, gupta2023convergence proposed an estimator in a reproducing kernel Hilbert space (RKHS) generated by a specified kernel and used conjugate gradient methods to regularize the solution. Although their method is similar to ours, the results are not directly comparable. Their main result concerns the convergence rate of an estimator in an RKHS, established under different assumptions. In particular, they assume a polynomial decay rate of certain operator eigenvalues, while we do not require any decay assumptions. Their source condition is also different, and it is unclear which one is weaker. Similarly, lin2021kernel studied conjugate gradient methods in the context of linear approximations to nonparametric regression with randomized sketches and Nyström sampling. While their results apply to general Hilbert spaces under stronger assumptions, they do not explicitly connect to the functional data analysis literature; for example, their simulations focus on linear approximations of $\ensuremath{\mathbf{E}}[Y|X]=f(X)$ with real-valued data $Y,X\in\mathbb{R}$ and Sobolev RKHS kernels. Neither paper shows that the proposed estimators are rate-adaptive or provide a practical early stopping rule for selecting the number of PLS components.

To illustrate the practical relevance of our results, we apply our method to climate science. Using a fine-grained county-level dataset of U.S. crop yields and temperatures recorded over 70 years, we estimate the impact of temperature on crop yields. We find that the critical temperature at which annual crop yields begin to decline is around 30${}^\circ$C, consistent with schlenker2009nonlinear, who relied on highly parameterized least-squares estimators. Our method provides additional insights by examining how the temperature effect curves have evolved over time. Interestingly, we find that the detrimental effects of high temperatures on corn and soybean yields have diminished over time, likely reflecting farmers’ adaptive actions, including the use of more resilient crops and improved irrigation systems. However, this finding is not conclusive when accounting for statistical uncertainty due to the small number of observed extreme temperatures.

The rest of the paper is organized as follows. Section (ref) introduces the functional regression model and the functional PLS estimator. Section (ref) establishes the theoretical properties, including the convergence rates for estimation and prediction errors, the minimax lower bound, and the adaptivity of the early stopping rule. Section (ref) develops inference for PLS. Section (ref) presents a Monte Carlo study. Section (ref) discusses an empirical application to nonlinear temperature effects in agriculture. Section (ref) concludes. Supplementary Material provides all proofs, comparisons to alternative estimators, and additional simulation results.

Functional Regression and PLS Estimator

Functional Linear Regression

Throughout the paper, we consider a generalized version of the functional linear regression model

where $(Y,X)\in \mathbb{R}\times \mathbb{H}$, $\beta \in \mathbb{H}$ is the unknown functional slope coefficient, and $(\mathbb{H} ,\langle.,.\rangle)$ is a separable Hilbert space with the induced norm $\|.\|=\sqrt{\langle.,.\rangle}$. The model in equation ((ref)) corresponds to the Hilbert space of square integrable functions, $\mathbb{H}=L_2[0,1]$, with the norm induced by the inner product $\langle f,g\rangle = \int_0^1f(s)g(s)$.

If $\mathbb{E}[X]=0$, the covariance restriction $\mathbf{E} [\varepsilon X]=0$ implies that the slope coefficient $\beta \in \mathbb{H}$ solves the moment condition

where $r\in \mathbb{H}$ and $K:\mathbb{H}\to \mathbb{H}$ is a compact covariance operator with summable eigenvalues whenever $\mathbf{E} \|X\|^2<\infty$. It is well-known that the inverse operator $K^{-1}$ is discontinuous and solving the equation $K\beta = r$ for $\beta$ is an ill-posed inverse problem; see carrasco2007linear, engl1996regularization, and hanke1995conjugate.

Roughly speaking, there are two popular strategies to regularize such problems:

Examples of (a) include the Tikhonov regularization when $R_\alpha(\lambda)=(\alpha+\lambda)^{-1}$, the spectral cut-off when $R_\alpha(\lambda)=\lambda^{-1}\mathbf{1}_{\lambda \geq \alpha}$ and the Landweber iterations. On the other hand, the estimators in group (b), often solve the empirical least-squares problem

where $\|v\|_n^2=v^\top v/n,v\in \mathbb{R}^n$ and we put $\mathbf{y} =(Y_1,\dots,Y_n)^\top$ and

for an i.i.d. sample $(Y_i,X_i)_{i=1}^n$. The basis $(h_j)_{j=1}^m$ spanning $\mathbb{H}_m$ can be either fixed (e.g. Fourier, polynomials, splines, wavelets) or adaptively constructed from the data.

The data-driven bases are especially attractive since they can adapt to the features of the population represented by the data and can approximate the slope parameter $\beta \in \mathbb{H}$ more efficiently; see delaigle2012methodology. The principal component analysis (PCA) \footnote{Using the PCA basis is also related to the spectral cut-off method described in (a).} and the partial least-squares (PLS) are two widely used methods to construct adaptive bases in practice. The PCA basis is constructed by identifying the directions in $\mathbb{H}$ where $X$ varies the most while the PLS basis is constructed in a supervised way taking into account the response variable as well. While the first $m$ elements of the PCA basis $ h_1,\dots,h_m$ usually capture most of the variation of $X$, these are not necessarily the most important vectors for approximating $\beta$ or predicting the response variable $Y$. It is easy to find empirical examples, where some of the last few low-variance components are important; see jolliffe1982note who documented the issue on datasets used in economics, climate science, chemical engineering, and meteorology.

PLS estimator

{ The PLS estimator constructs a data-driven basis iteratively maximizing the covariance with the response variable $Y$; see preda2005pls who introduced it in the functional data analysis setting. The iterative nature of the estimator makes it difficult to analyze its statistical properties. This prompted delaigle2012methodology to formulate an alternative functional PLS solving the problem in equation ((ref)) over the so-called Krylov subspace

where

are the estimators of $r$ and $K$; see also wold1984collinearity, helland1988structure, and phatak2002exploiting for the link between PLS and Krylov subspaces.

In this paper, we study a version of the PLS estimator with $m\geq 1$ components, denoted $\hat \beta_m$, characterized as a solution to the least-squares problem

over the Krylov subspace $\mathbb{H}_m$. The least-squares objective function is weighted by the adjoint operator of $T_n$

and corresponds to minimizing the first-order conditions to the problem in equation ((ref)), often called the normal equations. Equivalently, $\hat \beta_m$ fits the empirical counterpart to the equation ( (ref))

as it is easy to see that $\hat r = T_n^*\mathbf{y}$ and $\hat K = T_n^*T_n$. Importantly, the PLS estimator formalized in equation ((ref)) corresponds to the conjugate gradient method with a self-adjoint operator $\hat K$, cf. hestenes1952methods, known for its excellent regularization properties; see also hanke1995conjugate and nemirovski1986regularizing.\footnote{The method of conjugate gradients is one of the most efficient algorithms for solving high-dimensional systems of linear equations; see also nocedal1999numerical and references therein.} We provide a more detailed comparison between the two PLS estimators in the Supplementary Material, Section (ref). A related formulation of the PLS in reproducing kernel Hilbert spaces (RKHS) was recently studied in an independent work of gupta2023convergence who focus on the estimation error only and impose assumptions different from ours. Our work can be seen as using a kernel naturally adapted to the data which is unknown in practice. }

{ The estimator is uniquely defined for every $m\leq n_*$, where $ n_*$ is the number of distinct non-zero eigenvalues of $\hat K$; see Proposition (ref) in the Supplementary Material. It is also easy to see that for every $m\geq 1$, we have $\hat \beta_m = \hat P_m(\hat K)\hat r$ for a polynomial $\hat P_m(\hat K) = \sum_{j=1}^ma_j\hat K^{j-1}$ with coefficients $\mathbf{a} :=(a_1,\dots,a_m)^\top$, where $\hat P_0=0$ and $\hat \beta_0=0$. The coefficients vector solves the system of $m$ linear equations

where $\mathbf{K}:=\langle \hat K^j\hat r,\hat K^k\hat r \rangle_{1\leq j,k\leq m}$ and $\mathbf{r}:= \langle \hat K^j\hat r,\hat r\rangle_{1\leq j\leq m}$. From the practical point of view, it is more efficient to use an iterative conjugate gradient algorithm that bypasses the (potentially unstable) matrix inversion with an iterative multiplication by the operator $ \hat K$; see Algorithm (ref) in Section (ref). }

Adaptive Estimation

In this section, we will show that the functional PLS estimator achieves the (nearly) optimal convergence rate on a class of ellipsoids. We consider an early stopping rule to select the number of PLS components and show that it adapts to the complexity of the ellipsoid. Lastly, we study how rapidly, the number of selected components increases with the sample size and make some comparisons to the PCA estimator.

Optimal Convergence Rates

Since the operators $K:\mathbb{H}\to\mathbb{H}$ and $\hat K:\mathbb{H}\to \mathbb{H}$ are self-adjoint and compact, by the spectral theorem

where $\lambda_1\geq \lambda_2\geq \dots\geq 0$ and $\hat\lambda_1\geq \hat\lambda_2\geq \dots\geq \hat\lambda_{n}\geq0$ are the eigenvalues of $K$ and $\hat K$ and $(v_j)_{j=1}^\infty$ and $(\hat v_j)_{j=1}^{n}$ are the corresponding eigenvectors; see kress1999linear, Theorem 15.16. Note that the sample covariance operator $\hat K$ is a finite-rank operator with at most $n_*\leq n$ distinct non-zero eigenvalues.

For any bounded and measurable function $\phi:\mathbb{R}_+\to\mathbb{R}_+$, we define functions of operators through their spectral decompositions:

These definitions are commonly used in the inverse problems literature; see engl1996regularization.

The following inequalities for the operator norm will be often used:

where $\|A\|_{\rm op} =\sup_{\|x\|=1}\|Ax\|$.

We shall introduce several relatively mild assumptions on the distribution of the data next.

Assumption (ref) imposes mild restrictions on the data-generating process. Note that $\ensuremath{\mathbf{E}}\|X\|^4<\infty$ is satisfied when $X$ is a Gaussian process in $\mathbb{H}$. It implies that $K$ is a nuclear operator and, hence, compact.

Assumption (ref) ensures that the slope parameter $\beta$ is identified. If this assumption is violated, the focus would shift to the identified component of $\beta$ within the orthogonal complement of the null space of $K$; see babii2017completeness and engl1996regularization.

Assumption (ref) describes the complexity of the ill-posed inverse problem in terms of the smoothness of $\beta$ and the smoothing properties of the operator $K$. The parameter $\mu$ is known as the degree of ill-posedness. It restricts the rate of decline of the generalized Fourier coefficients $\langle \beta,v_j\rangle_{j\geq 1}$ relatively to the eigenvalues of $K$. A larger value of $\mu$ means that it is easier to estimate the slope coefficient $\beta$; see also carrasco2007linear. Recall also that the summability of eigenvalues holds whenever $\ensuremath{\mathbf{E}}\|X\|^2<\infty$. Note that Assumption (ref) is weaker than what is typically used to analyze the PCA estimators and does not require restricting the spacing between eigenvalues, cf. hall2007methodology.

Consider now the so-called residual polynomial $\hat Q_m(\lambda)=1-\lambda\hat P_m(\lambda)$, deriving its name from the identity $\hat r - \hat K\hat\beta_m = \hat Q_m(\hat K)\hat r$. It is known that the polynomial, $\hat Q_m$, has $m$ distinct real roots, denoted $\hat\theta_1>\hat\theta_2>\dots>\hat\theta_m>0$. The sum of inverse of these roots,

plays an important role in the analysis of the conjugate gradient regularization; see Lemma (ref) in the Supplementary Material.

Our first result characterizes the convergence rate of the estimation and prediction errors of the PLS estimator.

Note that the last condition in Theorem (ref) imposes that the number of components $m$ does not increase too fast with the sample size and is not binding. In fact, it is optimal to have $|\hat Q_m'(0)|\sim n^{\frac{1}{2(\mu+1)}}$, in which case we obtain the following convergence rate

When $s=0$, this shows that the convergence rate of PLS in the Hilbert space norm is of order $n^{-\frac{\mu}{\mu+1}}$. On the other hand, when $s=1/2$, we obtain the convergence rate of the out-of-sample prediction error, since

where $\ensuremath{\mathbf{E}}_X$ is taken with respect to $X$, independent of $(Y_i,X_i)_{i=1}^n$.

The following result shows that no estimator can achieve a faster than $n^{-\frac{\mu +s}{\mu +1}}\log^{-b} n$ rate on the class $\mathcal{S}(\mu,R,C)$.

Therefore, we conclude that the PLS estimator $\hat\beta_m$ achieves the (nearly) optimal convergence rate on $\mathcal{S}(\mu,R,C)$, simultaneously for the estimation $(s=0)$ and prediction $(s=1/2)$ errors.\footnote{It is possible to avoid the $1/\log n$ factor by considering the larger class of Hilbert--Schmidt operators.}

Adaptive PLS estimator

Next, we look at the adaptive PLS estimator, where the number of components is selected using the early stopping rule described in the following assumption.

Assumption (ref) states that the selected number of PLS components $\hat m$ equals the first non-negative integer $m$ such that the norm of the fitted “moment" is smaller than a certain threshold; see Supplementary Material, Section (ref) for a practical implementation of this early stopping rule. Note that the number of selected components is finite since $\hat m\leq n_*$, where $n_*$ is the number of distinct non-zero eigenvalues of $\hat K$; see Proposition (ref) in the Supplementary Material. In fact, the norm of “residual" is zero for $m\geq n_*$ in which case we have perfect overfitting.

The following result shows that the early stopping rule in Assumption (ref) is adaptive to the unknown degree of ill-posedness $\mu>0$.

Taking $\delta_n=1/\log n$ in Assumption (ref), we obtain from Theorem (ref) the convergence rate of the estimation and prediction errors of PLS with the number of components is selected iteratively with the early stopping rule:

Therefore, the adaptive PLS achieves the (nearly) optimal convergence rate simultaneously for the estimation and prediction errors without knowing the degree of ill-posedness $\mu>0$.

Number of Selected Components

In this section, we look at how rapidly the number of selected components in Assumption (ref) increases with the sample size. First, we consider a somewhat conservative bound that does not impose any assumptions on the spectrum of the operator $K$.

Taking $\delta = 1/\log n$, we obtain from Theorem (ref) that $\hat m = O_P\left((n/\log n)^\frac{1}{4(\mu+1)}\right)$. Next, we consider sharper estimates under additional assumptions imposed on the spectrum of the operator $K$.

Theorem (ref) shows that if the eigenvalues decline polynomially fast, then the selected number of components is $\hat m = O_P(n/\log n)^{\frac{1}{4(\kappa +1)(\mu+1)}}$ while in the case of the geometric decline, the number of selected components increases slowly with the sample size. Therefore, the adaptive early stopping rule will choose a smaller number of components if the eigenvalues of the operator $K$ decline faster and vice versa.

Inference

A Hypothesis Test

We aim to test the null hypothesis against the fixed alternative hypothesis:

for some known $b \in \mathbb{H}$. We are also interested in testing against a sequence of local alternative hypotheses

for some $\Delta\in\mathbb{H}$.

Recall that our main result in Theorem (ref) shows that

provided that the number of components $m$ is properly selected. The convergence rate is $O_P(n^{-1})$ for $s=1$ and it is slower than $O_P(n^{-1})$ for $s<1$. This suggests that the following statistic

may have a well-defined asymptotic distribution under $H_0$.\footnote{The statistic is similar to Wald's statistics which under homoskedasticity roughly corresponds to $s=-1/2$ and is not expected to converge at the parametric rate.}

Let $V = \mathbb{E}[\varepsilon^2 X \otimes X]$ be the variance operator with spectral decomposition

and let $Z_j \overset{\text{i.i.d.}}{\sim} {N}(0,1)$. The following theorem establishes the limiting distribution of $T_n$ under the null and alternative hypotheses.

Theorem (ref) implies that the test based on $T_n$ has correct size under the null hypothesis and is consistent under the fixed alternative hypotheses which we state as a trivial corollary below. To that end, let $z_{1-\alpha}$ be the $1-\alpha$ quantile of

The test rejects $H_0$ if $T_n>z_{1-\alpha}$. Then, the following result holds:

It is also easy to show that the critical value $z_{1-\alpha}$ can be consistently estimated with simple nonparametric bootstrap based on the i.i.d. draws from the sample $(Y_i,X_i)_{i=1}^n$. Alternatively, one can simulate the asymptotic critical values using the eigenvalues $\hat\omega_j$ of the sample variance operator $\hat V$.

A Confidence Set

Inverting the test, we obtain a confidence set

where $z_{1-\alpha}$ is the quantile of order $1-\alpha$ of $T=\sum_{j=1}^\infty\omega_jZ_j^2$. Indeed, since $T_n\xrightarrow{d}T$, the confidence set has the coverage probability $1-\alpha$:

Since $\beta$ is an infinite-dimensional object, computing this confidence set numerically requires finite-dimensional approximations which can be done as follows. Let $(h_j)_{j=1}^\infty$ be a basis of $\mathbb{H}$. Then we approximate $z_{1-\alpha}$ numerically by computing

For instance, in the functional data analysis, we often have $\mathbb{H}=L_2[0,1]$ in which case, we can use the Fourier basis or Legendre polynomials. If $\mathbb{H}=L_2(\mathbb{R})$, the Hermite polynomials are a natural choice; see Supplementary Material, Section (ref) for numerical implementation and simulation results.

Monte Carlo Experiments

In this section, we conduct several Monte Carlo experiments to evaluate the finite sample performance of the PLS estimator. We simulate the i.i.d. samples $(Y_i,X_i)_{i=1}^n$ from the functional linear model

where the predictors $X_i$ belong to the Hilbert space of square-integrable functions with respect to the Lebesgue measure, denoted $\mathbb{H}=L^{2}[0,1]$. The functional predictor is generated as

We specify the slope parameter $\beta\in L_2[0,1]$ and the spectrum $(\lambda_j,v_j)_{j\geq 1}$ correspond to one of the following three models:

All three models satisfy Assumption (ref) with the same complexity parameter. Model 2 emphasizes the importance of the first five coefficients of $\beta_j$, making the estimation and prediction tasks more challenging. Model 3 introduces repeated eigenvalues, causing the first five eigenvectors to be non-identifiable, although the slope parameter remains identifiable.

We compute the PLS estimator using Algorithm (ref), which is numerically equivalent to the estimator given by equation ((ref)); see hanke1995conjugate, Algorithm 2.1 and Proposition 2.1. This approach avoids direct inversion of the empirical operator $\hat K$ by utilizing iterative multiplication, making it suitable for solving high-dimensional linear systems of the form $\hat K\hat\beta = \hat r$ with a symmetric matrix $\hat K$.

The integrals in inner products and the operator $K$ are discretized using a simple Riemann sum approximation over a grid of $T=200$ equidistant points in $[0,1]$.\footnote{This provides a satisfactory approximation under weak assumption. Alternatively, one could employ quadrature rules with fewer points to reduce computational load.} The experiments feature $5,000$ replications, with each replication generating samples of size $n=100$. For each simulation experiment and an estimator $\hat\beta$, we compute:

We compare the performance of our functional PLS estimator with the number of components chosen using our early stopping rule, against several alternative methods:\footnote{See Supplementary Material, Section (ref), for more details on the practical implementation.}

Note that crambes2009smoothing, Proposition 2, demonstrates that the prediction error for the smoothing spline estimator is asymptotically equivalent to the infeasible optimal selection of tuning parameters. On the other hand, the kernel function employed in yuan2010reproducing is not adaptively selected from the data. We are not aware of any adaptivity results for other competing estimators, including PCA and the alternative PLS estimator with cross-validation. In contrast, our method for selecting the number of PLS components is adaptive for both estimation and prediction according to Theorem (ref).

Figure (ref) summarizes the distributions of estimation (ISE) and prediction (MSPE) errors from $5,000$ simulations. The early stopped PLS estimator achieves the lowest median estimation error, with consistently lower variability across all models. The penalized spline estimator performs effectively when coefficients decline rapidly but exhibits notably larger errors in Model 2, which involves high-frequency slope coefficients. Nevertheless, it shows excellent predictive performance consistent with crambes2009smoothing. The PLS estimator also demonstrates competitive predictive accuracy, although it does not consistently yield the lowest median prediction. The RKHS estimator achieves good prediction performance but inferior estimation across all models. The alternative PLS estimator generally performs comparably to our proposed PLS approach but exhibits notably larger prediction errors in Models 2 and 3. These higher errors may be attributed to numerical instability associated with finite-precision arithmetic, as discussed in delaigle2012methodology. Finally, PCA demonstrates relatively modest performance across all the considered models. Moreover, we find in simulations that the FPLS estimator outperforms FPCA with a smaller number of components; see frank1993.

Next, we evaluate the performance of the proposed test. We first assess whether the asymptotic distribution provides a good approximation to the finite-sample distribution under the null hypothesis $H_0$. Given that the PLS fitting error does not substantially decrease after approximately 10 components, we employ PLS with $m=70$ components to ensure that conditions of Theorem (ref) are satisfied.\footnote{We find that for these models the fitting error does not substantially decrease after $m=10$ components.} The results are presented in Figure (ref).

For each of the three models, the left columns in Figure (ref) display the finite-sample distribution of $T_n$ (in blue) under $H_0$, overlaid with the simulated $\sum_{j=1}^{100} \omega_j Z_j^2$ (in orange) obtained from Theorem (ref). The right columns present the corresponding QQ-plots, comparing the empirical quantiles of $T_n$ to the asymptotic quantiles of $T$. Overall, the simulated asymptotic distribution closely matches the finite-sample distribution, with only minor discrepancies observed in the extreme right tail.

To study the finite sample power of the test, we simulate the power curves corresponding to the local alternatives $\beta(s)+\delta s$ with deviations measured by the scale factor $\delta\in[-1,1]$. $\delta=0$ corresponds to $H_0$ while $|\delta|>0$ to $H_1$. We also consider doubling the sample size from $n=100$ to $n=200$. The results are displayed on Figure (ref), confirming that the test has more power once the null and the alternative hypotheses become sufficiently separated. The power also increases with the sample size as expected. Lastly, we provide additional simulation results in the Supplementary Material.

Nonlinear Temperature Effects in US Agriculture

The global surface temperature has increased by 1.1°C above pre-industrial levels and could increase up to 3.6°C to 4.5°C by the end of the century if current $CO_2$ emissions rise steadily according to the latest studies; see lee2023ipcc. The global warming will likely lead to more frequent and severe heatwaves, altered precipitation patterns, and intensified droughts. Of all major sectors, agriculture is arguably the most sensitive to climate change. While constituting a modest share of developed economies, it is vital for food security. Indeed, the intensified droughts could cause food shortages which in turn may potentially exacerbate mass migration and violent conflicts. Some have argued that the current climates are already warmer than is optimal for agriculture in many parts of Asia, Africa, and Latin America; see nordhaus2013climate.

Determining the precise functional form of the relationship between crop yields and temperature has recently attracted lots of attention; see schlenker2006nonlinear,schlenker2009nonlinear and cui2024model.\footnote{The influential study of schlenker2009nonlinear has more than 4,000 Google Scholar citations at the time of writing.} We argue that the methodology used to estimate such nonlinear temperature effects can be understood as a functional linear regression, where the outcome $Y_{it}$ is the log yield of a crop of a county $i$ in a year $t$, measured in bushels per acre, and the functional regressor $(X_{it}(s))_{s\in[0,40]}$ is a temperature curve, representing the crop exposure to temperatures between 0°C to 40°C during the growing season. Following schlenker2009nonlinear, the temperature curves are computed at discrete points $s_j\in\{0,1,\dots,40\}$ as $X_{it}(s_j) = \Phi_{it}(s_j+1) - \Phi_{it}(s_j)$, where $\Phi_{it}(s_j)$ is the length of time (measured in days) the crop was continuously exposed to temperature larger than $s_j$ for county $i$ in year $t$.

We focus on corn and soybeans which are the two major crops grown in the US. The dataset is comprised of fine-scale county-level crop yields and weather outcomes, spanning US counties east of the 100 degree meridian from 1950 to 2020.\footnote{The dataset is publicly available at the time of writing at \url{www.wolfram-schlenker.info/replicationFiles/SchlenkerRoberts2009.zip}.} We use the same set of controls as in schlenker2009nonlinear, namely: a constant, precipitation measured in mm from March through August, precipitation$^2$, county fixed effects, and a state-specific quadratic time trend to capture technological change. The crop yields $Y$ and the temperature curve $X$ are regressed on these controls to obtain the residuals which are subsequently used for the functional data analysis.

The slope coefficient is then estimated using: 1) our functional PLS estimator; and 2) a highly parameterized least-squares estimator with a step function approximation as in schlenker2009nonlinear. The latter fits a separate temperature effect for each 3°C bin from 0°C to 40°C, hence, it involves 13 parameters. On the other hand, our early stopping rule finds $\hat m = 4$ functional PLS components both for corn and soybeans; see Appendix Section (ref) for more details on the implementation.

Figure (ref) displays the estimated functional slope coefficient $\beta$ corresponding to our functional PLS (red cure) and step function approximation (black dash) for corn and soybeans. We find that the critical temperature after which the crop yields start declining is around 29-30°C which is similar to findings reported in schlenker2009nonlinear.

Lastly, we look at how the nonlinear temperature effects have changed over time. Figure (ref) reports the estimated functional slope coefficient splitting the data into three subsamples: 1950-1973 (blue dot), 1974-1997 (red dash), and 1998-2020 (green curve). The results indicate that the negative temperature effects were larger during 1950-1973 compared to the most recent 22 years, especially for the extreme temperatures. The mitigation of extreme temperature effects may come from two sources: the adaptation and the $CO_2$ fertilization. The $CO_2$ fertilization effects observed in our sample are likely to be small; see nordhaus2013climate who argues that doubling the atmospheric concentration of $CO_2$ would increase crop yields by $10$-$15\%$ only. In contrast, the adaptation effect is likely to dominate over time. It can be attributed to the actions taken by farmers, such as adjusting the sowing and harvesting dates to maximize yields, using more resilient crops, or building efficient irrigation systems. Our results, therefore, suggest some evidence of adaptation in US agriculture which has been reported in burke2016adaptation without properly accounting for nonlinearities. However, our confidence sets are too wide to have conclusive evidence when accounting for statistical uncertainty because there are only few observations of extreme temperatures observed in the data.

Conclusions

This paper proposes a new formulation of the functional PLS estimator related to the conjugate gradient method applied to an ill-posed inverse problem with a self-adjoint operator. We provide the first optimality result for functional PLS and consider a rate-adaptive early stopping rule to select the optimal number of functional components. The estimator has good estimation and prediction properties for a smaller number of principal components than PCA and the early stopping rule performs well in simulations. We find in an empirical application that the nonlinear temperature effects on crop yields have slightly decreased since 1950 which provides some evidence for adaptation of the US agriculture. However, this evidence is not conclusive due to statistical uncertainty arising from the limited number of observations of extreme temperatures in the data. Future studies need to develop more efficient methods to deal with this problem.

\spacingset{1.9}

\setcounter{section}{0}

Notation and Preliminary Results

In this section, we describe the notation and collect several propositions and lemmas.

\paragraph{Notation:} For two sequences $(a_n)_{n\geq 1}$ and $(b_n)_{n\geq 1}$, we will use $a_n\lesssim b_n$ if there exists a constant $c>0$ such that $a_n\leq cb_n$ for all $n\geq 1$. We will also use $a_n\sim b_n$ if $a_n\lesssim b_n$ and $b_n\lesssim a_n$. For two real numbers, we use $a\wedge b = \min\{a,b\}$ and $a\vee b = \max\{a,b\}$.

The following proposition states that the PLS estimator $\hat\beta_m$ is unique for every $m\leq n_*$ and that the tuning parameter selected in Assumption (ref) does not exceed the number of unique non-zero eigenvalues, $n_*$; see also blanchard2016convergence for a kernel regression model setting.

We will need the following tail inequality in Hilbert spaces.

Lemma (ref) allows us to control the tail probabilities for the PLS residual as well as the covariance operator errors on an event with probability at least $1-\gamma$.

We will also use the following two inequalities known in the perturbation theory.

As an immediate consequence of Lemmas (ref) and (ref), since $\|.\|_{\rm op}\leq \|.\|_{\rm HS}$, for every $\gamma\in(0,1)$, we have

on an event that holds with probability at least $1-\gamma$, where $\nu = \|\hat K\|_{\rm op}\vee \|K\|_{\rm op}$.

The following Lemma presents some useful results on the residual polynomials $\hat Q_m(\lambda)=1-\lambda\hat P_m(\lambda)$; see engl1996regularization and hanke1995conjugate. For the completeness of the presentation, we sketch proofs for the key results and refer to the aforementioned monographs for others.

In what follows, for $k\in\mathbb{Z}$, consider the following measure

where $\delta_{x}$ is the Dirac measure at $x\in\mathbb{R}$. For $\phi,\psi:[0,\hat\lambda_1]\to \mathbb{R}$, define

Lastly, let

be the orthogonal projection operators on the eigenspaces of $\hat K$ corresponding to eigenvalues smaller or equal to $a$.

The following lemma allows us to control the residuals of the PLS estimator.

The next lemma provides an upper bound for the derivative of the residual polynomial of degree selected by the stopping rule in Assumption (ref) with some fixed $\delta\in(0,1)$.

Proofs of Main Results

In this section, we provide detailed proofs of theorems.

Comparison to PCA

In this section, we shed some light on the behavior of functional PLS relative to PCA. We will show that for the same fixed number of components $m$, PLS fits the empirical moment better than PCA, hence, it may require a smaller number of components to obtain a comparable fit. We also show that the regularization bias part of the estimation and prediction risk of PLS is smaller than the one of the PCA. Therefore, the adaptive PLS basis is better suited for approximating the slope coefficient.

In what follows, we will use

to denote the functional PLS and PCA estimators. Note that the PLS estimator uses supervised regularization $\hat P_m$ while for the PCA estimator the regularization is fixed to select the terms related to the inverse of the largest $m$ eigenvalues of $\hat K$. We will also use

to denote the population counterparts.

The first part of Theorem (ref) shows that the PLS estimator fits the data better than PCA for the same number of components $1\leq m\leq n_*$. This is the functional version of a result of jong1993pls; see also phatak2002exploiting and blazere2014unified. For the second part of Theorem (ref), it is worth recalling that the estimation and prediction errors in Theorem (ref) can be decomposed as

where the second term is the so-called regularization bias. This shows that the PLS basis is more adapted for approximating the slope $\beta$ than the PCA basis.

Comparison to Standard PLS

We provide some comparison between our PLS which is a functional conjugate gradient method applied to the moment equation with the standard PLS. Note that the standard PLS solves

where $\mathbb{H}_{m}=\mathrm{span}\left \{ \hat{r},\hat{K}\hat{r},..., \hat{K}^{m-1}\hat{r}\right \}$. Equivalently, we solve \[ \min_{\alpha_1,\dots,\alpha_m\in\mathbb{R} }\left\| \mathbf{y} - \sum_{j=1}^{m}\alpha _{j}T_{n}\hat{K}^{j-1} \hat{r}\right\|_n^{2}. \] The first order conditions give

Remark that

where $T_{n}T_{n}^{\ast }$ is the $n\times n$ matrix with $\left( i,j\right) $ element $\left \langle X_{i},X_{j}\right \rangle$. Let $v_{1}$ be the $ m\times 1$ vector and $M_{1}$ the $m\times m$ matrix such that

Equation ((ref)) can be rewritten as $M_{1}\alpha =v_{1}$ where $\alpha = \left[ \alpha _{1},\alpha _{2},...,\alpha _{m}\right]^\top.$ So that $ \hat\alpha _{PLS}=M_{1}^{-1}v_{1},$ provided that the matrix $M_1$ is invertible.

Our version of PLS solves \[ \min_{b\in \mathbb{H}_{m}}\left \Vert T_{n}^{\ast }\mathbf{y}-T_{n}^{\ast }T_{n}b\right \Vert ^{2} \] where $\mathbb{H}_{m}=\mathrm{span}\left \{ \hat{r},\hat{K}\hat{r},..., \hat{K}^{m-1}\hat{r}\right \} .$ This can be understood as the PLS with respect to the weighted norm. Equivalently, we want to solve \[ \min_{\alpha_1,\dots,\alpha_m\in\mathbb{R} }\left \Vert T_{n}^{\ast }\mathbf{y}-\sum_{j=1}^{m}\alpha _{j}\hat{K} ^{j}\hat{r}\right \Vert ^{2}. \] Using the first order condition, we obtain the equation $v_{2}=M_{2}\alpha $ with

If we compare $v_{1}$ and $v_{2}$, $M_{1}$ and $M_{2},$ we see that the only difference is in the power of $\left( T_{n}T_{n}^{\ast }\right) $. Note also that $M_{1}$ and $M_{2}$ are Hankel matrices.

Additional Simulations: Early Stopping and Confidence Sets

In this section we report additional simulation results for our early stopping rule and confidence sets.

Early Stopping

Recall that our PLS amounts to fitting the norm of a “sample moment" (or a score) $$\left\| \hat{r}-\hat{K}\hat{\beta }_{m}\right\|.$$ The norm decreases monotonically to zero in $m$ and it becomes zero when the number of conjugate gradient steps reaches the number of non-zero eigenvalues of $\hat K$, i.e. $m=n_*$. To prevent overfitting, we select the number of PLS components as the first value of $m$ for which the norm of residual drops below a certain threshold:

where $\tau >1$ is a constant and $1-\delta$ is the confidence level of the rule; see Assumption (ref). To describe the practical implementation, we focus on the functional linear model

We compute the fitted moment as:

To compute the threshold in equation ((ref)), we estimate $\ensuremath{\mathbf{E}}\left \Vert X\right \Vert ^{2}$ by $\frac{1}{n}\sum_{i=1}^{n}\int_0^1 X_{i}^{2}\left( s\right) ds$ and $\sigma ^{2}$ by $\hat\sigma^2 = \frac{1}{n}\sum_{i=1}^{n}\hat{\varepsilon}_{i}^{2}$, where

and $\hat{\beta }$ is a preliminary estimator of $\beta $, described below. All integrals are discretized at the uniform grid of $200$ points.

We consider an iterative approach for estimating $\sigma ^{2}$ inspired by chernozhukov2024applied, Section 3.A. To that end, we first obtain a pilot estimator of $\beta$, denoted $\hat\beta_0$. Next, we compute the estimator of $\sigma ^{2}:$ \[ \hat{\sigma}_{0}^{2}=\frac{1}{n}\sum_{i=1}^{n}\left( Y_{i}-\langle X_i,\hat\beta_0\rangle\right) ^{2}. \] If $n>>T$, we can use the OLS estimator. Otherwise, any other regularized estimator can be used, e.g. PCA with generalized cross-validation. We set $k=0$ and specify a small constant $\xi \geq 0$ as a tolerance level and the maximum number of iterations $k_{\max }.$ The iterative procedure is described in Algorithm (ref).

We use the following values in all simulation designs: $\delta=0.1$, $k_{\max}=10$, and $\xi=0.01$. These values correspond to a $90\%$ confidence level, a sufficiently high number of iterations, and a sufficiently low tolerance level respectively. Our theory also tells us that we need $\tau>1$, so we set $\tau=1.01$. The same values are also used in the empirical application. We noticed in simulations that the algorithm tends to overestimate $\sigma^2$ which leads to underestimated values of $m$. Figure (ref) reports the population slope coefficient (solid black) and the average values of estimated slope coefficients $\hat\beta_{\hat m}$ (dashed red) when the selected number of PLS components $\hat m$ for the four simulation designs considered in the paper, Section (ref). We also report the median values of $\hat m$ selected by our data-driven adaptive rule in each case. Overall, we can see that the PLS estimator with the number of components selected using early stopping can successfully recover the global shape of the slope parameter across all simulation designs. It is worth noting that the early stopping rule selects the number of PLS components to achieve the minimax-optimal convergence rate for both the mean integrated squared error (MISE) and the mean squared prediction error (MSPE). In doing so, it simultaneously balances and minimizes the sum of the squared bias and the variance. This behavior is consistent with the results shown in Figure (ref), where a visible bias is observed as part of the optimal bias-variance trade-off.

Figure (ref) visualizes our early stopping rule. The left vertical axis corresponds to the values of MSPE (solid black) while the right vertical axis to the early stopping rule (blue). The number of PLS components is selected as the first $m$ (blue dot) when the norm of the fitted moment (dashed blue) drops below the value of the threshold (dotted blue). We also plot the minimum (black dot) of the simulated MSPE (solid black) together with $90\%$ confidence band (shaded gray). Given the uncertainty behind the lowest value of the MSPE as well as the norm and the threshold, the early stopping rule produces reasonable values of tuning parameters that are compatible with the global minimum of MSPE once the statistical uncertainty is accounted for. The simulation results also confirm that there is no reason to select higher values of $\tau$ since they would increase the threshold and lead to more convervative choices of $\hat m$. We conclude that our early stopping rule is prone to oversmoothing.

Confidence Sets

To compute the confidence sets using test inversion, we note that if $(h_j)_{j=1}^\infty$ is a basis of $L_2[0,1]$, then

with coefficients $b_j = \langle\beta,h_j\rangle$. Since $b_j\downarrow 0$, to reduce the computational cost, we search over functions $\beta$, where the sum is truncated to include only the first five basis functions. We use the cosine basis and create a uniform grid of $20$ points on $[0,4.5]^5$ corresponding to the first five coefficients $b_1,\dots,b_5$. We simulate the critical value using the approximation to the asymptotic distribution $\sum_{j=1}^{100}\lambda_jZ_j^2$, where $(\lambda_j)_{j=1}^\infty$ are the eigenvalues of $K$, and we use $50,000$ replications to compute the critical value $z_{0.95}$.

Figure (ref) displays the $95\%$ confidence sets. These confidence sets are quite informative about the global shape properties of the estimated slope parameter.